The Three AVHRR PCs

Roman Mureika perceptively observed about 10 days ago that the AWS reconstruction consisted of only 3 PCs. Jeff Id has extended this to observing that the AVHRR reconstruction consists of only 3 PCs. Jeff’s demonstration of this is correct, but a little awkward. I’ll show an alternate demonstration which shows a useful linear algebra property of singular value (eigenvalue) decompositions: the SVD of the transpose of a matrix is the same as the SVD of a matrix. If Jeff had simply transposed the short fat matrix to a long thin matrix, his algebra would have been simplified. I’ve also extended his analysis to plot the 3 eigenvectors (corresponding to the 3 PCs) on a geographic grid. Steig has archived something that’s been fitted. There is no possibility that the AVHRR archive represents original data.

If you’re wondering whether there might be some variation in this data over and above the three PC fit (notwithstanding the eigenvalues – say, you don’t exactly know what an eigenvalue is). Here’s a further simple proof that the AVHRR is fitted. Expand the 3 PC- 3 eigenvectors and compare the fit. The difference is negligible. As the others have noted, Steig’s AVHRR archive is a fitted product.There is NO possibility that this is original data.

A plot of the 3 PCs (identical to Jeff Id’s plot) can be produced as follows:

plot.ts(svd0$v[,1:3],main=”Stieg AVHRR PCs”)

This yields the following diagram – with the interesting change in amplitude of the PC3 in the early 1980s – about the time when AVHRR measurements actually began to exist.

Jeff did not plot the three corresponding eigenvectors. This is a bit trickier, but we have the grid info on the location of the AVHRR gridcells and, by color coding the eigenvector values, these can be illustrated in the same way as trend maps. The maps below do NOT show trends. They show the weights in each of the PCs. The following code produces the three maps. Using the method here, I’ve plotted each pixel as a round point. You can control the point size through cex; in each map using this method, I usually have to guess at the size and experiment a little to get a map that looks OK. (If nothing else, this shows that it doesn’t take much code in modern languages to produce pretty graphics.)

It is possible that this archive is derived from RegEM PTTLS, but, if so, it’s done in a different way than what we’ve seen for the AWS data by parsing Tapio Schneider Matlab code. In that code, real data where available is spliced with the fitted data. In the AVHRR case, it doesn’t look like any real data is spliced back in – if there was, you wouldn’t get such a perfect fit.

There’s another potentially interesting phenomenon here – in two cases, when RegEM is applied to truncated total least squares in which the truncation parameter is set at 3, the fitted output is made of 3 PCs. Is there a connection between the truncation parameter and the number of PCs required to decompose the solution? If so, is the RegEM methodology simply a long-winded way of converging to something that can be calculated directly? Sort of like a slow-converging infinite series where the answer can be calculated directly. Just a thought.

Steve, I’m not an Antarctic geographer by any means, but it looks to me like the maps you’re using include ice shelves as part of the antarctic land mass. In particular, the “bullseye” in PC2 appears to be over the Ross Ice Shelf.

I don’t know if this is significant or not, just thought I’d point it out.

The intent of my original post when I began was to rerun the satellite reconstruction from the 3 PC’s. To that end, I had already extracted the eigenvectors for the second set of plots in the original code, I posted a truncated version of the code on line. On another post I showed that the AWS reconstruction can be emulated from using only the 3 pc’s or using the data itself. So now if we have these 3 pcs, we can run Jeff C’s gridded data next to the satellite recon in a matrix and get a less peninsula weighted reconstruction.

When I plotted the 3 pc’s everything stopped because the variance in the 3rd pc in the reconstructed data doesn’t match the ‘real’ data very well confirming the point RomanM made about the RegEM process having a different variance in the reconstructed data vs the real. To me that puts a serious credibility burden on the reconstructed data and it is clearly visible in the 3rd pc.

#4. what the PC3 indicates to me is that during the period where we actually have AVHRR measurements, there is contrast between West Antarctica and East Antartica – perhaps this is related to the Southern Annular Mode, perhaps not. Sometimes East is warm and West is cold, sometimes opposite.

In the pre-instrumental period, the contrast is pretty much extinguished. They move more in lock step. So for all the huffing and puffing about complicated covariance matrices, they haven’t maintained the amplitude of this series. Does it “matter”? It seems pretty relevant when the entire purpose of the analysis is to reconstruct 1957-1980 values and you fail to preserve the amplitude.

