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 Bartholomew et al.: discussion

This is the 4th post in the Bartholomew et al. sequence in PCA/FA, but it’s an overview of what I did last time. Before we plunge ahead with another set of computations, let me talk about things.

I want to elaborate on the previous post. We discussed

• the choice of data corresponding to an eigendecomposition of the correlation matrix
• the pesky that shows up when we relate the of a correlation matrix to the principal values w of the standardized data
• the computation of scores as
• the computation of scores as where
• the computation of scores as projections of the data onto the reciprocal basis
• different factorings of the data matrix as scores times loadings

The three computations of scores, of course, are all the same, only looking different. I will show you an example where A is not invertible (but not in this post). Although it looks harder at first, it simplifies so considerably as to be easier in practice.

The choice of data is actually dictated by the SVD. If we do an eigendecomposition of the correlation matrix c,

then the data must be standardized or a constant multiple of the standardized data (e.g. small standard deviates), because that gives us that V from the eigendecomposition matches v from the SVD (to within signs of columns),

.

Similarly, if we do eigendecomposition of the covariance matrix, then the data must be centered or a constant multiple of the centered data. (We could define “small centered deviates” by dividing the data by .) Any multiple of centered data give us that V from the eigendecomposition matches v from the SVD (again, to within signs of columns).

Generally, then, if we do an eigendecomposition of XTX or XXT, then the data must be X or a multiple of it. (In that form it looks pretty obvious, huh?)

We had gotten used to not seeing that N-1 or . Both Davis and Malinowski used eigendecompositions in preference to SVDs, but they computed eigendecompositions of not of . I was able to replicate their results by using an SVD of X, instead of their eigendecompositions, and we had that : the principal values were equal to the square roots of the eigenvalues.

But Harman and jolliffe used eigendecompositions of the correlation matrix, and the square roots of the eigenvalues are proportional instead of equal to the principal values,

First remark: the same thing happens if we use the covariance matrix and centered data instead of the correlation matrix and standardized data. (It is not uncommon for people to refer to as a “covariance matrix”; I try to reserve that term for with X centered.)

Second remark: we saw that Harman’s model and Bartholomew’s is to factor the (preferably standardized, possibly centered) data matrix using a weighted matrix A:

and that in fact we could write that same factoring of X using the SVD:

Third remark: Jolliffe wrote his model with symbols similar to Harman’s but with two completely different meanings: he writes

Z = X A

but this time Z is the principal components instead of , and his loadings A are the orthogonal eigenvector matrix (i.e. my V), and X is still the data matrix (with variables in columns). That is, in my usual notation,

Y = X V,

where I have written Y for Z because the product X V is precsiely the new components of the data wrt V. His “principal components” are the new components of X.

Now, we can write that as a factoring of the data matrix, and because V is orthogonal we get

Funny thing, that looks just like Malinowski’s

It is.

Malinowski, of course, uses his decomposition for an arbitrary X: it may be standardized, or centered, or even raw data. Jolliffe uses his decomposition preferably for a correlation matrix, possibly for a covariance matrix. That difference aside, Jolliffe and Malinowski are using the same model.

Let me be explicit about something: we may start with v from the SVD and use that v in place of V in the eigendecomposition. Going the other way is harder: the SVD has consistent signs on u and v, but V may not be consistent with u. In the eigendecomposition, sign consistency is automatic between V and V^T, or between v and v^T. Ah, the analog of mixing u and V in the SVD would be to mix v and V in the eigendecomposition, but I can’t imagine that we would ever decompose the correlation matrix c as either or , using both v and V. That’s good, because it wouldn’t work in general.

Harman and Bartholomew are factoring the data as

which is similar to

.

while Jolliffe and Malinowski are factoring the data as

.

Note that I wrote singular values w not , hence “similar to” for Harman & Bartholomew et al.

The first are very nearly associating the singular values w with v, the second are associating the singular values with u. That strikes me as a fundamental difference.

