PCA / FA. Example 4 ! Davis, and almost everyone else

I would like to revisit the work we did in Davis (example 4). For one thing, I did a lot of calculations with that example, and despite the compare-and-contrast posts towards the end, I fear it may be difficult to sort out what I finally came to.

In addition, my notation has settled down a bit since then, and I would like to recast the work using my current notation.

The original (“raw”) data for example 4 was (p. 502, and columns are variables):

It would be prudent to compute the means…

and variances of the data…

It would also be prudent to calculate the row sums.

They are all the same (we have “constant row sums”).

Now, Davis chooses to analyze centered data. Let’s get it, by subtracting its mean from each column.

That is what he has on page 503. I remind us that we expect (HERE) that the row sums of the centered data are now zero. That in turn means that we have lost rank: the centered data should have two nonzero eigenvalues or singular values, instead of 3.

Next, Davis chooses to analyze instead of the covariance matrix . We recall that these two have the same eigenvectors, but eigenvalues differ by a factor of N-1. In fact, Davis also uses the singular value decomposition (SVD); that is a good reason for doing an eigendecomposition of .

Here is the SVD of X, . (Always remember that my u and v correspond to Davis’ V and U respectively, if you’re looking at both of us.)

Sure enough, there are only two nonzero singular values in the w matrix.

Let us now compute Davis’ loadings and scores. Using my notation for the SVD, he has

, but use the definition for computation: one less matrix product.)

(I should remark that we see that we could have done all of those calculations using the data X, the singular values w, and the matrix v: we do not need u.)

and we get

Davis does not have any of these extra columns of zeros, but all of our nonzero numbers agree (to within a sign for each column). At this point, I think the most significant thing to say is that these scores and loadings do not reproduce the data; in the sense that, the matrix product of corresponding scores S and loadings A is not X.

Let us see what Jolliffe would have done.

He writes his model as

Z = X A,

with Z as the new data, and A as an orthogonal eigenvector matrix; but I would write that as

Y = X v,

since I reserve A for the weighted eigenvector matrix.

(I have to point out that, as we saw in previous posts, .)

Now, Jolliffe would usually do an eigendecomposition of the covariance matrix; but we have enough experience to understand that we can do an eigendecomposition of . Here it is:

Here are its eigenvalues:

Here is an orthogonal eigenvector matrix (now I use my standard notation, V for an orthogonal matrix computed from an eigendecomposition, versus v from the SVD):

We can check that . (It does.)

We know that V should differ from v (from the SVD) by at most the signs of columns. Recall v from the SVD:

In fact, they agree completely.

Jolliffe would project the data X onto the eigenvector matrix V to get new data Y = XV:

Now we should recall , to confirm the equivalence:

We recall that Y = X v = X V is equivalent to the change of basis formula

x = v y,

where x is an old observation and y is the corresponding new one.

We learned – from Harman! – that the new data Y has maximally redistributed variance. It is crucial that we used an orthogonal eigenvector matrix v or V. Let’s confirm that by computing the variances of the columns of Y…

and their sum is 32.

Let’s recall the variances of the columns of X…

and their sum is also 32.

Now, what about Harman? He writes his model as

Z = A F,

where

,

(i.e. his variables are in rows, not in columns) and A is a -weighted eigenector matrix. And I generally favor his notation for Z, but I will follow Jolliffe for a few seconds later on.

Let us compute. From the eigenvalues of c,

we construct a diagonal matrix of their square roots:

We should expect that these are the nonzero nonsingular values w. Sure enough:

Note, however, that w is not square, hence is not equal to the diagonal matrix . This will lead to an irrelevant difference between A and , etc. In fact, I often ignore the distinction between and w conceptually; in practice, of course, whatever matrices we’re using have to be “the right size” for whatever matrix products we compute (conformable).

We should recall that the nonzero eigenvalues are the nonzero values of both and :

The nonzero eigenvalues were…

whereas

and

Resuming our thread, we compute the weighted eigenvector matrix A

and that should be Davis’ ; and it is, except for an extra – to my mind, irrelevant – column of zeros.