#6. Jeff, you’re simplifying things a little. Truncated TLS keeps 3 PCs in each of the 600 lines if the truncation parameter is 3. I don’t see that it is obvious from this that the RegEM fit after convergence likewise decomposes to 3 PCs. Indeed, we were all surprised when Roman noticed this.

Having said that, there’s reason to think that the 3 PC decomposition and truncation parameter of 3 are connected, but proving it seems non-trivial to me.

I looked into it further doing some single stepping through matlab. neigs 3 and regpar=2

rc = r(ir);
V11 = V(colA, 1:rc);
V21 = V(colB, 1:rc);

V21 and V11 is altered from V which in this case contains the 3 eigs down to two.

Xr(:,:,ir) = (V11 / (V11’*V11)) * V21′;

Xr is the return value from pttls and comes back as B. This is the only place Xr is used.

The only place B is used is in regem.m:

% missing value estimates
Xmis(j, kmisr{j}) = X(j, kavlr{j}) * B;

and finally

% update data matrix X
X(indmis) = Xmis(indmis);

The only other operations to X and Xmis I can find are scaling.

Since the regpar truncates the third vector and there’s no code anywhere else for it to have an effect, it looks like regpar is actually a truncation of neigs as Molan Labe stated.

Steve: regpar is the truncation parameter for the TTLS step. neigs enters at a different step – in the SVD decomposition of the correlation matrix C to give V entering the TTLS step. If you’ve already limited V through neigs, it may limit the scope of regpar. So the number of final PCs might be limited in two different ways. regpar probably has to be less than neigs to make sense, but I haven’t checked this.

I doubt that this will happen, but, in my judgment, a brief update or progress report on what has been elicited from the Steig et al. (2009) paper by one of the more technical participants in that endeavor would be helpful to the laypeople reading and posting here. I suspect that the full extent of the detective work done here could go unnoticed and under appreciated. As a novice using R and a layperson wanting to understand PCA and RegEM better these analyses has been a great learning experience. I have particularly appreciated Roman M’s patient and clear explanations of the methodologies involved – and help with R.

Without detail and from memory, my layperson’s view of the Steig paper analyses ongoing here at CA is the following.

The AWS and TIR reconstructions used the manned stations during the overlap period of measurements to calibrate and validate the reconstructions that then go back to pre-AWS and TIR (AVHRR) times to extract more fully covered data from that period than the manned stations would provide.

The AWS reconstruction used 3 AWS PCs for the entirety of the 1957-1980 period (or there abouts). For the period of actual AWS measurements (1980-2006) the measurements were used along with in-filled data from the 3 PCs. For the AWS reconstruction, the contributions of the PCs after the first 3 fall off perceptively. Could the contributions of more PCs affect the final reconstruction results? And is this less clear than in the case of the TIR reconstruction?

The TIR reconstruction used 3 TIR PCs for the entirety of the 1957-2006 time period in contrast to how the AWS reconstruction was handled. It would appear to me and from memory that the fall off from PC3 to PC4 for the TIR reconstruction is greater than in the AWS case.

Why would the Steig authors not use the TIR measurement data for the 1982-2006 time period? Could this change from the AWS case indicate that the authors had less faith in their adjusted TIR measurement data then the AWS measurement data? Do not the critical measurements in these reconstructions come from the manned station measurements during the 1980-2006 and 1982-2006 time periods? Evidently the manned station measurements, going back in time to the beginnings in 1957 when a number of these stations started measurements, are too sparse spatially (and temporally?) to be used alone to track trends from 1957-2006. That development appears to make the PCA and RegEM methods and resulting reconstructions provide more complete coverage during the pre-AWS and TIR period. But does it really?

The trends for the AWS reconstruction and the GISS temperature series for the zone 64S-90S for the period 1957-2006 (and time period segments within the 1957-2006 period) are in reasonably good agreement. I believe that GISS uses only manned station measurements going back from 1980 to 1957 and no TIR measurements for all of the 1957-2006 time period. I should determine whether how many of the AWS stations GISS uses.

In other words, the null hypothesis that there is no statistically significant difference between the satellite reconstruction and the ground data/AWS reconstruction is strongly rejected.

I think Ryan was showing this for individual stations. I did a difference comparison with annually reconstructed data between AWS and TIR using the trend of the difference over the period 1957-2006. I saw an increasing trend for the difference TIR minus AWS coming forward from 1957 to 1983. The significance with p less than 0.05 did not occur until I used the 1983-2006 time period. I judged that one could a prior select the 1982-2006 period for a comparison since that is the first year that TIR measurements were available. That period had a p = 0.07. I am not so sure one could find a good reason for the 1983-2006 period.