What is common is that they all are factoring the data matrix as scores times loadings.

Unfortunately, Davis is only sometimes in that camp, and sometimes not. From his definitions on p. 504 – the ones I showed – his scores times his loadings do not reproduce the data. He defines loadings

and then defines scores as projections of the data X onto the loadings:

.

We know that by projecting data onto non-orthonormal bases and , Davis is computing components of the data wrt the reciprocal bases instead of wrt the bases.

I don’t know why he uses scores and loadings that do not factor the data matrix. About the closest I can get to factoring the data is (making w square, and u, v conformable):

,

which just isn’t the same as

.

(So close and yet so far.)

Incidentally, if you are looking in Davis, remember that his uppercase U and V are my lowercase v and u respectively, not u and v; and his is my w (made square).

Those definitions notwithstanding, Davis defines different scores on pp. 536-537, and this new definition is a factoring of the data matrix; more to the point, by this new definition of the scores, they can be computed as .

That is, when he discusses “factor analysis” around p. 536, he agrees with Harman & Bartholomew et al.

Moving on.

Tell me again why we get that pesky ? Answer: because we want the eigenvalues of the correlation matrix, not the singular values of the data matrix.

And why do we want them? Because the eigenvalues are numerically equal to the variances of the new data Y wrt the basis V, i.e. the variances of the columns of Y,

Y = X V.

Why do we care that the eigenvalues are equal to the variances of the new data? Answer: I suspect that people wanted to make inferences about the new components of the data without computing it, i.e. without actually computing the scores. We can tell what the variances would be just from the eigenvalues.)

And yet, the new data does not depend on the eigenvalues or the singular values. Take the orthogonal eigenvector matrix V, premultiply it by X (“project the data onto the basis vectors”), and we’ve gotten data with the new variances, i.e. new data with the property that it has the same total variance as the original data X, but the variance has been redistributed maximally.

What’s important is that the new data does not depend on the matrix A. so why did they create A instead of using the orthogonal V? Answer: the F matrix.

What does depend on the matrix A is the F matrix. That pesky factor of gives us that

,

i.e. is not only standardized but also uncorrelated. That’s in complete contrast to Y = X V, which has the redistributed variances. Maybe that was the purpose, and the creation of A is the means to that end. (To look at it another way: u is orthogonal, the cut-down u0 is orthonormal, and is standardized, instead of orthonormal.)

Nevertheless, we struggle to have it both ways: either the new data is XV, with redistributed variances, or the new data is F^T with unit variances. I would go so far as to say that Harman & Bartholomew et al. are assessing one thing (the variances of the new data Y), while computing another thing (standardized data).

People talk about both. We’ll apply all this to an example.

PCA / FA Example 7: Bartholomew et al.: the scores

edited 17 Oct 2008 to round-off the first two columns of 5 u.

This is the 3rd post about the Bartholomew et al. book.

Introduction

It would be convenient if Bartholomew’s model were one we had seen before.

It is.

We got their scores and loadings just by following their instructions. Although they didn’t use matrix notation, their equations amounted to

loadings = the -weighted eigenvector matrix,

“component score coefficients” = reciprocal basis vectors,

scores = X cs.

where V is an orthogonal eigenvector matrix of the correlation matrix, is the diagonal matrix of (nonzero) , X is the standardized data.

Let’s recall Harman’s model. He had

where Z was the transpose of the standardized data X:

and A is the ]-weighted eigenvector matrix, so we have

.

Transposing, we get

.

There are two crucial facts:

• Harman has factored the data matrix X as scores times loadings.
• The right factor is not orthogonal: A^T rather than V^T.

Terminology: Harman says the elements of F are the factor values, and the “coefficients of the factors” are loadings. I take it, therefore, that A is the loadings.

So can we compute F or F^T? In this case, yes: all the eigenvalues are nonzero, so A is invertible, and then

.