Viewing A as a transition matrix, we know that Harman’s scores are the new data . We also know (HERE) that we may compute them as

,

that is, the “first few” columns of u, not all of u.

Ah, life gets interesting because A is of rank 2 instead of 3. The new data should be the first two columns of u, but if we want to see (Harman’s) Z = A F, then we must either cut down A or keep the first 3 columns of u, to have conformable matrices. I will cut down A.

We check that by computing A F and comparing it to ; we expect Z = AF. And we have it; here’s Z:

And here’s A F (using the two nonzero columns of A!):

I remind us that this new data does not have redistributed variance; in fact, this new data has constant variance:

There are a couple of (intertwined) ways to understand why the variances are not 1. Fundamentally, we did an eigendecomposition of rather than of the covariance matrix . Alternatively, our new data was given by columns of u rather than by columns of . (And that’s a consequence of the “fundamental” reason, but it may make sense by itself.)

That constant variance is in marked contrast to the maximally redistributed variances of the new data Y constructed from the orthogonal matrix v instead of from A.

It is my strong preference to compute the new data from the first few columns of u. To make sure I have the scaling correct, I then check the product A F.

Any number of people, including Davis himself in a subsequent example, show us how to compute the new data using a reciprocal basis instead; that is, one which weights the v matrix by instead of . (Whether or not an author uses the phrase reciprocal basis or dual basis, if he divides each column of v by , he is computing a reciprocal basis.)

Life is made interesting, but not difficult, by the presence of a zero eigenvalue. Frankly, the simplest way to deal with a zero eigenvalue – for finding a reciprocal basis – it is to set the zero eigenvalue to 1. Instead of the diagonal matrix

I use

and instead of , I compute the reciprocal basis (i.e. I divide each column by instead of multiplying by ).

Then the new data with respect to A can be found by projecting X onto B:

That had better be … and it is:

(Okay, that was the simplest way to get a reciprocal basis when we have a zero eigenvalue; but the simplest way to get the new data is to avoid the reciprocal basis entirely by taking the first few columns of u instead – assuming that you have the SVD available to you.)

For this example, Davis did not compute the new data with respect to A (i.e. wrt his ). As I have said before, I have seen him do it in another example; but his scores are not the new data , which is, as far I have seen in textbooks, what other people use for the scores corresponding to loadings A.

Moving on, we learned from Basilevsky (HERE) that the weighted eigenvector matrix was also the dispersion matrix between the new data and the old data X. I say “dispersion matrix” because A itself was not computed from a covariance matrix but from a more general “dispersion” matrix .

Let’s confirm that. For just a moment I want to denote the new data by (Jolliffe’s) Z instead of , so as not to get confused by transposes. A quick computation (on the side) shows me that I want to compare A with , i.e. that I want to compute

(finally reverting to my original notation). Fine, I compute that product, and I get

and I compare that to A:

Good. Our matrix A is, as we expect, a cross-dispersion matrix between old and new data. (Perhaps I should be less meticulous, and refer to A as a cross-covariance matrix; that term may be more familiar and hence clearer in sense, even though it is not the covariance matrix.)

Now let me show you something new. We had also discovered from Basilevsky (HERE) that the orthogonal eigenvector matrix v was not a dispersion matrix D between its old data X and its new data Y. The v matrix does not have this property of the A matrix. We had seen that v and D differed by eigenvalue column-weights:

.

But… we’ve just seen that weight here. Okay, it wasn’t , but it was either or . Well, it must have been , i.e. the one that multiplied v:

.

Of course, I should check that computationally. Here is :

and here is our :

Hang on. We need to explicitly see , the cross-dispersion matrix analagous to : we compute it…

I don’t know if it’s useful, but there it is anyway. I can hear my favorite Vulcan whisper, “Fascinating.”

Don’t misunderstand me: I am sure that the Q-mode scores and loadings are interesting in their own right. (I don’t understand the goal of a Q-mode analysis, but that’s another question.)

