Climate Audit

Basil

Let me get this straight. Each point in Figure 1 is an average, and not a single data point itself. Has that been taken into consideration in computed the regression SE’s?

• RomanM
  Posted Feb 26, 2009 at 7:15 PM | Permalink

  Re: Basil (#3),

  Basil, you are absolutely right! Each point is an average of 5509 grid point values. Wow!

  Only one problem… the values each happen to be a linear multiple of the same three PC series. With simple algebra, it isn’t difficult to show that that average is itself a linear combination of the same three series. In essence, it has the same amount of information in it as a single grid point series. The average of all these series is no better than any single one of the series because they are so completely interdependent with each other (becuase any three of them can predict every single other grid series).

  The averaging is a red herring. Steig et al. would have been just as well off by generating the continental “average” in the first place.

Hu McCulloch

RE Basil #3, RomanM #5,
The fact that each point in Fig 1 is an average of 5509 grid values is not in itself a problem. Nor is the fact that these 5509 values are just linear combinations of 3 PC series (aside from the fact that the 4 regions reported in my table in fact have only 3 independent degrees of freedom).

The real problem is how the 3 PCs were concocted from the 42 Antarctic + non-Antarctic station data. I suspect now that some of the stations must have had negative weights in their contribution to these 3 PCs, in order to get W. Antarctica to have such a strong trend in the late period. These negative weights would be the natural outcome of the negative empirical correlations that occur randomly when sufficiently many zero correlations are estimated with limited data, as documented by Steve on the recent Antarctic Correlations threads.

• Basil
  Posted Feb 26, 2009 at 9:35 PM | Permalink

  Re: Hu McCulloch (#6), Hu, it has been years since I was formally schooled in statistics, but I recall this being an issue. You’ve got to acknowledge the loss of degrees of freedom somehow. Now, I don’t know that it is a “big” issue, with that many observations (5509) for each data point, but it shouldn’t be dismissed. And then what about “errors in variables?”

  I’m not challenging your concerns about the PCA. Maybe that’s the most fatal issue here. But that doesn’t mean there are not other problems.

Willis Eschenbach

Steve, Quenouille gave a complex formula for calculating the loss in effective N when there is autocorrelation. This includes the autocorrelation, not just at lag 1, but at all lags. It is given by:

where r1, r2, etc. are the autocorrelation at lag 1, lag 2, and so on.

Is this the formula that you are simplifying above?

w.

Willis Eschenbach

Upon re-reading, I see that I should have addressed the above comment to Hu, who wrote the most interesting lead post.

w.

Manfred

#6,#12:
“The fact that each point in Fig 1 is an average of 5509 grid values is not in itself a problem. Nor is the fact that these 5509 values are just linear combinations of 3 PC series.”

the latter is quite an issue, because every other selection of the number of PC’s gives a lower temperature trend.

UC

Another important model of serial correlation, which UC has commented on here, but which goes beyond the scope of the present post, is that of long-memory or Fractionally Integrated errors.

A bit too short series for that kind of analysis. But for sure there must be independent white error process added to AR(1) to make the H0 model realistic. Try AR(1) with r=0.95 and variance 0.4 plus white noise with variance 1, for example. No more significant trends 😉

• UC
  Posted Feb 27, 2009 at 7:41 AM | Permalink

  Re: UC (#13),

  Try AR(1) with r=0.95 and variance 0.4 plus white noise with variance 1, for example.

  ..and to compare models it is always good to check what they predict, simple Kalman filter will do in this case:

  If I had time I would add trend uncertainties, but my wife said I must not audit any climate for a week, so be it 😉

Jean S

I suspect now that some of the stations must have had negative weights in their contribution to these 3 PCs

That’s easy to check. The following additions will do for Schneider’s RegEM code (similar additions to Steve’s R-code; Steve, I think this is something you also needed before in other contexts):
#1:

function [X, M, C, Xerr,BS] = regem(X, options)

#2:

BS = cell(1,n); % list for collecting all B’s
it = 0;
rdXmis = Inf;
while (it stagtol)

#3:

% missing value estimates
Xmis(j, kmisr{j}) = X(j, kavlr{j}) * B;
% collect all B’s
BS{j}=B;

% add up contribution from residual covariance matrices

Now call your RegEM as

[Xout, M, C, Xerr,BS] = regem(Xin, OPTIONS)

and you will have all regression coefficients for input values in a list BS (each index correspond a single record). My preliminary test with AWS reconstgruction indicates that there are plenty of negative station contributions (to imputed AWS values) as I also expected.

BTW, many of B’s are the same (all the records with the same “missing pattern”) as the algorithm does not account for the autocorrelation (only changes in “missing pattern”, which is relatively nonrandom in this application). That’s relatively funny in the light of Steig et al:

Unlike simple distance-weighting5,6 or similar7 calculations, application of RegEM takes into account temporal changes in the spatial covariance pattern, which depend on the relative importance of differing influences on Antarctic temperature at a given time.

Hu McCulloch

RE Willis, #5-7, I’ll get back to you later today. This will take a while to compose.

Kenneth Fritsch

Hu, thanks for formalizing this analysis. It appears that I am in reasonable agreement with your calculations for the entire Antarctica. The differences between what you and I calculate for the separate regions is probably caused by my arbitrary choices for divisions of the Antarctica into East, West minus the Peninsula and the Peninsula.

My calculations are here in Post #’s 164 and 341 at :http://www.climateaudit.org/?p=5151#comments

I am not certain how you finally use the Santer application for AR1 adjustments, but I obtain the following from the Santer Et al. (2008) text.

This persistence reduces the number of statistically independent time samples.
Following Santer et al. (2000a), we account for the on-independence of e(t) values by calculating an effective sample size ne: ne = nt (1 − r1)/(1 + r1)

By substituting ne – 2 for nt – 2 in Equation (5), the standard error can be adjusted for the effects of temporal autocorrelation (see Supporting Information). In the RSS example in Figure 2A, ne ≈ 16, and the adjusted standard error is over four times larger than the unadjusted standard error (Figure 2C).

Equation 5 gives the standard error squared as SE^2 = 1/(nt-2)*sum of [e(t)^2]

When I do this substitution the ratios of the unadjusted SE and adjusted SE becomes ((nt-2)/((nt/Factor)-2))^(1/2) where Factor = (1+r1)/(1-r1).

I have a question about using the AR1 correlation adjustment when that correlation cannot be shown to be statistically significant. There were a very few cases where AR1 was not significant in my calculations over the various time periods and regional divisions that I used. Also I see higher order ARn correlations that are statistically significant in my calculations. Should we be concerned about that when we attempt to adjust for autocorrelation?

I think it is most important to also keep the claims of Steig et al.(2009) in the perspective of the choice of starting time for the trends they report. There also should have been a discussion of the ice core reconstructions in Steig concerning the findings of West Antarctica that for the 20th century the warmest decade was from 1935-1945 – findings by Steig.

• bender
  Posted Feb 27, 2009 at 12:57 PM | Permalink

  Re: Kenneth Fritsch (#20),
  The higher order autocorrelations do matter. You can’t ignore them. See my comment to Steve on how to use R’s gls() function for specifying autocorrelated errors (all orders, low & high). I believe it was in “unthreaded”?

Layman Lurker

Re: Paul M #19

I don’t think Dr. Steig had any choice but to make this claim due to the decision to reconstruct with 3 PC’s. It imposed an artifact of .11C trend on the peninsula through reduced eigenvector weighting. As has been discussed progressively higher order eigenvectors would have eventually found something reasonable for the peninsula, but (if I understand this correctly) this would have come at the expense of the overall continental warming trend being claimed.

• David Jay
  Posted Feb 27, 2009 at 3:20 PM | Permalink

  Re: Layman Lurker (#21),

  It would likely mean that the trend doesn’t exceed the uncertainty – necessary to make the trend statisticaly “significant”.

  +.05 (+/-.07) per decade doesn’t get you the front cover of NATURE 😉

Kenneth Fritsch

I have compared the reconstructed TIR, reconstructed AWS and GISS trends over the periods 1957-2006, 1965-2007, 1975-2006 and 1980-2006 and one can readily see that those trends are statistically significant for the 1957-2006 period. Coming forward only a few years from 1957, the trends are no longer significant. AWS and GISS are in reasonably good agreement whereas the TIR trends diverge (warmer) from AWS and GISS coming forward in time.

Post #292
http://www.climateaudit.org/?p=5151#comments

I think the point in Steig et al. (2009) was to be able to show a warming trend were they could and then essentially ignore what occurs coming forwards (or backwards, for that matter), i.e. starting time sensitivity. They also do not detail differences between AWS and TIR over various time periods and for the 3 regions in the Antarctica.

Layman Lurker

#22 Keneth Fritsch

Hmm… let’s see we have the choice of start date, the autocorrelation with Chladni, the choice to truncate at K=3 rather than at a higher order, RegEm not being able to weight inputed stations by geography, a clearly un-natural artifact in PC3…

Layman Lurker

Continuing on #22

…Oh and CI issues (sorry Hu)

p.j. moran

I understand this is really a mathematical/statistics oriented site. In this instance, I do not see any comment about the fact that Steig acknowledges that from the 70’s to the present, the East Antarctic cooled at an approximate rate of .2 C per decade. Setting aside the artificial temperature creation for the East for all years prior to the beginning of satellite date ( 1979 if I recall correctly ), Steig admits that the East has cooled about .6 C and the clear conclusion is the East is cooling with no indication of abatement. 85% of the Antarctic ice is in the East. The East was not melting in 1970 when it was warmer. It is not melting now. He selects 50 years and extrapolates some temperatures from the West, applies them to the East, and, voila,”on average” it was warmer in (the Western extrapolation transplanted to) the East. Steig then states that ‘on average’ it has warmed. This is bunkem on its face. If the East has been cooling for at least three decades, we can all sleep easily for the next 500,000 years or so. Steig’s own report demonstrates the critical area of ice, the East, is not only not melting, but any suggestion that the Antarctic ice will melt and flood the planet is a bold faced fabrication. No stats, just his own words.

• Kenneth Fritsch
  Posted Feb 28, 2009 at 8:56 AM | Permalink

  Re: p.j. moran (#32),

  I understand this is really a mathematical/statistics oriented site. In this instance, I do not see any comment about the fact that Steig acknowledges that from the 70’s to the present, the East Antarctic cooled at an approximate rate of .2 C per decade. Setting aside the artificial temperature creation for the East for all years prior to the beginning of satellite date ( 1979 if I recall correctly ), Steig admits that the East has cooled about .6 C and the clear conclusion is the East is cooling with no indication of abatement. 85% of the Antarctic ice is in the East.

  Coming forward in time from 1957 into the early 1960s and beyond all one can say (using Steig’s TIR and AWS reconstructions) is that the data yield CIs sufficiently wide to say that the overall Antarctica temperature anomaly trend is not significantly different than zero, i.e. neither cooling nor warming.

  That is the issue of this thread: wide CIs (after adjustment) equal uncertainty in temperature trends.

  If one breaks down the Antarctica into regions of East, West minus the Peninsula and the Peninsula then one can say (using Steigs reconstruction) that the warming of the Peninsula is statistically different than zero and has been warming continuosly over the past 50 years. But then we already suspected that was happening.

VG

Will CA audit this?
http://nsidc.com/arcticseaicenews/
“near real-time data” adjustment still = AGW story thank you

• Geoff Sherrington
  Posted Feb 28, 2009 at 2:04 AM | Permalink

  Re: VG (#33),

  From the article you reference,

  While the AMSR-E satellite can provide a good check on SSM/I data, scientists prefer not to use it to extend the long-term sea ice record. Although AMSR-E has a lower absolute error, continuing to use SSM/I provides a lower relative error, which is more important for tracking long-term changes in conditions.

  I disagree. Absolute ice can melt and produce diverse effects. Relative ice is a figment of the imagination, a mathematical device of convenience.

  Maybe some climate scientists should seek the advice of terminologically didactic statisticians.

Mike Lorrey

The choice of starting points is important given there was a global cooling trend from the 50’s to the 70s. It means he gets to start his series at a time when it was known to be cooler, so it stands to reason it would be warmer now on average.

I’d like to see a chart of an 11 year running average temperature change and correlate that against the solar cycles…

• bender
  Posted Feb 27, 2009 at 10:23 PM | Permalink

  Re: Mike Lorrey (#34),
  I wouldn’t.

tty

Re 36

There are other problems with the “relative ice”. There is too much of it during the 1989-2003 period as even a cursory look at the time series at Cryosphere (http://igloo.atmos.uiuc.edu/cgi-bin/test/print.sh) shows. Take a look at the ice in the Baltic on July 1, a date when there is positively never ever any ice there. During 1989-2003 there is always a fair amount of ice shown in northern and eastern parts. Presumably the ice in the Arctic Ocean is similarly overestimated.

• AnonyMoose
  Posted Feb 28, 2009 at 10:43 AM | Permalink

  Re: tty (#37), I wouldn’t be surprised if there never is ice in the Baltic on July 1, but I’m not familiar with the region. Is there a non-CryosphereToday source of ice info for that period, or some temperature charts from Finland?

Hu McCulloch

Re Willis Eschenbach, #5-7,

Q was concerned with the correlation between two variables, say y and x. Let’s say
y = a + b x + ε
and that Cov(&epsilon) = σ^2 Γ, for a general stationary correlation matrix Γ = (γ(|i-j|), γ(0) = 1.

WLG assume mean(x) = 0, so that X’X is diagonal and inv(X’X) is also trivially diagonal. Set Sxx = ∑x(t)^2. Then if

Γ = I,
varOLS(b-hat) = σ^2 / Sxx.

However, when Γ ≠ I, the general formula for Cov(β-hat) becomes
var(b-hat) = σ^2/(Sxx^2) ∑∑ x(i)x(j) γ(|i-j|)
  = (varOLS(b-hat)/Sxx) (Sxx + 2γ(1) ∑ x(i)x(i-1) + 2γ(2) ∑ x(i)x(i-2) … + 2 γ(n-1)x(n)x(1))
  = varOLS(b-hat) (1 + 2γ(1) w(1) + 2γ(2) w(2) + … + 2γ(n-1) w(n-1)),
where
w(i) = ∑(i = j+1:n)x(i)x(i-j) / Sxx
is what I call the weak autocorrelation of x at lag i. It is “weak” because n-i terms are being divided by n terms, so that this is a downward biased (but still informative) estimator of corr(x(t), x(t -i)). For large n and/or small i, the bias is negligible, however.

Quenouille (Associated Measurements, 1952,p. 168) wrote,

The effective number of independent observations to be used in testing the correlation between the two series is roughly N/(1 + 2r(1)r(1)’ + 2r(2)r(2)’ + …). …. This formula provides an approximate correction for serial correlation provided the number of observations n each series is large. [Q in fact used subscripts rather than parens, but the sub tag doesn’t seem to work in WordPress comments.]

Following the source you copied back on 11/24/06, in Comment 143 of thread 903, you have the sum running to n, rather than n-1 as it should. In fact a series of length n can have only n-1 autocorrelations. Q himself was not explicit about where to end the sum, since he treated the all but the first few terms as negligible in is examples.

If the two series are y and x, r(i)’ should in fact be the weak autocorrelation of x, not an unbiased strong autocorrelation. Q himself probably had in mind corr(x(1;n-1), x(2:n)), which downweights x(1)^2 and x(n)^2, but that is a minor distinction. Also, r(i) should be the autocorrelation of the errors, not of y, though under the null b = 0 these are the same thing.

When estimating a time trend, the w(i) are nearly 1 for low values of i. As long as the r(i) die out quickly enough, var(b-hat becomes
var(b-hat) ≈ varOLS(b-hat) (1 + 2γ(1) + 2γ(2) … )
If in addition the errors are AR(1),
var(b-hat) ≈ varOLS(b-hat) (1 + 2 ρ ∑(i=0:inf)ρ^i)
  = varOLS(b-hat) (1 + 2 ρ / (1-ρ))
  = varOLS(b-hat) (1 + ρ) / (1 – ρ).
The ordinary estimate of ρ is the r(1) of the errors, whence the reduction of the effective sample size is what I called “Q” in the post,
Q = (1-r(1))/(1+r(1)).

The above “Q-factor” formula is equation (6) in Santer+Nychka+15 (2008). They in turn attribute it to a 2000 Santer+Nychka+others paper that I have been unable to obtain. It is also discussed (and a bias-reducing variant proposed) in a 2000 UCAR technical paper by Nychka+Santer+4 that is no longer online at UCAR. (See CA thread 4216, “If Nychka Standards applied to Mann.”) Conceivably it is due to Santer and/or Nychka, but I strongly suspect it goes much further back. We ought to find out, just to set the record straight. Meanwhile, I think of it as a special, limiting cae of Q’s concept, and therefore reasonably attributable to him.

Hu McCulloch

RE Bender, #12,

Using R’s gls() function you can pass it your AR1=0.318 and have it do the regression in the presence of this serial correlation. Look at the “correlation” parameter described in help(“gls”).

I’m not familiar with R, but GLS+AR(1) would re-estimate the coefficients, with an effect similar to the almost-as-old-as-Quenouille Cochrane Orcutt (aka CORC) method.

However, the GLS estimates are optimal only if you know the covariance matrix up to a multiplicative constant. If you in fact are estimating the shape of this matrix from the data, the results can be unstable. Hayashi’s recent text therefore suggests sticking with the imperfect but stable OLS estimates and just worrying about their precision. That’s what my CovAR does in the text, and what the Q factor approximately does.

• bender
  Posted Feb 28, 2009 at 8:42 AM | Permalink

  Re: Hu McCulloch (#40),
  Hu, my primary point is this: with modern software these days it is SO BLOODY EASY to try to make a correction for serial correlation that failure to do so should be considered to be an egregious error.
  .
  That the correction is approximate is a secondary issue. At least you’ve got to try!

bender

OT (sea-ice): #33, #36, #37

Re: Hu McCulloch (#40),

the GLS estimates are optimal only if you know the covariance matrix up to a multiplicative constant. If you in fact are estimating the shape of this matrix from the data, the results can be unstable

In my books that goes without saying. But I suppose it doesn’t hurt to say it. Any time you are estimating something and then assuming it to be correct you are propagating error – a fact which is often (i.e. almost always) conveniently overlooked. The reality is the estimated AR(1) = 0.318 in Steve’s example has a standard error associated with it. To measure the degree of instability pointed at by Hu you could run your regression with the top and low ends of the 95% confidence interval, i.e. AR(1) = 0.318 ± 0.2, for example.

• Kenneth Fritsch
  Posted Feb 28, 2009 at 9:21 AM | Permalink

  Re: bender (#41),

  The reality is the estimated AR(1) = 0.318 in Steve’s example has a standard error associated with it. To measure the degree of instability pointed at by Hu you could run your regression with the top and low ends of the 95% confidence interval, i.e. AR(1) = 0.318 ± 0.2, for example.

  I asked this very question in previous post on Santer et al. (2008) about how one would handle the CIs for the estimation of the AR1 when using them to adjust the standard error for AR1 and did not receive an answer.

  It would appear that the GLS method is iterative like the Cochrane Orcutt method applied by Lucia in determining the statistical significance of current temperature trends deviating from model predictions.

  I could do the R version, but it would be like a blackbox application for me. Could not hurt to compare adjustments, I suppose.

  Hu, we prefer to see R code here at CA, but that from Matlab will do in a pinch.

Hu McCulloch

Re Bender #41, 42,

The reality is the estimated AR(1) = 0.318 in Steve’s example has a standard error associated with it. To measure the degree of instability pointed at by Hu you could run your regression with the top and low ends of the 95% confidence interval, i.e. AR(1) = 0.318 ± 0.2, for example.

Good point — the uncertainty in ρ affects the distribution of β-hat in much the same way that the uncertainty of σ^2 does. But in the former case there is no simple Student t distribution that untegrates out the uncertainty. An ambitious researcher would try to do this numerically, but I’m feeling lazy at the moment. (BTW the .318 is my figure, not Steve’s. As post author, I get to have power-pink comments!)

Re #42, yes Steig committed an egregious error by not even applying the simple Santer-Nychka Q factor. It didn’t much affect the significance of the 1957-2006 trends he reported, or even those for the 1979-03 period that apparently was used for the cover diagram, but it makes trend much more fragile and sensitive to starting/ending period when the serial correlation is taken into account.

Ken Fritsch has already noted this sensitivity in his posts #164 and 341 on the Deconstructing Steig thread. (See #20 above)

Hu McCulloch

RE VG#33, Geoff S. #36, tty #37, AnonyMoose #46, Larry H #47,
Arctic sea ice is way off topic on this thread. Please move this discussion to Arctic Sea Ice – End of Game Analysis (#3819), or else to Unthreaded.

I’d use my Superpowers as poster to move these posts, but I’m afraid I might hit the wrong button and set off an irreversible melt-down of the Antarctic ice sheet instead if I tried!

• Geoff Sherrington
  Posted Feb 28, 2009 at 7:46 PM | Permalink

  Re: Hu McCulloch (#48),

  You are quite correct and I apologise. I have a general concern that errors, their calculations and representation are often poorly done and now I have made one. Geoff.

tty

Re 46

Here are daily AMSR-E ice maps for the Baltic from 2002 on:

http://www.iup.physik.uni-bremen.de:8084/baltic/baltic.html

If you check the maps for 2003 you can see that the last ice in the Bothnian Bay melted around May 20 (which is a few days earlier than average).

There is a vast amount of historical data on ice in the Baltic, going back to about 1700 AD, but very little orf it is online.

Kenneth Fritsch

Hu, I have a couple of questions:

Do you include the Peninsula in your West Antarctica calculations?

By error do you mean the residual error from the regression line and in terms of standard deviations?

Hu McCulloch

RE Ken, Roman, #53, 54,

I just followed the definitions Steig et al gave on p. 461 of their paper. See definitions and discussion in the post above, after Table 2. “West” does not include the Peninsula.

They define E Ant as 60W clockwise to 180E, W Ant as 180E to 60W, 72S to 90S , and the Pensinsula as “westerly longitudes” (presumably their 180E to 60W), above 72S.

Steig’s definitions aren’t the ones I would have used — I would perhaps have taken the Peninsula as 50W to 80W, down to say 78S, W Ant as 80W to 170W above 85S, and everything else in E Ant. This would sort the AWS as in Ryan O’s graph at #316 on thread #5151. Including or excluding an AWS site could make a big difference for the AWS recons, but since the TIR recons are averages of smooth functions of location, the definitions of the boundaries probably aren’t as critical. Note that Steig et al seem to place Butler Is and Limbert, I think it is, in E Ant, and leave out the tip of the peninsula altogether. This may have affected their AWS recon, but not the TIR discussed here.

My “errors” are in °C, just y – yhat as defined in the post. Since the trend lines are small numbers, the “errors” aren’t very different from the temperatures. The outliers are essentially the same either way.

• Kenneth Fritsch
  Posted Mar 1, 2009 at 8:25 PM | Permalink

  Re: Hu McCulloch (#56),

  If Steig cuts off the peninsula at latitude 72S, he is putting some of the peninsula, as I would define it,into West Antarctica. I think I like my partitioning better than Steig’s although mine was one of convenience. I would think a proper one would consider the mountain range running down the Antarctica (Trans Antarctica Range??) as a natural dividing line for East and West and separate (all of) the peninsula out from West Antarctica.

  I have all the data in a form to determine how much the cut at 72S or 78S would affect the West Antarctica trends.

Hu McCulloch

RE Ken, #55,
The code you give is for computing “HAC” standard errors for OLS coefficient estimates. I regard these as defective, since they truncate the correlations above a low threshold (3 for sample size 100, but slowing increasing with sample size), and downweight those above the threshold linearly with lag length. This makes the standard errors much too small, but researchers love it because it makes their t-stats almost as high as OLS!

For say ρ = .9, HAC does do an adequate job of estimating the standard error for samples greater than 1,000,000, hence its “consistency.” But for smaller sample sizes there could be a problem.

Hu McCulloch

RE Bryan Watkins, #59,

You make a valid point. However, I would hope that Antarctic data is at least collected consistently, so that errors would average out to some extent over many stations and days. This doesn’t work when the errors are systematic, such as occasional changes in notions about how to read values that are between the lines, round fractions, etc.

Hu McCulloch

In trying to reconstruct the regional AWS trends Steig et al would have obtained (but didn’t numerically report), bear in mind that if you take their definitions literally, Butler Island (72.2S) and Limbert (75.4S) are not in the Peninsula. Steig et al round both their longitudes to 300E, which is right on their boundary between E and W Ant, so which way they would go would depend on the tie-breaker rule. (Though with Mann on board, they may have put them in both regions…)

BAS in fact gives the longitude of Butler as 60.2W, which would place it unambiguously in W Ant, and that of Limbert as 59.9W, which would place it unambiguously in E. Ant, by the Steig definitions. However, since Steig takes the trouble to give the sometimes-rounded coordinates as he has actually used them, and puts these at exactly 300E, these stations should be taken as Steig-ambiguous with respect to E/W. But since he rounded their longitudes but not their latitudes, they are unambiguously not in Steig’s Peninsula.

But of course, the classification of the AWS sites does not affect the regional TIR trends given numerically in the article and approximately replicated in the post.

Willis Eschenbach

Hu, the reference I have is an article by Nir Shaviv, wherein he references Rahmstorff’s citing of

Quenouille, M.H., Associated Measurements, 242 pp., Butterworth Scientific Publications, London, 1952.

as the source of the equation I cited above. Kinda third hand, but there it is.

w.

Kenneth Fritsch

I have reported, in the first table listed below, the results of comparing the AWS and TIR 1957-2006 reconstruction trends using three different bounds for the Antarctica Peninsula from West Antarctica. My bounds for the East Antarctica were determined using clockwise 0 to 180 degrees longitude, while those for West Antarctica were clockwise from 180 to 360 degrees longitude. My West Antarctica minus Peninsula was then defined by subtracting out from West Antarctica the Peninsula using the longitude clockwise from 285 to 330 degrees and everything north of -76S. To approximate what Steig et al. used I made the bound for the Peninsula -72S and to approximate what Hu M suggested (I think) -78S.

As expected by Hu M, the bounds changed the West Antarctica TIR trends not at all and changed the Peninsula TIR trends slightly, with the middle bound of -76S showing the more positive trend.

With AWS the change from -76S to -72S reduced the Peninsula stations from 7 to 5 and increased the West Antarctica stations from 14 to 16, i.e. transferred 2 stations from the Peninsula to West Antarctica. On changing the bounds from -76S to -78S, no exchange of AWS stations occurred and therefore the trend changes from -72S to -78S are the same as those to -76S. For AWS, unlike TIR, the boundary change to include more of the Peninsula in the Peninsula does take some nominal warming trend from West Antarctica and gives it back to the Peninsula.

None of the boundary changes appear to affect the number of trends that are significantly different than zero for the various time periods used. Of course, East Antarctica remains unchanged and has no statistically significant trends for any of the time periods 1957-2006, 1965-2006, 1970-2006 and 1980-2006 in either the AWS or TIR reconstructions. West Antarctica has a significant positive trend for the time period 1957-2006 for the AWS and TIR reconstructions and also for 1965-2006 for the TIR reconstruction. I incorrectly reported on this thread that the West Antarctica trend in the AWS reconstruction was not statistically significant for the 1957-2006 time period. The Peninsula has significant and positive trends for all time periods in both reconstructions.

In reviewing my R code and tables for reporting the 3 Antarctica regional trends on the Deconstruction thread, I discovered that I had made some errors in transferring the R script and the AWS AR1 correlations to the post. The errors do not affect the results but I need to show the correct code and tables in this post.

download.file(“http://faculty.washington.edu/steig/nature09data/data/ant_recon.txt”,mode=”wb”,”TempReconTIR”)
download.file(“http://faculty.washington.edu/steig/nature09data/data/Tir_lats.txt”,”TempLat”)
download.file(“http://faculty.washington.edu/steig/nature09data/data/Tir_lons.txt”,”TempLon”)
tir=read.table(“TempReconTIR”)
lat=read.table(“TempLat”)
lon=read.table(“TempLon”)
Tir=t(tir)
Tirll=cbind(lat,lon,Tir)
EWAnt=ifelse(test=Tirll[,2]285,yes=ifelse(test=WestAnt[,3]-76,yes=”Pen”, no=”NoPen”), no=”NoPen”),no=”NoPen”)
PenWest=cbind(Pen1, WestAnt)
Penisula=PenWest[PenWest$Pen1==”Pen”,]
WestMinusPen=PenWest[PenWest$Pen1==”NoPen”,]
dim(Penisula)
[1] 221 604
dim(WestMinusPen)
[1] 1730 604
dim(EastAnt)
[1] 3558 603
WDates1957=1:600
lmW1957=lm(WTemps1957~WDates1957)
summary(lmW1957)
acf(residuals(lmW1957))$acf[2]

Hu McCulloch

RE #69, thanks, Ken.

On trying to plot maps of the data, it occurred to me that equal averaging of the 5509 cells in the AVHRR (TIR) recon, or any subset of them, gives too much weight to the cells far from the pole, since they in fact have smaller area than the ones near the pole. In order to correct for this, they should be weighted by sin(phi)/phi, where phi is the angle in radians from the pole, i.e. phi = (lats + 90) * pi/180.

I’ll rerun the tables in my post with these weights when I get a chance. Failure to do this weighting would give too much weight to the peninsula and coast in the continent-wide averages, and too much weigh to the coast in the E and W averages. The effect might not be large, but it will be interesting to see.

It’s not clear whether or not Steig et al used these factors. Climatologists routinely area-weight their 5° global grids, so perhaps it went without saying that they used them, or perhaps they just didn’t bother because the adjustment isn’t huge, or perhaps it just didn’t occur to them.

In the AWS recon, this issue doesn’t arise, since it is just the average over the 63 AWS reconstructions, without regard for the fact that their distribution does not represent area at all.

Also, I’m inclining now to computing r1 simply as
∑(e(t)e(t-1)) / ∑(e(t)^2) * n/(n-1),
rather than using s2 with its factor of n-k in the denominator, and then trying to use an analogous factor of n-k-1 on the numerator, as in the post. Unless rho = 0, s2 is biased anyway, so there is no point in using it to compute r1. This change will be barely noticeable with n = 600 or even 300, however.

Without the factor of n/(n-1), this becomes the “weak” autocorrelation w(1) that arises in #38-39 above. But I prefer to use the above “strong” version, which even with this factor is biased downwards as rho approaches 1.

• Kenneth Fritsch
  Posted Mar 4, 2009 at 12:33 PM | Permalink

  Re: Hu McCulloch (#71),

  Hu, I used your sin weightings to obtain an approximate area for my designations of the 3 regions of the Antarctica and then compared them to the TIR portions and the AWS stations for these same areas.

  The breakdown of AWS station assignments by way of my area designations was as follows: East Antarctica = 42; West Antarctica minus the Peninsula = 14 and the Peninsula = 7. My crude calculation to approximate percent of these 3 areas (per my designation) to the total of Antarctica is: East = 64%, West minus the Peninsula= 30% and the Peninsula = 6%. The percent breakdown for the AWS stations is: East = 67%, West minus the Peninsula= 22% and the Peninsula = 11%. The breakdown for grid areas for the TIR reconstruction is: East = 65%, West minus the Peninsula= 31% and the Peninsula = 4%.

  Posted #342 here at :

  http://www.climateaudit.org/?p=5151

Hu McCulloch

RE Ken F, #72,
In that post, you say that the Peninsula definition you are using places 221 of 5509 gridcells in the Peninsula, or 4.01% of the total, as in your quote. When you weight the gridcells by sin(phi)/phi, what is the share of these 221?

• Kenneth Fritsch
  Posted Mar 4, 2009 at 3:12 PM | Permalink

  Re: Hu McCulloch (#73),

  Hu, I used latitudinal zones cut by longitude to make my area calculations. For example, using the sin ratios I can determine the relative areas of a zone going around the world from say -65S to -76S and apportioned the amount of this zone that is bounded by the longitudinal lines. My areas were approximate because I eye-balled the latitudes and longitudes. My area for the Peninsula was 6% of the Antarctica.

  I think my calcualtions are correct but I am not so sure about by eyes.

  A = 2*pi*R^2 [sin(lat1)-sin(lat2)][lon1-lon2]/360
  = (pi/180)R^2 [sin(lat1)-sin(lat2)] [lon1-lon2]

  http://mathforum.org/library/drmath/view/60748.html

Hu McCulloch

Ken —
It sounds like you are including a lot of open water in your zone-areas. Since the Peninsula zone contains more water than the other zones, you are therefore overstating the relative area of the Peninsula itself.

If your definition of the Peninsula places 221 of 5509 gridcells in it, then since sin(phi)/phi is generally smaller for these cells than for the non-Peninsula cells, the Pensinsula should end up with distinctly less than 221/5509 = 4.01% of the total area.

• Kenneth Fritsch
  Posted Mar 4, 2009 at 7:41 PM | Permalink

  Re: Hu McCulloch (#75),

  Hu, not to beat a dead horse here, but my eyeballing and estimating took into account the open water. There has to be a more direct way to determine whether the Steig grid cell sizes are of equal area. I would guess that they are and that would mean that the ratio of grid cell numbers would be equal to the ratio of areas.

  My estimates were just for putting the AWS and TIR weights for my designations of the three Antarctica regions in the ball park. I should have known that was not going to be good enough for someone pointing to the n-1 out of 600 differences. Seriously though if I understand what you attempting to discover here, I think we can quickly determine it. Did not the Steig paper mention the area of each grid cell?

Kenneth Fritsch

Antarctica (ice covered) = 13,720,000 km²

http://en.wikipedia.org/wiki/Antarctica

From Steig et al. (2009)

Monthly average surface temperature anomalies were obtained from the TIR data for the domain covering all land areas and ice shelves on the Antarctic continent, at 50kmby50km resolution.

13,720,000/50*50=5488. Is that sufficiently close to 5509?

Hu McCulloch

Re Ken F #78,
I ran the counts with and without sin(phi)/phi area weights, using the regional definitions given by Steig et al p. 461 (see above post, after table 2). Out of 5509 TIR coordinates, I get:
West — 1223/5509 = 22.20% of cells, 22.30% of area.
East — 4146/5509 = 75.26% of cells, 75.19% of area.
Pen. — 131/5509 = 2.38% of cells, 2.35% of area.
Other — 9/5509 = 0.16% of cells, 0.16% of area.

As I discuss in the post, Steig’s “peninsula” is rather restrictive. Also, there are 9 “orphan” cells in the tail of the Peninsula with lon > 300 and lat > -65 that don’t belong in any region, if Steig’s definitions are taken literally.

As it turns out, the area weighting doesn’t make much difference, even for the Peninsula, so I didn’t bother rerunning the trends with area weighting.

Although the Steig grid are of (approximately) equal area on an azimuthal equidistant projection, they are not of equal area on the ground. The sin(phi)/phi weighting (where phi is the angle in radians from the pole) takes this into account, but as it happens doesn’t matter much with this data.

When a polar plot (azimuthal equi-distant projection) of the cells is converted to “x” and “y” Cartesian coordinates, they turn out to be about 50.5 nominal km square, not 50. Also there is a slight tilt to the grid. However, in the “x” direction, round(x/50.5) (in Matlabian) gives an integer index, since essentially 0 is one of the x-values used. For some reason, in the y-direction, the first positive value is about +28, while the first negative value is about -28, so that this band is a little wider than the others, and the whole y-grid is offset by about 1/2 cell. However, round(y/50.5 + .5) gives a nice integer index. These integers can then be used to construct a nice matrix grid that can be input to surf( ) for color plots.

Although an azimuthal equi-distant plot measures true distances from the point of projection, distances are not true in other directions. Thus, the grid “squares” have equal “height” along the 0-180, but slightly shrinking true widths as you move away from the pole. Likewise, “squares” along the 90-270 line have equal “widths”, but shrinking heights as you move away from the pole.

Hu McCulloch

RE UC # 81,

With small sample, can we distinguish between fractional and ‘normal’ noise ? I’ve played with data sets with more than 200 000 samples, and then AR(1) + WN etc simple models can be easily thrown out if fractional noise is present.

With a “small” sample (like Steig’s 600 or 300), I’d guess you couldn’t tell them apart, but then an Information Criterion like Akaike’s would select FI, since it has only 2 parameters (d + one variance), whereas AR(1) + WN has 3 (rho + 2 variances, or equivalently, 2 ARMA parameters plus one variance). Either way, the long covariances are going to be jack the standard errors way up.

Squinting at Steig’s TIR data for all-Ant (Figure 1 in the above post, I see why you like AR(1) plus WN — it looks like WN on top of a series with perhaps only 6 swings from positive to negative and back. Do you have actual estimates of the ARMA parameters? I’m not setup to do this in Matlab at the moment.

Hu McCulloch

Re #38, 66, I have now obtained Mitchell’s 1966 WMO technical report, with the Santer-Nychka factor (1-r1)/(1+r1) that I am calling “Q”.

It’s there on p. 71, but he attributes it to Brooks and Carruthers, 1953, Handbook of Statistical Methods in Meteorology, Air Ministry M.O. 538, HM Stationery Office, London., and also to a 1960 article by R. Sneyers in French that appeared in a German journal that OSU doesn’t have.

Also, Mitchell is only using it to adjust the se of an estimate of the mean, not of a trend slope line. When estimating a mean, Quenouille’s r(i)’ are all exactly 1, which simplifies the derivation. But when estimating a trend or the slope on any other highly trending variable (like GDP), r(i)’ are very near 1 for value of i for which rho^(i-1) has not yet vanished. Therefore it works for a trend slope se as well, if not for a less highly correlated regressor.

Perhaps Nychka’s student and co-author on Nychka et al 2000, Rachel Buchberger, was the first to apply it to the trend slope se.

Anyway, I’ve sent for our library’s copy of Brooks and Carruthers to see what they have to say.

Hu McCulloch

RE UC #81,

The easiest way to estimate an AR(1)+WN process (aka ARMA(1,1)) is to match its moments to r(1) and r(2) (which is basically all OLS does to estimate an AR(2) model). But using the residuals of the all-Antarctica full sample TIR regression, r(1) = .3183 and r(2) = .3188. This would imply the AR(1) rho = r(1)/r(2) is slightly > 1!

It sounds like an AR(2) would raise fewer red flags. Or AR(p), where p is a fixed increasing function of sample size that guarantees consistent estimation of the covariation function.

(Above I used the “strong” sample autocorrelations I define in #71 above, rather than the “weak” autocorrelations w(i) defined in #38. For present purposes, it doesn’t much matter which are used. In any event, when estimated from regression residuals rather than the errors themselves, each is the ratio of two inconsistent estimates of error moments. For present purposes we can just live with that bias.)

• UC
  Posted Mar 9, 2009 at 4:08 AM | Permalink

  Re: Hu McCulloch (#88),

  But using the residuals of the all-Antarctica full sample TIR regression, r(1) = .3183 and r(2) = .3188. This would imply the AR(1) rho = r(1)/r(2) is slightly > 1!

  Something wrong here, using Stoffer’s notation, for h more than 0, the ACF of the observations is

  

  Should’n we start with r(2)/r(1) ? And is r(2)=0.3188 correct ?

Hu McCulloch

RE #23,

Bender, even though I’m not using R, I’ve downloaded the 2722 pp R Ref Index PDF, which presumably contains the same info as online Help. However, the entry for gls( ), on p. 2274, is pretty cryptic. Am I looking the wrong place?

Thanks!

Hu McCulloch

Re #90,

Should’n we start with r(2)/r(1) ? And is r(2)=0.3188 correct ?

Right on both counts, UC — The AR(1) ρ (Stoffer’s φ) would correspond to r(2)/r(1), not r(1)/r(2). Also, r(2) is only 0.1529 (I don’t know where I got 0.3188 from). This is more than r(1)^2, and would make ρ = 0.4804.

Using Quenouille’s general formula in #38, this would jack the variance of the trend coefficient up by a factor of approximately (1 + 2 r(1) ∑(0:inf)ρ^i) = (1 + 2r(1) – ρ)/(1-ρ), rather than (1 + ρ)/(1 – ρ) as in the pure AR(1) case.

Since r(1) = .3183, this makes the variance factor 2.22 and the se factor 1.49, instead of 1.93 and 1.39, resp., as without the WN. This whittles away at the already reduced significance of the trends, but still does not exactly overwhelm it.

Incidentally, in Stoffer’s formula, if y is the observed series, I think the ratio in the denominator should be var(noise)/var(AR innovation), rather than var(y)/var(ω), regardless of what ω is.

• UC
  Posted Mar 14, 2009 at 11:50 AM | Permalink

  Re: Hu McCulloch (#92),

  This is more than r(1)^2, and would make ? = 0.4804.

  Not sure what would be the best command in Matlab to try this (and one should always prefer the Ritson method ), but

  armax(data,[1 1])
  Discrete-time IDPOLY model: A(q)y(t) = C(q)e(t)
  A(q) = 1 – 0.5354 q^-1

  C(q) = 1 – 0.2208 q^-1

  Estimated using ARMAX from data set data Loss function 1.2284 and
  FPE 1.23665
  Sampling interval: 1

  suggests 0.5354, quite close to your estimate.

Hu McCulloch

It seems that the image that appears on the cover of Nature is posted at NASA at: http://earthobservatory.nasa.gov/IOTD/view.php?id=36736. It is calibrated to trend rate, and is supposed to represent the entire period 1957-2006:

(Click on image for 2.3MB blowup).

I had assumed from the text that it was based on Figure 4e, and pertained only to 1979-2003.

I first found this image incorporated into the Wikipedia Antarctica article, giving Steig et al as its source.

• curious
  Posted Mar 13, 2009 at 4:20 PM | Permalink

  Re: Hu McCulloch (#93), Hu – Just to clarify – the trend is degC/decade. That pic. has appeared all over without the accompanying scale. Here is the 1st para of the NASA release you link:

  “For a long time, it seemed that Antarctica was immune to global warming. Most of the icy southern continent, where temperatures can plummet to minus 80 degrees Celsius (-112 degrees Fahrenheit), seemed to be holding steady or even cooling as the rest of the planet warmed. But a new analysis of satellite and weather station data has shown that Antarctica has warmed at a rate of about 0.12 degrees Celsius (0.22 degrees F) per decade since 1957, for a total average temperature rise of 0.5 degrees Celsius (1 degree F).”

  Last para:

  “This allowed the group to map out the relationship between ground measurements and satellite measurements so that they knew roughly what the satellite temperature would be when the thermometer at a weather station registered -5 degrees Celsius, for example. The team used this relationship to extrapolate what the satellite would have recorded over the whole continent had it been in orbit when the weather station record began in 1957. Once the group reached the period when the satellite was in orbit, they checked the extrapolated values against the actual measurements to confirm that the method was sound. In the end, they generated a 50-year record of temperatures across Antarctica. Their work was published in the January 22, 2009, issue of Nature.”

  Sorry if this repeats what loads here know but the shot is so dramatic I think the words that go with it are important.

• UC
  Posted Mar 14, 2009 at 10:04 AM | Permalink

  Re: Hu McCulloch (#93),

  That seems to be the new hockey stick,

  http://yle.fi/uutiset/luonto_ja_ymparisto/2009/03/ilmasto_lampenee_pahimpiakin_ennusteita_nopeammin_613048.html

  ( in Finnish, the key message is that the worse than worst-case IPCC scenario trajectories are being realised )

Hu McCulloch

Re #95,
Thanks, Curious. Yes, the page is worth reading, and you have highlighted the most interesting parts. I had to load the image twice in order to get a version with the scale.

Given that co-authors Comiso and Shindell are both with NASA, I think this can be taken as a semi-official release from the authors.

Perhaps Comiso can be prevailed upon to post his pre- and post-processing TIR data, and cloud-screening algorithm, on a subpage of this page. This would not require imposing on Steig to put it on his site.

Hu McCulloch

RE UC, #98,

The image certainly sounds twice as scary when I see it described in Finnish! 😉

The impression one gets is that 1/4 of the frying pan is already on fire, and it’s just a matter of time until the conflagration spreads to the other 3/4.

I must admit that TIR EOF2 does an excellent job of delineating the Transantarctic Range, and therefore the boundary of W. Ant, just from AVHRR correlations and without using any topographical information. Kriging-like distance weighting is ordinarily more than adequate, but could never be quite so precise.

But until we see how the calibration to surface temperatures was done, and why PC3 unnaturally flatlines before 1982, there is no way of knowing how valid this recon is.

• UC
  Posted Mar 14, 2009 at 11:42 AM | Permalink

  Re: Hu McCulloch (#99),

  The impression one gets is that 1/4 of the frying pan is already on fire, and it’s just a matter of time until the conflagration spreads to the other 3/4.

  • Kenneth Fritsch
    Posted Mar 14, 2009 at 11:55 AM | Permalink

    Re: UC (#100),

    How do you put CIs on a picture?

    Is that picture in your link part of the new campaign I have heard that implores the climate scientists to better show the urgency of the AGW crisis? Can I assume that that ring of fire is being stoked with coal?

    • UC
      Posted Mar 14, 2009 at 11:59 AM | Permalink

      Re: Kenneth Fritsch (#102),

      new campaign

      It has been in Mann’s homepage for very very long time!

    • Kenneth Fritsch
      Posted Mar 15, 2009 at 9:36 AM | Permalink

      Re: UC (#104),

      You must have access to Mann’s site that I do not. Anyway I have never gone beyond the 10 pictures of the Mann himself that can be shown in 10 Megabyte renditions where I have attempted to estimate his age by way of a wrinkle proxy.

    • bender
      Posted Mar 14, 2009 at 1:43 PM | Permalink

      Re: Kenneth Fritsch (#102),
      No uncertainty whatosever when it comes to religious iconery. Like the original hockey stick of MBH98.

Hu McCulloch

RE Bender, #12, #97,

Thanks for the help. I’ve now looked up nlme (Non-Linear Mixed Effects), gls, and corAR1 in the 13.0 MB R PDF manual, and have an idea what they do.

As I understand it, you are feeding corAR1 the .3188 that was estimated from the OLS residuals, and then corAR1 constructs the correlation matrix from this, and gls( ) re-estimates the parameters with this correlation matrix. If the errors were truly AR1 with this correlation, the resulting standard errors will be valid, and the new parameter estimates will be efficient.

GLS is the efficient estimator if you know the correlation matrix, but if you’ve estimated it from the data, and so are doing “Feasible GLS”, authors like Hayashi point out that its properties are not so clear. GLS will give different coefficient estimates than OLS. A full GLS AR(1) ML estimation (which R’s gls( ) doesn’t seem to do, but EViews’ AR(1) feature does) would actually go one step further and re-estimate the AR1 coefficient simultaneously with the slope coefficients.

The approach I follow in the post is to just take the OLS estimates as given, and try to correct their se’s. This is the same approach Santer, Nychka et al (2008) and Nychka, Santer et al (2000) take. (The Newey-West HAC procedure, which is very popular in econometrics, likewise takes the OLS coefficients as given, and just tries to correct the se’s, but I argue that it under-adjusts, and that this procedure is therefore inadequate.)

When I was a young perfectionist, Feasible GLS seemed to me like the only way to go, but now that I’m older and mellower, I’m agreeing more with Hayashi that OLS point estimates have a nice solidity to them, even if the OLS se’s need major correcting.

• bender
  Posted Mar 14, 2009 at 1:39 PM | Permalink

  Re: Hu McCulloch (#103),

  A full GLS AR(1) ML estimation (which R’s gls( ) doesn’t seem to do, but EViews’ AR(1) feature does) would actually go one step further and re-estimate the AR1 coefficient simultaneously with the slope coefficients.

  Interesting. Can you provide the full Hayashi ref. and any pointers to software documentation? One of the things I like about R is nothing is proprietary, not even the source. Other software might be more powerful, but if it lacks the transparency, then that’s a definite downside. (Imagine using a package called, I dunno, “TeamStat” – where none of the decades-old Fortran source could ever be examined.)

Hu McCulloch

RE UC, #101,

Is ARMAX in an add-on MATLAB toolbox? I don’t find it in the online help with my configuration, so maybe I need to put in a requisition.

• UC
  Posted Mar 14, 2009 at 1:57 PM | Permalink

  Re: Hu McCulloch (#105),

  In System Identification Toolbox, I used Version 7.0 (R2007a)

Hu McCulloch

RE #5, 38, 66,

It turns out that the “Q-factor” (1-r1)/(1+r1) for approximately adjusting the effective sample size does go back to Maurice S. Bartlett (1910-2002), “Some Aspects of the Time-Correlation Problem in Regard to Tests of Significance,” J. Royal Statistical Society, 1935, pp. 536-543, at p. 537.

Bartlett in fact gave an even more general version that does not require the regressor to be perfectly or almost perfectly correlated (as is the case with a trend line regression). If r1′ is the first order serial correlation of the regressor, Bartlett gives the formula as (1 – r1×r1′)/(1 + r1×r1′). This is the correct limit of Quenouille’s even more general 1952 formula quoted in #38 above, in the case when both autocorrelations are AR(1). Bartlett in fact seems to take it for granted that all serial correlation is AR(1), although Quenouille does not.

Both Bartlett and Quenouille are interpreting r1 etc as the serial correlation of y rather than of the error terms in a regression of y on x, but under the null that the slope is 0,this are the same thing.

The 1966 WMO report by Mitchell et al, cited in #66 above, does give the “Q” formula (which perhaps should now be the “B” factor, for Bartlett?) on p. 71, as a factor for adjusting a test for equality of two means. Its use in climatology therefore goes far back before Santer et al (2008) or Nychka et al (2000). In fact, I’ve received a new e-mail from Doug N., who reports, “Well gosh, I always thought this was a standard formula everyone in times series knew!”

I’m guessing that the famous Bartlett Window (aka Filter, Kernel) is named for electrical engineer Albert C Bartlett (who wrote in 1920s), and not for this MS Bartlett.

Speaking of vintage literature on serial correlation, it turns out that exact p-values for the Durbin-Watson stat can easily be obtained with MATLAB function DWTEST. A routine by the same name exists in R. In large samples, r1 is approximately N(0, 1/n), but in small samples the exact distribution depends on the number of regressors and their own serial correlation. The traditional DW tables just give upper and lower bounds for the 5% 1-tailed critical values, corresponding to the best and worst case X matrix. DWTEST instead requires you to give it the actual X matrix, and then it computes the exact critical value, with no need to worry about the traditional “dl” and “du”.

Note, however, that the default in DWTEST is a 2-tailed test for either positive or negative SC, rather than the traditional 1-tailed test for positive SC. I think DW originally just used a 1-tailed test, because the felt that positive SC should be corrected for since it can overstate results, whereas negative SC can only understate results and therefore is not as big a concern. Today, some practitioners would instead jump at any opportunity to overstate results!

Unfortunately, DWTEST requires you to input the residuals and X, rather than just the DW-stat and X, so it can’t be used to construct a graph of critical values for a trendline regression as a function of sample size. If anyone knows of such a routine, please let me know!

I have some results with using DWTEST on annual vs monthly AWS and TIR trend regressions that I’ll be reporting on the “It’s a Mystery” thread.

• UC
  Posted Mar 15, 2009 at 9:11 AM | Permalink

  Re: Hu McCulloch (#109),

  I have Bartlett’s book Bartlett, M. S. An introduction to stochastic processes with special reference to methods and applications (1960) somewhere, used it before in another thread here. Will check if there’s something useful.

• Kenneth Fritsch
  Posted Mar 15, 2009 at 10:00 AM | Permalink

  Re: Hu McCulloch (#109),

  I have some results with using DWTEST on annual vs monthly AWS and TIR trend regressions that I’ll be reporting on the “It’s a Mystery” thread.

  I posted some DW test results at http://www.climateaudit.org/?p=5151 on the Deconstructing the Steig AWS Reconstruction thread at posts # 281 and #292 (for a GISS Antarctica comparison).

Hu McCulloch

RE Mann-Made Global Warming self-parody, #115:
Despite the error message, the upload went through and works OK. The archived URL is http://www.climateaudit.org/wp-content/uploads/2009/03/global_warm.gif, but be sure to give credit to http://www.meteo.psu.edu/~mann/Mann/research/research.html as its source if you use it.

• Ryan O
  Posted Mar 23, 2009 at 11:04 AM | Permalink

  Re: Hu McCulloch (#116), I note with irony that they have used an inverted cryosphere for their graphic.

Kenneth Fritsch

Hu, I did the dwtest(lmtest) in R for the TIR annual reconstruction 1957-2006 (trend slope = 0.11818 degrees C per decade) and obtained a DW stat = 1.5933 and a 2-tailed p-value = 0.05331.

The p-value differences between R (my rendition) and MATLAB are not trivial.

Hu McCulloch

RE #119 —
This is as good a place as any for this question, but I’m not sure I’ll have the answer. Perhaps others who have figured out CE and RE can help out.

My impression of CE and RE is that they were well-intentioned efforts by Fritts in his c. 1972 book to check the out-of-sample fit of a model. He invented these ad-hoc statistics by analogy to the in-sample R^2, but they don’t have a well-defined distribution and aren’t used outside of dendro so far as I know. It’s not clear to me what they mean even if you simulate their distribution with artificial data whose characteristics sort of match the real data. Even the verification R^2 (or r^2 ?) doesn’t have the clear interpretation of the calibration R^2.

So I’d generally put primary credence in the standard full-sample calibration statistics instead, such as t and F (adjusted, of course, for at least AR(1) serial correlation), rather than these ad hoc “verification” statistics. The calibration R^2 is just a transform of the “Regression F” test, that tests that all coefficients other than the intercept are 0.

One could then break the sample in half or thirds and do a Chow (switching regression) test if desired, just to double check that there is no conspicuous variation in the coefficients. This has a clear distribution (for pre-specified break points). Breaking in thirds is in the spirit of what MBH do, with their “verifications” on the first and last third of their data. Or, more ambitiously, one could do something like k-fold cross-validation, for which there is also an established theory.

I’ve done some playing with the Steig manned station data, together with the unrestricted covariances from his AVHRR matrix, but without PCA or RegEM, and have and have been getting an even stronger point estimate of the full-period slope for W. Ant. than he does. However, this is dominated by Byrd station in the heart of W. Ant, which was only active for the first 15 years or so. After it stops reporting reliably, the W. Ant. record interpolated from more distant stations like Adelaide-Rothera and McMurdo becomes pretty flat, even when the up-trending extreme Peninsula stations are included. I haven’t done any significance tests yet.

One reason I’m getting stronger point estimates of the early trend is that I’m using what I call “trending anomalies” rather than the usual zero-sum anomalies Steig presumably used. If every station had a full or almost full record, it wouldn’t make any difference. But when most have only fragmentary records, computing zero-sum anomalies for each could erase any trend signal. If you’re looking for a 50-year trend in the continent, you have to concede that there might be a 50-year trend in each of the individual station records. A station like Byrd, then, should be anomalized about its 50-year trendline (passing through 0 at mid-sample), rather than about its own mean. This creates a lot of extra measurement error in the anomalies, to be sure, but I haven’t tried to take this into account in the se’s yet.

Unlike Steig’s primary AVHRR-guided recon which takes no account of geography except to classify gridcells by region, I’m taking the relatiavae coordinates of the stations and the gridcells into account. I identify each station with the 50-km AVHRR gridcell it lies in. Then, using the comprehensive AVHRR covariance matrix, I can use the standard MV Gaussian conditional expectation formula to construct the expected temperature at each AVHRR gridcell, using the stations that are active at each point in time.

Doing this, I couldn’t help but notice that 6 of the extreme “Peninsula” stations turn out to be on the S. Sandwich Islands (most or all on K. George Island), which are separated from the Peninsula per se and the nearest AVHRR gridcell by a c. 100 km channel. I’m including these as representing their nearest gridcell (#1, with O’Higgins) despite the channel, but am excluding the 5 oceanic stations, which are at least 700 km from the nearest AVHRR gridcell and so can’t reasonably be identified with one.

Be all this as it may, you’ve done a great job replicating Steig with his own arcane methods, so by all means write up what you found, complete with CE and RE! I’ll think some about your question, but don’t guarantee any insights.

• Ryan O
  Posted Jul 26, 2009 at 8:11 AM | Permalink

  Re: Hu McCulloch (#120), I’ve mentioned before that I don’t like RE or CE – especially RE. Really they seem simply to be a test for stationarity. Be that as it may, the reason for using RE and CE in the response to Steig et al. is for comparison purposes. So I think RE and CE must be included, though I am certainly not opposed to showing more informative and definitive statistics.
  .
  By the way, I did a similar thing (but using the eigenvector coefficients rather than the raw covariance matrix) and calculated a reconstruction using regularized least squares to estimate the PC loadings. To do this, however, I first have to use the EM algorithms to estimate missing station data, so I’m very interested in what you get using Gaussian conditional expectation. I assume that you just use the station data as-is (without infilling) and estimate grid temperatures using whatever stations happen to be available for a particular time.
  .
  Along with that, I also directly used the AVHRR covariance matrix to obtain estimates for missing station data rather than using an EM algorithm. The difficulty I had was that this method was subject to extreme overfitting during periods where few stations were available. The results from experiments where I deliberately withheld data and tried to predict it were much poorer than the results from the EM algorithms. Because of that, I did not pursue this more direct method much further.
  .
  I should also mention that including or excluding the off-grid Peninsula stations does not have any material impact on the reconstruction using either an EM-estimation of the PCs or regularized least squares – and this goes for both trend magnitude and geographical location.
  .
  In the meantime, I would be appreciative of any insights you might have on the way Steig calculated his verification statistics (regardless of the particular statistic chosen). To me, it makes no sense to substitute red noise for the original data and compare that back to the original data. What would make sense is either to:
  .
  1. Substitute red noise for the original data and perform a reconstruction, and then compare the resulting reconstruction to the original data; or,
  .
  2. Substitute red noise for the PC loadings and compare the resulting reconstruction to the original data.
  .
  All Steig’s method accomplishes is to show that red noise doesn’t fit the original data very well. It cannot make any statement about the accuracy of the reconstruction because the reconstruction is neither used nor compared to the result.

Ryan O

BTW, Hu, of the two statistics (RE and CE), my opinion is that CE is the least crappy. It’s the same as average explained variance (R^2) with calibration period information censored. What I’ve been doing lately is simply withholding all the data from a particular station and calculating r, R^2 (though calling it CE), and rms error from the resulting reconstruction to the original data.
.
I forgot to mention that the regularized least squares method is identical to using the full AVHRR covariance if all AVHRR eigenvectors are included. I’ve done this, and the answer doesn’t change much once you get past ~13 eigenvectors. The difference, then, is how much the final answer depends on the EM algorithm infilling. That’s why I’m very interested in what you find.

Ryan O

Hu,
.
On the autocorrelation issue, the maps provided in the SI are not corrected, either. This may have been pointed out before. Anyway, just using the 95% CIs from the OLS standard errors, you get this:
.

.
Which exactly replicates the pattern shown in the SI.
.
After correction for serial correlation, the picture becomes this:
.

.
Script for calculating the CIs:

get.trends=function(dat, st, en, decade=TRUE, se=FALSE) {

vec=FALSE

### Set window
dat=window(dat, st, en)
f=frequency(dat)
dat=as.matrix(dat)
N=nrow(dat)
C=ncol(dat)
sum.f=colSums
if(C==1) {vec=TRUE; sum.f=sum}
dat=ts(dat, st, freq=f)

### Per decade?
dec=ifelse(decade==TRUE, 10, 1)

### Initialize storage
temp=matrix(nrow=1, ncol=ncol(dat))

### Get linear regression
fm=lm(dat~I(time(dat)))

if(se==TRUE) {

### Initialize variables
temp=matrix(nrow=3, ncol=ncol(dat))
I=seq(1:N)/frequency(dat)
rmI=ssq(I-mean(I))

### Calculate sum of squared errors
SSE=sum.f(fm[[2]]^2)/(N-2)

### Calculate OLS standard errors
seOLS=sqrt(SSE/rmI)*dec

### Calculate two-tailed 95% CI
ci95=seOLS*1.964

### Get ACFs
r1=vector()
for(i in 1:ncol(dat)) {r1[i]=acf(dat[, i], lag.max=1, plot=FALSE)$acf[[2]]}

### Calculate CIs with Quenouille (Santer) adjustment for autocorrelation
Q=(1-r1)/(1+r1)
ciQ=ci95/sqrt(Q)

### Save results
temp[2, ]=ci95
temp[3, ]=ciQ
}

### Save trends
if(vec) {temp[1, ]=fm[[1]][2]*dec} else {temp[1, ]=fm[[1]][2, ]*dec}
temp
}

Hu McCulloch

Several comments about the Steig Corrigendum that were formerly here have now been moved to a new thread on that topic.

Ryan O

Hu, Nic Lewis and I have been talking about DoF issues with respect to anomalies. Conversion to anomalies results in the loss of 12 DoF. Should not the trend CIs also be adjusted for these lost DoF?

• Pat Frank
  Posted Sep 14, 2009 at 6:28 PM | Permalink

  Re: Ryan O (#144), Ryan, I’ve been wondering whether temperature anomalies should be treated as statistical residuals, in terms of loss of DoF. After all, anomalies are supposed to reflect some deterministic process, such as trends in net temperature change. Therefore they should reflect deterministic processes to some degree and, according to theory, therefore show some autocorrelation. Reducing the DoF according to a statistical criterion implies that the autocorrelation is random; perhaps showing red noise. So applying a noise model to temperature anomalies implies the physical theory that temperatures have no deterministic trends.

  Of course, some of the autocorrelation may be due to random persistence, but some of it (an unknown amount) should be deterministic. So, it seems to me that the loss of DoF in temperature anomalies should be smaller than the prescription calculated for autocorrelation of strictly random residuals.

Hu McCulloch

RE Ryan O #144,

Conversion to anomalies results in the loss of 12 DoF. Should not the trend CIs also be adjusted for these lost DoF?

You’re quite right. In fact, 12 seasonals have been removed from the data, not just 1 mean, so the OLS se’s for the slope should be computed with n-13 DOF, not just n-2.

But with n = 600, the difference between these values is pretty small, so I just ignored it. If this is taken into account, the true se’s are a little larger than I indicated. Serial correlation makes a much bigger difference.

• Ryan O
  Posted Sep 15, 2009 at 7:02 AM | Permalink

  Re: Hu McCulloch (#146), Cool. While from a standpoint of 600 it may not make much difference, after correction for serial correlation, you end up with ~350 (even fewer for some regions with higher autocorrelation), it does affect things a bit. But mostly, Nic and I just want to make sure that the math is right. It’s not a hard thing to correct for the DoF.
  .
  The other thing is that this DoF correction should also be incorporated into the imputation algorithms. There it might make more of a difference, as you go through 100-1000 or more iterations. Small differences may compound.

Hu McCulloch

In the original post, I mentioned that I was working on a new approach to estimating ρ, which I now call the “Moment Ratio” estimator, since it matches r1, which is the ratio of two sample moments, to the corresponding ratio of population moments. The AR(1) case (which is all I have worked out so far) is treated in a working paper at http://www.econ.ohio-state.edu/jhm/papers/MomentRatioEstimator.pdf.

As I mentioned in the post, this refinement doesn’t make much difference for Steig’s trend regressions, because we are not very close to a unit root and the sample is pretty big. Here are the “AR” se’s from the post and the “MR” counterparts:


***** seAR seMR
All. .0458 .0461
Pen. .0259 .0260
West .0429 .0432
East .0504 .0507


So although in my mind “seMR” is the way to go, for purposes of discussing Steig, there is no point in introducing this complication. In fact, as I showed in the paper, “seQ” gives almost identical results to “seAR”, so if you just want to use the “Q” adjustment, you may as well go ahead. The climate people won’t give you any argument, in any event.

As noted in the post, I obtained seQ just by multiplying seOLS by sqrt(1+r1)/(1-r1)). This is equivalent to adjusting the DOF (n-k) by this factor, rather than n itself and then subtracting k. Quenouille and Bartlett actually did this differently, but Quenouille pointed out that since it is only a large sample approximation anyway, it doesn’t matter which interpretation you give the factor.

One thing B and Q didn’t discuss was whether the adjusted DOF should also be used to find the t critical value. This would make sense, but there isn’t any theory that I know of behind it. Here, even the adjusted DOF are so high that you’re off the end of the standard t-table, so you’re essentially using the normal critical value anyway.

Back to Ryan’s question in #144, then, to properly take into account that 12 seasonals have been removed from the data and not just 1 mean, the trendline should ideally be estimated in a regression with 12 seasonal dummies plus a time trend. OLS will then divide e’e by n-13 instead of n-2 when estimating the regression variance, and may also give a slightly different slope estimate than when the deseasonalized temperatures are used.

But the difference in slopes, if there is any, would be very slight, so it would be reasonable to just do OLS on a constant plus a trend with deseasonalized data, and then multiply those standard errors by sqrt((n-2)/(n-13)) to allow for the seasonals before multiplying again by sqrt((1+r1)/1-r1)) to allow for the serial correlation.

Hu McCulloch

RE Ryan #149,

It would seem appropriate, then, when estimating the sample covariance matrices that n-13 should be used instead of n-1.

The adjustment for a standard trendline regression (constant plus slope) would be n-2, not n-1. Often we don’t think about the constant, because with demeaned data and demeaned trend variable, its estimate always ends up being zero. Nevertheless, this is an uncertain “zero” because we didn’t really know what the population mean of the data was. Since we effectively estimated this with the sample mean, we have really estimated two parameters (constant and slope), not just one (slope). OLS won’t know this if you just give it the demeaned variables and ask for only the slope, but it will get it right if you ask for a constant and a slope.

With 12 seasonal dummies and a slope, OLS will correctly count 13 parameters and use n-13.

• Ryan O
  Posted Sep 16, 2009 at 10:32 AM | Permalink

  Re: Hu McCulloch (#150), Thanks. But I wasn’t referring to the regression, I was referring to calculation of sample covariance matrices – such as that in RegEM. Basically,

  cov(X) = [t(X) * X] / sqrt(n-1)

  If you have 12 missing DoF due to anomalies, then it would seem appropriate to estimate covariance by:

  cov(X) = [t(X) * X] / sqrt(n-13)

Hu McCulloch

RE Ryan #151,
But there are just 12 monthly means and no slope being estimated, the correction would be (n-12), not (n-13).

• Ryan O
  Posted Nov 12, 2009 at 11:14 AM | Permalink

  Re: Hu McCulloch (#152), Thanks. Duh. Nic Lewis corrected me on this as well.
  .
  BTW, I was wondering if you might be interested in reviewing and commenting on the paper once I have a draft completed. Right now, I’m at about 70% complete on the main paper and 20% complete on the Supplemental Information.

Dr Michael Koch

To Hu et multi al.,

thank you all for an extremely interesting and convincing reading. I want so sum it up (for a forthcoming book for Sweden-Norway-Finland-Denmark-Iceland and Germany-Switzerland-Austria on the incredible climate hoax) after telling the ‘McKitrick/McIntyre vs Mann’, the ‘Peiser vs Oreskes’ and the ‘Hu vs Steig’ story (because it is, again, so convincingly typical for deliberate and biased ‘junk science’ in this field) like this:

(Re: Antarctica HU vs Steig)
“Ironically, just this situation could by now be changing due to the paradox effect of the higher reflection from the large white Antarctica Polar ice and snow cap, as explained by Svensmark.* It has no corresponding area in the North Pole region, as Greenland is so much smaller. It is to consider if not this phenomenon alone could be responsible for what has been vaguely noted as the ‘North-South swinging’.

(*As Henrik Svensmark pointed out, Antarctic temperatures tend to go the opposite way to the rest of the world. His explanation is that clouds are less reflective than the Antarctic snow and ice, hence when albedo decreases globally and thus heats the planet, Antarctica becomes more reflective so cools down. And vice versa. Since global cloud cover has started increasing and we are now in a global cooling phase, we can expect Antarctica to become warmer.)

…It is also clear by now, that the time period the authors (Steig et al.) have chosen for their study includes many years of the cooling period 1940-1976, before the onset of the recent moderate Global Warming period, that of comprehensible reasons (Svensmark again) may have spared the Antarctica. What the authors have shown with low significance is thus, that it was a little colder in the 1950s than in 1990-2006 – which is a fact already well known from every temperature record on the surface of the earth. Their own data show, that – chosen a starting point between 1965 and 1980 – there was no significant warming of the Antarctica whatsoever: Kenneth Fritsch (27.2. 2009) checked the trends coming out of the Steig et al. data set over the periods 1957-2006, 1965-2007, 1975-2006 and 1980-2006 and could readily show that the published trends are statistically significant only for the 1957-2006 period. That means that coming forward only 8 years from 1957, there was no significant signs of Global warming any more. A truly sensational finding that the authors failed to mention.

“There also should have been a discussion of the ice core reconstructions in Steig et al. concerning the findings of West Antarctica, that for the 20th century the warmest decade was from 1935-1945…” (Kenneth Fritsch, 27.2. 2009), which neither were mentioned in their final conclusions, which thus omitted two more serious blows to the AGW theory from their summarizing comments.

That is exactly the one-eyed style climate ‘research’ that for decades has confused and bluffed the world. It was just a ‘search’ for some biasing verifications. (K R Popper would rotate in his grave.)

No newspaper in the world has ever (aside from the Hu McCulloch critic) made these findings known, in spite of their truly belonging to the ongoing AGW discussions that according to the IPCC and other recently praised Nobel laureates, such as Al Gore, Rajendra Pachaury and Barack Obama, are “over, finished, closed, the facts clear, the science settled.” They could have given this Nobel Prize even to Bert Bolin, Michael Mann and Jim Hansen, because they are the prior culprits by spreading their biased disinformation for so many years. They turned referenced science into main stream ‘referended’ bias.

‘Climate research’ was no real science under their auspices, but a politicized Great Carneval, for which the mentioned Nobellaureates acted Cheer Claque leaders.

But – how come?

My psychological explanation for how this could happen, put down as when I signed the ‘Minority Report’ of the US Senate:

“I am physician, epidemiologist and book author, especially interested in the history and the psychology of collective errors, which is also my motive for signing the Minority Report of the US Senate on the AGW climate hoax:

The whole incredible AGW climate issue (in its core a HGW issue (H for Heliogenic) and since 1998 turning into a HGC issue (C for Cooling) is the deplorable outcome of classical ‘group think’. Like so many other eclectic groups before them (such as the famous historical Academy of Science in Paris that found plague and tuberculosis “not infectious”), the scientists of the IPCC haven’t understood the clues of Karl R Popper, the philosopher who has shown that no hypothesis ever can be ‘proven’ by ‘verification’ alone. His convincing cloue: Thousands of positive confirmations cannot really confirm any truth. However, one single falsification can disprove the whole theory cogently.

Aware of this crucial insight the IPCC scientists should have actively searched for counterarguments, instead of disposing all those scientists that held contrarian views. Only by due diligence in incessantly disproving every single objection and counter-argument that may be found, their theory could have become more and more plausible and finally convincing. They chose the inverse strategy – a road astry: to scold and sever the sceptics and ‘deniers’. Thus, it was an indeed fatal blunder to publish a sea-level chapter unchanged in the IPCC report that by the world’s leading expert on sea-level measurement (NA Mörner) had been judged pure rubbish. It doesn’t help then to put his name on the report to assure scientific solidity – that’s just chutzpah.

‘Group think’ always tends to – in the course of time – slide on its own bias by selective omissions, arrogance, ignorance and intriguance into collective errors ad absurdum. To revise own convictions is, for sure, no favorite human pastime, but science has to resist the insideously comfortable trait of self-assurance and petitio principii often prevailing in selected groups. I’m sure, the AGW hoax will serve for generations after us as a warning example for this unlucky human shortcoming.”

For us witnesses of an important scientific debate it is, of course, crucial to hold it stringent and scientifically detailed see the oustanding, sagacious analysis of McCullough & McKitrick: Check the Numbers: The Case for Due Diligence in Policy Formation, Fraser Institute 2009, 43 pp (ISBN 1918-8223).

Who said: “The object of life is not to be necessarily on the side of the majority, but to escape from finding oneself in the ranks of the insane.”? Already Marcus Aurelius. [I chose this as the motto for my book, working title: The Power of Ideas over Reality: the ‘Climate Catastrophy’] So the problem is not a new one…

As I assume you all know McKitrick’s Taken by Storm, I only want to add two similarly urgendt recommendations: Nigel Calder: The Manic Sun (1997), Nigel Calder & Henrik Svensmark: The Chilling Stars.

But I would also like to mail to the leader of your discussion group three extremely informative texts of some of the most brilliant and experienced solarists in the world – if you are interested. In so case let me know, please – to which e-mail address?

Yours sincerely Michael

Hu McCulloch

RE Ryan O #157,
I’m not quite sure what you mean. Computing the full 5509×5509 covariance matrix and then extracing the 42×5509 submatrix you need should be equivalent to just computing that submatrix in the first place, without having to use up .2 GB or so of memory for the full matrix. AVHRR covariances between unoccupied gridsquares can easily be computed, but aren’t necessary for inferring unoccupied square temperature anomalies from occupied square anomalies.

In any event, Steig’s results are obviously wrong — given this data, there has to be more trend in the peninsula than in W. Ant., but he has it the other way around.

I ended up modifying Steig’s literal regional definitions a little to make them more reasonable — I included the tip of the peninsula and its entire west coast in “Peninsula” rather than cutting it off at 60 dW and above some latitude as Steig said he did. I also included a few degrees in all directions from the S pole in “East Ant”, rather than allowing “W Ant” to come right up to the pole as Steig said he did. I think I used 3 degrees, but maybe 5.

See the BAS map of “Antarctic Peninsula” at http://wattsupwiththat.com/2009/11/14/every-cloud-has-a-silver-lining-antarctica-glacier-retreat-creates-new-carbon-dioxide-store/ for a conventional definition. They go all the way down to 75S (not 72S like Steig), and then impose an 80W boundary. I just used Steig’s 72S, plus an eastern boundary somewhere to exclude the tip of Thurston Island.

If Steig really did impose the 65dS north boundary that he describes. that would have eliminated the strongly trending King George Island stations, and might account somehow for his weak Peninsula trend. (Just a guess.)

• Kenneth Fritsch
  Posted Nov 16, 2009 at 11:58 AM | Permalink

  Re: Hu McCulloch (#158),

  As I recall Steig admits that he gets the Peninsula trend wrong, but waves that off by stating it is small price to pay for getting the rest of Antarctica correct and besides we know that the Peninsula is heating up due to all the measurements we have of that area.

  The main thesis of the Steig paper, in my judgment, was that counter to treating the Peninsula of Antarctica as a special and isolated climate area, one could show that the warming of the Peninsula had spread into West (especially) and East Antarctica.

  Showing that the Peninsula is warming at a much faster rate than the remaining 95% of the Antarctica area and to even the exclusion of much warming in that larger area would not have made the cover of Nature. If the authors artifically spread some of that Peninsula warming to West Antarctica their conclusion as captured by Nature is misleading.

• Ryan O
  Posted Nov 16, 2009 at 12:23 PM | Permalink

  Re: Hu McCulloch (#158), My bad. If you do a truncated SVD on just the 42 locations you get different spatial weights than if you do a truncated SVD on the full 5509 locations and extract the 42 . . . but you’re not doing any truncation because you’re not doing SVD in the first place. I’m just so used to the SVD methods now . . . 🙂
  .
  I did masks that follow the traditional delineations if you find it helpful.
  .
  Peninsula: A line from Cape Adams to the mainland south of the Eklund Islands ( http://en.wikipedia.org/wiki/Antarctic_peninsula ).
  .
  West: West Antarctica is that portion of Antarctica lying on the Pacific Ocean side of the Transantarctic Mountains. ( http://en.wikipedia.org/wiki/West_Antarctica )
  .
  East is everything else. I also did a mask for just the Ross Sea. These are refined from what I did earlier. The full West Antarctic mask is simply “mask.west.a” plus “mask.ross” and east would be all grid cells not allocated in these masks.
  .
  mask.pen=c(1:114, 117:128, 140:152, 166:181, 195:208, 223:235, 250:260, 278:287, 310:317, 345:350, 383:386, 420:422, 457:459)
  .
  mask.west.a=c(115, 116, 129:139, 153:165, 182:194, 209:222, 236:249, 261:277, 288:309, 318:344, 351:382, 387:419, 423:456, 460:495, 497:535, 538:575, 578:615, 620:657, 663:700, 707:743, 752:787, 797:832, 843:878, 891:924, 938:970, 988:1020, 1043:1074, 1101:1131, 1163:1184, 1187:1192, 1227:1246, 1252:1256, 1292:1310, 1357:1374, 1424:1442, 1493:1505, 1560:1567, 1628:1634, 1696:1701, 1765:1770, 1835:1840, 1907:1912, 1978:1982, 2049:2053, 2120:2123, 2191:2194, 2262:2265, 2335:2337, 2407:2410, 2482:2483)
  .
  mask.ross=c(1185:1186, 1247:1251, 1311:1319, 1375:1383, 1443:1449, 1506:1516, 1568:1582, 1635:1649, 1702:1717, 1771:1785, 1841:1854, 1913:1925, 1983:1996, 2054:2066, 2124:2136, 2195:2206, 2266:2277, 2338:2348, 2411:2416, 2418:2419, 2484:2486)
  .
  mask.west=c(mask.west.a, mask.ross)