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.
First download the AVHRR recon data:
download.file(“http://faculty.washington.edu/steig/nature09data/ant_recon.txt”,”temp.dat”,mode=”wb”)
grid=scan(“temp.dat”,n= -1) # 3305400
#600*5509 #[1] 3305400
recon=ts(t(array(grid,dim=c(5509,600))),start=c(1957,1),freq=12)
dimnames(recon)[[2]]=1:5509
Now download the AVHRR grid info (changing to -180/180 long for plotting):
grid.info=1:5509;grid.info=data.frame(grid.info);names(grid.info)=”id”
grid.info$lat=scan(“http://faculty.washington.edu/steig/nature09data/data/Tir_lats.txt”)
grid.info$long=scan(“http://faculty.washington.edu/steig/nature09data/data/Tir_lons.txt”)
temp=grid.info$long>180
grid.info$long[temp]= -( 360- grid.info$long[temp])
Now transpose and do the singular value decomposition:
X=t(recon);dim(X); rm(recon) #5509 600
svd0=svd(X)
Next plot the eigenvalues. This proves that only 3 PCs contribute to the recon.
plot(svd0$d,xlim=c(0,100),main=”Steig AVHRR Eigenvalues”)
This yields the following plot:
For good order’s sake, the first six eigenvalues are:
round(svd0$d[1:6],7)
#[1] 2272.6863885 993.2727255 708.5284449 0.0000024 0.0000024 0.0000024
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= svd0$u[,1:3]%*% diag(svd0$d[1:3]) %*% t(svd0$v[,1:3])
range(A-X)
#[1] -6.12266e-08 6.07958e-08
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.)
#library(GDD)
library(mapproj)
col1=c(“#0000BB”,”#0014FF”,”#006BFF”,”#00C2FF”,”lightblue2″,”lightblue1″,”grey85″,
“lightgoldenrodyellow”, “lightgoldenrod”,”#FF8700″, “#FF5B00″,”#D70000″,”#800000”)
K=length(col1)
M=nrow(grid.info)
for (eigen in 1:3) {
grid.info$eof=svd0$u[,eigen];
ylim0=range(grid.info$eof,na.rm=T);ylim0 # 0.007698207 0.200473573
y0=ceiling ( max(abs(ylim0))/.005) *.005;y0 #.025
breaks0=seq(-y0,y0,2*y0/K); (N=length(breaks0))
h=function(x) { if (is.na(x)) h=NA else {
if(x< =max(breaks0)) h= min( (1:length(breaks0)) [ x<breaks0 ]) -1 else h= length(breaks0);
h[h==0]=1}
h}
grid.info$class=rep(NA,M); for(i in 1:M) grid.info$class[i]=h(grid.info$eof[i])##PLOT
# GDD(file=paste("d:/climate/images/2009/steig/avhdd.eof",eigen,".gif",sep=""),type="gif",h=420,w=420)
par(mar=c(3,3,2,1))
map('world',proj='azequidistant',orient=c(-90,0,0),ylim=c(-90,-60),fill=TRUE,col="white") #irrelevant
a=seq(-180,150,30); for(i in 1:length(a)) lines(mapproject(list(x=rep(a[i],181),y= c(-90:90) )),col=4,lty=2)
a=c(-180,0); for(i in 1:length(a)) lines(mapproject(list(x=rep(a[i],181),y= c(-90:90) )),col=4,lty=1)
a=seq(-50,-80,-10);
for(i in 1:4) lines(mapproject(list(y= rep(a[i],72),x= seq(-177.5,177.5,5) )),col=4,lty=2 )
points( mapproject(list(y= grid.info$lat,x= grid.info$long) ),pch=19,col=col1[grid.info$class],cex=.5 )
title(paste("Steig AVHRR Eigenvector ",eigen) )
# dev.off()
Here are the three eigenvectors:
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.
Climate Audit 
98 Comments
Steve, is it my lyin’ eyes, or is that very heavily weighted area in PC 2 suspiciously close to Harry?
Steve: I wondered that too. Harry is not on the bulls-eye, but it’s near.
Re: Mike B (#1),
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.
Steve:
A different color scheme that does not connote temperature might avoid confusion.
I’m used to red and blue for positive and negative. Plus color schemes take a while to work out. I’ve got a usable one right now. So no plans to change.
It’s also important to understand that the PC2 and PC3 are contrasts – so that something is being flipped over at some point in the calcs. Things flipping over are always worth watching out for.
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.
I believe the answer to this is that the truncation parameter is equivalent to the number of PC’s.
#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.
Re: Steve McIntyre (#7),
You may be right but I ran a whole bunch of different parameters when I started.
If you look at neigs=3 and regpar=3 and neigs=4 and regpar=3, the output is visually the same. Molan Labe pointed this out to me in the comment thread. I didn’t confirm it any other way.
Re: Jeff Id (#8),
Also, it is intersting to see the loss of amplitude pre-1980 when neigs=3 and regpar=3 vs neigs=3 and regpar=2. PC3 effects?
– Cliff
Re: Jeff Id (#8),
Jeff, regpar is the same in the two results. So that doesn’t test regpar.
Re: Steve McIntyre (#11),
I looked into it further doing some single stepping through matlab. neigs 3 and regpar=2
V21 and V11 is altered from V which in this case contains the 3 eigs down to two.
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:
and finally
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.
Re: Jeff Id (#17),
Theres a section of the code which limits regpar to less than neigs.
Yup, here it is pttls.
SteveM,
Is it time to ask if the AVHRR reconstruction file matches the plot on the cover of Nature? Has anyone checked?
Did Steig post an intermedary file instead of the reconstruction? It would not be the first time for such an event.
– Cliff
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.
Re: Kenneth Fritsch (#12),
Ryan O over at Lucia has an interesting write up…
http://rankexploits.com/musings/2009/steigs-antarctica-part-three-creative-mathemagic/#more-3409
Re: MJT (#16),
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.
Steve,
Do you have a reference for this property?
Thanks
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.
Re: Steve McIntyre (#21),
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 .
So good luck with your “experiments” Steve and keep us informed .
Re: Ross McKitrick (#19),
LOL Excellent.
Re: Ross McKitrick (#19), “the entropy production vector keeps vanishing into the matrix kernel”–you’re pulling our leg, right?
Steve,
I get an ‘Unexpected symbol error in function h’ error message at this line: if(x <= max(breaks0)) h= min( (1:length(breaks0)) [ x h[h==0]=1}.
I removed the extra space between the less than symbol and the equals sign and placed it in front of the max(breaks0) call. This didn’t fix the error.
Keep up the good work.
Steve: it’s something to do with wordpress.this came up before. I’ll post a text version some time.
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: MC (#23),
In this case a different instrument was used which actually measures surface temperature by IR emission rather than air temp done by microwave absorption.
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 (#24),
And on the ground, apparently they use a PIR… which people like to lie under. Can we see the station log entries where they remove the contaminated data?
I’m not sure exactly if you’re agreeing or not but here’s my take:
neigs is the number of decompositions calculated
regpar is a truncation of the number of eigs used. This happens by default in the TTLS analysis.
I’m not sure why anyone would ever have these variables different unless TTLS is used in a more general form sometimes.
Anyway, I’m going to work on other stuff.
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’.
Patrick M,
I don’t have a reference, but it follows from the definition of eigenvalue:
For matrix A, a value l is said to be an eigenvalue of A if there exists a non-zero eigenvector K such that AK = iK.
If you substitute the transpose of matrix A, the eigenvector is the transpose of K, but the eigenvalue remains the same.
Re: Dishman (#28),
Math correction:
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
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.
Re: jc-at-play (#37),
In the previous message, all the T’s are supposed to be superscripts (for transpose).
Re: jc-at-play (#36),
In response to Steve’s in-line reply:
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.”
Re: Mark T (#40),
Yes, I realized belatedly that I had been thinking only of the case where A is a square matrix.
Re: jc-at-play (#63),
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!
^RomanM: agreed, PC1 does resemble the elevation.
Mark
Steve,
Do you have a more accurate date for the change in amplitude of series 3?
Which series went from spliced to real data on that date?
Re: MarkR (#31),
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).
Re: MarkR (#31),
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
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).
Let me see if I understand this correctly:
TTLS produced a best-fit of 3 new constructed series from the raw series.
Prior to 1981, one or more of the raw series was filled with reconstructed data from other series.
Some time around 1981, one or raw series switched from reconstructed data to (presumably) actual data.
The switch was different enough from previous data that it warranted its own constructed series (series 3).
Prior to the switch, series 3 had only slight differences from the reconstructed data. After the switch, it reflected the newly introduced data.
Series 3 may be one or more stations. It seems likely to me that it is actually composed of several stations and is actually a composite of all stations which began recording real data after 1981.
This is just my stab at understanding, and I’m not prepared to defend it, but others can probably confirm or reject it.
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.
Re: Mark T (#40),
Oops, D is still diagonal, you just can’t say DT = D unless D is square.
Mark
Re: Mark T (#40),
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.
There ought to be a rule on this blog that after, say, 40 comments, Steve has to write a translation of this stuff for us folks in the cheap seats. My head hurts- I’m going to bed now.
Has ANYONE told Nature (journal) about ALL this?
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: VG (#46),
It is WAY to early for anybody to be:
1. Telling nature about this
2. Assuming that the ultimate result will render Steig et al useless or
3. Gloating
I don’t understand the impatience.
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.
Re: bender (#50),
The folks at real climate almost immediately admitted that East Antarctica has been cooling recently.
In the best case, this analysis will show that the Steig et al calculation of a long term continent wide trend was flawed.
In the more likely case, this analysis will show that the Steig et al trend is merely not statistically significant.
Steig et al wasn’t terribly earth shattering. Whether Antarctica is warming or cooling, Real Climate has told us that it is consistent with the model predictions.
The main consequence of Steig et al may be what it reveals about the utility (or lack there of) of Nature’s peer review process.
I doubt very much that they are eager to print such a analysis of their own failures.
Any accusation or broken review processes must be withheld until the evidence supporting it is absolutely rock solid, lest the scientific community dismiss it with a wave of its collective hand.
Re: Jason (#52),
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.
Re: bender (#53),
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.
Re: Jason (#54), (Steve Mc’s inline comment)
.
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.
Note to frustrated commenters (Kenneth F, Ryan O, Craig L, Geoff S, etc 🙂 )…
AFAIK, you are having to fill out a “captcha” when you comment, because Javascript is disabled in your browser.
I sometimes disable JS while writing comments, because the (usually handy) preview thing can make everything glacially slow… one solution: temporarily enable JS just before pressing “Submit Comment”.
Perhaps there’s a different comment preview plugin, that puts the preview in another tab so it is not always active.
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:
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.
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?
Re: RomanM (#55),
Roman, you observed:
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.
Re: RomanM (#55),
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).
Added note:
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.
Darn! Learn to read properly!
west = 18.66 PC1 +9.25 PC2 +10.07 PC3
In the methods section of the paper the authors state:
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.
Re: Jeff C. (#60),
I am more down to earth (ice?). I have been looking at topographic relief maps and to my eyes, PC1 sort of looks like this to me:
(From: http://terraweb.wr.usgs.gov/projects/Antarctica/accdemsr.html )
Jeff, you quoted the following from Steig:
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.
As for PC1 they say:
and PC2, they say:
Or it could just be what you get from doing principal components on spatially autocorrelated data – i.e. Chladni patterns.
Steve, this line:
doesn’t seem to make any sense. 😦
Steve: Download the ASCII version. It’s a wordpress problem that I don’t know how to fix.
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.
Re: Jean S (#70),
with regpar=1 (which has much to recommend it), there’s a slight negative trend from 1957 to the present. regpar=2 looks like regpar=3.
Re: Steve McIntyre (#71),
Yes, and running the analysis with standard RegEM parameters (ridge regression) seemed to reduce the positive trend considerably in my preliminary tests.
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)
Re: Steve McIntyre (#73)
It’s Prunuspicker’s rule. Use the number of PCs required to get the desired trend. No, wait, Preisendorfer.
Re: bender (#74), Not a comment to read drinking hot tea.
I can almost see the snark in your face every time you post one of these replies. Maybe it’s because Bender is my 2nd favorite Futurama character (my son calls him Vender). 🙂
Mark
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)
Regpar Neigs Slope (deg/decade)
1 3 -0.073
2 3 +0.164
3 3 +0.157
4 4 +0.168
5 5 +0.113
6 6 +0.091
Steve – I’m asuming your R version allows you to have a regpar larger than neigs? Was your second plot regpar = 6, neigs =3?
Re: Jeff C. (#76),
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.
Did someone see something different?
Re: Jeff Id (#77),
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.
Re: Jeff C. (#78), Jeff, I didn’t set up neigs as I couldn’t see any purpose for it.
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.
Mark
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.
A general question on PCA . . . if your processing power allows, why would you limit the number of PCs?
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).
Mark
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. 🙂
Re: Ryan O (#84),
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
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.
Mark
Re: Mark T (#86),
.
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. 🙂
Re: Ryan O (#87),
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).
Mark
Thanks, Mark. 🙂
.
If you want to see the eignenvectors in action based on how they predict temperature, here’s a little movie in R. The two text files are 1) with the AWS recon and manned stations centered and spliced; and, 2) with just the manned stations centered and spliced. They’re easy to pick out . . . as are the eigenvectors.
.
To use, set sat_recon to the matrix you want (either spliced.manned or spliced.all). The rest of the instructions are commented in the last 4 lines of the script.
.
http://www.mediafire.com/file/imkymyqz4fw/MeansPlot.R
.
http://www.mediafire.com/file/mnyw15yy51m/spliceall.txt
.
http://www.mediafire.com/file/5qdzcuxynut/splicemanned.txt
Re: Ryan O (#89),
The R script download link wasn’t working. Here it is again.
http://www.mediafire.com/?imkymyqz4fw
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.
Mark
Re: Mark T (#90),
.
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.
Mark
Hopefully not too OT . . . I started looking at the PCA reconstruction and eigenvector 2 looks to be the inverse of the AVHRR eigenvector 2. I didn’t know where else to put this . . .
.
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: Jeff C. (#95), After thing for a few minutes (should do that before I type)the fact that it flipped probably isn’t that significant. Ryan – can you put up the plots of the 3 PCs for this recon?
Re: Jeff C. (#96), Here’s the PCs:
.
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.