But I do find it… fascinating… that Davis’ Q-mode scores are the cross-dispersion matrix between the old data X and the new data Y (in addition to our earlier finding that his Q-mode loadings are the new data Y with respect to the transition matrix v).

To put that another way…. The orthogonal eigenvector matrix v could have been computed from an eigendecomposition of the dispersion matrix , and defines new data Y; and is the cross-dispersion between old data X and new data Y.

I had asked a question: is the orthogonal matrix v the dispersion matrix D between the old data X and the new data Y? The answer was, “No” (HERE).

Now we know the rest of the story: not v, but the Q-mode scores are that dispersion matrix D.

PCA / FA example 4: davis. Davis & harman 3.

Hey, since we have the reciprocal basis, we can project onto it to get the components wrt the basis. After all that work to show that the are the components wrt the reciprocal basis, we ought to find the components wrt the basis. And now we know that that’s just a projection onto the reciprocal basis: we want to compute X B.

Recall X:

Recall B:

The product is:

What are the column variances?

Does that surprise you? the 1/3 comes from our using eigenvalues of instead of eigenvalues of the covariance matrix. The new variables have a common variance, but it isn’t 1. this corresponds to our initial work in harman, when we discovered that using the weighted eigenvector matrix gave us new variables with variance 1.

at first glance, it seems not worthwhile to work that out using the covariance matrix, but i suppose it might be worthwhile to come at the covariance matrix from the SVD. (oh, yes! it is!)

Let’s do it.

The covariance matrix of X is

We should just replace X by and find the SVD of X2. but there’s a catch.

Do it anyway. here’s our X2:

Compute the variances…

and their sum … 10.6667 .

Hang on. The data X2 has 1/3 the variance of the original data. We’re trying to find new forms of X (not of X2) with variance 1. but let’s just keep going.

Get the SVD of X2:

here are new u…

new v…

and new w…

Then i compute new A’s and S’s from

.

.

.

.

(using X2, u2, v2, and w2, of course)

are the components of the data X2 wrt the orthogonal eigenvector matrix v2.

Compute the variances…

and their sum is 10.6667 . Yes, the new orthogonal eigenvector matrix v has redistributed the total variance of X2 among two new variables.

What about the R-mode scores?

Compute the variances:

i say again, yikes!

For the record, here are the new

and

From the orthonormal basis v2 , i construct a reciprocal basis for the -weighted matrix . i’ll just use the diagonal of w2 (instead of the since i haven’t computed them) to scale two of the vectors. Recall the new w matrix:

i construct a diagonal matrix using 1/w when possible, 1 otherwise.

then i premultiply that by the new v, and keep only the first two columns; the new reciprocal basis is:

Having the reciprocal basis, we project the data onto it get the new data wrt the basis. But we project the original data. That’s the catch. As i said, it’s the variance of X we’re trying to scale to 1, not the variance of X2. we compute X B2:

And the variances?

.

Indeed. this corresponds even more closely to what we first saw in harman. (if instead we project X2 onto B2, we get a common variance of 1/3 again.)

PCA / FA example 4: davis. Davis & harman 2.

What about redistributing the variance?

We learned from harman that if we use the orthogonal eigenvector matrix v to compute new data, the result has redistributed the original variance according to the eigenvalues; and we also learned that if we use the weighted eigenvector matrix instead, we get new variables of unit variance.

(i think i thought this required us to start with the correlation matrix; it does not.)

let’s see this again for davis’ centered data X:

We have variances…

or…

and the sum of the variances is…

32.

Now is a good time to recall :

The sum of the diagonal elements is 96, and if we divide by 3 (= N-1, for the covariance matrix instead of for ), gee, we get 32.

stop for a moment. We have two, admittedly proportional, sets of eigenvalues, one set from the covariance matrix, the other from . this gives us two different possible sets of weights to apply to the eigenvector matrix v. (recall that v is an eigenvector matrix both for and for the covariance matrix of X, as well as coming out of the SVD of X.)

Let’s stay with what i know: our new data, using the orthogonal eigenvector matrix v, is .

