Steig’s Verification Statistics – Part 2 (A little more majic)

A guest post from Ryan O.  The second part of Ryan’s Steig et al. verification statistics.  It’s amazing what you can find by replicating a paper.  The title above is my own.

————————–

REPLICATING STEIG’S VERIFICATION STATISTICS – PART TWO

Fig. 1-Geographical assignments for ground data comparison.

.
This is kind of a misnomer since I won’t actually be doing any replication, but as there was a Part One, there needed to be a Part Two.  Originally I had intended to replicate the r, RE, and CE statistics in Table S2 of the Supplemental Information, but doing so would require that I first replicate the restricted 15-predictor TIR reconstruction, an activity that I consider has little value at this point.

Instead, this post will examine how well the main TIR reconstruction matches the ground data.  Note that this is something not done anywhere in Steig’s paper.  To recap, the verification statistics that Steig presents are:

.
1. r2, RE, and CE for the main TIR reconstruction (though the displayed images were of the PCA reconstruction) where the AVHRR data was compared to the reconstructed values.
2. r, RE, and CE for the AWS reconstruction (Table S1).
3. r, RE, and CE for the restricted 15-predictor TIR reconstruction.

.

There is a noticeable absence of verification statistics for the PCA reconstruction (except by accident) and pre-satellite verification statistics for the main TIR reconstruction.  Steig does not tell us how well the ground data matches the reconstruction.

.

Why is this important?

.

Remember that the methodology for the main TIR reconstruction is the following:

.
1. Perform cloudmasking/removal of outliers for the raw AVHRR data.
2. Perform single value decomposition of the cloudmasked data and extract the first 3 principle components.
3. Place the first 3 principle components alongside the data from 42 surface stations and use RegEM to impute the 1957-1982 values for the 3 principle components.
4. Generate the 1957-2006 reconstruction from the 3 principle components, which are now populated back to 1957.

.

I placed step 3 in italics because it is a step that requires some thought.  Aside from problems with RegEM, a critical prerequisite for this step to work properly is for the AVHRR data and the ground data to be measuring the same quantity.  If they are not measuring the same quantity – i.e., not properly calibrated to each other – the output of this step might very well be gibberish.

.
A check one could use to determine if the output is gibberish is to calculate verification statistics of the ground data vs. the reconstruction.  By itself this will not be able to distinguish between calibration problems, imputation problems, and resolution problems (insufficient number of PCs), but it will at least provide some insight as to whether there are problems that should be concerning.  Since Steig did not include such an analysis with the paper, we shall take it upon ourselves to provide them.

.
The procedure for doing this is a bit different than Part One because the ground stations are incomplete during different timeframes.  To handle this, I take the first 50% of data that actually exists and designate that as the calibration period.  The last 50% of data that actually exists is then the verification period.  I do this for each station, individually, and then calculate r, RE, and CE.  I then reverse the order of the calibration and verification periods, recalculate, and store the minimum values obtained.  This makes the procedure exactly analogous to the method used in the SI for the AVHRR data.  Additionally, to ensure that there are sufficient degrees of freedom for the comparison to be meaningful, I exclude any station without at least 4 years of history.  I apply the above procedure to all AWS stations and all manned stations except Campbell, Macquarie, Grytviken, Orcadas, and Signy.  The reason for excluding them is because the nearest gridcell is 100+ km away, so there are no equivalent reconstruction values to which I can compare them.

.
The plots are color-coded with the geographical assignments of Fig. 1.

.
First up, correlation coefficient:

Fig. 2 – Correlation coefficient (ground station to reconstruction).

.
Take special note of the relatively high correlation coefficients for the red group, which is West Antarctica, and the low coefficients for the blue group, which is East Antarctica.  This may be unexpected, and it will become important.

Now coefficient of determination:

Fig. 3 – Coefficient of determination (ground station to reconstruction).

.
From the above 2 plots, we might be tempted to draw the following two conclusions:  1) the correlation coefficients and coefficients of determination are fairly low; and, 2) the reconstruction performs the best in West Antarctica and the Ross Ice Shelf.  This would seem to be counter-intuitive, since the plots of the eigenvectors and the AVHRR-to-reconstruction verification statistics would have us assume that the best performance would be in East Antarctica.  However, before we draw too many conclusions, let’s look at RE and CE:

Fig. 4 – RE (ground station to reconstruction).

.
Hm.   In this one, the Peninsula and West Antarctica that look good.  But isn’t the Peninsula supposed to be the least well-reconstructed region?  Maybe CE will tell us something:

Fig. 5 – CE (ground station to reconstruction).

.
Ah.  On the CE plot, the Peninsula stations look suitably poor.  However, the poor performance of East Antarctica seems as unexpected as the better performance of West Antarctica and the Ross Ice Shelf.

Now let’s step back and think about this for a moment.  After a review of all of these plots, we would conclude:

.

1. The reconstruction does the worst job in East Antarctica on all plots – r, r2, RE, and CE.
2. The reconstruction does the best job in West Antarctica on all plots, with the Ross Ice Shelf running a close second.
3. The overall performance of the reconstruction with respect to the ground data is poor, and in the Peninsula and East Antarctica, simply taking the mean of the station data explains more variance than fitting the reconstruction to the ground data.
4. The trends from the reconstruction should therefore be the most accurate in West Antarctica and the least accurate in East Antarctica.

.
So let’s see if our conclusion #4 is on the mark.  The minimum number of data points required for the following plot is 48 months.

Fig. 6 – Difference in slope for common points (reconstruction minus ground data).

.
Funny.  The most faithful replication of slope is in East Antarctica!

.

Here are the region means and standard deviations:

.

1. Peninsula:  -0.4282 (mean);  0.5059 (standard deviation) – 20 stations
2. West Antarctica:  -0.3576 (mean);  1.1581 (standard deviation) – 9 stations
3. Ross Ice Shelf:  0.1782 (mean);  0.5595 (standard deviation) – 26 stations
4. East Antarctica:  -0.0358 (mean);  0.5215 (standard deviation) – 40 stations

.
Why is this so?

First, the sample sizes for the mean/sd calculations are low.  In all cases, zero is well within the 95% confidence intervals.  So we cannot say that any of the means are different from zero at any reasonable significance levels – except for perhaps the Peninsula.  It could be that East Antarctica performs the best simply because we have the largest number of samples in that area.  Indeed, the distribution of peaks as we go from left to right almost looks random.

.
Second, the signal-to-noise ratio is very low.  In a 10-year sample, the signal we are looking for is on the order of 0.5 deg C or less, while the noise exceeds 10 deg C peak-to-peak.  Because of this, r, r2, RE, and CE penalize missing the noise to a much greater extent than they penalize missing the signal.  This leads to preferentially high scores for matching the high-frequency noise over the low-frequency signal.  The problem is made worse by the potential for spurious correlations between noise of physically unrelated stations because of the short record lengths and spotty coverage.

.
In short, I believe had the authors presented plots of the relevant verification statistics for the ground stations – which show negative REs for about half of them and negative CEs for about 3/5ths – they might have had a more difficult time with publishing.  The reconstructed values simply do not match the ground records.

.
With that being said, I personally do not think the situation is intractable.  I think that a reasonable reconstruction can be done using both the satellite and ground data.  Hopefully over the next few weeks I will be able to present a method that does a fairer job of representing actual Antarctic temperatures.  It will be based in a large part on the work that Jeff Id has already done.

.
Additionally, I intend to run some experiments using manufactured temperature fields to perform sensitivity studies with RegEM and this particular implementation of PCA.  There should be some fairly generic rules that result.  There are also a few different statistical measures I intend to test to determine a set of tests that provide more statistical power than r, r2, RE, and CE alone.

.
For the interested, here are the station numbers as used in the above plots:

Fig. 7 – Ground station numbers and associated names for the plots.

.

Lastly, the R script for calculating the above quantities:

ssq=function(x) {sum(x^2)}
get.stats.matrix2=function(orig, est) {
###  Set up the placeholder matrix for intermediate calculations
stats.matrix.early=matrix(ncol=ncol(orig), nrow=14)
stats.matrix.late=matrix(ncol=ncol(orig), nrow=14)
###  Loop through each column
for(i in 1:ncol(orig)) {
###  Extract station i
orig.all=orig[, i]
est.all=est[, i]
###  Fill in NAs for all points that are not shared
map=is.na(est.all); orig.all[map]=NA
map=is.na(orig.all); est.all[map]=NA
###  Find linear trends for all common points
tr.orig=lm(orig.all~I(1:(length(orig.all))/120))
tr.est=lm(est.all~I(1:(length(est.all))/120))
###  Remove NAs for r, RE, and CE calculations
orig.all=na.omit(as.vector(orig.all))
est.all=na.omit(as.vector(est.all))
###  Split vectors into calibration and verification periods
orig.c=orig.all[1:trunc(length(orig.all)/2)]
orig.v=orig.all[(1+trunc(length(orig.all)/2)):length(orig.all)]
est.c=est.all[1:trunc(length(orig.all)/2)]
est.v=est.all[(1+trunc(length(orig.all)/2)):length(orig.all)]
###  Set up temporary storage for calculation results
stats.matrix.early[, i]=rep(9999, 14)
stats.matrix.late[, i]=rep(9999, 14)
###  Ensure at least 24 months of data in each period before calculating r, RE, and CE
if(length(orig.c)>23 & length(orig.v)>23 & length(est.c)>23 & length(est.v)>23) {
###  Get stats for early calibration (first 50%)
stats.matrix.early[, i]=c(get.stats(orig.c, orig.v, est.c, est.v), tr.orig$coef[2], tr.est$coef[2])
###  Get stats for late calibration (last 50%)
stats.matrix.late[, i]=c(get.stats(orig.v, orig.c, est.v, est.c), tr.orig$coef[2], tr.est$coef[2])
}
}
###  Find the minimums between early and late
stats.matrix.min=stats.matrix.early
map=stats.matrix.early>stats.matrix.late
stats.matrix.min[map]=stats.matrix.late[map]
###  Replace 9999’s with NAs
map=stats.matrix.early==9999
stats.matrix.early[map]=NA
map=stats.matrix.late==9999
stats.matrix.late[map]=NA
map=stats.matrix.min==9999
stats.matrix.min[map]=NA
###  Return a list with minimums, early cal stats, and late cal stats
stats.list=list(stats.matrix.min, stats.matrix.early, stats.matrix.late)
}
get.stats=function(orig.cal, orig.ver, est.cal, est.ver) {
###  Get calibration/verification period means
orig.mean.c=mean(orig.cal)
orig.mean.v=mean(orig.ver)
est.mean.c=mean(est.cal)
est.mean.v=mean(est.ver)
###  Get residuals
cal.r=orig.cal-est.cal
ver.r=orig.ver-est.ver
###  Calculate Hurst parameter
hurst=HurstK(c(orig.cal, orig.ver))
###  Calculate average explained variance [Cook et al (1999)]
R2c=1-ssq(cal.r)/ssq(orig.cal-est.mean.c)
###  Calculate Pearson correlation for verification period
###  [Cook et al (1999)]
pearson=cor.test(orig.ver, est.ver, method=”pearson”)
r.ver=pearson$est
r.ver.p=pearson$p.value
###  Calculate Pearson correlation for all available data
###  [Cook et al (1999)]
pearson=cor.test(c(orig.ver, orig.cal), c(est.ver, est.cal), method=”pearson”)
r.all=pearson$est
r.all.p=pearson$p.value
###  Calculate Kendall’s Tau for the verification period
tau=cor.test(orig.ver, est.ver, method=”kendall”)
tau.ver=tau$est
tau.ver.p=tau$p.value
###  Calculate Kendall’s Tau for all available data
tau=cor.test(c(orig.ver, orig.cal), c(est.ver, est.cal), method=”kendall”)
tau.all=tau$est
tau.all.p=tau$p.value
###  Calculate RE [Cook et al (1999)]
RE=1-ssq(orig.ver-est.ver)/ssq(orig.ver-orig.mean.c)
###  Calculate CE [Cook et al (1999)]
CE=1-ssq(orig.ver-est.ver)/ssq(orig.ver-orig.mean.v)
###  Return vector
stats=c(hurst, R2c, r.ver, r.ver.p, r.all, r.all.p, tau.ver, tau.ver.p, tau.all, tau.all.p, RE, CE)
stats
}
Note:  Ground data was taken from the MET READER site on 4/4/2009.  Recon data taken from Steig’s site (ant_recon.txt)