So let’s compute F^T. (You may have recognized . It is in fact the reciprocal basis to A, and that equation says that F^T is found by projecting the data onto the reciprocal basis. That’s exactly what Bartholomew had us do, without telling us that we were computing along the way.)

Before we do that, I want to be very careful about the signs of columns. In the previous post, for example, the A matrix which I computed differed from the A matrix which Bartholomew et al. computed. We know that’s ok, because eigenvectors of length 1 can be multiplied by -1 and still be eigenvectors of length 1. But I was not consistent about how I changed signs to confirm that our graphs matched: in fact while I had a different A matrix from them, I defined my “component score coefficients” to match theirs. That was not consistent.)

But now, if I want to check my computations by multiplying the scores times the loadings, I need a consistent set of signs for the columns.

Rewriting Bartholomew’s scores

My personal preference for achieving that is to start with the orthogonal v matrix from the SVD (singular value decomposition) and construct A from it. You’ll see why, below.

What follows seems to be my adopted recipe for doing PCA / FA, although it’s cluttered here by comparisons to Bartholomew’s calculations. You’ll see it again.

First, decide what “the data” is to be: raw, centered, standardized, etc. We know that Bartholomew used standardized data, so let X be the standardized data, and let

be its SVD.

I remind us that the dimensions are

X: 26 x 9

u: 26 x 26

w: 26 x 9

v: 9 x 9

That is, u and v are square (as they must be: they are orthogonal); and w is the same shape as the data X.

Right now, I wish to use Harman’s & Bartholomew’s model, i.e. the correlation matrix. Rather than do algebra in my head, I’m going to compute the correlation matrix of X… get just its eigenvalues …

That’s all I want the correlation matrix for: its eigenvalues. Instead of its eigenvector matrix V, with whatever signs Mathematica® returns, I use the v matrix from the SVD.

Then create a diagonal matrix whose entries are the square roots of those eigenvalues; that is, the diagonal entries are

Let me remind us that w0, the cut-down w, is the same shape as . That is, let w0 be the square 9×9 submatrix of w:

I promised in the previous post that there would be one caveat, and here it is: w0 is proportional to , but not equal to it.

The proportionality constant is ; in this case, 5 because we have 26 observations. We have

It’s that proportionality that I know I will mess up unless I explicitly start with the eigenvalues of the correlation matrix.

then the -weighted matrix A is

and its first 2 columns, rounded off, are…

Those are loadings. They happen to match Bartholomew’s choice of signs.

The “component score coefficients” are computed as -weighted matrix v… that inverse square root tells me, and you, I hope, that the “component score coefficients” are the reciprocal basis to A.

The first two columns of the “component score coefficients” are…

They, too, match Bartholomew et al., as we would expect, since the first two columns of A matched. (I have labeled those “cs” but they are only the first two columns of the cs matrix, which is 9×9. I will consistently abuse that notation, because I’m going to be displaying matrices which are equal but computed differently. I really am displaying the results of computations, not just copying the thing that I’m supposed to get.)

The scores are computed by projecting the data onto the “component score coefficients”. That told me that Bartholomew’s scores were the new components of the data wrt the A basis. Here they are, as Bartholomew would compute them:

I need to plot these, to confirm that they are the same as Bartholomew’s.

Yes.

We have reproduced the calculations in the previous post, but whereas I was changing signs all over the place in the previous post – in order to match their pictures – this time I know that I have a consistent set of signs for whatever calculations I choose to do.

Computing F transpose instead

Now I’m ready to get . Let’s compute and look at its first two columns: they should be the “component score coefficients”.

They are. Recall the “component score coefficients”, as computed according to Bartholomew:

They match.

Then F^T is the projection of X onto the reciprocal basis, and must agree with the scores as we already computed them:

Yes, the same.

Let’s check things. We can compute AF, compare it with Z = X^T. They should be equal…. I get differences of about 10^–15 or less.

Good.