The variances are…

and we don’t even need to consciously add those: the sum is 32. We have redistributed the original variance among our two new variables.

Now, everyone seems to say that the -weighted eigenvectors have variances equal to the eigenvalues. What those eigenvectors really have is lengths equal to their .

here indeed be dragons.

as linear algebra is my touchstone in general for all of this, that redistributed variance using the orthogonal eigenvector matrix is my vorpal blade (dragon, jabberwock, no big deal – snicker-snack either way). any different choice of basis must lead to different results.

we looked at that redistribution of variance in harman. We learned then that the correct way to get new data with redistributed variance is to use the orthogonal eigenvector matrix. We just reconfirmed it for davis’ example.

Instead of actually kicking one of the dragons awake, let’s try counting its scales while it sleeps. If the R-mode scores are some form of data (and they are), what are their variances?

my is

the variances of the first two columns are

yikes! if nothing else, that should convince you that the R-mode scores are weird.Life is so much simpler if we use the orthogonal eigenvector matrix. The new data is trivial to compute (X v), and we can see that it has redistributed the variance of the original data.Trying to use any weighted eigenvector matrix is dicey if there is a zero eigenvalue (because our transition matrix is not invertible); and confusing if we used instead of the covariance matrix (because they have different eigenvalues).i will eventually calculate the data wrt both of the weighted eigenvector matrices. The dragons will be awake and aloft and pretty to look at.

PCA / FA example 4: davis. Davis & harman 1.

Let us now recall what we saw harman do. Does he have anything to add to what we saw davis do? if we’re building a toolbox for PCA / FA, we need to see if harman gave us any additional tools.

here’s what harman did or got us to do:

• a plot of old variables in terms of the first two new ones
• redistribute the variance of the original variables
• his model, written Z = A F, was (old data, variables in rows) = (some eigenvector matrix) x (new data)

We should remind ourselves, however, that harman used an eigendecomposition of the correlation matrix; and the eigenvectors of the correlation matrix are not linearly related to the eigenvectors of the covariance matrix (or, equivalently, of .) having said that, i’m going to go ahead and apply harman’s tools to the davis results. i will work from the correlation matrix when we get to jolliffe.

What exactly was that graph of harman’s? he took the first two columns of the weighted eigenvector matrix, then plotted each row as a point. He had 5 original variables. His plot

shows that the old 1st and 3rd variables, for example, are defined very similarly in terms of the first two new variables. let me emphasize that this is a description of defining variables, not a display of the resulting data. Simple counting can help with the distinction, e.g. 3 variables but 4 observations.

from davis we only have 3 variables to start, and our weighted eigenvector matrix has only two nonzero columns; and it’s one form of what we called . we have to be careful, because there are sign differences among the various forms of ; not because the definitions are different, but because eigenvectors are not unique. In particular, my u and v from the SVD have some different signs from davis’ U and V from eigendecompositions.

most of the time, i will use my SVD results, but right now let’s compare davis with harman. Davis had as…

Davis showed us a plot of these numbers viewed as three points.

So the R-mode loadings tell us how the 3 old variables are described by the new ones; we see that a plot of corresponds to harman’s plot.

i want to emphasize that even the orthogonal eigenvector matrix shows us that the data is really only 2D. Don’t think that we had to scale the eigenvectors in order to get 2D.

my orthogonal eigenvector matrix v is

but wrt that basis, the new data is :

i.e. the new data is described by the first two eigenvectors, hence it lies in a plane. but then so does the original data, because the transformation between them is linear (a plane cannot be mapped to 3D, although it could be mapped to a line). we could plot the data in 3D and rotate the point of view until we see the plane edge on.

Do davis and harman’s models resemble each other?

Let’s recall the SVD form of the A’s and S’s.

.

.

.

.

Then

.

in particular,

and then .

is harman’s model (Z = A F, where , A is a weighted eigenvector matrix, F is the new data), while is jolliffe’s (Z = X A, where Z is the new data, X is the original data, and A is the orthogonal eigenvector matrix). i think it would be unfair to say that davis is using either of those models, especially since harman and jolliffe (usually) work from the correlation matrix.