What I find in these reconstructions is a good deal of uncertainty in warming and cooling trends and differences between AWS and TIR.

#12. There sure is a lot of material and I’m sympathetic to the request, but I don’t think that things are ready to be pulled together. Yes, I will make a go at it at some point, but first we have to figure out what Steig did. We’re still guessing.

I’ll show an alternate demonstration which shows a useful linear algebra property of singular value (eigenvalue) decompositions: the SVD of the transpose of a matrix is the same as the SVD of a matrix.

There’s an alternative explanation: maybe the actual Antarctic surface and atmosphere is not of full rank. Outside the circumpolar vortex the climate system has hundreds of degrees of freedom, whereas inside, everything is a linear combination of only 3 vectors. This might also explain why it’s so cold there: the entropy production vector keeps vanishing into the matrix kernel.
Geez I bet nobody’s ever thought of this before, but it would definitely explain your PC results.

Re: Ross McKitrick (#19), remember the IMO interesting posts last year on spatial eigenfunctions (arising out of the humdrum Stahle SWM network)? I produced things that pretty clearly equated to Chladni figures. I’ve been experimenting with similar spatial covariance things in the Antarctic. (Gavin, get busy.)

But I’m 99.999% sure that the 3 PC thing has nothing to do with anything physical, but is a property of Regem truncated TTLS with regpar=3.

remember the IMO interesting posts last year on spatial eigenfunctions (arising out of the humdrum Stahle SWM network)? I produced things that pretty clearly equated to Chladni figures. I’ve been experimenting with similar spatial covariance things in the Antarctic. (Gavin, get busy.)

I am waiting with impatience on the results of this experiment .
I remember well the thread about spatial autocorrelations and consider that this issue is largely misunderstood and understudied .
If I had not a job or if I was lucky to be paid for doing “climatology” , I would focus on spatial patterns , stability of chaotic attractors and relations between both .

Normally when a physicist finds that a process A (temperatures) analysed in a creative manner presents surprising similitudes to apparently unrelated processes B (Chladni figures) and C (Eigenfunctions for rectangular potential pits in QM) , he is either on something interesting or has incredible bad luck .

This is interesting. I’m beginning to see some sense in what they were trying to do. I need to get this paper and the data but it appears to be that by generating the PCs they get the inherent variance over time in the measurement set (or at least the set presented). They then fit linear trends to this and extrapolate into the past and then fill in the best estimates (from linear trend) at each grid points using each PC trend and eigenvalue weights.
So that if you undid the PC matrix transformation you would come back to a set of data that was a best estimate over the extrapolation period. That’s actually quite clever.
The problem is does your PC generation method catch all the data i.e. do 3 PC’s catch all the AVHRR data. It looks like they catch a large portion but you don’t know the values of PCs above 4. It may be a truncation.
Another interesting titbit was Eli Rabbett’s answer to a post I made at RC. Apparently the satelite data measures temps a few kms from the surface so you have to include a lapse rate estimate from AVHRR to actual surface before comparing apples to apples. They just considered trends and variation instead of actual estimates. Fascinating.

Re: Jeff Id (#24), So Jeff it is apples for apples then. Yet the stuff you guys show is that there is a disparity between AVHRR and AWS trends for each station in the calibration period. Sounds like we need a standback moment (and data) so we can have a proper look. What are the chances though?

Re: Jeff Id (#26), Yeah Jeff I was agreeing. The standback moment was a more general thing as in you guys have uncovered a very interesting property, enough for everyone else to say, ‘wait hold the machines a bit’.

It is NOT generally true that the transpose of an eigenvector is an eigenvector of the transpose. In fact, “transpose of a vector” doesn’t really make sense, in this context.

It IS true, however, that each eigenvalue of a matrix also is an eigenvalue of its transpose. This follows from the fact that the determinant of a matrix is equal to the determinant of its transpose, since each eigenvalue of A is a root of the equation det(rI – A) = 0.

Re: Patrick M. (#16),
I believe that strictly speaking, the SVD of a matrix is NOT the same as the SVD of its transpose, but something similar is true. Here’s a derivation from basic principles of linear algebra. The SVD of a matrix A takes the form A = U*D*VT, with D a diagonal matrix having some special properties. Since the transpose of a product is the product of the transposes, multiplied in reverse order, it follows that

AT = (U*D*VT)T = V*DT*UT= V*D*UT

using the fact that DT=D, since D is a diagonal matrix.

This is not exactly the same as the SVD of the original matrix A, although the D part remains unchanged.

Steve: C’mon. The SVD is identical for the relevant purpose. Of course, the U and V matrices are exchanged and transposed. You can make the thing wide and fat if you want by transposing again. You’re being a bit pedantic here. You know the point.

You may have mistaken me for someone who actually understands applied linear algebra. I trust your statement that “The SVD is identical for the relevant purpose;” but personally, I don’t really “know the point.”

Not sure if it is the point, but since SVD is a linear transformation, then the transpose of the SVD is the same as the SVD of the transpose (which is what Steve meant, but did not specifically state, I think). That implies that the eigenvalues (singular values) are identical, which is really the point, since I guess it is computationally easier to work more rows than columns (which seems to make sense from a brute-force computation standpoint).

Btw, I figured you simply mis-stated, as I did in my first follow-up, hehe. I probably just caught my boo-boo quicker because I’m not exactly busy at the office right now!

UPDATE: An additional question has been brought up related to why the data seems to be missing from the poles. Dr. Christy also responded:

As the spacecraft rolls over the pole it does so at an inclined orbit so
that the highest nadir latitude is about 82 deg with the scanner looking
out a bit closer to the pole. Since we apply the scan line data mostly to
the nadir area directly below the satellite, the actual data only go to
about 83 deg. In the gridded data I interpolate over the pole, but I
wouldn’t trust the data too much beyond 85 deg.

This is a different satellite instrument. The instrument you referenced is the microwave sounder which measures temperature through microwave absorption in a thickness of air. The data in the antarctic paper is actually physical surface temperature measurement by infrared emission (not surface air temp).

As the spacecraft rolls over the pole it does so at an inclined orbit so
that the highest nadir latitude is about 82 deg with the scanner looking
out a bit closer to the pole. Since we apply the scan line data mostly to
the nadir area directly below the satellite, the actual data only go to
about 83 deg.

It is always interesting to note that the ratio of the area south of 82S to the area south of 60s is about 7% and for the ratio of 83s to 60s is about 5% and for the ratio 85S to 60s is 3%.

Yes, this is the case here. The imputed values are obtained through the linear regression model (eq. (1) in Schneider’s paper) in (Reg)EM. Now the regression matrix (denoted by B) is calculated in RegEM-TTLS by truncated total least squares, which means that B is obtained by TLS truncated to the given rank. Now the truncation rank is given in Schneider’s code as

trunc = min([options.regpar, n-1, p]);

Since there are much more rows and columns (n and p) in input here, trunc is set to regpar (=3). In other words, for each time instant, the matrix B is of rank regpar. Hence, if you a looking only the final result, the data appears to be a linear combination of three vectors (PCs) although it is actually a rank three linear combination of higher number of vectors (take also into account that B is different for each time instant).

I have commented on this before over in the comments of Roman’s post. I apologize for banging the same drum, but I think it is relevant to this discussion.

Dr. Steig actually does give us a peek behind the curtain as to what the real cloud masked data looked like prior to any RegEM or reconstruction manipulation. On page 1 of the SI, he discussed various cloud masking methods explaning why his results differ from Comiso 2000. He included this graphic.

Look at figure C and figure D. Steig says these are the trends from the cloud masked data using his new methodology. This is supposed to be actual measured results, with date/locations where clouds covered the continent removed from the data set (i.e. cloud masked). Figure C is the trend for 1982 to 1999 (same as Comiso) and figure D is the trend for 1982 to 2006.

Now compare figure C to this plot. These are the trends for the 5509 gridcells of the satellite reconstruction windowed to the 1982 to 1999 time frame of figure C. I’ve tried to use the same color scheme as above.

Note that the huge area of cooling in the interior is gone. Also note that the spatial detail seems less. The reconstructed trends look to have been smoothed (as a result of being boiled down to 3 PCs presumably).

There is also the issue of why figure C looks so different from figure A, but that is a discussion for a future thread.

So – PC3 looks completely goofy. A few questions from a layman: Would Steig likely have graphed the PCs as a normal part of doing the research/calcs? And has anyone from the Hockey Team commented on this? I haven’t seen anything from them posted here, the Air Vent, or RC. Just curious.

D is only diagonal if A is square, though even if not, the diagonal elements of D are unchanged after the transpose. You cannot make the generalization in your last equation, i.e., DT is not generally equal to D, but their nonzero values are equal (obviously).

Mark

Steve: Nope. You get a diagonal D even from rectangular matrices. Please, I don’t want to spend space on this as there is nothing to disagree about.

See my comment in (#42), I corrected myself before you posted this. Your statement in the original post above, however, should be corrected, as the SVD of a matrix is not the same as the SVD of its transpose, though you do get the same singular values.

Mark

Steve: As I said above, you get fat and wide back simply by transposing the U and V matrices. we’re talking about the same point.

So Steig used only 3 of over 100 eigenvectors characterising the measured data.

Then, why didn’t he publish his results on the back of an envelope ?

Steve: I don’t object to 3 eigenvectors per se. IMO, there might be a case for only using one eigenvector, tho I’m still toying with this. Right now we’re just observing and analysing this sort of thing. The 3 eigenvector thing is interesting mathematically, that’s all that can be said for now.

Steve: We’re just exploring this data right now. Maybe something will come of it, maybe not. Right now, we don’t even know for sure what Steig did. Nor do we have any AVHRR data. Tho perhaps it’s time to ask Nature to deal with the AVHRR data unavailability problem.

Re: Jason (#47),
It is not too early to tell Nature the real story about Antarctica: no continental-scale warming in decades and any warming is restricted to the peninsula. Impatience because the issue is topical. Wait, and you’ve lost your audience. Steig’s study is not “useless”. For some special interests it’s very useful. That’s not to say it isn’t a major distortion of fact.

I doubt very much that they are eager to print such a analysis of their own failures.

When you submit a correspondence they have no choice but to hand it to an Editor to deal with it. Whether they’re eager to do so is irrelevant. As scientists they are compelled to at least pay lip service to the search for truth. Which means they can’t just dismiss a correspondence addressing substantive aspects of the science. They would look too ridiculous.

If what has been uncovered thus far were submitted to Nature, I think it is very safe say that the submission would be ignored, and Nature would not look terribly ridiculous in the process.

See the Strumpf vs. Liebowitz discussion in Ross’s recent paper. It was only bad luck that resulted in the editor’s decision making process even seeing the light of day.

Steve: I totally agree. Most readers are far too quick to see gotchas. My sense of this thing is there are probably a bazillion calculations going around in circles, but all you really have to work with from 1957-1980 are a dozen or so surface records and, try as you may, you’re going to end up with linear combinations of these records in the recon and there’s a limited amount of harm that you can do to this sort of data. The big interest for me (and hence the focus) is that I’ve been nibbling around the edges of RegEM for a couple of years, but never really dug into it. This data set gives an excellent foothold for analysing the method without all the complications of bristlecones and upside-down Tiljanders. And the extra input from the Jeffs and Roman has really moved the analysis along. I expect that the main dividend will be in the proxy studies as we’ll end up with usable RegEM tools.

I totally agree. Most readers are far too quick to see gotchas. My sense of this thing is there are probably a bazillion calculations going around in circles, but all you really have to work with from 1957-1980 are a dozen or so surface records and, try as you may, you’re going to end up with linear combinations of these records in the recon and there’s a limited amount of harm that you can do to this sort of data.

.
AFA investigating how the AVHRR reconstruction was done, yep, there’s a lot more that needs to be done before it could see the light of day. But there is no question that the reconstruction does not match the instrumental record for those records that are present, and the differences are systematic, not random.

What are the trend stats on eigenvalues 1 and 2? They appear to be temperature trends with slightly different phasing schedules; 1=continental, 2=regional. Red/Blue loadings therefore implying warming/cooling.
.
Eigenvalue 3 is bizarre.

I am not sure whether everyone understands just how pervasive the 3 PCs in the AVHRR reconstruction really are.

In Figure 2 of the paper, the authors compare the AWS and AVHRR reconstructions through a visual comparison (no statistics) of the annual averages of the series:

Solid black lines show results from reconstruction using infrared satellite data, averaged over all grid points for each region. Dashed lines show the average of reconstructed AWS data in each region.

Since this involves averaging 5509 grid points in one case and 63 AWS in the other, this has to produce results with a extremely high level of certainty, right? Not exactly.

We have seen in the AVHRR case that each of the grid points is a simple weighted sum of the same three PCs. What happens if we average them? A little high school algebra easily shows that the result must itself be a simple weighted sum of the same three PCs where the weights are the average of the weights of the series you are averaging. Thus, the East Antarctic, West Antarctic and the continental monthly averages are nothing more than simple linear combinations of the same three PCs. The uncertainty level is not reduced by the averaging process because of the high correlation of the series being averaged. It is easy to see that the annual averages will also be just the weighted sum of the annual averages of the three PCs as well.

For the AWS situation, it is a little better. For the time period before 1980, where the sequences are pure reconstructions, the same considerations hold for the monthly and annual averages. Post 1980, there is real data mixed in and this may offer some reduction in uncertainty levels. Note that any comparison of AWS and AVHRR prior to 1980 is basically controlled by the manned stations which are used to construct both sequences.

So what does the AVHRR monthly mean sequence look like?

We can continue from Steve’s script above to express that sequence in terms of the PCs.

pcs = svd0$v[,1:3]

#calculate monthly means
# east is defined as long greater than zero
antarc = rowMeans(recon)
east = rowMeans(recon[,grid.info$long>0])
west = rowMeans(recon[,!grid.info$long>0])

Since I am not sure of the paper’s definition of east, I have taken it to be longitude between 0 and 180. It is easy to change the above script if you wish to have a different definition. the results:

antarctic = 28.74 PC1 + 2.06 PC2 + 2.15 PC3

east = 34. 34 PC1 – 1.95 PC2 – 2.26 PC3

west = 18.66 PC1 + 8.25 PC2 +1.01 PC3

The difference between east and west becomes pretty evident in the influence of PC2. I think that the interpretations made in previous comments on the “meaning” of the PCs were pretty accurate in light of these results.

By the way, did you notice the fairly strong spatial coherence of the PCs evident in Steve Mc’s maps despite the fact that there was no specific spatial weighting involved in their calculation?

By the way, did you notice the fairly strong spatial coherence of the PCs evident in Steve Mc’s maps despite the fact that there was no specific spatial weighting involved in their calculation?

as I’ve mentioned, I’ve done some experiments to see what the eigenvectors would look like for a process in which you have exponentially decaying spatial autocorrelation and sites distributed as above. (I took a random sample of 500 of 5509 sites to make the calculation more manageable on my computer.) You get intriguing patterns and I can sort of convince myself that I see Chladni patterns (a la the posts on Stahle spatial covariance that we chatted about). If Antarctic were a perfect circle, my guess is that the match would be closer, but the asymmetry appears to emphasize some “tones” more than others.

If you have highly asymmetric sampling e.g. the AWS-surface network, the eigenvectors appear to be affected by the site pattern. The concentration of sites in some locations promotes those patterns.

By the way, did you notice the fairly strong spatial coherence of the PCs evident in Steve Mc’s maps despite the fact that there was no specific spatial weighting involved in their calculation?

I did notice that and plotted out the three eigenvector of the AWS reconstruction for comparison. While they look vaguely similar, there are significant differences, particularly in the West/East divide of #2 and #3. It makes sense that there should be a spatial pattern as there was a pattern in the original information (the AWS and satellite data). Since the pattern in the eigenvectors is different for the AWS and satellite recons, I think that means the station data was used differently in the infilling (or infilling and fitting for the satellite recon).

I forgot to show a verification that in fact the calculated results really are expressible as a combination of the three PCs. To do that, e.g. type summary(antarc.reg) in R. Notice that the R-squared is exactly 1 indicating a perfect linear fit. The same holds for the other two regressions.

The first three principal components are statistically separable and can be meaningfully related to important dynamical features of high-latitude Southern Hemisphere atmospheric circulation, as defined independently by extrapolar instrumental data. The first principal component is significantly correlated with the SAM index (the first principal component of sea-level-pressure or 500-hPa geopotential heights for 20u S–90u S), and the second principal component reflects the zonal wave-3 pattern, which contributes to the Antarctic dipole pattern of sea-ice anomalies in the Ross Sea and Weddell Sea sectors.

I did some googling of the SAM index and the zonal wave-3/Antarctic dipole and found some charts of these over Antarctica. They do look similar to Steve’s maps #1 and #2 above. Makes sense as the Antarctic temperature pattern should be related to the known climactic phenomena. This probably gave the authors confidence they were on the correct path.

The first three principal components are statistically separable and can be meaningfully related to important dynamical features of high-latitude Southern Hemisphere atmospheric circulation, as defined independently by extrapolar instrumental data.

Leaving PC1 and PC2 along for a moment, I wonder what “meaningful relationship” the PC3 has to an “important dynamical features of high-latitude Southern Hemisphere atmospheric circulation”. They don’t discuss it – I wonder why 🙂

Taking their assertion at face value, there was sure a remarkable increase in its amplitude around 1980-82.

Here are the three eigenvectors for Roman’s reconstruction of the AWS (the AWS recon from only the 3 PCs, not measured data). I ran these a few days ago but had set them aside. I’ll put them up for comparison to Steve’s maps above. I would think these should be comparable to the satellite recon eigenvectors since these are entirely fitted with no measured data.

Although there are some vague similarities to the satellite maps, they are notably different in magnitude and even sign over large sections of the continent. If these are well-correlated with “important dynamical features of high-latitude Southern Hemisphere atmospheric circulation”, I would think they’d agree better.

Re: Jeff C. (#70),
It seems to me that the “robustness” purported in Steig’s paper boils down to the fact that RegEM-TTLS with regpar=3 and PCA-based method (with 3 PCs) are ultimately very similar. So instead of thinking more of wonders of RegEM-TTLS, I think a more fruitful approach would be running data with different RegEM parameters. My preliminary investigations with default RegEM parameters (just remember to increase the maximum iteration value) and AWS data indicate that the final result is not so “robust” they indicate in the paper. Unfortunately, I’m very busy right now, and can’t investigate that any further for awhile. Maybe you, or someone who has already data and (Matlab) code ready to run, would like to continue to that direction.

Jean S, here are AWS recons for regpar=1 and regpar=6. A slight negative trend overall for regpar=1; and no trend since 1965 for regpar=6 (or 5). Looks like the regpar for maximum trend is 3. What a surprise. (This is all with old Harry etc)

This looks like it was covered earlier in the thread, but I’ll put it up FWIW. The Matlab version limits you to having regpar less than or equal to neigs. If regpar is larger than neigs, Matlab throws an error.

With this limitation, here are the Matlab slopes increasing neigs for regpar greater than 3 (AWS, old Harry, 1957-2006 anomaly calculation window)

Jeff, unless I missed something the neigs represents the number of eigenvectors calculated and the regpar is a truncation of the matrix. I didn’t see where the neigs were used anywhere else in the R code.

I’m not quite sure what is going on, and I wasn’t very clear in my comment. I was having difficulty replicating Steve’s second plot in #73 using Matlab. Does your comment in #22 mean that Matlab RegEM will force regpar to a value equal to neigs if you try and set regpar to a number larger than neigs?

I don’t understand why Matlab RegEM would have that limitation, but I was speculating that maybe Steve’s R version doesn’t, thus the discrepancy.

My guess would be that using neigs larger than regpar results in a situation, in the MATLAB code, in which neigs is used to index into the array after it has been truncated by regpar. MATLAB would, by design, generate an “index out of bounds” error since the dimension would no longer exist. Whether R does this or not I cannot say (in a compiled C program, the index would simply to into subsequent memory which may nor may not create a seg-fault), but (#77) seems to indicate that the R version compiled by Steve avoids this. Perhaps there needs to be some check to choose the smaller of the two, i.e., the smaller of the array dimension or the index.

Re: Mark T (#80), I think your close as that is the error I was getting (wording different but same intent), but it was when regpar was larger than neigs, not visa versa. I’m away from my computer, but I’ll experiment with neigs when I get back to my office.

Re: Ryan O (#82), It’s not necessarily an argument about processing power, but one of dimensionality (or rank). You could have, for example, 500 series, but the entire space is only described by 3 or 4 vectors. There’s no reason to maintain the full data set, and it may actually be detrimental to do so for other reasons (rounding, overflow, etc.). Also, part of the goal is to find the dominant signals in the underlying data, with the assumption that the smaller signals are either noise, or not relevant because they contribute so little.

Btw, there is a counterpart to PCA called minor component analysis which seeks to understand the information contained in the smaller components (at least, that’s my take, I’ve never directly studied MCA).

Re: Mark T (#83), Wasn’t quite what I was getting at . . . my terminology here is poor. Speaking from the point of someone who hasn’t done any PCA, I would continue including PCs until the results stopped changing. Above, we see significant differences (in terms of the magnitude of the effect we are investigating) by including #4, 5, and 6. This would indicate to me that PCs 1-3 are insufficient to fully capture the variation in the data set and PCs 4-6 (and maybe more) are relevant.
.
That is the context of the question. In a case where including additional PCs makes a difference in the results, what justification would there be to exclude them?
.
I apologize for the poor wording of my initial question. 🙂

Speaking from the point of someone who hasn’t done any PCA, I would continue including PCs until the results stopped changing. Above, we see significant differences (in terms of the magnitude of the effect we are investigating) by including #4, 5, and 6. This would indicate to me that PCs 1-3 are insufficient to fully capture the variation in the data set and PCs 4-6 (and maybe more) are relevant.

Well, that is what I was saying you should do, but that may not be the clear case here (I haven’t run the numbers myself). What you are referring to as “results stop changing” may also be more extreme than is required. If the MSE only changes by 1%, even though it may look to the eyes as significant, it doesn’t really mean anything in the end and is therefore not interesting. Now, if PCs 4-6 contain, say, 10% of the variance, than maybe it makes sense to include them, but even with that, they will only change the outcome no more than about 1 dB, which is pretty small. In other words, your question

In a case where including additional PCs makes a difference in the results, what justification would there be to exclude them?

needs to be put in terms of actual increase in error without them, which is, as I’ve discussed, rapidly diminishing below 10%.

Re: Jeff C. (#85), OK. Indexing either way, which is what it sounded like. I’d have to dig into the code to better understand. Maybe this weekend if you want me to.

What you are referring to as “results stop changing” may also be more extreme than is required.

.
What I meant was that it’s dependent on the effect you are looking for. 🙂
.
When it comes to the recons, we’re looking for 1/10 degree C changes in information with a noise amplitude of +/-5 to 10 degrees . . . in other words, our signal is very small compared to the noise and the MSE matters more than if our signal was on the order of 1 degree, or 2, or 5.
.
Not only that, we’re looking for trends in specific geographic areas. As Steve’s plots show, the eigenvectors have geographical significance. So while including an additional PC might not change the overall error much, it would seem imprudent to exclude them on that basis alone as they may reduce the error significantly in local regions.
.
Or am I missing something? Sorry if these seem overly pedantic or basic. 🙂

MSE matters more than if our signal was on the order of 1 degree, or 2, or 5.

MSE as I mentioned is a relative value, i.e., in terms of percentage of the total signal. I commented on PCs 4-6 in the other thread, btw, particularly noting why they might matter in this case (which I think is an unphysical representation).

Hehe, sounds like I have a “learn R” notation on my agenda now. 🙂 Thanks.

Maybe I’ll look at this stuff more this weekend. I’d be willing to bet, btw, that if you look at the impacts the various PCs have, PCs 3-6 will be unnoticeable pre-1980, and much larger post-1980, i.e., the percentage of variance they describe pre-1980 is near zero, but much larger post-1980, which averages out to the values used in the recon.

Maybe I’ll look at this stuff more this weekend. I’d be willing to bet, btw, that if you look at the impacts the various PCs have, PCs 3-6 will be unnoticeable pre-1980, and much larger post-1980, i.e., the percentage of variance they describe pre-1980 is near zero, but much larger post-1980, which averages out to the values used in the recon.

.
This one is a different situation than the AWS recon. In this case, there are only the 3 PCs because the entire recon is generated exclusively from them. Without the processed satellite data, we can’t see what 4-6 would be.
.
And I would be willing to bet that PCs 4-6 are significant, because the trends of the grid points corresponding to the manned station locations are very different than the trends from the actual ground observations.

Uh, I was speaking in general terms, but that’s cool, too, since I didn’t notice the distinction anyway – good to know when I actually delve into it. It’s getting confusing jumping between threads, I think.

Re: Ryan O (#94), Don’t forget your day job! Seriously, this is good stuff. The PCA is supposed to use similar methodolgy to the RegEM except they replace missing values with the monthly mean instead of letting RegEM do it. Why would the eigenvector flip? Doesn’t make any sense, particularly when they say “The first three principal components are statistically separable and can be meaningfully related to important dynamical features of high-latitude Southern Hemisphere atmospheric circulation” as they do in the paper.

There is also a 15 station recon using RegEM that may have some interesting features.

re: neigs
In the RegEm code neigs controls the number of eigenvalues actually calculated in the TTLS regression (peigs.m). Of course, it should be larger than “the number of eigenvalues actually used” (regpar). By default, neigs is set to the maximum (i.e., all eigenvalues are calculated), so you don’t have to worry about that. The purpose of neigs is purely computational: there is no point of calculating all eigenvalues of a large matrix if you only plan to keep a few of those.