## Data Science Done Right (the Kitchen Style) #12

Example of using PCA for signal separation problem

Before going into more details of more sophisticated multidimensional scaling techniques, and criteria for selecting principal components/factors, let’s take a look at data that has the same measurement units in all their dimensions, and relatively easy to guess what criteria we may use to distinguish the real-life principal components. The dataset reference (https://homes.esat.kuleuven.be/~smc/daisy/daisydata.html, De Moor B.L.R. (ed.), DaISy: Database for the Identification of Systems, Department of Electrical Engineering, ESAT/STADIUS, KU Leuven, Belgium), we’ll be using, is borrowed from (A.J.Izenman 2013). These are the 8 channel ECG readings of a pregnant woman. What we may want to do with these data is to separate those mixed “crazy statistical” variables, we were talking before, into proper analytical or real-life tangible variables, like reading of mother’s and child’s heartbeats, perspiration rhythms, etc.

Fig. 12.1. 8 channel ECG of a pregnant woman with mixed in child’s heartbeat.

This is a typical “cocktail party problem” in which we want to separate voices of the party participants out of recordings made from differently placed microphones in the party-room while all invited people babble simultaneously. Of course, more serious applications of the problem may include submarine or drone detection out of sonar or radar readings, electric grid invasions, sources of rocket engine explosions out of sensor data, etc…, you name it.

dim <- c(“V2″,”V3″,”V4”, “V5″,”V6”, “V7”, “V8”, “V9”)
ekg2 <- ekg[, dim]
ds <- ks_eigen_rotate_cov(ekg2)
names(ds) <- c(“V2″,”V3″,”V4”, “V5″,”V6”, “V7”, “V8”, “V9″)
ds[,”V1″] <- ekg[,”V1”]
describe(ds)

vars n mean       sd median trimmed mad min max range                    skew kurtosis
V2 1 2500 1.57   215.17 17.69 19.92 53.02 -1366.55 474.04 1840.59   -3.69 18.36
V3 2 2500 0.03 44.47 -3.04 -2.88 24.02 -172.36 286.30 458.66         2.51 13.30
V4 3 2500 0.07 19.66 -0.85 -0.60 14.62 -104.64 147.41 252.05           0.74 6.58
V5 4 2500 0.18  6.13 0.46 0.48 4.60 -34.43 26.13 60.57                      -0.65 2.75
V6 5 2500 0.23  5.36 0.86 0.28 5.04 -25.50 21.63 47.13                      -0.04 1.04
V7 6 2500 0.10  3.32 0.01 -0.02 3.03 -11.75 16.15 27.90                         0.53 1.44
V8 7 2500 0.04  2.23 0.06 0.06 2.11 -8.35 7.98 16.33                            -0.12 0.39
V9 8 2500 0.02  2.01 0.07 0.04 2.03 -6.83 7.17 14.01                            -0.11 -0.03

Fig. 12.2. Covariance transformation matrix based eigenbasis, with V2 (blue) – mother’s and V6 (green) – child’s separated heartbeats (and perspiration), and V9 (red) noise.

As one can easily see, that even simple eigenbasis rotation around the covariance matrix transformation allows us to see clearly mother’s heartbeat (14 of them) (dark blue V2, also V3 and V4 with additional signals), child’s heartbeat (22 of them) and perspiration rhythm (dark green V6), and noise V9, with V5, V7, V8 as intermediate mixes of various proportion of these three signals.

While the formal parametric descriptive statistics may be still hard to use as a guiding indicators (except obvious V2 with really high variance, skewness and fat tails), non-gaussian shape of the distributions hints more decisively which variables may be not  noise and we may be interesting in.

Fig. 12.3. Pdf’s of the synthetic eigenbasis variables.

Appendix 12.1

As usual, working files are in: https://github.com/NSelitskaya/kitchen-style-r

## Data Science Done Right (the Kitchen Style) #11

Meanings of Covariance Metrics

The Covariance Family metrics (including coefficient of correlation and r squared) may be interpreted in various ways, which indicates that may be poorly understood, hence, and inappropriately used. For example, the usual interpretation of covariance/correlation is a measure of the variables’ dependency may be quite confusing – on the one hand we may want to reduce it, to deal with independent variables, on the other – we do want to find that dependency. Of course we may say we want independent explanatory variables, and dependent response, but that division is quite subjective, therefore that dependency interpretation of covariance/correlation is, at least, subjective, too.

Other interpretations, such as relation to regression model or paying attention to their moment functions nature, may be more useful. We may envision the covariance family (as well as, for comparison, rss/sse or variance/std. deviation) mechanistically. For example, rss (residual sum of squares) could be envisioned as little springs/ rubber bands attached on one side to the free floating regression bar, and on another – to data points. Equilibrium of such a configuration will give as regression bar position we seek – the one with minimal “internal tension” of the model. Using the same pattern, we may envision the same little springs/rubber bands attached to data points on one their end, but, now, the other end would be attached to a coordinate line with a rotation point fixed in the mean point, such way creating a rotation moment.

It’s quite interesting what will be the equilibrium point for the system. If we let go our coordinate system it will rotate to the eigenbasis coordinate system (relative to our covariance matrix of course) – such a system that sum of all rotation moments from the individual data points will be 0, i.e. intuitively we can see (we’ll demonstrate that more rigorously later) that all the springs/bends will be in a minimal stressed state, which means that one of our new axes will be parallel to the regression line (just with the different baseline of the “minimal tension”, which is defined by the mean rotation point constraint), while another (2D case would be easier to intuitively imagine) will be orthogonal (for orthogonal basis, which is usually the case with the bases of abstract statistical spaces).

When Covariance Analysis may be useful

When we model real-life processes we may encounter situations, especially in “soft” sciences like psychology or medicine, or pseudosciences 🙂 like sociology or economy, when we are able to measure some parameters, but suspect that there exist more primal “hidden” or “latent” parameters which we can’t measure directly, such as Intelligence Level, or Health Rank, or Social Stability or Economic “Health”. Of course in hard sciences we measure many parameters indirectly, for example already mentioned forces, accelerations, or masses are usually measured via deformations (of a spring) or stresses (of a piezo crystal) with the following electromagnetic field generation, which, in its turn, can be measured via charged particles displacement/rotation (of inductor), etc… But in hard sciences we usually have analytical models developed which could be rigorously verified and reproduced, and we can get reliable indirect inferences of one parameter via direct measurement of others, so such inference can be seen as a mere inconvenience (or even convenience when (more) direct measurement is possible, but more complicated/costly).

If we don’t have a proven analytical model, which is usual from “soft” sciences, or in the beginning of research in hard science, or our culture/civilization (of statisticians and data scientists :)) not yet developed tradition of advanced analytical models, even choice of the observable parameters may be incomplete or redundant. For instance, in our previous example of the “crazy statistical” 2nd Newton’s law of motion linear modelling (applied to throwing books out of a window) b + c = d + e, we created the model based on 1..i-th..n observations of bi, ci, di, ei parameters made in respect to an abstract statistical (orthogonal) basis i,j,k,l in a statistical space A, where bi = bi*i, ci = ci*j, etc… Let’s imagine these parameters were measured by weird instruments (for example by cameras with some video image processing) made by those who had no physics insights on what’s going on.

If we take a 4-dimensional space, linear mode of the type b + c = d + e, or b*i + c*j – d*k – e*l = 0 (or a typical generic form y = b1*x1 + b2*x2 + … + b0*y), would mean that out of all possible members of the space we select only those that satisfy the above equations, which, actually mean that for the selected data points the 4-dimensional basis is excessive, because linear model equations are also conditions for the linear dependence of the (basis) vectors. Which means linear models describe subspaces of, at least, 1 less dimensionality than the original space. Which subspaces we select to project our original data to is a separate question. We may use methods of the least deltas (residues) between our data and their projections into a candidate subspace, or least angles between original data hull and subspace hyperplane, or various weighted methods, but for understanding PCA analysis we want to work with the least RSS algorithm in which residues are measured in orthogonal to the subspace hyperplane direction (say, b + c = d + e model/subspace/hyperplane was chosen using that method).