i would summarize davis with all four A’s and S’s.

i can’t, however, say too many times that and can be constructed using u or v:

,

and

.

i see no grounds for associating with just u or just v. and no grounds for associating either u or v primarily with R-mode or Q-mode; they’re equally important to the SVD.

i would probably focus on v and : one defines new variables, the other gives me the associated new data.

PCA / FA example 1: harman. discussion 2

recall that we had discovered that harman’s new F variables were uncorrelated and standardized: they had mean 0 and variance 1. he had implied, however, that the variances of the new variables would be the eigenvalues of the correlation matrix of Z.

 

he didn’t get that. can we?

 

instead of his model

 

(where Z is the standardized data, A is a weighted eigenvector matrix, and F is the transformed, new, data), which implies 

, 

we take

,

(where P is an orthogonal eigenvector matrix, and Y is the transformed, new, data), i.e.

hence

.

(in words, we fix his messy by using an orthogonal eigenvector matrix P, for which .) here’s the transpose of the Y matrix:

and the variances of the Y variables?

ta da! those are so close to the eigenvalues…

that we all are quite sure they are equal in principle. right?

 

for the record, here is the covariance matrix of Y:

like the Fs, the Ys are uncorrelated; unlike the Fs, the Ys have different variances, and in fact, the first Y variable, , has the largest possible variance among linear combinations of the Xs, subject to the constraint that the total variance of the five variables is constant (= 5.) and so on: it’s no accident that has the second largest variance.

we have somethig new here, but it does not come from the F variables in harman’s model. using an orthogonal eigenvector matrix as a change-of-basis, we have redistributed the variances of the Z data. harman did not accomplish this with the weighted eigenvector matrix.

it gets worse. if we were to recreate that wonderful picture using the first two coumns of the orthogonal eigenvector matrix -P… 

(yes, it is convenient to change the sign of P; we did it to A, too.) we plot those pairs of numbers and get the left-hand image; the right-hand image is the previous plot…

     

we see exactly the same implications (2 & 5, 1 & 3 together, 4 alone), but the scales are different. maybe we should have stayed with P all along. 

next? let’s see how jolliffe does PCA in “principal component analysis”.

PCA / FA Example 1: Harman. discussion 1.

harman gave us two outputs as a result of his analysis: a table and a picture.

let’s consider the picture first. it clearly shows z2 and z5 similar, z1 and z3 similar, and z4 roughly in the middle between them.

 

i don’t know about you, but this surprised me. should it have? well, let’s take a look at the correlation matrix, from which we got our results. i’m going to round it off just a little bit, so i can refer to 3-digit numbers instead of 5.

we confirm that there is a high correlation (.972) between z1 and z3, and and one almost as large (.863) between z2 and z5. z4 is about equally related to z2 (.691) and to z5 (.778). the same information as the picture seems to be there, but i didn’t see it.

let’s learn from that blindness. how might we have expected it? this time, round off the correlation matrix severely (to the nearest .5).

there we are: on the off-diagonal, we have 1s for z1 & z3, z2 & z5 (good), and also for z4 & z5 (not so good). at the other extreme, we have 0s for z1 & z2 and z1 & z5, so z2 & z5 are far from z1. finally, we have either .5 or 1 for z4 vs everything else. 

there was nothing sacred about .5. if we round to the nearest .25, instead, we get … 

we see similar relationships, but now the only off-diagonal 1 is between z1 & z3; z2 seems to be equally related to z4 & z5. but z4 is still the only variable related to everything else. bear in mind that we – i at least – still think the third variable, P3 (cf. F3) might be important.

conclusion: the picture we got showed us some strong relationships that are not quite so clearly seen by playing with the correlation matrix. 

i find it reassuring that the relationships can be somewhat seen in the correlation matrix; i begin to trust that the relationships were not produced by our computations.

conclusion: severe rounding can be informative.

before we leave the picture for his table, let’s go to 3D. we take the first three columns of the weighted eigenvector matrix A…