Bartholomew’s “component score coefficients” are the columns of the reciprocal basis , and his scores are the matrix. Bartholomew et al. have factored the data matrix the same way Harman did. (Good. it’s nice that they agree with a classic text.)

There is another way

But there is another approach to all this. After all, we know how to factor a data matrix. In fact, I started this process by factoring it!

There’s just that one little caveat: the singular values w are not quite the square roots of the eigenvalues – because there’s a pesky factor of .

To be quite explicit, suppose we do an SVD of the standardized data X:

.

Now, compute the correlation matrix c..

but the eigendecomposition gave us

and we may choose, as I did, to set V = v. (That’s pretty major: we are requiring that u and V have consistent signs.) Then we see that

Hence (if we cut down w to w0 so that it is square)

or

.

(as we confirmed earlier.)

Now go back and rewrite the SVD (cutting down u to u0 that it is conformable, i.e. the right size) :

.

Collect terms:

Now recognize (since ) and we have

.

For A invertible, we may then infer that

.

in this case, with N = 26, we have .

In other words, the scores are 5 u0.

(I have to tell you that I found that his scores were 5 u0 before I understood why. I was a little annoyed that they were so much simpler than projections onto a reciprocal basis.)

So let’s compute 5 u, and take its first two columns. (In fact, we could compute the full matrix 5u rather than 5 u0; I would never take more than 9 columns of it, because that’s all the columns I have of data X, and the scores are the new components of the data.)

Here are the first two columns of 5 u:

edit 17 oct 2008: while it isn’t essential, it would help if i rounded those off, and showed F^T:

Those are, indeed, first two columns of the scores, and the first two columns of F^T.

It’s the computation of 5u that requires the signs of A to be consistent with the signs of u. That’s why I will use v from the SVD instead of V from the eigendecomposition of a correlation or covariance matrix, even though I will calculate their eigenvalues. (If, and it’s a big if, I choose to use the correlation or covariance matrices at all.)

I want to confirm the factoring of the data, that the product of the scores and the loadings is the standardized data. Let’s compute (5 u0) A^T. For this I need the cut-down u. Since A is 9×9, u0 must have 9 columns.

So I compute it. the difference between the product and the data, , is the zero matrix (to within 10^-15 on each entry).

Let me emphasize that the showed up because we used instead of w. If we’d never seen the eigenvalues of the correlation matrix, computing scores would have been simpler, because we’d have used w instead of the square roots .

But, of course, the scores and loadings would have been different: we would have gotten a different matrix A. Instead of

we would have had

i.e. .

Anyway, we know 3 equivalent ways to compute the scores for Harman’s or Bartholomew’s model:

• project onto the reciprocal basis: either as or as .
• compute .

Ah, they’re not quite equivalent: what we just worked is a case where A was invertible.

What if A is not invertible? Equivalently, what if there’s a zero eigenvalue?

That is a very real possibility, even for real data: if the rows sum to exactly 1, then, as we saw, both the correlation matrix and the covariance matrix have a zero eigenvalue.

In such a case, is probably the way to go, although we can make the other computations work. Let’s look at such a case next time..

PCA / FA Example 7: Bartholomew et al. Calculations

The familiar part

This is the second post about Example 7. We confirm their analysis, but we work with the computed correlation matrix rather than their published rounded-off correlation matrix.

I am pretty well standardizing my notation. V is an orthogonal eigenvector matrix from an eigendecomposition; is the associated eigenvalues, possibly as a list, possibly as a square diagonal matrix. is the square roots of , possibly as a list, possibly as a square matrix, and possibly as a matrix with rows of zeroes appended to make it the same shape as w (below).

Ah, X is a data matrix with observations in rows. Its transpose is .

The (full) singular value decomposition (SVD) is the product , with u and v orthogonal, w is generally rectangular, as it must be the same shape as X.