Had we have some insights on which hidden physics parameters are involved into 2nd Newton’s law of motion (and what it is), we could come up with mappings of our “crazy statistical” variables  to the following hidden physics parameters: Fm/2*w = F*m/2*w; c = 2*a = 2*a*u; d = m*a = m*a*u; e = – m/2*w + 2*a – m/2*w + 2*a*u,  (1) (relative to some spacial basis, say u,v,w of a physical space B, where u is a basis vector of Lagrangian type attached to a book and collinear to the acceleration of a book, v – to the force applied to it (say, we don’t know that u and v are linearly dependent), and w is, say, an orientation of our mass-measuring instrument (yeah, imagine we measure a mass vector, but we are “crazy statisticians”, aren’t we?).

Knowing these hidden physics variables are related via some law of motion (which we pretend we don’t know yet), we would perhaps wanted to shrink our initial abstract statistical basis i,j,k,l of A-spaceand rotate it to a smaller latent statistical basis f,a of abstract statistical space C which is more conveniently mapped to the basis of hidden physics variables u,v,w of B-space. Here, let’s call the unknown physics variables “hidden”, and corresponding unknown abstract statistic variables “latent” (of course these are interchangeable terms in literature). Say, we want mapping g:B->C by a simple diagonal unit matrix, such that f=1*v, a=1*u and 0*w. That may seem as a redundant transformation, but we still want to have B and C spaces distinct, for example because abstract statistical basis is orthogonal, but spacial u,v,w is not necessarily, and maybe not really a basis (linearly independent) after all.

Knowing that g*f = h, where f:A->B and h:A->C, i.e. f(b+c-e) = F = F*v, g(F) = F*f; f(d) = m*a*u, g(m*a*u) = m*a*a; and mappings are linear (f(a)+f(b)=f(a+b)), we would expect that h(b+c-d-e) = H*(b+c-d-e) = F*f – m*a*a=0.

If we transform our dataset of bi,..ei tuples into C space where they become fi,ai pairs, first, and then model them by a linear model that minimizes squared sums of epsilon orthogonal to the regression line, we’ll get the same model (as soon all transformations are homomorphic, or in vector spaces linear, which they are): = C*a. Although we may think that C=m stands for the mass of the book in our predictive model, however that won’t be a parameter of the linear model, but rather its coefficient for the parameter a. If we use standard RSS minimization regression, it will be an averaged mass ma that allows our model to go through the mean values of F and a of our dataset: ma=muF/mua (see below).

Regression line we just found would be parallel to one of the eigenbasis vectors (while epsilon-line would be parallel to another eigenvector), therefore rotating our bais f,a to the one spawning regression end epsilon lines, via, say, transformation e (represented by matrix E), would be transformation to the eigen basis v1,v2 if we also shift origin of our coordinates. Therefore, again, because of the linearity of transformations, our regression would be like e(F-m*a) = E*(F-m*a) =E*(F) – E(m*a), reg = (e11*F – e12*m*a)*v1, eps = (e21*F – e22*m*a)*v2 = 0 (therefore e21=e22=0);   and, because regression in eigen coordinates would be var1 = b1*var2 + b0, where b1=0, then e11*F – e12*m*a – C’1 = 0.

On the other hand, we could have rotated basis of A space, in which our data represented by bn..en tuples, into eigenbasis of D space, in which n-th tuple would be var1n, var2n,…var4n, where var1n = (a11*bn + a12*cn + a13*dn + a14*en)*v1, etc…, or

(a11   a12  a13   a14)  (b)  i       var1   v1
(a21  a22  a23  a24)  (c)  j   =  var2  v2
(a31  a32  a33  a34)  (d)  k      var3  v3
(a31  a32  a33  a34)  (e)  l       var4  v4

Now, again we could have tried to find a linear model of our data in the eigenspace D with eigenbasis v1,v2,v3,v4, but because we know that all regression lines/planes/hyperplanes in eigenspace will be parallel to some eigenvectors (and orthogonal to another), one of the regressions vari = b1*varj + b0, where b1=0, therefore vari – C = 0 will be the same regression we found above in space C with f,a basis, and converted also into subspace of eigenspace D, we can write ai1*b + ai2*c + ai3*d + ai4*e = reg = e11*F – e12*m*a – C”1.

So, by performing PCA, i.e. finding eigenbasis relative to covariance matrix, we implicitly find the least lossy linear models of the hidden variables that describe our data, and use them as an eigenbasis for a “linear model subspace”, while the rest of the eigenvectors compose an orthogonal “error subspace”, which direct sum produces whole space our data reside in. In the case above we have only one set of hidden parameters that are part of one linear model equation, therefore we have only one-dimensional eigenbasis of the “linear model subspace”.

Haven’t we known nothing about hidden variables, and would rotate basis i,j,k,l around covariance matrix of our dataset into the same eigenbasis e1 of the “linear model subspace”, and would get additional e3,e4 eigenvectors to already known e2 eigenbasis vector of the error subspace. However, we would not know which eigenbasis vector is part of the “linear model subspace”, and which – “error subspace”. Although, having some guesses about distributions of the linear experimental parameters (for example that we through our books out of window with steadily increased force, and expecting uniform distribution of the data in linear model dimension), also of the non-linear parameters (for example we know that there are standard book sizes and weights, around which modes of book weight distribution would congregate), or the noise, coming from numerous minor sources we don’t know and care, and therefore (according to central limit theorem) Gaussian. The latter dimensions we, perhaps, want either drop, or set data in them to 0 or other baseline, and data in dimensions with mixed Gaussian and non-linear error – “de-Gaussianize”, and sometimes rotate back into original basis.

More rigorous on Regression and Covariance connections

A matrix can always be seen as a particular transformation of an object from one space (with one particular basis) into another (with other basis). Particular, covariance matrix S, can be seen as a transformation, defined by our dataset, of, for example the unit vector x of the basis our dataset expressed in ijk.., into object y of another space – space where basis vectors uvw… are defined via sums of ij, jk, ki etc,… covariances of our data set:

y = Sx, s.t.:

(sii^2  sij  sik …)   (x1) i     (vi * x1 + sij * x2 + sik * x3 + …*…)     u
(sji  sjj^2  sjk …)   (x2) j = (sji * x1 + vj * x2 + sjk * x3 + …*…)     v
(ski skj skk^2 …) (x3) k    (ski * x1 + skj * x2 + vk * x3 + …*…w

However, for a given matrix-transformation S we may find such a basis i’j’k’, s.t. S-transformation of  x will be equal to its scalar (lambda) multiplication, and, hence, asymmetric covariances in such a basis will be nulls:

y = lambdax, s.t.:

(lambda  0  0 …)   (x1′) i’     (lambda * x1′)   u’
(0  lambda  0 …)   (x2′) j’ = (lambda* x2′)   v’
(0  0  lambda …)  (x3′k’    (lambda * x3′w’

Actually, it’s not the space with basis u’v’w’ we are interested in, but the basis i’j’k’, in which asymmetric pairwise covariances are equal 0’s, or the basis in which (synthetic (or latent)) variables are independent, or, as we noted above, regression line(s) are parallel to basis vectors. Though, it’s not every regression, but quite particular (least RSS) ones. Let’s see for which regression bij = sij / sii^2, or in a more usual for 2D regression notation: b1 = cov(x, y) / var(x). Let’s express residual square sum as:

RSS = Sum i=1..n (e^2) = Sum (y – b1*x – b0)^2

@RSS/@b1 = -2*Sum(x*y -b1*x^2 – b0*x), which we want to set to 0 to find out RSS’s minimum:
Sum(x*y)/n – b1*Sum(x^2)/n – b0*Sum(x)/n= 0

if cov(x, y) = E((x – E(x))(y – E(y))) = E(x*y) – E(x)*E(y) – E(x)*E(y) + E(x)*E(y) = E(x*y) – E(x)*E(y) = Sum (x*y)/n – mux*muy

and var(x) = E(x*x) – E(x)*E(x) = Sum (x*x)/n – mux^2

then:
c0v(x, y) + mux*muy – b1*var(x) – b1*mux^2 – bo*mux= 0

i.e. b1 = cov(x, y)/var(x) (hence b1=0 when cov(x,y)=0) if muy = b1*mux – b0 or mux=0,

which, with bo inference:
@RSS/@b0 = -2*Sum(y – b1*x – b0)/n = 0 = muy – b1*mux – bo,
is exactly the case.

If we generalize derivations above onto multivariable (but still univariate) regression yr = b1*x1 – b2*x2 – … b0, and express dataset’s y as xn with bn coefficient (which we always set to -1 to end up with the “traditional” regression notation), then:

RSS = Sum i=1..n (exn^2) = Sum (bn*xn + b1*x1 + b2*x2 + … + b0)^2

cov(xi, xj) = Sum(xi*xj)/n – muxi*muxj
var(xi) = Sum(xi^2)/n – muxi^2

partial derivatives of RSS would be:

@RSS/@b1 = -2*Sum(x1*(bn*xn + b1*x1 + b2*x2 + … + b0)) = 0 = bn*Sum(x1*xn)/n + b1*Sum(x1^2)/n + b2*Sum(x1*x2)/n + … + b0*Sum(x1)/n
@RSS/@b2 = -2*Sum(x2*(bn*xn + b1*x1 + b2*x2 + … + b0)) = 0 = bn*Sum(x2*xn)/n + b1*Sum(x2*x1)/n + b2*Sum(x2^2)/n + … + b0*Sum(x2)/n

@RSS/@bn = -2*Sum(xn*(bn*xn + b1*x1 + b2*x2 + … + b0)) = 0 = bn*Sum(xn^2)/n + b1*Sum(xn*x1)/n + b2*Sum(xn*x2)/n + … + b0*Sum(xn)/n

@RSS/@b0 = -2*Sum(bn*xn + b1*x1 + b2*x2 + … + b0) = 0 = bn*Sum(xn)/n + b1*Sum(x1)/n + b2*Sum(x2)/n + … + b0

or

b1*var(x1) + b2*cov(x1,x2) + … + bn*cov(x1,xn) + mux1*(b1*mux1^2 + b2*mux2 + … + bn*muxn + bo) = 0
b1*cov(x2,x1) + b2*var(x1) + … + bn*cov(x2,xn) + mux2*(b1*mux2*mux1 + b2*mux1 + … + bn*muxn + bo) = 0

b1*cov(xn,x1) + b2*cov(xn,x1) + … + bn*var(xn) + muxn*(b1*mux1 + b2*mux2 + … bn*muxn + bo) = 0

b1*mux1 + b2*mux2 + … + bn*muxn + bo = 0

then matrix form will be:

cov(x, x) * b = – mux * (muxT * b + bo) = 0

where x = (x1, x2, … xn); b = (b1, b2, … bn); mux = (mux1, mux2, … muxn)

cov(x, x) * b = 0 means that for covariance matrix cov(x, x) = S of the eigenbasis i’j’k’… only non-zero members will be variances, and each bi*var(xi) = 0 means bi must be null, i.e. making least RSS regression (S*b=0, or multivariate S*B=0, where B is NxM matrix where M is number of dimensions for simultaneous regression, and N is the dataset size) lines/planes parallel to the remaining, after regression, eigenbasis vectors (as we can see on Fig 11.1). Which means that the least RSS regression in eigenbasis may be done just by dropping dimensions of our choice.

Orthogonality of the eigenbasis of the symmetric matrix, which covariance matrix is, comes from here:

A=AT – for symmetric matrix

(A*v)T = vT*AT, which is obvious from example:

( (a11 a12)   (b1) )T  =   (b1*a11+b2*a12, b1*a21+b2*a22)
( (a21 a22)  (b2) )

(b1, b2) (a11 a21)     =   (b1*a11+b2*a12, b1*a21+b2*a22)
.              (a12 a22)

A*v1 = lambda1*v1, which is an eigenvector definition, we multiply by another distinct (with different lambda2) transposed eigenvector v2:

lambda1*v2T*v1 = v2T*lambda1*v1 <= v2T*A*v1 => v2T*AT*v1 = (A*v2)T*v1 = (lambda2*v2)T*v1lambda2*v2T*v1, or

lambda1*v2T*v1 – lambda2*v2T*v1 = 0, or (lambda1-lambda2)*v2T*v1 = 0

Fig. 11.1. SPLOM (Scatter Plot Matrix in 10 synthetic variable eigenbasis of the Diabetes dataset

However, now, the question we want to answer would be what dimensions we may want to collapse, or, in the extreme projection just to one dimension, – which one it would be?  As we discussed before, RSS is a measure of the original dataset Structure preservation (or, rather, not preservation, or level of granularity our transformation (regression) becomes non-homeomorphic). In our case of the least RSS regression in eigenbasis, variances and RSSs will be the same. Therefore, it would be reasonable for starters, to collapse those dimensions which variances are smallest – i.e., when doing regression, we lose minimal Structure of the original dataset. However, that approach depends on the variables’ measurement scale, use of normalization, and we may not be after that aspect of the dataset Structure which gives the highest variance. Because some aspects of the data Structure may be inherently random, we may be more interested in those synthetic variables which distribution is not gaussian (assuming that multisource randomness comes in a gaussian form), using either parametric metrics (like skewness or kurtosis), or nonparametric methods like Kolmogorov-Smirnov tests, etc…

If we run an eigenbasis analysis based on covariance matrix for our Diabetes dataset, using all but glucose related variables, and not normalizing data:

dim3 <- c(“chol”, “hdl”, “ratio”, “bp.1s”, “bp.1d”, “age”, “height”, “weight”, “waist”, “hip”)
diabetes3 <- diabetes[complete.cases(diabetes[, dim]), dim3]
ds <- ks_eigen_rotate_cov(diabetes3)

Fig. 11.2. Histograms of 10 synthetic eigenbasis variables in comparison to Glycosylated hemoglobin histograms

We’ll able to see that distributions of the synthetic eigenbasis variables somewhat differ (in central tendencies and dispersion metrics) from the Glycosylated hemoglobin distribution (v1, v3, v4, v5 (especially v3 with obvious opposite skewness)) not only for variables with the high variance, and actually not necessarily do for the high variance variables (v2) (Fig 11.2). Because of the distributions’ differences, we can use those variables for clustering purposes of the high- and low-risk patients.

Those eigenbasis variables, again, show no surprises: high cholesterol and high blood pressure are very bad, age and extra weight are moderately bad, hdl ratio is rather good (maybe in some aspects). Bold indicates more significant, and cursive – “good” variables.

vars n     mean            sd          median trimmed mad min max       range      skew kurtosis
v1 1 377 280.76           45.82 276.67 278.32 39.69 136.48 520.84   384.36    0.81   2.38
v2 2 377 118.88           41.31 115.47 116.17 36.65 22.77 265.49            242.72    0.69  0.75
v3 3 377 -108.09         24.92 -105.24 -106.43 22.85 -216.86 -41.18 175.68   -0.73   1.14
v4 4 377 -60.25           16.49 -58.30 -59.45 14.79 -132.42 -10.14       122.28  -0.52   0.73
v5 5 377 7.80               14.43     6.88 7.19 14.25 -32.48 56.33              88.80    0.45   0.41
v6 6 377 -11.56              8.72  -12.07 -11.74 8.85 -35.94 19.96              55.90    0.27    0.11
v7 7 377 -32.98             4.55  -32.83 -32.93 5.04 -43.54 -19.71            23.82   -0.06 -0.46
v8 8 377 47.39              2.50    47.35 47.38 2.30 37.25 57.28               20.04    0.09   0.99
v9 9 377 35.88              2.03    36.00 35.88 1.76 27.95 44.25               16.30    0.04   1.44
v10 10 377 -2.99           0.67 -2.86 -2.91 0.34 -9.71 -1.64                      8.07     -4.17    32.36

where:

v1 = 0.9399*chol +0.0475*hdl +0.01945*ratio +0.1650*bp.1s +0.0821*bp.1d +0.1083*age  -0.0014*height +0.2552*weight +0.0417*waist +0.03362*hip

v3 = 0.2021*chol -0.0128*hdl +0.0052*ratio -0.8554*bp.1s -0.3213*bp.1d -0.3448*age +0.0089*height +0.0653*weight -0.0221*waist -0.0070*hip

v4 = 0.0565*chol -0.8737*hdl +0.0732*ratio -0.0435*bp.1s -0.2519*bp.1d +0.3873*age -0.0085*height -0.1089*weight +0.0253*waist -0.0106*hip

v5 = -0.0800*chol +0.4424*hdl -0.0345*ratio -0.1739*bp.1s -0.4054*bp.1d +0.7636*age -0.0105*height +0.1208*weight +0.0591*waist +0.0261*hip

Fig. 11.3. Dataset is initially normalized. SPLOM (Scatter Plot Matrix in 10 synthetic variable eigenbasis of the Diabetes dataset

Fig. 11.4. Dataset is initially normalized. Histograms of 10 synthetic eigenbasis variables in comparison to Glycosylated hemoglobin histograms

or, for initially normalized data, high variance (0f v1 and v2 variables) may be a more obvious indicator of the dataset Structure we may interested in:

vars n mean         sd median trimmed mad min max range skew kurtosis se
v1 1 377 2.31        0.31 2.30 2.29 0.32 1.67 3.21 1.54             0.39 -0.16 0.02
v2 2 377 -0.53     0.25 -0.51 -0.52 0.27 -1.45 -0.06 1.39    -0.42 -0.19 0.01
v3 3 377 -1.11       0.19 -1.10 -1.11 0.19 -1.65 -0.47 1.18          -0.15 0.17 0.01
v4 4 377 -1.77      0.17 -1.75 -1.77 0.18 -2.25 -1.30 0.95        -0.07 -0.48 0.01
v5 5 377 1.71        0.16 1.70 1.70 0.15 1.24 2.24 0.99              0.28 0.28 0.01
v6 6 377 0.32      0.13 0.30 0.31 0.11 0.03 1.15 1.12                1.34 5.58 0.01
v7 7 377 0.26       0.08 0.25 0.25 0.08 0.04 0.62 0.59          0.57 1.07 0.00
v8 8 377 -0.21     0.07 -0.21 -0.21 0.07 -0.42 0.05 0.47       0.19 0.33 0.00
v9 9 377 -1.10      0.05 -1.10 -1.10 0.05 -1.30 -0.89 0.41       0.00 1.64 0.00
v10 10 377 -0.14  0.03 -0.14 -0.14 0.02 -0.42 -0.07 0.35    -3.43 23.53 0.00

v1 = 0.065*chol_n -0.1978*hdl_n +0.1322*ratio_n +0.1549*bp.1s_n +0.1814*bp.1d_n +0.1371*age_n +0.0486*height_n +0.5314*weight_n +0.5925*waist_n +0.4767*hip_n

v2 = -0.1774*chol_n -0.1444*hdl_n -0.0270*ratio_n -0.3777*bp.1s_n -0.2264*bp.1d_n -0.8219*age _n+0.1369*height_n +0.2067*weight_n +0.0446*waist_n +0.1175*hip_n

Appendix 11.1

As usual, working files are in: https://github.com/NSelitskaya/kitchen-style-r

## Data Science Done Right (the Kitchen Style) #10

Adding Principal Component Analysis into our cooking mix…

If we stop our stepwise regression from the previous chapter on last two dimensions (don’t look at the eigen parameter yet :)), and draw a scatter plot (Fig.10.1):

dim <- c(“glyhb”, “chol”, “hdl”, “ratio”, “stab.glu”, “bp.1s”, “bp.1d”, “age”, “height”, “weight”, “waist”, “hip”)
diabetes2 <- diabetes[complete.cases(diabetes[, dim]), dim]

ds <- ks_lm_dim_red(diabetes2, dim, n_dim=2, eigen=FALSE)

names(ds)[1:2] <- c (“v1”, “v2”)
ds\$glyhb <- diabetes2\$glyhb
ds_nd <- ds[ds\$glyhb<7,]
ds_pd <- ds[ds\$glyhb>=7,]

ggplot(ds, aes(x=v1, y=v2, colour=glyhb))+
geom_point(data=ds, alpha=0.5)+
geom_point(data=ds_nd, alpha=0.5)+
geom_point(data=ds_pd, alpha=0.5)+

Fig.10.1. Scatter plot of two last synthetic variables for the Diabetes dataset without eigenbasis rotation

We’ll easily see where (at about the bissectrice of the first quadrant) the regression line would be, and what kind of pdf’s of the population and “non-diabetes and “diabetes” groups would be after projection of the data to that line in normal direction (instead of collapsing dimension v2). In such a regression (the multivariate one, which we’ll talk about later) difference in samples’ distributions would be more noticeable than in the univariate regression case (Fig.10.2).

Fig.10.2. Histogram of the Diabetes dataset projection on the last synthetic variable (no eigenbasis rotation)

The reason we may want to prefer normal projection to the being projected regression line is because, although all projections are not a homeomorphic transformations and, therefore, do not preserve structure of the domain space, they may not preserve it on different degrees of granularity. For example (Fig.10.3), if we project point a1 to the line r, while the projection p is a continuous transformation (neighbours in the domain space remain neighbours in the range space, i.e. if distance between two points in the range space is less than the chosen epsilon – radius of the open ball there, we always may chose such a radius of the open ball in the domain space – delta – that the inverse images of those range points will be inside that ball).

Fig.10.4. Scatter plot of two last synthetic variables for the Diabetes dataset with their eigenbasis rotation

But the inverse image p-1 is not continuous (if we take a point a2 that lays on the projection line and the very line r, then its projection a2’=a2=a1′, and if we chose epsilon’ < d(a1,a2), where d is a distance (usual Euclidean) function, then there is no way to find such a small delta’ that would bring images a1=p-1(a1′) and a2=p-1(a2′) into neighbourhood < epsilon’, because delta’ is already 0, and d(a1,a2) >epsilon’). The same applies to the projection p”(a1)=a1″=a3″=a3 to the axis x, with epsilon” < d(a1, a3). However epsilon”= epsilon’ / cos (alpha), alpha being an angle between r and x. Therefore, with epsilon” > epsilon’, we lose structural homeomorphity with the normal projection to the regression line r on a smaller granularity level.

Fig.10.5. Histogram of the Diabetes dataset projection on the last synthetic variable (with eigenbasis rotation of last 2 synthetic variables)

It’s not only the multivariate regression that can help with that – we may also rotate the basis to make one of the basis vectors to spawn the regression line (Fig.10.4), and collapse (with our ready routines) another dimension:

ds <- ks_lm_dim_red(diabetes2, dim, n_dim=2, eigen=FALSE)
ds <- ks_eigen_rotate(ds, std=TRUE)
ds <- ks_lm_dim_red(ds, eigen=FALSE)

What we’ve just done (rotation to eigenbasis) is considered a part of Principal Component analysis. The need for it comes from the very nature of the statistical modelling. I love how Noam Chomsky put it in this lection at Google (33:50):

We can, indeed, throw a book out of a window, and, instead of using analytical mathematical models (based on the 2nd Newton’s law of motion F=ma), make a bunch of video recordings and process them with statistical modeling methods, and even predictive models of the Department of Statistics and Data Analysis may be better than the ones of the Department of Physics. However, those statistical methods, because they don’t bother building causality models of the driving forces and reasons for development of the processes, are prone mistakenly take a deterministic process for a random one because of not including all necessary variables/dimensions in the model; or including the same, “real” analytical one, partially, into number of “phony” statistical variables; or/and multiplying number of such variable above necessary dimensionality.

For example, we may build our model by training it by throwing out books of the approximately same mass, and get F=Ca model formula, where C is some constant coefficient. Then, lacking dimension m, we’ll take a deterministic process for a random one, because, having the same F, we’ll get multiple a due to the mass difference of the books being thrown out.

Or, by some caprice of the statistical mind, we may come up with the following model formula: F – m/2*i + 2a = ma – m/2*i + 2a; where Fm/2, c = 2a , e = – m/2*i + 2a, md, where i is some unit vector of our basis.

b + c  = d + e

with the obviously excessive and dependent variables/dimensions. Principal Component analysis, or, actually, covariance analysis is meant to address exactly last case…

Appendix 10.1

ks_eigen_rotate <- function(df, std=FALSE){

ei <- eigen(cov(df))
#print(ei\$values)

ds <- as.data.frame(as.matrix(df) %*% ei\$vectors)
colnames(ds) <- matrix_symvect_mult(t(ei\$vectors), names(df))

if(std){
dim <- colnames(ds)
n_ds <-lapply(dim, std_norm_ds, ds)
names(n_ds) <- dim
ds <- as.data.frame(n_ds)
}

ds
}

## Data Science Done Right (the Kitchen Style) #9

Picking up loose ends…

When previously we were talking about the transformation matrix N-1 – the matrix of changing basis from the original to a basis of the element of the quotient space we project our dataset to – we picked the simplest model just to keep computational ball rolling, planning to come back and address a more reasonable choice later. Before starting playing with word vector regressions let’s improve our transformation to a bit more useful form.

Our N-1 matrix was defined as:

(a 0 .. b)
(0 c .. d)
(.. .. .. ..)
(0 0 .. 0)

and we decided to go with coefficients a=b=c=d=..=1  for simplicity.

Let’s see what that transformation that matrix would define. And let’s do it in 3D space for easier understanding, then, when we come up with a better matrix coefficients, generalizing it to a higher dimension case again).

Those matrices that convert bases define how unit vectors of these bases relate (expressed via linear combination) to each other. Therefore the following matrix:
N-1         x                  y
(1 0 1)   (1)     (1*1 + 1*1u
(0 1 1)  (0) j =  (          1*1v
(0 0 0)  (1) k

basically says that if we take a unit vector’s x=1*i+1*j+1*k projection in a ik plane, then the length (in respect to base vector u (the same will be with v)) of the transformed (rotated to) vector y will be additively combined from the lengths of vector x in respect to base vectors i and k (k is the dimension we are collapsing). And the norm of the resulting vector in respect to u will be ||y||^2 = uT*u = 2*2 = 4, or ||y|| = 2, therefore the unit vector in respect to basis u will be as twice as smaller than the transformation of the vector that has unit components in respect to the basis ik. Or, ||x||^2= xT*x = 1*1 + 1*1 = 2, or||x|| = sqrt(2), or ||y||/||x|| = sqrt(2), which would lead to a “dimension inflation”, or shrinking unit vector length’s of each synthetic dimension (or growing coefficients of the synthetic bases) after each regression step (see Fig. 9.1). Which is obviously may be inconvenient if we want to compare distributions between regression steps.

We could have set a,b,c,d,etc… coefficients in N-1 matrix to something else than 1, for example keeping a=b and c=d, but setting a!=c, also counting in account that transformation R (regression’s projection) will shrink unit coefficient of k to t31 (in this context we are looking at the transformation of the unit vector’s x projection in ik plane), and asking for ||y||=1:

N-1            x                  y
(a 0 a)      (1)     (a*1 + a*t31u
(0 c c)      (0) j =   (          c*t31v
(0 0 0)  (t31) k

then: ||y||^2=a^2*(1+t31)^2=1; a=sqrt(1/(1+t31)^2)

However, we still assign the same weight a to components of i and k when composing length of u (even we take only projected part t31 of the original unit vector k), which, if we don’t have particular reason for doing that, will distort input of the original dimensions i and k into the synthetic u. Therefore it makes sense to use projections of the unit vector components in respect to i and k on our projected dimension u to compose length of the synthetic unit vector:

N-1                          x                  y
(d*a’ 0 d*b’*t31) (1)     (d*a’*1 + d*b’*t31u
(                          ) (0) j =
(                          ) (1) k

Where: ||y||^2=d^2*(a’+b’*t31)^2=1;
a’=1*cos(alpha); b’*t31=sin(alpha)*t31;
sin(alpha)=t31/sqrt(1+t31^2); cos(alpha)=1/sqrt(1+t31^2);
then: d=1/sqrt(1+t31^2); d*a’=1/(1+t31^2); d*b’*t31=t31^2/(1+t31^2)

Now, with new N-1 matrix, we’ll get consistent, normalized range distributions of the synthetic variables on each regression step, and multi-modality of the synthetic eigenbasis distribution is even more pronounced:

Fig. 9.3. Histogram/pdf of diabetes dataset in normalized synthetic 1D eigenbasis

with the synthetic eigenbasis vector v1 = -4.46389046671904*bp.1s – 0.000335763837256707*hip + 0.0197350621146874*waist + 0.0266855254863463*bp.1d + 0.00308836195333824*height + 0.866182912400952*age + 0.00535941186779068*stab.glu – 0.00021277667540017*hdl + 0.000649570489160293*chol + 0.00146620009377012*glyhb

where coefficients less than 1e-04 were rounded down to 0 (for weight and ratio).

Fig. 9.4. Histogram/pdf of diabetes dataset in normalized synthetic 1D eigenbasis for groups with glyhb < 7 (green), and glyhb > 7 (red)

That is our current model (though without rotating synthetic basis to get rid of the variable dependency), where: glyhb – Glycosylated Hemoglobin, chol – Total Cholesterol, hdl – High Density Lipoprotein, ratio – Cholesterol/HDL ratio, stab.glu – Stabilized Glucose, bp.1s – Systolic Blood pressure, bp.1d – Diastolic Blood pressure, while other parameters’ names are self explanatory (height, hip and waist circumferences are measured in inches, weight in pounds). Standardized value of a variable is computed as sv = (v-min)/range, where:

Fig. 9.3. Histogram/pdf of diabetes dataset in normalized synthetic 1D eigenbasis for groups with bp.1s < 150 (green), and bp.1s > 150 (red)

vars        min      max  range
glyhb       2.68     16.11    13.43
chol        78.00  443.00  365.00
hdl          12.00  120.00  108.00
ratio         1.50     19.30    17.80
stab.glu  48.00  385.00  337.00
bp.1s       90.00  250.00  160.00
bp.1d      48.00  124.00    76.00
age          19.00    92.00    73.00
height     52.00    76.00    24.00
weight    99.00  325.00  226.00
waist       26.00    56.00    30.00
hip           30.00    64.00    34.00

Fig. 9.3. Histogram/pdf of diabetes dataset in normalized synthetic 1D eigenbasis for groups with chol < 250 (green), and chol > 250 (red)

V1 value greater than 0.35 seams is a threshold between first three modes of the distribution that cluster data with low parameters in the cholesterol, glucose, blood pressure groups, and last two-three modes that cluster data with high such values, which we can see at Figs. 9.4-6 with overlaying histograms for groups with low/high glyhb (>7 Fig. 9.4), bp.1s (>150, Fig. 9.5), and chol (>250, Fig. 9.6).

Appendix 9.1

# Linear univariate regression that minimizes rss
# of variable ‘yname’ of ‘df’ dataset
ks_lm <- function (df, yname, xnames, full=1, error=0){
all.names = c(unlist(xnames), unlist(yname))
clean.df <- df[complete.cases(df[, all.names]), all.names]

nr <- nrow(clean.df)
nc <- ncol(clean.df)

model <- list()
model\$X <- matrix(do.call(rbind,clean.df), nc)
rownames(model\$X) <- all.names

x <- clean.df
x[, yname] <- NULL
x[, nc] <- rep_len(1, nr)

y <- clean.df[, yname]

X <- matrix( do.call(cbind,x), nrow(x) )

b1 <- solve(t(X) %*% X) %*% t(X) %*% y

if(full | error){
e <- X %*% b1 – y
model\$e2 <- t(e) %*% e
model\$e2n <- model\$e2/nr
model\$en <- sqrt(model\$e2)/nr
}

if(full){
#Projection onto 0 member of quotient space X/Y
model\$R <- diag(rep_len(1, nc))
model\$R[nc,1:nc] <- b1
model\$R[nc,nc] <- 0

da1 <- 1/(1+as.vector(b1)**2)
model\$Nm1 <- diag(da1)
dbt1 <- 1/(1/as.vector(b1)**2+1)
model\$Nm1[,nc] <- dbt1

model\$Nm1 <- model\$Nm1[-nc,]

model\$h0 <- rep_len(0, nc)
model\$h0[nc] <- b1[ncol(b1)]

model\$X0 <- model\$R %*% model\$X
rownames(model\$X0) <- all.names
model\$Xl <- model\$X0[, 1:nr] + model\$h0
rownames(model\$Xl) <- all.names

model\$Nm1R <- model\$Nm1 %*% model\$R
model\$DimNames <- matrix_symvect_mult(model\$Nm1R, all.names)

model\$Y <- model\$Nm1 %*% model\$X0
rownames(model\$Y) <- model\$DimNames
}

model
}

sym_var_mult <- function(item, m=1){
item <- item_conv(item)

mult_res <- as.double(item[1])*m
if(abs(mult_res) < 1e-04)
mult_res <- 0

c(as.character(mult_res), item[2])
}

## Data Science Done Right (the Kitchen Style) #8

Using linear regression for dimensionality reduction

In our linear model, designed and developed (in R) in the previous chapters, we calculate not only new images of the data, being projected on the chosen element of the quotient space that satisfy particular requirements (we use minimal SSE), in the original space basis, but also in a basis of the quotient subspace. However, we will need one more function that would calculate in the symbolic form what would be new synthetic basis vectors composed by the linear combination out of the original space basis vectors. Basically we need to implement a symbolic multiplication of the matrix N-1R of our transformation into the quotient space element by symbolic vector of the variable names (see those functions in Appendix 8.1).

Having symbolic matrix by vector multiplication functions ready, we may start playing with our dataset, trying to reduce its higher dimensionality into a more intelligibly manageable and comprehensible lower dimension image. As we mentioned before, regression is a way to reduce dimensionality of a set, of course, by a price of losing complexity of the data Structure. We try to minimize that loss of the continuity of the transformation function (projection we use in regression) using some criteria (loss function) – in our case SSE. So far we have an univariate multivariable regression function ready, therefor we can collapse dimensions consecutively, one by one. With this approach we may need to decide how many dimensions we want to collapse and/or in what order.

There could be multiple approaches and criteria for that, and we may come back to this discussion later, but for now let’s start with simple ones. For obvious reasons of saving as much Structural complexity (i.e. continuity of transformation) we may want to collapse those dimensions which projections give the minimal loss function of our choice (of course, first, we normalize those dimensions to their ranges). On the other hand, it may be not so important if we go to the ultimate dimensionality reduction to the 1D space, basically making its basis vector an eigenvector.

Let’s take the Diabetes dataset we used for experimentation in the previous chapter, which has 12 dimensions (obviously not independent), and run consecutive univariate regressions until 1 dimension left:

# Dimensionality reduction
diabetes2 <- diabetes
dim <- c(“glyhb”, “chol”, “hdl”, “ratio”, “stab.glu”,
“bp.1s”, “bp.1d”, “age”, “height”, “weight”, “waist”, “hip”)
ds <- ks_lm_dim_red(diabetes2, dim, n_dim=1, reorder=1)

Now we can see what will be our synthetic eigenvector variable:

resp_name <- names(ds)[1]
resp_name
[1] “573.865492043754*chol+25.8440157742691*weight-507.345171188756*hdl+455.427566923276*hip+251.259082865673*waist-12.0580460009216*bp.1s+17.2494522189963*glyhb+142.022683378986*stab.glu+233.150806921495*height+59.1169607418513*bp.1d+29.0382173062296*age”

On the first look, not much surprises – cholesterol is really bad, HDL (“good cholesterol”) is good, hip circumference is bad, as well, while waist circumference is as twice as less bad, surprisingly, height is noticeable worse than weight, and age and blood pressure have less importance than other factors, but high values are still bad. But, anyway, that’s only a formal statistical model.

Now, let’s see how the data distribution looks like in our new synthetic 1D eigenbasis v1:

names(ds)[1] <- “v1″
ggplot(ds, aes(x=v1))+
geom_histogram(data=ds,aes(y=..density..),fill=”blue”, bins=40, alpha=0.5)+

Histogram/pdf of diabetes dataset observation in synthetic 1D eigenbasis

We see it’s noticeably multi (quatro or quinto) modal. We don’t know, yet, what that clusterization mean, and correlates with, so let’s see how it relates to the usual criterion for diabetes glycosylated hemoglobin > 7. Let’s overlay histograms of two groups of observations divided by that criterion:

clean.diabetes <- diabetes2[complete.cases(diabetes2[, dim]), dim]
ds\$glyhb <- clean.diabetes\$glyhb
ds\$glyhb_c <- 0
ds[ds\$glyhb>=7, “glyhb_c”] <- 1

ds_nd <- ds[ds\$glyhb<7,]
ds_pd <- ds[ds\$glyhb>=7,]

ggplot(ds, aes(x=v1))+
geom_histogram(data=ds, fill=”blue”, bins=40, alpha=0.6)+
geom_histogram(data=ds_nd, fill=”green”, bins=40, alpha=0.6)+
geom_histogram(data=ds_pd, fill=”red”, bins=40, alpha=0.6)

Histogram of the full dataset (blue), and datasets w/o (green), and w/ glyhb >= 7 in v1 eigenbasis

Also we may look at these two subdatasets and the whole dataset on scatterplot:

ggplot(ds, aes(x=v1, colour=glyhb))+
geom_point(data=ds, y=.5, alpha=0.5)+
geom_point(data=ds_nd, y=.75, alpha=0.5)+
geom_point(data=ds_pd, y=.25, alpha=0.5)+

Whole diabetes dataset (in the middle), and subsets w/o (above), and w/ glyhb >= 7 (below) in v1 eigenbasis

As well as descriptive statistics of these datasets:

> describe(ds[[1]])
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 377 508.01 192.12 504.76 505.21 192.78 92.88 1215.27 1122.39 0.2 -0.14 9.89

> describe(ds_nd[[1]])
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 319 478.84 177.23 482.3 476.71 182.29 92.88 1041.15 948.27 0.12 -0.31 9.92

> describe(ds_pd[[1]])
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 58 668.48 193.08 694.31 671.4 169.63 149.86 1215.27 1065.41 -0.13 0.27 25.35

Difference in the means between two subsets is just around one sigma, which has not very impressive confidence level of statistical significance (like 60-70%), but that’s not exactly the point. Glycosylated hemoglobin and it’s level 7 itself is not the ultimate and definitive diabetes type II criterion, but just one useful estimation. Diabetes type II is not just about “blood sugar”, it’s a general metabolic condition, not necessarily disorder per se. Some researchers hypothesise it’s an evolutionary adaptation (to cold climate), and when it was developed the “inconveniences” it causes in the old age and “normal”, comfortable  conditions of contemporary life was not an evolutionary priority.

Anyway, that clustering of observations with the positive glyhb criterion at last two modes (v1 > 650) may be used as an indicator of such adaptation, or likelihood of diabetes type II condition. Let’s see how this clustering look in the original (cholesterol-blood pressure-glycosylated hemoglobin) space:

clean.diabetes\$syn <- ds[[1]]
clean.diabetes\$syn_c <- 2
clean.diabetes[clean.diabetes\$syn<320, “syn_c”] <- 1
clean.diabetes[clean.diabetes\$syn>=650 , “syn_c”] <- 3
cloud(chol ~ glyhb * bp.1s, clean.diabetes, groups=syn_c)

Observations grouped by modes on histograms above: blue – first mode, magenta – next two modes, blac – the last one(s)

What we can see on this scatterplot is that the last mode clusters around observations with largest variabilities of at least one of the displayed parameters, which is a well known risk factor for diabetes. Therefore, indeed, it may make sense to investigate more the above described criterion in association with diabetes type II.

It’s easy to notice that out of the original 12 variables many of them are dependent, and could be grouped into 5 more or less independent backets (glucose related, cholesterol related, blood pressure, physical characteristics and age). Collapsing one of the dimensions in the group may lead to collapsing others, therefore it may be enough to do much less projections than the number of dependent dimensions to convert synthetic dimensions into the eigenbasis. For our particular case 5 projections (and respective basis conversions) are enough to put the data in a relatively straight line which will be the eigenbasis vector we are looking for.

From the computational point of view we are looking for a state when median sigma of the coordinates on all the remaining dimensions are less than a particular value, say, 5 %:

ds <- ks_lm_dim_red(diabetes2, dim, sd_dim=0.05, reorder=1)

or setting that we want 7 dimensions left:

ds <- ks_lm_dim_red(diabetes2, dim, n_dim=7, reorder=1)

And plotting the transformed data in the first 3 of the synthetic dimensions (in the others it will be the same), we’ll get a pretty straight line:

Our Diabetes dataset in the synthetic basis one step away from the eigenbasis transform

names(ds)[1:3] <- c (“v1”, “v2”, “v3”)
clean.diabetes <- diabetes2[complete.cases(diabetes2[, dim]), dim]
ds\$glyhb <- clean.diabetes\$glyhb
ds\$glyhb_c <- 0
ds[ds\$glyhb>=7, “glyhb_c”] <- 1
cloud(v1 ~ v2 * v3, ds, groups=glyhb_c, pretty=TRUE)

To roughly converting our dataset into 1D eigenbasis, and preserving the basis vector length (which doesn’t really matter for the eigenbasis), there would be enough to multiply one of the coordinates by sqrt(n) (where n is number of remaining dimensions (assuming the vector length preservation in two bases: A = bu = ai+aj+ak+…+an, AT*A = b*2 = n(a*2)).

If we stop just one step before:

ds <- ks_lm_dim_red(diabetes2, dim, n_dim=8, reorder=1)

We’ll still have some variability of data points:

Our Diabetes dataset in the synthetic basis two steps away from the eigenbasis transform

Whether we run all projections till the last dimension, or stop by the mean sigma criteria, we’ll have practically the same distribution along the eigen basis vector. In our case, if we stop earlier, numerical values of the synthetic variable will be different (we’ve been lazy in previous chapters defining the N-1 transform as a diagonal and last column unit matrix, therefore after each projection the length of the basis vectors shrinks), but for eigenbasis that’s not important. Eigenbasis coordinate will be just proportionally adjusted:

v1 = 23.70*chol+1.09*weight-18.89*hdl+18.57*hip+10.66*waist-0.51*bp.1s+0.73*glyhb+5.88*stab.glu+9.42*height+2.211*bp.1d+0.917*age

Histogram/pdf of diabetes dataset observation in synthetic 1D eigenbasis with the mean sigma stop parameter

We’ll play with a more intelligent Principal Component Analysis of the synthetic basis on projection later, and for now take look at another applications of the linear regression for the eigenbasis transform, in particular such, it seems, unlikely area as Natural Language Understanding…

Appendix 8.1 – Symbolic matrix by vector multiplication

#Examples:

str1 <- “2.*a-3.1*b-0.4*c-d”
str2 <- “.01*c+35*b+3.24*e-5.*d”
# Symbolic multiplication of string ‘str1’ by number 2.5
sym_lc_mult(str1, 2.5)
# Symbolic sum of 2 strings ‘str1’ and ‘st22’

str_v <- c(“a+.5*b”,”b-2*c”,”3*c+d”)
str_u <- c(“”,””,””)
m <- matrix(c(1,2,3,4,5,6,7,8,9),3,byrow=TRUE)
# Symbolic multiplication of numeric matrix ‘m’ by string vector ‘str_v’,
# result is a string vector ‘str_u’
str_u <- matrix_symvect_mult(m, str_v)
str_u

# Symbolic multiplication of numeric matrix ‘m’ by string vector ‘str_v’,
# result is a string vector
matrix_symvect_mult <- function(m, str_v){
nr <- nrow(m)
nc <- ncol(m)
str_u <- rep_len(“”, nr)
for(i in 1:nr){
for(j in 1:nc){
}
}

str_u
}

# Symbolic sum of 2 strings
str_out = NULL
str_t1 <- gsub(“-“,”+-“, str1)
l_tup1 <- strsplit(strsplit(str_t1, “[+]”)[[1]],”[*]”)
str_t2 <- gsub(“-“,”+-“, str2)
l_tup2 <- strsplit(strsplit(str_t2, “[+]”)[[1]],”[*]”)

for(item1 in l_tup1){
item1 <- item_conv(item1)
found = FALSE

for(item2 in l_tup2){
item2 <- item_conv(item2)

# item in both strings
if(identical(item1[2], item2[2])){
found = TRUE
mult = as.double(item1[1]) + as.double(item2[1])

if(mult != 0){
if(!is.null(str_out))
str_out = paste(str_out,”+”, sep=””)

str_out = paste(str_out,
as.character(mult),
“*”, item1[2], sep=””)
}
}

}

# item in first string but not in second
if(!found){
if(as.double(item1[1]) != 0){
if(!is.null(str_out))
str_out = paste(str_out,”+”, sep=””)

str_out = paste(str_out, item1[1],
“*”, item1[2], sep=””)
}
}
}

# item in second string but not in first
for(item2 in l_tup2){
item2 <- item_conv(item2)
found = FALSE

for(item1 in l_tup1){
item1 <- item_conv(item1)

# item in both strings – not processing
if(identical(item1[2], item2[2])){
found = TRUE
}
}

if(!found){
if(as.double(item2[1]) != 0){
if(!is.null(str_out))
str_out = paste(str_out,”+”, sep=””)

str_out = paste(str_out, item2[1],
“*”, item2[2], sep=””)
}
}
}

if(is.null(str_out))
str_out = “0”

str_out2 <- gsub(“[+]-“,”-“, str_out)
str_out2
}

# Symbolic multiplication of string ‘str’ by number m
sym_lc_mult <- function(str, m=1){
str_t <- gsub(“-“,”+-“, str)
l_tup <- strsplit(strsplit(str_t, “[+]”)[[1]],”[*]”)

res <- lapply(l_tup, sym_var_mult, m)

str_out = NULL
for(item in res){
if(!identical(as.double(item[1]), 0)){
if(!is.null(str_out))
str_out = paste(str_out,”+”, sep=””)

str_out = paste(str_out, item[1], “*”, item[2], sep=””)
}
}

if(is.null(str_out))
str_out = “0”

str_out2 <- gsub(“[+]-“,”-“, str_out)
str_out2
}

item_conv <- function(item){
if(length(item)==1){
if(is.na(as.double(item[1]))){
if(identical(substr(item[1],1,1),”-“))
item = c(“-1”, sub(“-“,””,item[1]))
else
item = c(“1”, item[1])
}
else
item = c(“0”, “0”)
}
else if(length(item)==0){
item = c(“0”, “0”)
}

item
}

sym_var_mult <- function(item, m=1){
item <- item_conv(item)
c(as.character(as.double(item[1])*m), item[2])
}

Appendix 8.2 – Linear regression model

# Linear univariate regression that minimizes rss
# of variable ‘yname’ of ‘df’ dataset
ks_lm <- function (df, yname, xnames, full=1, error=0){
all.names = c(unlist(xnames), unlist(yname))
clean.df <- df[complete.cases(df[, all.names]), all.names]

nr <- nrow(clean.df)
nc <- ncol(clean.df)

model <- list()
model\$X <- matrix(do.call(rbind,clean.df), nc)
rownames(model\$X) <- all.names

x <- clean.df
x[, yname] <- NULL
x[, nc] <- rep_len(1, nr)

y <- clean.df[, yname]

X <- matrix( do.call(cbind,x), nrow(x) )

b1 <- solve(t(X) %*% X) %*% t(X) %*% y

if(full | error){
e <- X %*% b1 – y
model\$e2 <- t(e) %*% e
model\$e2n <- model\$e2/nr
model\$en <- sqrt(model\$e2)/nr
}

if(full){
#Projection onto 0 member of quotient space X/Y
model\$R <- diag(rep_len(1, nc))
model\$R[nc,1:nc] <- b1
model\$R[nc,nc] <- 0

#X to Y Basis rotation
model\$Nm1 <- diag(rep_len(1, nc))
model\$Nm1[,nc] <- rep_len(1, nc)
model\$Nm1 <- model\$Nm1[-nc,]

model\$h0 <- rep_len(0, nc)
model\$h0[nc] <- b1[ncol(b1)]

model\$X0 <- model\$R %*% model\$X
rownames(model\$X0) <- all.names
model\$Xl <- model\$X0[, 1:nr] + model\$h0
rownames(model\$Xl) <- all.names

model\$Nm1R <- model\$Nm1 %*% model\$R
model\$DimNames <- matrix_symvect_mult(model\$Nm1R, all.names)

model\$Y <- model\$Nm1 %*% model\$X0
rownames(model\$Y) <- model\$DimNames
}

model
}

Appendix 8.3 – Dimensionality reduction via linear regression

# Reduce dimensionality of dataframe ‘df’ on variables listed in ‘dim’
# until ‘n_dim” dimansions left
ks_lm_dim_red <- function(df, dim, n_dim=1, reorder=1, sd_dim=0, hb=1){
clean.df <- df[complete.cases(df[, dim]), dim]

# normalize variables by their range
n_ds <-lapply(dim, norm_ds, clean.df)
names(n_ds) <- dim
norm.df <- as.data.frame(n_ds)

#order variables by rss of their univariate regressions
i_dim <- order(v_dim[2,])
order.df <- norm.df[, i_dim]

# run univariate regressions for variables in rss order
# until n_dim dimensions left
ds <- order.df
n_max <- length(dim)-n_dim
i <- 0
while(i < n_max){

if(hb)
print(i)

resp_name <- colnames(ds)[1]
pred_names <- colnames(ds)[-1]

mod <- ks_lm(ds, resp_name, pred_names)
ds <- as.data.frame(t(mod\$Y))

i <- i+1
if(i >= n_max)
break

# check if we already have effectively eigenvector (all coordinates are the same)
if(sd_dim > 0){
m_sd <- mean(sapply(as.data.frame(t(ds)), norm_sd))
print(m_sd)
if(m_sd < sd_dim){
# convert length to 1D eigenbasis
col_name <- names(ds)[1]
ds <- as.data.frame((ds[,1] * sqrt(length(dim)-i)))
names(ds)[1] <- sym_lc_mult(col_name, sqrt(length(dim)-i))
break
}
}

#reorder synthetic variables by rss on each iteration
if(reorder){
tmp_dim <- colnames(ds)
i_dim <- order(v_dim[2,])
ds <- ds[, i_dim]
}

}

ds
}

# Create list of rss for univariate regression on each normalized variable
i_range <-range(norm.diabetes[,item])
norm.diabetes[,item] <- norm.diabetes[,item]/(i_range[2]-i_range[1])

predict_dim <- colnames(norm.diabetes)
predict_dim <- predict_dim[which(predict_dim!=item)]

mod <- ks_lm( norm.diabetes, item, predict_dim, full=0, error=1 )

}

# Normalize a variable by its range
norm_ds <- function(item, norm.diabetes){
i_range <-range(norm.diabetes[,item])

norm_col <- (norm.diabetes[,item]-i_range[1])/(i_range[2]-i_range[1])
norm_col
}

# Normalize st.dev. by mean
norm_sd <- function(item){
sd(item)/mean(item)
}

## Data Science Done Right (the Kitchen Style) #7

We just created, in the previous post, a simple univariate, multivariable regression function ks_lm0 that calculates regression slopes and intercept.

However, it would be interesting not only to find out what our quotient (or factor) space X/Y would be (which is defined by the bi (i !=0 ) regression coefficients), and what particular equivalence relation R would be of our dataset A with its image B on Yl we project it to, and what the projection function (T or R) we have to Yl and Y0 spaces respectively, but also look at what our projected dataset B becomes in the new space of the reduced dimensionality Yl (or the same in Y0, because we found out before that in the linear space d=f). Actually that is the main goal of the whole exercise we just done – we projected our data from the space with the less “convenient” structure into one with a more “convenient” (with the reduced dimensionality) structure.  BTW, that thing is a part of Factor Analysis (no surprise we are dealing with factor spaces).

What we want to find (see previous chapter) is the transform M. We just found (in ks_lm0) the above mentioned transform R (which is the diagonal matrix with one row made from bi coefficients (!= 0)), and we know that NM=R, i.e. M=N-1R. In general case, if we want to rotate uv basis relative to the ij one (BTW, which is a part of Principal Component Analysis) N-1 transform would be all rows, except one 0 row matrix, but we won’t do it now, and save it for later, and just declare that basis for our Space of new “synthetic variables” Yl would be u=ai+bk, v=cj+dk, and even a=b=c=d=1. I.e. we do not rotate new basis (yet), and do not significantly scale the unit vectors length (only do that trigonometrically). Therefore N-1:
(1 0 … 1)
(0 1 … 1)
(.  .  …  .)
(0 0 … 0)

So, let us stop for a moment at a more visually “comprehendable” 3D representation of our data, where we look only at glucosylated hemoglobin (glyhb), cholesterol/hdl (rati)o and systolic blood pressure (bp.1s) dimensions, and see how our data would look in 3D and 2D representation, after we have collapsed glyhb dimension, and created 2 synthetic variables of the linear combination of glyhb & ratio and glyhb & bp.1s. To do that we will enhance our new ks_lm function to not only calculate slopes and interception as in ks_lm0, but also transformation matrices R and N-1, and the very transformed datasets Xl and Yl.

First, we use lattice package for 3D scatter plot of the original (blue) data X in coordinates V3=glyhb, V1=ration, V2=bp.1s, and projected (magenta) data Xl  into quotient space member Yl (regression plane) and then converted back to the X space basis:

Scatter plot of X (original) and Xl (regressed) data in X space

ks_model <- ks_lm( diabetes, c(“glyhb”), c(“ratio”,”bp.1s”) )

xdf <- as.data.frame(t(ks_model\$X))
xdf[, “V4”] <- rep_len(1, nrow(ydf))

xldf <- as.data.frame(t(ks_model\$Xl))
xldf[, “V4”] <- rep_len(2, nrow(ydf))

xxldf <- rbind(xdf, xldf)

cloud(V3 ~ V1*V2, xxldf, groups=V4)

Scatter plot of X (original) and Xl (regressed) data in X space (view from another direction)

Scatter plot of X (original) and Xl (regressed) data in X space (and yet another direction)

And then, we use ggplot2 package to draw a 2D scatter plot of our new synthetic variables in the Yl space:

Projection of our data to Glyhb-bp.1s (black), and to Yl synthetic basis

ydf <- as.data.frame(t(ks_model\$Y))
ydf[, “V3”] <- rep_len(2, nrow(ydf))

xdf <- as.data.frame(t(ks_model\$X))
xdf[, “V3”] <- rep_len(1, nrow(ydf))

xydf <- rbind(ydf, xdf)

ggplot(xydf, aes(x=V1, y=V2, colour=V3))+
geom_point(alpha=0.5)

Now, having visually seen how our linear regression function work on 3D data, and how it may be used for dimensionality/factor reduction, and where in it is a place for principle component manipulation, let us proceed playing with it on a larger dimension data…

Appendix

ks_lm <- function (df, yname, xnames){
all.names = c( unlist(xnames), unlist(yname) )
clean.df <- df[complete.cases( df[, all.names] ), all.names]

nr <- nrow(clean.df)
nc <- ncol(clean.df)

model <- list()
model\$X <- matrix(do.call(rbind,clean.df), nc)

x <- clean.df
x[, yname] <- NULL
x[, nc] <- rep_len(1, nr)

y <- clean.df[, yname]

X <- matrix( do.call(cbind,x), nrow(x) )

b1 <- solve(t(X) %*% X) %*% t(X) %*% y

#Projection onto 0 member of quotient space X/Y
model\$R <- diag(rep_len(1, nc))
model\$R[nc,1:nc] <- b1
model\$R[nc,nc] <- 0

#X to Y Basis rotation
model\$Nm1 <- diag(rep_len(1, nc))
model\$Nm1[,nc] <- rep_len(1, nc)
model\$Nm1 <- model\$Nm1[-nc,]

model\$h0 <- rep_len(0, nc)
model\$h0[nc] <- b1[ncol(b1)]

model\$X0 <- model\$R %*% model\$X
model\$Xl <- model\$X0[, 1:nr] + model\$h0

model\$Y <- model\$Nm1 %*% model\$X0

model
}

## Data Science Done Right (the Kitchen Style) #6

Let’s implement univariate multivariable linear regression (see derivations in previous chapter) in R the way it is usually used, i.e. mapping our dataset A not just to an element of the Quotient space Yl, but, after that, mapping it back into original space X. Our ks_lm0 R function is presented in Appendix 1.

Let’s run it against a dataset we used in one of the Statistical Methods classes (Statistical hypotheses testing for Diabetes Type II risk factors), also, for comparison, running these data through the standard lm function, making sure that results do match:

> ks_lm0( diabetes, c(“chol”), c(“glyhb”,”ratio”,”bp.1s”,”bp.1d”) )
[,1]
[1,] 1.5868593
[2,] 11.2724797
[3,] 0.1593560
[4,] 0.3243256
[5,] 98.7703565

> cols <- c(“chol”,”glyhb”,”ratio”,”bp.1s”,”bp.1d”)
> clean2.diabetes <- diabetes[complete.cases(diabetes[,cols]),cols]
> View(clean2.diabetes)
> lmodel <- lm(“chol ~ glyhb + ratio + bp.1s + bp.1d”, clean2.diabetes)
> lmodel

Call:
lm(formula = “chol ~ glyhb + ratio + bp.1s + bp.1d”, data = clean2.diabetes)

Coefficients:
(Intercept) glyhb ratio bp.1s bp.1d
98.7704 1.5869 11.2725 0.1594 0.3243

Which, of course (no rocket science there, just simple matrix arithmetic), they do 🙂

However, let’s not stop here, and continue visualizing and playing with our data further…

Appendix 1

ks_lm0 <- function (df, yname, xnames){
all.names = c( unlist(yname), unlist(xnames) )
clean.df <- df[complete.cases( df[, all.names] ), all.names]

x <- clean.df
x[, yname] <- NULL
x[, ncol(clean.df)] <- rep_len(1, nrow(clean.df))

y <- clean.df[, yname]

#X <- matrix(unlist(x), nrow(x))
X <- matrix(do.call(cbind,x), nrow(x))

solve(t(X) %*% X) %*% t(X) %*% y
}