… and we plot them as points in 3D:

 

 

that’s disappointing: z1 & z3 are still close, but z2 & z5 are not. ok, we have something to look at down the road. what we’re nibbling on the edges of is reduction of dimensionality: can we replace 5 variables by fewer? we are not done with this, not by a long shot.

open: investigate reduction of dimensionality.

now let’s go look at his table. recall it:

there are a few things which i find confusing at best, misleading at worst. (in fact, i think the table is inconsistent.) 

i do not count his use of instead of , although you might. maybe i should not have used his notation. he used instead of simply because his is a principal components example, but he will later do a factor analysis of the same data. eventually, he will obtain some other coefficients which he will multiply by . 

this is just like my distinquishing the orthogonal eigenvector matrix P from the weighted eigenvector matrix A. you haven’t seen me use P since we got A, but i want both available; we won’t see him compute the other coefficients, but he wanted both labels available. i do plan to discuss his factor analysis; unfortunately, the example i need to use is not based on this data.

the first thing i object to is the column headings. whether or , they are misleading. the column headed P1 contains the first column of the weighted eigenvector matrix A. P1 (cf. F1) should denote the first column of the F matrix. instead, he is using the column heading as a reminder that in the equation…

we will multiply, for example, the 3rd column number .319 by P3. numerically, that just says that the second row of Z is the matrix product of the second row of A with (each column F1, F2, … of) the F matrix. his headings really stand for the columns of F. 

his fundamental model is

and for the second row of Z, that equation would be written with instead of :

as i said, maybe i made a mstake in switching to his notation; i’ll keep it in mind for the future. and no, he did not write that model as Z = A P; although he changed to , he did not change the matrix F to P. good thing, or i’d have had to use a symbol different from P for the orthogonal eigenvector matrix.

the second thing i object to is both of his “variance” labels. in fact, the numbers in the second-to-last row are the squared lengths of the columns of the weighted eigenvector matrix A. and the numbers in the last column are the squared lengths of the rows of the weighted eigenvector matrix A. 

i shouldn’t be hasty. maybe those squared lengths are the variances. 

ok, let’s compute F. hmm. we need Z. harman never computed Z. as i said, what he wanted was that drawing showing the relationship between the Z variables and the F varables.

ok, let’s figure out what Z might be. we can certainly construct a Z matrix with variances equal to 1, which is what he says they are. 

(the only reason i found this confusing is that for his theoretical work he has taken Z to be “small standard deviates” – where you divide each datum by – instead of “standardized”.) 

in order for the variances of the Z variables to be 1, all we need to do is standardize the data. we recall our data matrix D…

we standardize D and call it X:

i have changed the name; this let’s me distinquish the data matrix D from a design matrix X. further, we will let 

that may seem like a silly step to you, but i have enough trouble sorting out what’s going on without having to translate things in my head: i really want both X and Z. oh, why do i want both? because Z in harman’s model has variables in rows, not in columns. back to our model:

(BTW, A inverse exists because none of the eigenvalues were zero. A is an eigenvector matrix weighted by the (square roots of the) eigenvalues. it came from an orthogonal matrix P of unit eigenvectors. P is trivial to invert; just take the transpose: A will be invertible so long as the weights are nonzero, because the corresponding column of A will be zero if and only if a weight is zero.)

next, bizarre as you may think me, i’m going to take the transpose of the last equation. and i’m going to abuse notation (so i’ve been told), writing for the inverse transpose of A…

but that’s

why am computing , the transpose of F? for one thing, both Z and F have 12 columns and don’t fit in the page! for another, Mathematica wants variables in columns: its Mean, Variance, Standardize, Covariance, and Correlation commands each want , not F. (did you think i was computing the correlation matrix step-by-step?)

here is

ok, like the Z variables, the F variables have zero means:

that is to say,

and, like the Z variables, the F variables have unit variances:

oops. this is not ok. the table is wrong: harman said the variances were the eigenvalues. (and if you happen to know, vaguely or precisely, that PCA was supposed give us new variabes that redistribute the variance, you were really, really expecting to see the eigenvalues. welcome to the club.) 

what’s going on? it’s time to do a little algebra. 

i need one key relationship. if we have a matrix M with variables in columns, with N rows, and with each variable having zero mean (called centered data), then the covariance matrix of M is given by

further, if the variance of each variable is 1, then the covariance matrix c is also the correlation matrix r.

let’s lay out all the equations we have. 

one, X contains columns of standardized variables, so the correlation matrix r (= the covariance matrix) of X (or of Z) is given by

two, the covariance (not necessarily correlation) matrix of F is given by

(i’m going to take it for granted that the means of the F values are zero; after all, the Fs a linear combination of Xs, which do have zero mean.)

three, the eigendecomposition was

with P orthogonal and diagonal.

four, the weighted eigenvector matrix A is

where denotes a diagonal matrix whose entries are square roots of the entries of . (all of the eigenvalues are non-negative in general; and in particular, my changing all the signs of columns of A amounts to choosing the negative square roots! eigenvectors are not unique.)

perhaps i should remind you that multiplying diagonal matrices is almost trivial: we multiply corresponding diagonal elements. that means, for example, that and have diagonal elements which are the inverses of each other. that’s also why it makes sense to take square roots of the elements of and call the result .

finally, five, the fundamental model is

assuming that A is invertible (i.e. assuming all eigenvalues of r were positive, not just non-negative).

let’s compute the covariance (not correlation) matrix c of F. we have

.

first we substitute and :

;

then we substitute :

;

 

then we substitute :

whoa! an identity matrix? boy, do we need to compute the covariance matrix of F for this data!

here it is, truncated so i don’t have to reformat a bunch of numbers of the form a x 10^-b:

an identity matrix, conceptually if not numerically. yes, the algebra was right. not only are the F variables of unit variance, they are completely uncorrelated with each other. yessssss. now i remember seeing that statement somewhere.

conclusion: the F variables are uncorrelated variables of unit variance.

conclusion: the eigenvalues are not the variances of the F variables. not for Z = A F, with A a weighted eigenvector matrix. tsk, tsk.

conclusion: harman’s table can be fixed most simply by changing the first-column label “Variance” to “Eigenvalue”.

nevertheless, his “variance” label suggests that we should be able to construct combinations of the Zs which have the eigenvalues for variances.

open: how do we find new variables whose variance is the eigenvalues? 

but we’ve done enough for one post.

ok, ok, i won’t quite leave you hanging. what about using the orthogonal eigenvector matrix P instead of the weighted eigenvector matrix A?

PCA / FA Example 1: Harman. what he did.

Less is more. And more is huge. It is easy for me to end up with huge posts to put out here, but I’d rather go with smaller.

Let’s get started with PCA / FA, principal components analysis and factor analysis.

In case it matters, I am using Mathematica to do these computations.

Here is an example, the first of several. This comes from Harman’s “factor analysis”. In order to emphasize the distinction between PCA and FA, he has one example of principal component analysis, and this is it.

Let me tell you up front what he did:

• Get some data;
• Compute its correlation matrix;
• Find the eigenstructure of the correlation matrix;
• Weight each eigenvector by the square root of its eigenvalue;
• Tabulate the results;
• Plot the original variables in the space of the two largest principal components.

I also need to say that his conceptual model is written

Z = A F,

And from the dimensions of the matrices it is clear that

A is square, k by k

Z and F are the same shape, with k rows.

We infer from its size that A will be derived from the eigenvector matrix, and that Z is derived from the given data matrix. From the shapes, we conclude that Z has observations in columns, rather than in rows. (If you’re used to econometrics or regression, you expect the transpose, observations in rows.)

But this is a fine thing, because we recognize that Z = A F is a change-of-basis equation for corresponding columns of Z and F; A is a transition matrix mapping new components (any one column of F) to old components (a column of Z).