The cut-down SVD makes a square matrix w0 out of w, and I write . I only do it in order to make an invertible form of w. For the product to work, we must drop some columns of u, and possibly columns of v.

By the same token, if any eigenvalues are zero, we would drop them, and make the diagonal matrices and invertible. And, again, conformability requires that we drop columns of V. (There is no U, unless I also do an eigendecomposition of , which is not all that likely.)

A is a weighted matrix derived from V. It is usually defined as

,

which means that each column of V has been weighted by a common .

You will note that I am being very sloppy about the sizes of , , and V. I do distinguish things derived from the eigendecomposition (, , V, A) from things derived from the SVD (u, v, w); in addition, I usually distinguish the cut-down SVD (uo, vo, wo) from the full one.

Here are the eigenvalues of the computed correlation matrix:

I can get a scree plot and confirm their fig 5.6. That is, we simply plot the each eigenvalue against the order in which it occurred. Personally, i would have used percentages rather than the actual eigenvalues. Regardless, we would probably keep 3 or 4 of the eigenvalues. (We discussed both of these way back in the Jolliffe posts, which followed Harman, the beginning of the PCA / FA category.)

I start to say that their definition of loadings is, in fact, Davis’

but it’s not so. edited: = instead of a 2nd

Harman and Jolliffe both computed eigendecompositions of a correlation matrix; Davis and Malinowski both computed eigendecompositions of . There’s quite a difference. Do you recall it? (I will, later and in some detail.)

Their definition of loadings is the weighted eigenvector matrix: best to write it as .

First we get as a diagonal matrix with the following diagonal elements, which are :

For reference, here are the first three columns of the orthogonal eigenvector matrix V:

Get loadings as and display only the first 3 columns, rounded.

That’s their table 5.7. (Except for the sign of the 2nd column, but we know that’s ok: any eigenvector can be multiplied by any number, and any unit eigenvector can be multiplied by –1.)

Now plot the first two. But change the sign of the second (only to match them for the drawing).

That’s their fig. 5.7.

Next, we get “component score coefficients”. They define them as eigenvectors weighted by . (Doesn’t that division remind you of a reciprocal basis, namely the reciprocal basis of A? That’s exactly what it is.) For now, anyway, I simply compute what they tell me, and get, for the first two columns:

That’s their table 5.9. and those two columns are two of the reciprocal basis vectors to A.

The unfamiliar part

To compute scores, we need to introduce the data.

Well, what is the data?

Please don’t curse me out. That wasn’t as stupid a question as it sounded.

First off, we started with “data”. We had to: we couldn’t use the published correlation matrix because it was not positive definite. But what we started with was the “raw data”, and what we did to it was compute its correlation matrix.

So what’s the appropriate data for the analysis? After all, over the past several months we’ve discussed four possible forms of the data

1. raw
2. centered
3. standardized
4. small standard deviates

and every one of them has exactly the same correlation matrix. How do we choose?

There are only two sensible possibilities, but the fact is there are two, not just one. We can use standardized data, or we can use small standard deviates.

(In fact, we could use any multiple of the standardized data but why would we? Small standard deviates are one particular multiple, all divided by , where N is the number of observations.)

But why must we use some multiple (possibly 1) of the standardized data? Because then the eigenvector matrix V from the eigendecomposition of the correlation matrix matches the v matrix from the singular value decomposition of the standardized data (to within signs of columns).

Multiplying the standardized data matrix by any constant simply multiplies the w’s or the ’s, it has no effect on the orthogonal eigenvector matrices u, v, V.

Oh, if we had done an eigendecomposition of the covariance matrix instead of the correlation matrix, we would have to now choose our “data” as the centered data, not the standardized. The reasoning is the same: the eigenvector matrix V from the eigendecomposition of the covariance matrix matches the v matrix from the singular value decomposition of the centered data (or any constant multiple of it).

So now we do what most people take completely for granted, and we standardize the data, because we worked from the correlation matrix. But we could choose small standard deviates. (More to follow.)