The following data comes from Harman, p. 14. We would customarily say it has 5 variables (k = 5) and 12 observations (n=12); to be precise, we mean it has 12 observations per variable, since the total number of data points is 60. I have chosen to use regression notation, and denote the number of variables by k.

Here’s the data.

I have called it D, for data matrix; and I have displayed it with observations in rows. This is what I’m used to from regression analysis. This is also how Harman displayed it. The point of the discussion about Z is that Z is a transposed matrix, relative to D: D has variables in columns, Z has variables in rows. In addition, if Harman were to to compute Z – which he does not – it would almost certainly be standardized data. (That is, subtract the mean of each variable, and then divide each by its sample standard deviation.)

He printed means and standard deviations. I checked his printed means…

and, more importanty, his standard deviations…

Our numbers agree: despite using N in his discussion, he correctly used N-1 in the computation of the sample variance. I would have been stunned to disagree with his means, but he led me to expect different computed standard deviations. I am pleasantly surprised.

We compute the correlation matrix r:

 

 

This also agrees with his tabulated correlation matrix on p. 14. (I have rounded the display to match his.)

Now (jumping from his p. 14 to p. 135) we get the eigenstructure of the correlation matrix r. Here is an orthogonal eigenvector matrix:

(I remind you that eigenvectors – even orthonormal ones – are not unique; as it happens, every one I got is the negative of his. that doesn’t matter, but for subsequent computations, I changed the sign of my matrix.)

(We have “diagonalized” the correlation matrix r. That is, we have computed a matrix P whose columns are unit eigenvectors, and  a diagonal matrix whose diagonal elements  are the eigenvalues, such that .)

We also agree on the eigenvalues (which are unique):

 

Now he weights each eigenvector by the square root of its eigenvalue; e.g., the first eigenvector will have length while the fifth will have length .

We end up with the following (rounded) matrix of weighted eigenvectors (where I have also multiplied all of mine by -1):

He presents the following summary table, whose center is precisely our weighted eigenvectors:

I’ll have a lot to say about that table, but not yet. Oh, I should point out that the row labeled “variance” consists of the five eigenvalues and their sum (yes, 5.).

There’s one more thing he did. Back on p. 16, as a preview of things to come, he wrote

(where he used rather than to emphasize that these were principal components from PCA rather than factors from FA; you should read them as , but I don’t want to misrepresent exactly what he wrote).

But that is precisely the equation Z = A F written for the second row of the Z matrix: multiply the second row of A by each of the columns of F. I am happy to see this. It explains why people never display Z, and never compute F: they just want to describe the old variable names in terms of the new variable names.

Let us select all the coefficients of P1 and P2, not because only two new variables are important, but because it’s easy to plot things in 2D. By weighting the eigenvectors by , he has emphasized the first eigenvector over the second, and the second over the third, etc. So let’s see what the first two tell us.

Here are the first two columns of the weighted eigenvector matrix A, i.e. the first two weighted eigenvectors. 

 

Let’s plot those pairs of points. Recall

 

 

but retain only the first two terms: 

 

 

Now plot the point (.767, -.545) and label it “2″. Do the same for the other four equations for z1 and z3 thru z5. (Please don’t be offended at my being explicit: it took me forever to get my mind off the data itself and to grok this relationship between variable names.)

 

Verily, points 1 and 3 (i.e. variables z1 and z3) are similar, and variables z2 and z5 are even more similar, and variable z4 is about equally different from the two clusters, but closer to 2-5.

Let me summarize Harman’s analysis to this point. For this example, he has:

• given us data: 5 variables, 12 observations;
• computed the eigenstructure of the correlation matrix;
• Weighted each eigenvector by the square root of its eigenvalue;
• Tabulated the results;
• Plotted the original variables in the space of the first two principal components.

What has he not done?

• He did not actually compute Z or F;
• He did not explain why his table looks the way it does;
• He did not explain the graph he drew;
• I wonder if the fact that z4 is different from the two clusters suggests that there really ought to be three new variables P1, P2, P3 instead of just two;
• Maybe we should do a 3D graph.

Enough for now. Next, I’ll talk about what he did.