I would check the variances; they are each 1, as they should be.

Henceforth, X is the matrix of standardized data. Here are its first two columns:

Finally, having decided upon the standardized data X, they compute the scores by projecting X onto the “component score coefficients”,

Hey, that’s why they computed the reciprocal basis: so they could get the components wrt the A basis by doing a projection. Here are the first two columns of the scores.

Let’s plot those numbers, as xy-coordinates.

That matches their figure 5.8.

To summarize: the familiar part was the eigendecomposition of the correlation matrix, yielding eigenvalues and their square roots , and the orthogonal eigenvector matrix V and the -weighted eigenvector matrix . We got a scree plot of the eigenvalues, and plotted the first two columns of A.

For the unfamiliar part, we computed the scores by computing the reciprocal basis (“component score coefficients”, i.e. weighted by the inverse square roots, or ), and then projecting the standardized data onto the reciprocal basis.

There are a couple of ways to get the scores without the reciprocal basis, but there is one wrinkle. Let’s look at it next time.

Review: reciprocal basis

if the reciprocal basis confuses you, look for it in the linear algebra category, and in the PCA / FA category…. or you could pick a month, and search for “reciprocal”.

Or, for this special case, when we only change the sizes but not the directions of an orthonormal basis, think of it as follows.

If we have old components of data, as a column vector z, and an orthogonal transition matrix V, then the new components of the data, as a column vector y, are described by

.

Transpose:

Let , since :

Post-multply by V (edit: not by ):

since V is orthogonal, i.e. .

What we’ve just seen is that the usual change-of-basis equation

or

which represents old and new data components as columns is consistent with the projection

,

which represents old and new data components as rows, if V is orthogonal.

Now suppose we were to make all the vectors in V larger, say by a factor of 2; then the new components y would have to be smaller by a factor of 2:

.

If we were to project X onto the new basis (2 V), however, we’d get numbers that were larger by a factor of 2 instead of smaller. That’s just wrong.

This what the reciprocal basis takes care of. If the new basis were 2 V, then the reciprocal basis would be V / 2. If the new basis has columns , then the reciprocal basis has columns .

And then we can get new components of the data by projecting onto the reciprocal basis. If you insist on projecting, then you need the reciprocal basis. (If you want the scores to be the new components of the data wrt the A basis. Not everyone does, apparently. Harman’s or Jolliffe’s definitions of the scores are not universal: Davis is an exception. More to follow.)

As I said, however, there are a couple of ways to get the scores without the reciprocal basis, but there is one wrinkle.

Let’s look at it next time.

PCA / FA Example 7: Bartholomew et al. Correlation matrix

edit 5 Oct 2008: I had omitted the word “constant”. see edit.

The following example comes from Bartholomew et al. “The Analysis and Interpretation of Multivariate Data for Social Scientists.”

It is an excellent example with which to wrap up PCA / FA. (There’s a lot we haven’t done, but it’s almost time for me to move on.)

The example is “employment in 26 European countries”, “eurojob” for short, from chapter 5 (either 1st or 2nd edition), and data for both editions is available at http://www.cmm.bris.ac.uk/team/amssd.shtml . Please note that I am using the 1st edition of the book, and the 1st edition data.

When I first worked this example, I knew that something interesting happened, but not why; and, there was one thing I didn’t understand at all, back then.

One reason why this example is so good is that they provide both the data and a correlation matrix. In fact, most of their analyses in chapter 5 seem to be based on the correlation matrix, except in those few cases when they compute scores (for which they need the data). We get to return to our starting point, using the correlation matrix.

Here’s the raw data.

I remark that the matrix is 26×9. What we have is: for 26 countries, the % of people employed in 9 categories (agriculture, mining, etc.).

This would be a good, rather than excellent, example to end with simply because it takes us back to our roots, using the correlation matrix for its PCA analysis. I compute the correlation matrix, and I save it; but I also round it to 2 places, and show it to you:

I do not quite match them. In fact, we differ in one number, at the end of the first row (and at the end of the first column, of course). To five places, I have the number as -0.56492. They show -.57 where I round to -.56.

That’s irrelevant. What is not irrelevant is the result of computing the eigenvalues of the rounded correlation matrix.

It could be a big deal, but it isn’t. Here be dragons, in either case. The point is that their apparent starting point is their printed (hence rounded) correlation matrix. (To be more precise: they clearly computed an eigendecomposition of the correlation matrix. I used their printed correlation matrix.) Whether I use theirs or mine, so long as it’s rounded to 2 places, life is interesting.

Let’s just do our thing: get the eigendecomposition of the rounded correlation matrix… and look at the eigenvalues.

I hope you took a deep breath right there. The smallest eigenvalue is negative. But the correlation matrix is supposed to be positive semi-definite. (Because it’s of the form .)

Alarm bells should be ringing in your head: eigenvalues of a positive semi-definite matrix are non-negative. Where did that negative eigenvalue come from?

Presumably the smallest eigenvalue is very close to 0, and the negative number is numerical error.

Well, sort of, but we can say more. Here are the eigenvalues of the unrounded correlation matrix:

We see that the last one is very small but positive. The computed correlation matrix is positive definite. The problem comes from rounding off the correlation matrix.

But that must be positive definite, too, you cry?

Why?

Symmetry is not enough to guarantee positive definite.

This is just the tip of an iceberg. There was nothing special about rounding to 2 digits, nothing special about a correlation matrix. It is not trivial to manipulate a positive (semi-) definite matrix and guarantee that subsequent matrices are still positive (semi-) definite. But that’s about all I’ll say about the numerical hassles.

As I said, the computed correlation matrix is, in fact, positive definite: its determinant and all its eigenvalues are positive. Mathematica® does just fine with it.

But the rounded correlation matrix – which is what they presented – is not positive definite: its determinant and one eigenvalue are negative.

So we now know two good reasons to provide data instead of just providing a correlation matrix: one, we can do more with the data; two, a rounded-off correlation matrix may not be positive semi-definite. And, bless them, they did provide the data.

Just this insight alone upgrades this to a very good example. And this negative eigenvalue was the “something interesting”. We’ll see very shortly why it happened. (Perhaps you can guess.)

To summarize. The smallest eigenvalue of the correlation matrix is extremely close to zero, but non-negative, as it should be. The real problem is that the correlation matrix is almost singular (not of full rank); then the rounded correlation matrix is almost singular, too, but the smallest eigenvalue went to the other side of zero.

The lesson we have learned is: if the correlation matrix is barely of full rank, working from a rounded version of it could work out badly.

Now, why is the correlation matrix almost singular? I.e. why is the smallest eigenvalue so close to zero?

I swear I had no idea that this example would tie in so closely with recent posts. I knew that inadvertent reduction of rank was important, but I didn’t have a clear reason why.

What might cause the smallest eigenvalue to go to zero? All we did was compute a correlation matrix. That is, all we did was implicitly center the columns.

Ah ha! Just what are the row sums of the raw data?

Here they are:

The min, mean, and max are 99.7, 99.9923, and 100.4 respectively.

They’re not constant, but they’re awfully close to constant. Hence, the smallest eigenvalue of the correlation matrix isn’t zero, but it’s awfully close to zero. I explained here that making (edit: constant) non-zero row sums followed by centering the columns leads to constant zero row sums.

So: we just got a real-life example of the inadvertent loss of rank caused by using the correlation matrix on constant-row-sum data. The computed correlation matrix is positive definite but almost singular (determinant almost zero but positive, smallest eigenvalue almost zero but positive); and then because we were nearly singular, the rounded correlation matrix could fail to be positive definite, and in fact failed to be positive semi-definite.

Next, we will confirm that part of their analysis based on the correlation matrix.