Climate Audit

Hu McCulloch

Nice job, Roman! Steig et al should take you on as a co-author on the revision, since they apparently will need you to explain to them what they did the first time around!

I’m afraid I don’t understand your first figure. If this is the difference between your reconstructed reconstruction and theirs, shouldn’t it just be rounding error? Why is this a “good fit”?

Do you have any idea yet where the 3 EOFs used before 1980 came from, or how the satellite data was used?

I’m betting the +20°C anomalies at the end of the first figure are for Racer Rock, #52. Am I right?

Nicolas Nierenberg

Based on answers from Dr. Schmidt I believe that the AWS data was used in the same manner as the satellite data. The idea was to fill in the spatial gaps between manned stations. I believe the two were used separately as a check, although they may have been combined in the final result. The point is that the trends should be driven from the manned station data with each manned station being weighted differently as to the amount of area it represents.

• RomanM
  Posted Feb 6, 2009 at 4:52 PM | Permalink

  Re: Nicolas Nierenberg (#5),

  Part of the reason I looked at pre-1980 was that neither the satellites nor the AWS could influence the variability of the reconstruction PCs. Although the level (the constant term in the reconstruction) would be influenced by the AWS and the relative weights for combined would be possibly influenced by the satellite data, the ups and downs would be solely due to manned surface stations.

Jeff Id

Very nice Roman, didn’t read the name. This is excellent.

JTaylor

I have not seen any mention of UHI adjustments for any stations in Antarctica. If they can detect UHI effects in Barrow Alaska with a population of 4000, there should be some measurable affect of 5000 researchers in Antarctica. Research vessels, cruise ships, aircraft, research villages and eco tourists have to throw off some heat.

“The number of people conducting and supporting scientific research and other work on the continent and its nearby islands varies from about 1,000 in winter to about 5,000 in the summer. Many of the stations are staffed year-round.”

Which ones are staffed year round, what is the poulation nearby? Any parking lots? 🙂

Scott Gibson

Re: ianl (#18),

Au Contraire (sp?), there is a scale bar :). It is scaled in polar coords, not linear.

jeff Id

Here we go.

Principal component analysis of the weather station data produces results similar to those of the satellite data analysis, yielding three separable principal components. We therefore used the RegEM algorithm with a cut-off parameter k=3. A disadvantage of excluding higher-order terms (k>3) is that this fails to fully capture the variance in the Antarctic Peninsula region.We accept this tradeoff because the Peninsula is already the best-observed region of the Antarctic.

John Baltutis

According to Tapio Schneider, his method* (Steig’s code) requires that the temperature data be Gaussian and that the missing temperature values are missing at random: “the EM algorithm and the methods that will be derived from it in subsequent sections are only applicable to datasets in which the missing values are missing at random.”. Has anyone verified and proven that the data from AWS is Gaussian or is it just assumed so that the mathematics/statistics described can be used in an analysis?

He also says that the data missing at random only “if correlations between anthropogenic temperature changes and the availability of data are negligible.” How are these temperature changes determined, how does one tack on the anthropogenic attribute to those changes, and how does one correlate those changes with data availability and what criteria determines if the correlations are negligible?
*

• RomanM
  Posted Feb 6, 2009 at 8:33 PM | Permalink

  Re: John Baltutis (#25),

  According to Tapio Schneider, his method* (Steig’s code) requires that the temperature data be Gaussian and that the missing temperature values are missing at random

  The mathematics of the methods can be applied to anything, even a sow’s ear. What fails, when assumptions such as normality or other such distributional considerations are not satisfied, is the interpretational aspect.

  This could lead to an introduction of bias, a systematic over- or under-estimation of parameters in those cases when a missing value is missing for reasons related to that value, e.g. a temperature sensor stuck in a deep snow drift for several months or longer. When the distribution of the data is substantially different from the assumed shape, it could lead to large changes in the interpretation of test statistics, as in deciding that a trend is (statistically) significant when in fact it is business as usual for data of that type.

  • John Baltutis
    Posted Feb 6, 2009 at 9:20 PM | Permalink

    Re: RomanM (#28),

    The mathematics of the methods can be applied to anything, even a sow’s ear. What fails, when assumptions such as normality or other such distributional considerations are not satisfied, is the interpretational aspect.

    Understood. So, has anyone, to your knowledge, proven the normality assumption? If not, then I consider that the exercise is fruitless.
    
    Steve: the exercise is not necessarily fruitless if the data is non-normal. I’d view that as secondary for now.

    • RomanM
      Posted Feb 6, 2009 at 10:02 PM | Permalink

      Re: John Baltutis (#34),

      You generally don’t “prove” that a variable has a particular distribution. Rather, you try to show that the data behaves in a manner consistent with having that distribution. In most cases, however, these types of assumptions don’t necessarily apply to the raw sample, but instead might be a property of a particular aspect of it.

      For example, in a regression you don’t expect the response, Y, variable to have a particular distribution since there is a dependance on the predictor x which may or may not be randomly chosen. Rather the “error term”, i.e. the distance from a linear relationship is assumed to be normal. In order to check that, one needs to estimate the regression line first thereby complicating things.

      Re: Bob North (#39),

      In a word, no. The complexity of the whole process in multivariate techniques in general often makes it impossible to devise suitable nonparametric techniques. As you mentioned, these are less powerful, but I see it as paying insurance, you give up a little of your information in order to allow you to test in situations where the parametric test is not tenable.

    • John Baltutis
      Posted Feb 6, 2009 at 10:31 PM | Permalink

      Re: RomanM (#40),

      Thanks! So where does Tapio Schneider’s method take us, if we make all of his assumptions? Is Steig responsible for determining the distribution, showing that it’s normal, and validating the randomness of missing data? AFAICT, that’s not addressed in his paper.

    • RomanM
      Posted Feb 7, 2009 at 7:39 AM | Permalink

      Re: John Baltutis (#33),

      I agree with Steve’s assessment of the situation with regard to normality. At this stage, I would be more concerned with systematic biases. Without knowing the details of how RegEM was applied to the data by the authors, it is difficult to be specific. For example, a positive trend could be induced into the reconstruction be simple underestimation by some constant valueat times where there is no AWS data. Since most of these times are in the early record, a regression would show a positive slope. Whether this might be the case here is speculation at this point.

      Re: Jeff C. (#41),

      Just to clarify a point, it is not necessary to think of pre and post-1979 for the PCs. I separated the data this way only for the purposes of the actual values of the principal components. Since the weights for combining for a particular AWS turned out to be the same in both periods, the PCS should be thought of as applying to the entire 50 year period.

      Did the actual measured AWS data compliment the existing trends or was it just so lightly weighted that it had minimal effect?

      This is not at all clear. At this stage key parts of the analysis are missing so we are limited to looking at what the reconstruction itself can tell us.

      Re: Geoff Sherrington (#48),

      Unfortunately, your graph is not particularly clear and lacks time scale so I can’t say anything specific about possibly added trends. However, it was striking similarities between reconstructions that led me to look at the data in the way I did. Since the manned surface stations were major players and since spatial distribution was important in determining the weighting of the PCs your observation is likely correct that the similarities are not coincidental is correct.

      I am having some trouble sympathizing with the “problem” of your 46 degree weather. This morning the thermometer outside my window read -21 C and last month we broke a 130 for the coldest day locally. On the average, however, for the two of us it is a balmier 12.5 C day. 😉

    • Jeff C.
      Posted Feb 7, 2009 at 12:05 PM | Permalink

      Re: RomanM (#50),

      Roman – thanks for your response. I really don’t understand how RegEM infills the data, but this is what’s bothering me.

      Since there is no AWS data before 1980, the assumption is the AWS reconstruction is based on manned station data up to this point (i.e. the AWS measured data isn’t extrapolated back in time). Your three PCs using the AWS recon data up to 1979 show they all have common roots.

      In the AWS-era, we have sporadic measured data points that are added to the mix. Steve has shown that the measured data is added directly to AWS recons without modification. When infilling missing data between measured data points, how heavily does RegEM value the measured AWS data?

      Let me give an extreme example to make a point. Suppose the manned station data used in the recon for an individual AWS station had a slope of +0.2 deg C/decade. Then we had measured data from an individual AWS from 1980 to 2003 but only for a single month of each year (say March) that had a slope of -0.2 deg C/decade. Since it was for the same month each year, it is probably a somewhat reliable indicator of the true temperature trend for that period. However, when infilling, the algorithm doesn’t recoginize the importance of all the data being from the same month in establishing a trend. It just recognizes that 11/12ths of the data is missing (months other than March) and discounts the importance of the measured data. The infilled data continues with the slope of around +0.2 deg C/decade.

      This is all speculation but it seems reasonable. I need to learn R to play around with it.

    • RomanM
      Posted Feb 7, 2009 at 12:56 PM | Permalink

      Re: Jeff C. (#61),

      I don’t claim to understand everything that they do but, by accident, I noticed that there are seasonal adjustments made to the data. When I was doing the reconstruction for 1980 forward, I first calculated the reconstruction values and then substituted the AWS data (not anomalies). When I plotted the difference between that result and Steig’s recon, I got graphs that looked like this:

      You will note the periodic differences indicating the seasonal adjustments. How this was done for sporadic single month observations is beyond me at this point.

    • Ryan O
      Posted Feb 7, 2009 at 1:16 PM | Permalink

      Re: RomanM (#64),
      .
      I thought that RegEM doesn’t substitute when actual values are present. Is that effect, then, the difference between the actual data as anomalies (Steig’s) vs. what the reconstruction would have predicted had no data been available? I may be behaving densely, but I’m not sure I follow why the difference is assumed to be a data adjustment???

    • RomanM
      Posted Feb 7, 2009 at 1:29 PM | Permalink

      Re: Ryan O (#66),

      I thought that RegEM doesn’t substitute when actual values are present.

      My graph doesn’t contradict that. My point with that graph is that it appears that some sort of anomaly calculation was done to the data before the substitution was made. In response to Jeff’s statement: “However, when infilling, the algorithm doesn’t recoginize the importance of all the data being from the same month in establishing a trend”, it would seem that the infilling was done using anomalies which inherently would be a recognition of the month in which the data was measured. What I don’t know is how the anomaly calculation would be done for extremely sparse data.

    • Ryan O
      Posted Feb 7, 2009 at 1:46 PM | Permalink

      Re: RomanM (#68),
      Ah! Got it. Thanks!

    • Jeff C.
      Posted Feb 7, 2009 at 1:33 PM | Permalink

      Re: RomanM (#64),

      This plot seems reasonable as the infilling source (presumably the manned station data) would have a seasonal variation contained in it.

      A different tack on the same theme, i.e. the measured data does not have a significant effect on the infilling in the recons. If you took the AWS recons over the full period of 1957 to 2006 and deleted the actual measured data points, would they still boil down to three PCs? If so, wouldn’t this indicate the measured AWS data has little influence on the infilling?

      I’ll play with this today. I’m getting the impression that AWS recons have been pre-determined by manned station data and the measured AWS data is window dressing.

    • Nic L
      Posted Feb 10, 2009 at 7:30 AM | Permalink

      Re: RomanM (#64),
      “You will note the periodic differences indicating the seasonal adjustments. How this was done for sporadic single month observations is beyond me at this point.”

      The seasonal adjustments to the AWS data, to convert them into temperature anomolies, seem to be done by subtracting from each month’s actual temperature the mean of all observations for that month. However, as per an earlier post I made, for a number of the AWS the figures in Steig’s AWS reconstruction differ from the thus calculated temperature anomolies in the case of one or two particular months each year (generally by the same amount for each year in which there is an observation for the affected month), and in some cases there is disagreement for all or most months. December is by far the most common month for differences to appear. In the case of some stations, the figures in the reconstruction differ from the calculated temperature anomolies for all months, usually by the same amount every year for each particular month. Very odd.

      What is also very strange is that using observations up to December 2006 to derive mean monthly temperatures does not always give the fewest differences between the temperature anomolies calculated from the observations and the figures in Steig’s reconstruction. In the case of 89266 Butler Island graphed by Roman, the reconstruction figures differ from the calculated anomolies for all months where the mean observed temperatures are calculated up to December 2006. However, they agree exactly for all but July and December data if the mean observed temperatures are calculated up to July 2008. For Limbert and Terra Nova Bay also, calculating mean observed temperatures through to July 2008 results in fewer months showing a difference between the reconstruction and the calculated temperature anomolies (using means to June 2008 is even better for Terra Nova Bay). This seems to suggest that Steig has, in some cases, derived figures in his AWS reconstruction making use of post 2006 data.

      I have also found some oddities in the figures not derived from observations. As Roman says, the actual AWS readings are simply spliced into the prediction series calculated as per the Principal Components Analysis. However, on my calculations there are a small number of points where there is no actual data but the figure in Steig’s reconstruction differs from that calculated using the PCA. These are (+ve being Steig’s figure higher):

      89266 Butler Island Dec 1986 +0.2237
      89324 Byrd Dec 1999 -1.1484
      89376 Gill Jul 2003 -0.8622
      89758 LGB10 Oct 1993 +0.1070 Nov 1993 -1.9814
      89568 LGB35 May 2003 -1.6159 Jun 2003 -1.1941 Jul 2003 +1.0424
      89811 Law Dome Summit Feb 2005 +0.8768 Apr 2005 -0.9409 Jun 2005 +0.1890
      Jul 2005 -3.1503 Aug 2005 -1.3090 Dec 2005 +0.0408
      89864 Manuela Jun 2003 -0.2159 Jul 2003 -1.0821
      89866 Marble Point Jul 2003 +0.9882
      89667 Pegasus North Jan 1991 -0.3045 Feb 1993 +0.1019

      It seems very odd that (assuming no error in my calculations) Steig’s AWS reconstruction should differ from that calculated using PCA in 19 of some 30,381 cases.

      PS Congratulations to Roman for his work: I had only managed to work rather clumsily through the PCA up to December 1994 before he cracked the whole period.

    • Dave Dardinger
      Posted Feb 10, 2009 at 11:29 AM | Permalink

      Re: Nic L (#168),

      If you’ll let me give you a Nic L’s worth of advice, start using the B-Quote quicktag when using quotes. This would turn your first paragraph into:

      “You will note the periodic differences indicating the seasonal adjustments. How this was done for sporadic single month observations is beyond me at this point.”

      instead of just:

      “You will note the periodic differences indicating the seasonal adjustments. How this was done for sporadic single month observations is beyond me at this point.”

      You just need to highlight what you want as a quote and then press the B-quote button. The B-quotes are much easier to pick out and this avoids confusion.

    • Jeff Id
      Posted Feb 10, 2009 at 11:47 AM | Permalink

      Re: Dave Dardinger (#172),

      The controls on this blog are great. It really makes posting easy.

      Re: Kenneth Fritsch (#171),
      Thanks, I learn more when lurking here than I can anywhere else on the internet. Roman’s post definitely connected a few PCA synapses for me.

      I’ve read that same paragraph a dozen times, there’s a lot of conclusion there and it seems to be based on the unseen satellite data processed by the unknown code.

    • Nic L
      Posted Feb 12, 2009 at 1:08 PM | Permalink

      Re: Nic L (#168),
      I have now discovered that all of the apparent differences between Steig’s reconstruction and values calculated using PCA that I reported in the above post appear to be due to actual data rather than infilled being used by Steig. Although data values for none these points are included in the current READER files (which is what I was working from), in Steve’s data file (downloaded earlier, and therefore more in line with what Steig must have used) they all have values. See my post #142 in the Carnage thread. When I changed to using temperature anomoliescalcualted from actual readings for these data points, the differences disappeared.

    • Jim Owen
      Posted Feb 7, 2009 at 9:24 PM | Permalink

      Re: Jeff C. (#61),

      Thank you – I’m almost ready to read what he wrote and chase some data. I spent some time today looking at background info (mainly history, photos, etc) – and thinking (remembering) the early 60’s as well as 40+ years of spacecraft work and how it related to this. And wondering why I suddenly smelled a rat in the woodpile earlier today. Or meybe – why I didn’t smell it earlier. We’ll see where it leads me. Usually into trouble. 🙂

  • Nic L
    Posted Feb 7, 2009 at 4:43 PM | Permalink

    Re: RomanM (#27),
    “According to Tapio Schneider, his method* (Steig’s code) requires that the temperature data be Gaussian and that the missing temperature values are missing at random”
    I done some very rudimentary investigation of the effect of missing data on the AWS temperature records, by comparing trends using all the data in the HTML pages with trends using only the data where over 90% of the possible observations were actually available (as per the text pages, and used by Steig). The simple average over all AWS of the temperature anomoly slope (in degrees C per decade) that I get is 0.15 with 90%+ observation data and 0.03 with all observation data, using data up to 2006 as used by Steig.
    I have also looked at the effect of including data up to January 2009 (where available) rather than just to December 2006, using 90%+ observation data. The overall average AWS temperature anomoly slope falls marginally, from 0.15 to 0.10. This equates to an average reduction in slope of 0.23 for the 14 AWS which had post-2006 data. Limbert was particularly affected, with its slope changing from +0.81 to -1.21 upon inclusion of post-2006 data. However, for the manned stations including post 2006 data made virtually no difference to their average slope.

• Gary A.
  Posted Feb 6, 2009 at 8:49 PM | Permalink

  Re: John Baltutis (#25),

  and that the missing temperature values are missing at random

  I’m not sure missing data is random. Seems like a lot of missing data would occur in seasons of extreme cold or large amounts of snow.

  • Lucy Skywalker
    Posted Feb 7, 2009 at 5:30 AM | Permalink

    Re: Gary A. (#29),

    Re: John Baltutis (#25),
    “and that the missing temperature values are missing at random”
    .
    I’m not sure missing data is random. Seems like a lot of missing data would occur in seasons of extreme cold or large amounts of snow.

    I’d second that concern, following my observations that in the records I’ve seen that show seasonal values, winter temperatures spike, whereas summer temperatures are far smoother; thus cold records over 10 deg C anomaly could get ejected where summer records don’t.

John Baltutis

Arrrrgggghh! Link doesn’t work properly. One more try:

* http://web.gps.caltech.edu/~tapio/papers/imputation.pdf

Vernon

RomanM, nice work! There was one study that found that RegEM had a warming bias.

Alan Wilkinson

Global warming getting too hot here and I can’t type anymore. Please ignore all my attempts to spell Schneider properly.

Jeff C.

I’m a grunt engineer, not a statistics whiz, but this really seems odd. If you compare Steve’s plot above (which I think uses Roman’s PCs from 1957 to 1979) to Steve’s plot in the next thread of Steig’s AWS reconstruction from 1979 to 2003, it appears that only five points out of the 63 change their sign. Based on the dot size, those that do change sign appear to be so close to zero that you can’t really say the trend changed much.

Despite all the infilling of actual measured AWS data from the last 20+ years into the reconstructions, there was no significant change in the pre-existing trends that existed prior to 1979? Did the actual measured AWS data compliment the existing trends or was it just so lightly weighted that it had minimal effect? When infilling missing data during the AWS-period, is the algorithm still dominated by pre-1979 trends? It probably depends on how much AWS data is missing.

Seems to me the AWS reconstructions are just an extension of the pre-1979 trends apparently derived from the manned station data. How this verifies the satellite reconstruction beats me.

Steve McIntyre

It’s twice daily information on different grids than used in Steig. Steig processed the data conditioning it for cloud cover in a different way than their predecessors. For the statistic analysis in Steig, they would have prepared a data set of monthly AVHRR results. Steig says that he provides this information to “legitimate” researchers, but refused to provide it to me.

• Jeff Alberts
  Posted Feb 7, 2009 at 3:11 PM | Permalink

  Re: Steve McIntyre (#44),

  You’re a published researcher. How can he claim you’re not “legitimate”?

Fabius Maximus

Thank you for the explanation.

This might be one cause of confusion. If I may restate your words, Steig has processed the AVHRR data — and it is this intermediate process that you would like to know. Or at least get the intermediate work result. The NOAA data he used as a starting point is available to the public. I believe this complaint has to be clearly expressed.

As an outsider to the debate, it is perhaps the key point — far more important IMO than the debate (where you were largely a by-stander) over precedence in discovery of the BAS data errors, and that unfortunately seems to have obscured more important issues.

Steig’s definition here of who has “legitimate” access to “data, including the intermediate steps” deserves discussion. He has defined who gets access to key information on which public policy will be made — policies that might determine the fate of the world (as Al Gore has told us). Apparently much of the climate science community agrees with him that only a few should be allowed to participate in this science.

This deserves its own thread, and might have far greater results Steig 2009 itself.

RomanM

Oops, it should read 130 year old record for the coldest day…

Jim Owen

I’m likely to reveal my ignorance here, but – so be it. I have a different question that, at least from my viewpoint, needs an answer. To wit – how the hell did he do it? Meaning specifically – Steig claims an AVHRR reconstruction starting in 1957 – BUT – the AVHRR data set doesn’t start until 1981. Where did those 14 years from 1957 to 1981 come from?

Understand that I’m not exactly a novice where satellite data is concerned. In 1964 I was doing video data analysis for the Nimbus series spaceraft – and there simply was NOT any higher resolution data available from any satellite program. And certainly nothing comparable to AVHRR in terms of resolution or accuracy.

So where did his AVHRR reconstruction data come from?

Steve: It’s extrapolated using surface data to “infill” prior AVHRR data.

Jim Owen

Steve said: It’s extrapolated using surface data to “infill” prior AVHRR data.

IOW – – snip: language not allowed

… I’m more than a little irritated.

In 1957 there was nothing even resembling AVHRR. The first available satellite data was from the Nimbus 1 HRIR instrument in late 1964 and there was less than 28 days of data (the spacecraft “died”). The next possible data set would be the Nimbus 2 through 4 HRIR and MRIR sets, for which – interestingly, the NSSDC has only a truncated data sets. NOAA has data for the period prior to 1979 but until the advent of the MSU in 1979 none of it had anywhere close to the resolution of the Nimbus data.

OK – I’m rambling – but what I’m leading up to is this – If I remember correctly, Steig said specifically that the AVHRR reconstruction used NO AWS data. But IMO, extension (extrapolation?) of the AVHRR algorithms 14 years or more into the past when only AWS (or other pre-AWS)data would have been available, fails to meet any (OK – “my”) standard of reasonability. And belies his statement that he used no AWS data.

Hmmm – I didn’t check on when the AWS data became available, did I? Shoulda done that. Was it available as early as 1957. If not – then what did he use for ground truth for the 1957-(AWS data start) time period? He apparently used the AWS (or other) stations to provide ground truth to calibrate the AVHRR data, so what was the reliability/accuracy/etc of those ground truth sources? Again – satellite data HAS TO BE tied to some source of ground truth or it’s simply not valid (like tree ring squences). So Stieg HAD to use the AWS (or other) station data as a basis for anything he did. Otherwise, the entire structure has no basis.

Not to mention – why did he pick 1957? Why not 1981 when the AVHRR dataset began? Or 1979 when my questions wouldn’t be a factor? Is the paper available anyplace online? The PDF link apparently doesn’t work anymore. – snip = Too many questions, too few answers.

Anyway – thank you. I think I’ll go look for some answers.

• Hugo M
  Posted Feb 7, 2009 at 11:55 AM | Permalink

  Re: Jim Owen (#58),

  Jim, the Steig paper is here: http://www.nature.com/nature/journal/v457/n7228/full/nature07669.html

  By the way, this “letter” to Nature appears very unusal, as it is accompanied by a “making of”, just as it were Hollywood script. It must have been “Nature”, who actually determined the headline: “Making the paper: Eric Steig”

  http://www.nature.com/nature/journal/v457/n7228/full/7228356a.html

  […] In 1999, Steig realized he could use the surface-temperature measurements gathered by satellites. But at that point these had collected only 18 years’ worth of data, so he faced a wait of almost a decade before he would have enough to produce meaningful results. In 2001, David Schneider, then a PhD student in Steig’s lab, began to dig into the satellite data that was available and analyse it more thoroughly. Satellites do not make direct measurements of surface temperatures, instead calculating them from measurements of the radiation emitted from snow, so the data come with uncertainties. “David’s work showed us the limitations of the satellite data, but also that it was more reliable than most people had assumed,” says Steig.

  By 2007, Steig finally had the 25 years’ worth of data he needed, and teamed up with Michael Mann at Pennsylvania State University, University Park. They applied an algorithm that allowed them to fill in missing data across the Antarctic continent as far back as 1957, revealing the warming in West Antarctica. Careful examination of climate modelling results suggested that this warming relates to the regional loss of sea ice over at least the past 25 years (see page 459). Their results were initially met with scepticism. “Everyone I talked to said ‘There’s no way that’s right.'” says Steig. But, with detailed explanations, he managed to change opinions. “Some then realized that they had data that supported us pretty strongly,” he says. […]

  And one of the goals of his trendy 2009 “letter” seems to be this:

  […] Now Steig is heading back out into the field as a member of a team aiming to drill an ice core to the bottom of the West Antarctic ice sheet. Other researchers will take temperature measurements in the borehole, providing a definitive test of Steig’s results — the cooler temperatures of the 1950s should still be detectable at depth.

  http://www.nature.com/nature/journal/v457/n7228/full/7228356a.html

  • Bill W
    Posted Feb 7, 2009 at 2:25 PM | Permalink

    Re: Hugo M (#60),

    David = Tapio (?)

• Kenneth Fritsch
  Posted Feb 7, 2009 at 1:05 PM | Permalink

  Re: Jim Owen (#58),

  OK – I’m rambling – but what I’m leading up to is this – If I remember correctly, Steig said specifically that the AVHRR reconstruction used NO AWS data. But IMO, extension (extrapolation?) of the AVHRR algorithms 14 years or more into the past when only AWS (or other pre-AWS)data would have been available, fails to meet any (OK – “my”) standard of reasonability. And belies his statement that he used no AWS data.

  Hmmm – I didn’t check on when the AWS data became available, did I? Shoulda done that. Was it available as early as 1957. If not – then what did he use for ground truth for the 1957-(AWS data start) time period? He apparently used the AWS (or other) stations to provide ground truth to calibrate the AVHRR data, so what was the reliability/accuracy/etc of those ground truth sources? Again – satellite data HAS TO BE tied to some source of ground truth or it’s simply not valid (like tree ring squences). So Stieg HAD to use the AWS (or other) station data as a basis for anything he did. Otherwise, the entire structure has no basis.

  It is my understanding that the TIR and AWS measurement data goes back to the early 1980s, but the manned station measurement data goes back (or some of the stations) to 1957. The separate TIR and AWS reconstructions use the overlap time period to separately calibrate and validate the AWS and TIR against the manned stations. The calibration includes relating the manned station data spatially to the AWS and TIR measurements. Using the calibration, including the spatial correlations, the manned stations measurements going back to 1957 are used to “fill in” the AWS stations and TIR grids back to 1957.

• davidc
  Posted Feb 7, 2009 at 3:20 PM | Permalink

  Re: Jim Owen (#58),

  Thanks for going over those dates. When I looked at the plot of their “recontructions” (sorry, I’ve lost the link, maybe someone can provide it) it appeared to me that any positive trend you might establish relied on data before 1980. If you take only data after 1980 it looks to me like a small negative trend.

  Whatever methodology they have used has to be different before and after the satellite data became available. It seems that they have made two very brave extrapolations, first spatially after 1980 then back in time to before 1980. Whatever the statistical details (even, say, simple linear regression for both steps) I think we need more than the author’s assurance that all is well.

  In my view the only evidence we have at this point suggests that the method is flawed: the trend changes when the methodolgy changes.

  • Jim Owen
    Posted Feb 7, 2009 at 10:01 PM | Permalink

    Re: davidc (#81),

    David –
    The first AWS was apparently installed in 1980. I believe the AVHRR data begins in the same general time frame.

    I’m not looking at this from the statistical POV, but rather from the historical/operational POV because there’s something here that just doesn’t fit. My next step is to look at what Steig & Co are saying, at what ‘s being batted around here and at what I know of ground truth testing, satellite data collection and processing, etc. If I come up with anything (and I usually do in this kind of situation), I’ll let y’all know.

    In the meantime, I wonder if Steig knows about kadabatic winds?

    I’ll leave you with this quote from the Wegman Report that apparently applies to the situation:

    Especially when massive amounts of public monies and human lives are at stake, academic work should have a more intense level of scrutiny and review. It is especially the case that authors of policy-related documents like the IPCC report, Climate Change 2001: The Scientific Basis, should not be the same people as those that constructed the academic papers.

    But then – I’m preaching to the choir. 🙂

theduke

Roman: it’s as if they are purposefully making the science impenetrable so that the findings can’t be challenged. Is that a fair summation?

• RomanM
  Posted Feb 7, 2009 at 1:17 PM | Permalink

  Re: theduke (#62),

  No, I would not ascribe such motives to the users of the methology. My personal experience in academia is that a lot of users of more complicated statistical methods don’t always understand the methods but have been taught by (non-statistician) teachers and colleagues that these are correct and appropriate without being taught enough to know when they shouldn’t be employed. The impenetrability is just a side-effect – a convenient protection against criticism of the results.

  • theduke
    Posted Feb 7, 2009 at 2:04 PM | Permalink

    Re: RomanM (#67),

    Thank you for the response. You’re right: I shouldn’t ascribe motives.

    But I was struck by what you wrote in the paragraph which now comes after my original post but at one point preceded it:

    Sounds simple, but there are a lot of pitfalls. RegEM is not particularly transparent in its workings. They are estimating values at 5509 locations for 600 monthly periods – more than 3,300,000 values, and those values are not only interrelated, but may require ad hoc fine tuning to make the procedure work. Calculating reliable error bars becomes a problem. Quantifying the relationship involves assumptions such as linearity which may or may not be true. The real question is just how much faith can one put into the validity of the results.

    What some climate scientists haven’t paid enough attention to is the KISS principle (no – not the one people usually state, but the one statistical consultants tell their clients: Keep It Statistically Simple). It is a much easier matter to convince someone who understands the situation that you have done the right thing.

    For myself and the other laymen out there with limited scientific and statistical knowledge, allow me to phrase the question this way: Rather than ask the pure scientific question, “What does the data say?,” does the methodology in the Steig paper and others like it ask “What does the data say if we assume this adjustment number is reliable and then massage the data with it?”

    One more question: was it a simple reading of the temperature data that led to the generally held belief that Antartica was cooling in contrast to the ROW?

    • RomanM
      Posted Feb 7, 2009 at 2:38 PM | Permalink

      Re: theduke (#73),
      The electrons in my original post #67 accidently disintegrated, but I managed to collect enough of them to repost it 😉

      does the methodology in the Steig paper and others like it ask “What does the data say if we assume this adjustment number is reliable and then massage the data with it?”

      One more question: was it a simple reading of the temperature data that led to the generally held belief that Antartica was cooling in contrast to the ROW?

      Don’t think of it as an “adjustment number”. Suppose you have a (simplistic) problem such as “Have you noticed, over a period of years,Site A’s temperatures have almost always been a lot closer (than site B’s) to the temperature at Site C. Site C’s thermometer is broken, but we know A and B. How can we guess C’c temperature?” I think it is obvious we can do better than a simple average of A and B as a guess. Multivariate methods have been designed to deal with the type of information available in the problem I describe to produce improved estimators of C’s temperature. it’s not really “adjustment numbers” but rather a technical implementation of the relationship between A, B and C.

      Is it a “simple reading of the temperature”? No. Have the statistical methods been done properly and appropriately? That’s what we would like to be convinced of by having the research materials placed into the public forum.

curious

re: Roman at 63 above. This is how it seemed to me but what is the point of drawing the satellite data set into the argument? Have I understood that the calibration of the satellite data set in the “present period” (ie where satellite records exist) is being done from a (continuous/gridded?) spatial representation derived from the physical stations? If so then to go back in time would this spatial representation not serve to provide the backcasting? Or are the satellite records being used to provide a description of how the spatial representation of the extrapoloated/interpolated physical data has changed over time and thus provide an additional correction factor for the backcasting? Or is the satellite data set in the present period being used to evaluate and validate/reject anomalous physical station records? This latter seems unlikey given the problems recently identified.

I hope that makes sense and apologies if this is explained elsewhere and I haven’t caught it – I haven’t been following blow by blow – please point to a ref. if available.

Re: East vs West Antarctica – I saw in one of the posts a top down (bottom up?) view of the continent showing all the weather station locations. What struck me was the physical landmass imbalance accross the approx 175-355deg (ie East/West) longitude. The albedo would seem to be skewed to be more reflective on the West half of the continent and I’d expect this to have an effect on the relative still surface air temps? Is there any point looking at the temp. records and investigating those with low wind speed?

As above, sorry if this is already covered elsewhere. Certainly support the KISS principle and fully acknowledge my need to study more!

• RomanM
  Posted Feb 7, 2009 at 2:16 PM | Permalink

  Re: curious (#69),

  There are several reconstructions in the paper. One is for the temeperatures at the locations occupied by the AWS. The other is at 5509 gridpoints across all of the continent. The continental reconstruction must use satellite data in the RegEM procedure. I don’t think any of us here know exactly the details of either of these reconstructions beyond the generalities I stated earlier since we don’t have any of the satellite data. I can’t answer in more detail as to how so I am concentrating on trying to take the result apart to see what makes it tick.

  Re: Jeff C. (#70),

  Jeff, I am not sure what you mean by “delete the measured data”. My reconstruction of the three PCs used only reconstructed values – no measured values were used directly. However, when I looked at the “corrected version of the data set yesterday, I was a bit surprised to see that there were some changes in the new PCs (which translates into changes in all of the other AWS besides the two series whose data was corrected).

  • Jeff C.
    Posted Feb 7, 2009 at 3:04 PM | Permalink

    Re: RomanM (#74),

    Roman, I apologize as I realize I’m not making myself clear. Let me try another explanation. I may be completely off base, but I’m suspecting the AWS reconstructions were generated using the following method:

    1) Manned station data alone from 1957 to 2006 used to generate 3 PCs (methodology TBD).
    2) 3 PCs used to generate AWS reconstructions from 1957 to 2006 for each AWS location.
    3) Where available, measured AWS data is plugged into recons replacing data generated using the 3 PCs.

    Steig’s seems to imply this on page 3 of the SI.

    This would explain why the AWS measured data is unchanged in the AWS reconstructions as it was added after the fact and was not manipulated by the algorithm. This would also explain why the Harry/Gill mashup had a significant effect on the reconstruction for that site as it would have had a large percentage of measured data plugged in due to the large number of measured data points. Most of the other sites have a small number of measured AWS data points and have little effect on the overall trend generated by the 3 PCs.

    By “delete the measured data” I mean the following:

    1) Start with the AWS recons for all 63 sites from 1957 to 2006.
    2) Wherever measured data exists for a given AWS location (as shown in the BAS records) replace the value in the recon with (null)
    3) Calculate the PCs on the AWS recons with the real, measured data missing (as produced in step 2)
    4) If only 3 PCs are non-zero, then the measured data had no significant effect on the reconstuction except where the measured data points were plugged in. If more than 3 PCs are non-zero, then the measured data (or something else) was affecting the reconstruction.

    I realize it is a lot of work. Steve mentioned he scraped the BAS data awhile back. Does anyone know if it is available in a format comparable to the AWS recons?

    • Jeff C.
      Posted Feb 7, 2009 at 3:22 PM | Permalink

      Re: Jeff C. (#77),

      “If only 3 PCs are non-zero, then the measured data had no significant effect on the reconstuction except where the measured data points were plugged in.”

      I meant if only three eigenvalues are non-zero. College stats was 25 plus years ago for me. Senior moment.

Steve McIntyre

#77, Jeff, I’ve been experimenting a little with this. Roman’s technique, as you observed, is really very brilliant as, among other things, it yields the “fitted” values the RegEM output writes over.

By recovering the “fitted” values that RegEM overwrites, it’s possible to develop some comprehensible statistics on the fitted portion. The idea behind RegEM for fitting things isn’t crazy; but the lack of proper exposition of why things are doing what they do is disquieting. I know that all they’re interested in is the “answer”, but some people actually care about how you get the answer.

curious

Re: Roman at 87. Sorry if this is nitpicking and/or showing ignorance but I don’t think this works:

In an extreme case, for illustration, suppose A and B stations always read the same at the same time and I can only see A. Predict B at time x in the future. No amount of looking at A is of any use. However, if I want to know what B did x amount of time in the past, I can just check A’s records for an exact result.

You still need to know about the relationship between A and B in order to draw conclusions about past B from looking at past A. The fact they have tracked each other for a period of time does not mean they always did, or as you say, always will. Lets say A was installed 10 years before B (which could be more or less remote), then B is installed and then they track each other for ten years – would n’t it be stretching it a bit to say you can give an exact value for B say 15 years ago? Surely there would still be a qualifying level of confidence associated to this value – this is what you appear to be saying in reply 84? And wouldn’t this level of confidence be influenced by the amount of knowledge of the physical situation/relation as well as the numerical record? The ten year equivalence period could just be a transitory phenomenon? This is what I was trying to understand with the satellite record in comment 69 (thanks re: response) and seems to be what Douglas in 85 is querying?

Or are saying, for example, that they were both installed 20 years ago and they have tracked each other since, then you lose contact with B and check A’s record to find out what B said 15years ago? This is ok but only possible because you’d already had the evidence of direct equivalence?

Sorry if this is factor 10 dim – feel as if I’m showing a lack of stats. knowledge!

curious

re: me at 92 above. Sorry Roman – just reread your last line of 87. I think we are in agreement and my query is is whether this relationship over time can be/is being clarified by use of the satellite records to calibrate the surface station records.

Geoff Sherrington

You might have these figures, but Australian data gathering in the Antarctic can be tabulated approximately as:

ID (Aust) NAME START END LAT LONG

300000 DAVIS Jan-57 -68.577 77.9725
300001 MAWSON Jan-54 -67.6014 62.8731
300002 TAYLOR ROOKERY Jan-57 Dec-58 -67.45 60.8667
300003 WILKES Jan-59 Dec-69 -66.25 110.5833
300006 CASEY (TUNNEL) Jan-69 Jan-89 -66.2833 110.5333
300007 MOUNT CRESSWELL Jan-72 Feb-03 -72.7333 64.3833
300008 MAWSON (MOORE) Jan-72 Dec-84 -70.3 65.1
300009 MAWSON (KNUCKEY) Jan-75 Dec-84 -67.8 53.5
300010 MAWSON (MT KING) Jan-75 Dec-84 -67.1 52.5
300011 CASEY (LANYON) Jan-83 Dec-87 -66.3 110.8667
300015 LAW BASE Jan-87 Feb-92 -69.4167 76.5
300016 PRINCE CHARLES Jan-88 Dec-92 -70.2333 65.85
300017 CASEY Feb-89 -66.2792 110.5356

Kenneth Fritsch

I much enjoy these analyses/audits/deconstructions of papers that are performed here at CA and that has been particularly so with Steig et al. (2009) even when realizing that the analyses are far from exhaustive. Having said that, I judge the claims in Steig to be an exercise in overreaching and under-detailing based on the following results of a simple-minded analysis of my own – and with R script.

Below I provide R script for regressions of the monthly temperature anomalies from the Steig reconstruction with AWS (downloaded at approximately 10:00 AM February 7, 2009) for the time periods 1957-2006 and 1965-2006. Below I report the summary of those time period result along with the regressions from 1970-2006 and 1980-2006 for which the R script is not posted since the procedure can be readily adapted from the first two scripts. 1957-2006 regression is adjusted for lag 1 auto correlation of the regression residuals while the other time period regressions are not. The p value below is the probability of the trend not being 0 by chance. TS is the trend slope (in degress C per decade) and CI is the 95% confidence intervals for the trend slope.

1957-2006: Adj R^2= 0.02; p= 0.00; TS = 0.138, CI = 0.046 to 0.230
1965-2006: Adj R^2= 0.00; p= 0.30; TS = 0.052, CI = -0.046 to 0.149
1970-2006: Adj R^2= 0.00; p= 0.93; TS = 0.006 CI = -0.120 to 0.131
1980-2006: Adj R^2= 0.00; p= 0.54; TS = -0.057, CI = -0.239 to 0.126

What these calculations show me is that given accurate measurements, adjustments and RegEM procedures in Steig, the conclusion of a warming Antarctica can be made only when going back to 1957 and that going forward in the starting time of the series (and not very far) the warming/cooling trends can be shown to be statistically very uncertain.

I think the combining these calculations with the following statement from RC about 1935-1945 being the warmest decade of the 20th century makes any conclusions about long term trends in the Antarctica and potential effects of AGW even more tenuous.

“In our own published work to date (Schneider and Steig, PNAS), we find that the 1940s [edit for clarity: the 1935-1945 decade] were the warmest decade of the 20th century in West Antarctica, due to an exceptionally large warming of the tropical Pacific at that time.”

http://www.realclimate.org/index.php/archives/2009/01/state-of-antarctica-red-or-blue/

Download file from Steig site and calculate 1957-2006 monthly trends AWS:

>download.file(“http://faculty.washington.edu/steig/nature09data/data/ant_recon_aws.txt”,”temp.dat”)
> aws=read.table(“temp.dat”) # Probably not needed.
> mataws=as.matrix(aws)
> y=rowMeans(mataws) #Calculate the monthly averages for all 63 stations.
> x=1:600
> MonthAveAws=cbind(x,y) #Not needed for regression but for later use.
> lmaws=lm(y~x)
> summary(lmaws) # Gives regression results in degrees C per month – see conversion below to degrees C per decade.

Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-6.43293 -0.76924 -0.01838 0.82107 4.73797
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.5174155 0.1101518 -4.697 3.27e-06 ***
x 0.0011498 0.0003176 3.620 0.000319 ***
—
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1
Residual standard error: 1.347 on 598 degrees of freedom
Multiple R-squared: 0.02145, Adjusted R-squared: 0.01981
F-statistic: 13.11 on 1 and 598 DF, p-value: 0.0003191

> ResidualAws=resid(lmaws) # Calculating the regression residuals and offsetting by one month to do lag1 correlation for determining adjustment to slope Stdev CIs.
> MatResidAws=as.matrix(ResidualAws)
> Lag1=MatResidAws[c(2:600),]
>UnLagged=MatResidAws[c(1:599),]
> CorrMat=as.matrix(cbind(Lag1,UnLagged))
> Lag1Corr=cor(CorrMat)
> Lag1Corr
Lag1 UnLagged
Lag1 1.0000000 0.2022246
UnLagged 0.2022246 1.0000000
> Factor=(1+0.2022)/(1-0.2022) # Calculating the adjustment to the CIs for the trend slope using the procedure from Nychka et al., 2000 and as applied in Santer et al. (2008).
> Factor
[1] 1.506894
> n=600
> AdjStdevSlope=((n-2)/((n/Factor)-2))^(1/2)
> AdjStdevSlope
[1] 1.228598
>SlopeCperDecade=0.0011498*12*10
>SlopeCperDecade=0.13798
>StdevSlope=0.0003176*12*10
> StdevSlope
[1] 0.038112
AdjLower95=SlopeCperDecade-1.96*AdjStdevSlope*StdevSlope
> AdjLower95
[1] 0.04620032
AdjHigher95= SlopeCperDecade+1.96*AdjStdevSlope*StdevSlope
> AdjHigher95
[1] 0.2297517

Calculate 1965-2006 monthly trends AWS:

> x=96:600
> mataws1965=mataws[c(96:600),]
> y=rowMeans(mataws1965)
> MonthAve1965=cbind(x,y)
>lmaws1965=lm(y~x)
> summary(lmaws1965)
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-6.54779 -0.77952 -0.01415 0.80131 4.77842
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.2156203 0.1562245 -1.38 0.168
x 0.0004308 0.0004141 1.04 0.299
Residual standard error: 1.356 on 503 degrees of freedom
Multiple R-squared: 0.002147, Adjusted R-squared: 0.0001635
F-statistic: 1.082 on 1 and 503 DF, p-value: 0.2987
> Trend = 0.0004308*12*10 #Trend in degrees C per decade
> Trend
[1] 0.051696
StdevTrend =0.0004141*12*10
> StdevTrend
[1] 0.049692
> Lower95Trend=0.051696-1.96*0.049692
> Lower95Trend
[1] -0.04570032
>Upper95Trend=0.051696+1.96*.049692
Upper95Trend
[1] 0.14909

• RomanM
  Posted Feb 8, 2009 at 12:52 PM | Permalink

  Re: Kenneth Fritsch (#105),

  I hesitated to offer improvements for fear of discouraging Kenneth on his first effort, but hey…

  You can avoid all of this:

  > ResidualAws=resid(lmaws) # Calculating the regression residuals and offsetting by one month to do lag1 correlation for determining adjustment to slope Stdev CIs.
  > MatResidAws=as.matrix(ResidualAws)
  > Lag1=MatResidAws[c(2:600),]
  >UnLagged=MatResidAws[c(1:599),]
  > CorrMat=as.matrix(cbind(Lag1,UnLagged))
  > Lag1Corr=cor(CorrMat)
  > Lag1Corr
  Lag1 UnLagged
  Lag1 1.0000000 0.2022246
  UnLagged 0.2022246 1.0000000

  with this

  acf(residuals(lmaws))$acf[2] #[1] 0.2018417

  I’m not sure why it is a bit different. Maybe the data has changed on the website again.

  The value can then be incorporated into the script that follows without having to retype it.

  Check out ?lag and ?acf for more info on autocorrelation and lagging of sequences.

  • Kenneth Fritsch
    Posted Feb 8, 2009 at 1:13 PM | Permalink

    Re: RomanM (#107),

    Thanks, Roman for your help – again. I know that I am often too hesitant to ask for help, but solicited or unsolicited help with R is always welcome. I can already see where R is much more efficient than Excel (and I think I know Excel well)- even with my clumsy programming.

    That ?command asking R for help has worked well for me as has searching the internet. Your example points to just how efficient R can be when one learns most of the best commands.

    Thanks, again.

  • Alan Wilkinson
    Posted Feb 8, 2009 at 9:47 PM | Permalink

    Re: RomanM (#107),

    The two acf calculations do give different results and not because the data has changed on the website. I reproduced both Kenneth and Roman’s values so there is either an algorithmic or precision difference in the two methods.

    As a complete R novice that’s all I can contribute.

    • Kenneth Fritsch
      Posted Feb 9, 2009 at 8:45 AM | Permalink

      Re: Alan Wilkinson (#137),

      The two acf calculations do give different results and not because the data has changed on the website. I reproduced both Kenneth and Roman’s values so there is either an algorithmic or precision difference in the two methods.

      As a complete R novice that’s all I can contribute.

      Alan, as an incomplete R novice, I agree that the difference is not due to the input data changing, but not to worry, as I am going with Roman’s procedure in the future.

    • RomanM
      Posted Feb 9, 2009 at 2:57 PM | Permalink

      Re: Alan Wilkinson (#137),

      I did some experimentation by exporting Ken’s data to a statistics program called Minitab. Interestingly enough, it gave me the identical two results that R did. There must be a difference in the way the correlations are calculated for ACF. Offhand, I don’t know what that difference is.

      Re: bernie (#150),

      Because of the complexity of the procedure, who knows what happens with Regem. Rotating the PCs does not affect the reconstructions, but it changes the three (in this case) pieces that are added together in the reconstructions. If a rotation would cause two of the coefficients for say, for a given AWS, to be close to zero, then the remaining PC has the same pattern of variation as the AWS. If other AWS’s have a similar configuration of coefficients for that rotation, one might consider that PC as identifying the group (e.g. close in distance or from a particular region). I doubt that it would pick out particular individual AWS’s.

    • Alan Wilkinson
      Posted Feb 9, 2009 at 6:24 PM | Permalink

      Re: RomanM (#155),

      I wondered (lying awake in bed – shouldn’t do that!) if the difference is because the acf calculation is doing a multivariate correlation whereas Kenneth was looking at a single lag? So the single variable regression soaks up some of the extra trend that gets allocated to higher order lags in the acf calculation?

    • RomanM
      Posted Feb 9, 2009 at 7:47 PM | Permalink

      Re: Alan Wilkinson (#156),
      Re: Bernie (#157),

      Hey! Don’t badmouth Minitab. It’s in version 15 and pretty good now – a lot easier to teach in an elementary course.

      Alan, take a sleeping pill. No, a correlation is a correlation. Something is different in the calculation, but I’m not sure what since I haven’t worried about how modern programs do these things. Other people have already done the worrying. There might be some sort of adjustment because the “sample sizes” are different for the various lags: n-1 for lag 1, n-2 for lag 2, etc.

    • RomanM
      Posted Feb 9, 2009 at 8:38 PM | Permalink

      Re: RomanM (#158),

      Final comment on the acf correlation. Got it figured out. When Ken did the correlation by hand, for the denominator, his formula used the square root of the product of the variances of the shortened (599 point)original sequence and the variance of the (599 point) lagged sequence. The formula for acf replaces each of these with variance of the complete 600 point sequence. That is the difference.

    • bender
      Posted Feb 9, 2009 at 8:54 PM | Permalink

      Re: RomanM (#161),
      Well done!

    • RomanM
      Posted Feb 9, 2009 at 9:05 PM | Permalink

      Re: bender (#162),

      Got the clue from here. 🙂

Kenneth Fritsch

acf(residuals(lmaws))$acf[2] #[1] 0.2018417

And when one leaves off the $acf[n]in the above command one obtains a neat/cool graph showing all the lag correlations from 0 through 25. Lag12 comes in second highest to lag1 at something like 0.14, as I recall.

• jeff Id
  Posted Feb 8, 2009 at 3:00 PM | Permalink

  Re: Kenneth Fritsch (#109),

  Cool..

Jeff C.

Recall that above Roman generated his 63 estimated 1980-2006 AWS reconstructions using a few series that had overlaps of no measured station data. When he subtracted his estimated reconstruction from the Steig aws reconstruction, he got the first chart in his post. He speculated that the spikes were from the actual measured data that was spliced in the reconstructions.

I took Roman’s estimated recon and replaced elements in the time series with the corresponding element from Steig’s recon ONLY at points where the BAS data shows that real, measured data exists (from Steve’s data.tab file). I then subtracted the two and got this, the difference is virtually zero.

[IMG]http://i404.photobucket.com/albums/pp127/jeffc1728/diff_plot.jpg[/IMG]

Note the scale on the left. The biggest difference between the two is 5.189498e-07.

What does this mean? Here is my take, some of this I said earlier in the thread, but now I am much more confident.

1) Roman’s methodolgy nailed the basis of the aws recons.
2) The AWS recons aren’t really reconstructions (to my way of thinking at least). They are products of the 3 PCs with the real, measured data sprinkled in almost as an afterthought. In my opinion, a true AWS recon would have used to measured AWS data as part of the infilling algorithm. It doesn’t use it at all. Regardless of the trend in the measured AWS data, the infilling is solely based on the 3 PCs.
3) On page of the SI, Steig says the same methodology was used for the AWS data as the satellite data. My strong suspicion is that he used these same 3 PCs for the satellite recons, and then spliced in satellite measurements that passed some cloud masking test. Presumably the AWS recon coeficients for each AWS site are determined from on lat/long or the distance/direction to some fixed position. The satellite recon would use the same method to determine their coeficients for each gridcell.
4) Steig talks of how the satellite recons and AWS recons “agree very well”. Not surpising if the roots of both are the same 3 PCs. Gavin stated that fixing the Harry/Gill mixup made the satellite and AWS agreement improve. Makes sense. Although Harry/Gill was corrupted, it had more independent data (data from other than the 3 PCs). Once Harry was fixed, more of its data came from the 3 PCs than before, hence better agreement.
5) It all comes down to the source of the 3 PCs. Steig seems to say it comes from the manned station data on page 3 of the SI. In comment 14 above, Steve notes that PC2 causes the “contrast” between West and East Antarctica. (Note, if running Steve’s code in #14, comment out line 4 (x=rec.svd$v[,3]) or it gives you PC3 instead.) If we can figure out how the manned station data was used to get these, this thing breaks wide open.

I’ll put the code in a follow-up comment.

Steve McIntyre

Roman’s decomposition of the AWS recon is very cool. Maybe some of the people who are complaining about the blog could spend some time with this analysis.

Here are a few supplementary notes on what the RegEM “is” in terms understood off the Island.

Roman has shown that the RegEM reconstruction can be written as follows”

where k =3.

is orthogonal in the pre-AWS period (and “near-orthogonal” overall. isn’t. That the fit has 3 degrees of freedom is highly relevant.

Then they (in effect) substitute actual values for fitted values, making a splice of fitted and real data.

By overwriting fitted data with actual data, they get a “perfect” fit to observations and thereby make it impossible to evaluate things.

But Roman’s methods enable the fitted data to be reverse engineered (in this case) and the calculation of things that you expect to see in statistical analyses – like residuals. Things that are never mentioned in Mannian studies.

I’ve also applied Roman’s method to pre- and post-Harry results and will report on this tomorrow probably (unless, of course, Gavin beats me to it.)

A question here: over at RC, an AWS reconstruction is in this figure:

I haven’t seen any information on how they go from the AWS reconstruction of the 63 stations to an Antarctic recon. Has anyone seen any info?

• RomanM
  Posted Feb 8, 2009 at 6:50 PM | Permalink

  Re: Steve McIntyre (#113),

  Yes, being able to get “under” the replaced values was in fact the part I was happiest with in the analysis. The can of course be orthogonalized in the entire period which is what I did in #55 although that is not crucial to the reconstruction. It can help to do so when looking at averages of trebds from the different reconstructed series.

  Re: Jeff C. (#111),

  I agree with most of your comments hoowever, I have the feeling that the aws reconstruction itself does use the observed data beyond the simple replacement. I may not have dug deep enough, but the all the reconstructions appeared to change somewhat in the “corrected” data set posted on Steig’s website.

  I spent time today looking at ant_recon.txt. I am happy to report that in fact the 5509 series in that file are ALL reconstructed from a single set of 3 PCs. With only 600 values in each series there is a heck of a lot of colinearity in these series. I ended up having to use svd since princomp (which gives more statistical content in the output) refuses to do the calculations when the number of variables exceeds the length of the series. I am not quite finished with the analysis and I will post relevant scripts, probably tomorrow afternoon. The data set is large. With all of the 5509 complete regressions. etc, I ended up saving a 250 MB RData workspace! I’d like to see Ken do that in Excel.

  By the way, the few places where I compared the sat recon to the AWS recon did not look all that close to me. Of course, I didn’t smooth everything first with a 40 year filter.

  More tomorrow on this.

  • Jeff C.
    Posted Feb 8, 2009 at 8:37 PM | Permalink

    Re: RomanM (#114),

    Roman – yes it does seem as something is going on there. In the SI on page 3 Steig says:

    We apply the same method as for the TIR-based reconstruction to the AWS-data, using RegEM with k = 3 and using the READER occupied weather station temperature data to produce a reconstruction of temperature at each AWS site for the period 1957-2006

    This is an interesting sentence as it seems to say that the satellite and AWS reconstructions use the same methodology and that the basis for the recons (your 3 PCs) are the “READER occupied weather station temperature data” (i.e. the manned station data, not the AWS data).

    However, up in comment #95 Hugo M notes:

    […] It’s not like the AWS data is all that informative, it only contains wind speed, wind direction, temperature and pressure, but for the last 20 years, it’s the only recording of winter weather that we have, so it’s important to continue it so we can compare it with the new acquisition systems that will start with the opening of Concordia. […]

    The thinking being that many of the manned stations are deserted during the Winter. Perhaps what Steig refers to above as “READER occupied weather station temperature data” actually contains some AWS data to supplement the winter months.

    Thus the 3 PCs are affected by changes in some AWS data series resulting in a small change in ALL of the recons. However, changes in a specific AWS series does not affect ONLY its infilling as that is governed strictly by the 3 PCs.

    Hope that makes sense. When I replaced the the elements in your recons with the actual measured data, the difference between it and Steig’s recons was virtually zero (5.18e-07 worst case for ANY point in any series.) If the measured data from individual AWS sites was used to in its reconstruction (other than affecting the 3 PCs), there would have been non-zero outliers.

    Things sure move fast around here.

    • RomanM
      Posted Feb 8, 2009 at 9:10 PM | Permalink

      Re: Jeff C. (#127),

      Thus the 3 PCs are affected by changes in some AWS data series resulting in a small change in ALL of the recons.

      Yes, what surprised me was that corrections to Harry (and presumably Racer) produced differences in all of the AWS recons before 1980 when the only existing information was the surface stations. In RegEM every variable (surface or aws or satellite) can have an effect on every other variable at every time in the infilling process.

      Re: Steve McIntyre (#128),

      That’s odd. The large blue dot overlapping with a somewhat smaller red dot in the upper left of PC1 is strange. I realize that the trends for the two could still both have the same sign, but the pattern involving the three PCs is different. If the PCs were physically meaningful, there would have to be some explanation for this.

    • Kenneth Fritsch
      Posted Feb 8, 2009 at 9:51 PM | Permalink

      Re: RomanM (#133),

      In RegEM every variable (surface or aws or satellite) can have an effect on every other variable at every time in the infilling

      The authors of Steig state that the AWS and TIR reconstructions are totally independent and therefore the interactions should only be AWS stations with manned stations and TIR grid data with manned stations and never AWS stations with TIR grid data.

      Is that your understanding also.

    • RomanM
      Posted Feb 9, 2009 at 6:31 AM | Permalink

      Re: Kenneth Fritsch (#139),

      My statement applies to all of the variables that are used in a particular analysis. From the Steig SI pdf:

      We apply the same method as for the TIR-based reconstruction to the AWS-data, using RegEM with k = 3 and using the READER occupied weather station temperature data to produce a reconstruction of temperature at each AWS site for the period 1957-2006.

      I would interpret that as saying that the AWS and satellites were not usen in the same analysis.

      Re: Kenneth Fritsch (#136),

      I could do it in one line even incorporating unequal weights for the months (which I didn’t do in my quick analysis of the graph that Steve’s post contained). When I programm in R, I generally try to avoid giving names to variables if I won’t use that variable anywhere other than that specific calculations. It reduces clutter in the workspace. To calculate the variable y in the script in your post, (assuming no missing values in matws) I might do something like

      y =rowMeans(mataws)*rep(c(31,28,31,30,31,30,31,31,30,31,30,31)/365,50)
      y = colSums(matrix(y,nrow=12,ncol=50))
      x=1957:2006
      …

      The first two lines could be combined, but there is no gain in doing that and you might end up spending 10 minutes figuring out why you parentheses didn’t match up. In the process only one variable y was formed.

    • jeff Id
      Posted Feb 9, 2009 at 7:48 AM | Permalink

      Re: RomanM (#143),

      I would interpret that as saying that the AWS and satellites were not usen in the same analysis.

      From the methods section.

      Additionally, we use a truncated total-leastsquares (TTLS) calculation30 that minimizes both the vector b and the matrix A in the linear Regression model Ax5b. (In this case A is the space-time data matrix, b is the principal component time series to be reconstructed and x represents the statistical weights.)

      Aside from Einstein’s interpretation, I believe the space-time matrix A may be the 3 pc’s from satellite which divide the space according to the trend, Steve’s plot of trend #14 smakes me think that. I’m really not sure what’s going on exactly but I don’t believe location of individual stations was incorporated into the analysis except by correlation to the 3 pc’s in the RegEm analysis.

      You probably have a much better handle on RegEm than myself but it looks like the TTLS method incorporates an additional input matrix, I haven’t quite figured out yet how to implement it in matlab. I’ve spent a few hours this weekend pushing matrices at it, I’m getting some errors still. If A is another input (not even sure where that goes in matlab), it would be how the satellite trend was incorporated into the AWS measurement.

      If someone has thoughts on this I’m all ears.

      Wish we had the code.

    • RomanM
      Posted Feb 9, 2009 at 8:29 AM | Permalink

      Re: jeff Id (#144),

      From the Methods Summary section on the last page of the paper (bold mine)

      We use the RegEM algorithm to combine the data from occupied weather stations with the TIR and AWS data in separate reconstructions of the near-surface Antarctic temperature field.

      In my interpretation, this is using stations with TIR in one analysis and stations with AWS in the other.

      I interpret the “space-time” matrix differently. I think that space-time is used gratuitously to be impressive here when a more straightforward description is more informative. From T. Schneider’s 2001 paper:

      For the following formal description of the EM algorithm, let (X) be a data matrix with n records consisting of p variables, with the values of some of the variables missing in some records. The p variables might represent a geophysical field at p different locations, and the n records might represent incomplete measurements of the p variables at n different times. In the conventional EM algorithm, the number n of records is assumed to be much greater than the number p of variables, so that the sample covariance matrix of the dataset completed with imputed values is positive definite.

      Without going through the matlab programs, my take is that in the Tir and station recon, p would be 5509 + the number of stations used (I believe they did both the case of all and stations and also using only some smaller selected number) and n =600 months. RegEM infills all the missing values (in this case close to half of the values in the matrix are missing since there are no satellite values prior to 1982). In the AWS construction, p = number of manned stations + number of AWS, and again, n = 600. In this case, n is larger than p, but the AWS data is extremely sparse. Presumably the data could be either raw data or anomalies although, I suspect that the former might lead to some problems with the assumptions of the analysis.

    • Jeff Id
      Posted Feb 9, 2009 at 9:24 AM | Permalink

      Re: RomanM (#146),

      Thanks Roman, I think that will help.

    • Kenneth Fritsch
      Posted Feb 9, 2009 at 8:39 AM | Permalink

      Re: RomanM (#143),

      I could do it in one line even incorporating unequal weights for the months (which I didn’t do in my quick analysis of the graph that Steve’s post contained).

      I checked your quick analysis with and without using the weights for the months and the trend slope and stdev of the trend slope differ very little. Actually to be more concise, we have to adjust for leap years (do that in one line).

      Seriously this has been a good learning opportunity/experience for me. I have actually been keeping notes.

• RomanM
  Posted Feb 8, 2009 at 7:29 PM | Permalink

  Re: Steve McIntyre (#113),
  It appears that this is a plot of the simple annual means of all aws.

  plot(1957:2006,colMeans( matrix(rowMeans(recon_aws),12,50)), type=”l”, main = “Annual Means”, ylab =”Anomaly C”, xlab = “Year”)

  (a one line program, great shades of APL!) creates:

  which is a perfect match… except for the fact that it is about .25 degrees lower than the one in the RC graph. Hmmm…

  • Kenneth Fritsch
    Posted Feb 8, 2009 at 9:43 PM | Permalink

    Re: RomanM (#120),

    Roman, I just finished my script for doing the AWS annual trends. Can you really do most of what I did below in one line? Does it weight the monthly means for the different number of days in the months. I am going to try your script tonight or tomorrow AM.

    Anyway the summary of what I found was the following with the 1957-2006 AWS trend in degrees c per decade.

    Adj R^2= 0.09, Trend Slope= 0.139, p= 0.02, Lag1Corr= 0.27, Adj CI of Trend Slope= -0.009 to 0.287

    What I find of interest here is that the lag1 auto correlation of the regression residuals is larger than the corresponding correlation was for the regression using monthly data. Not that I have done that many time series regressions but it has been my experience that the reverse is usually true, i.e. the monthly larger than the corresponding annual.

    Also when the CI limits are adjusted for the lag1 auto correlation the lower limit is less than 0 indicating that we cannot be confident that the trend slope is not 0.

    > MonthAve=rowMeans(mataws)
    > MonthWts=c(31,28,31,30,31,30,31,31,30,31,30,31)/365
    > MonthWtsAll=rep(MonthWts,50)
    > MonthAve2=MonthAve*MonthsWtsAll
    > PreAnnualAve=matrix(MonthAve2,nrow=50,ncol=12, byrow=TRUE)
    >x=1957:2006
    >y=rowSums(PreAnnualAve)
    > lmaws=lm(y~x)
    > summary(lmaws)
    Call:
    lm(formula = y ~ x)

    Residuals:
    Min 1Q Median 3Q Max
    -1.1956 -0.3538 -0.0507 0.3487 1.2978

    Coefficients:
    Estimate Std. Error t value Pr(>|t|)
    (Intercept) -27.784614 11.137698 -2.495 0.0161 *
    x 0.013935 0.005621 2.479 0.0167 *
    —
    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1

    Residual standard error: 0.5735 on 48 degrees of freedom
    Multiple R-squared: 0.1135, Adjusted R-squared: 0.09505
    F-statistic: 6.146 on 1 and 48 DF, p-value: 0.01673

    > corr=acf(residuals(lmaws))$acf[2]
    >Factor=(1+corr)/(1-Corr)
    >n=50
    > AdjStdevSlope=((n-2)/((n/Factor)-2))^(1/2)
    > AdjStdevSlope
    [1] 1.346737
    >StdevSlope= 0.005621
    >TrendSlope=0.0139
    > LowCI=TrendSlope-1.96*AdjStdevSlope*StdevSlope
    [1] -0.0009372217
    >HighCI= TrendSlope+1.96*AdjStdevSlope*StdevSlope
    [1] 0.02873722

Steve McIntyre

Roman, to simplify the comparison of pre- and post-Harry, I’ve made your script into a function. For Ken and others, the function shows a nice feature of R – you can make a list of odds and ends done in the calculation so that when you want to compare pre- and post-, you have easy nomenclature. I’ll show below.

Here’s the function implementing Roman’s decomposition of Steig’s recon_aws

romanm=function(regem=recon_aws,target=Data$aws){ #convert roman script into function
#A. scrape early portion
#A1. recon
early.rec = window(regem,start=1957,end = c(1979,12))
svd.early=svd(early.rec)
M=ncol(regem);M
rec.pc = princomp(early.rec);
#plot(rec.pc$sdev)
rec.score = rec.pc$scores[,1:3];dim(rec.score)
#matplot(time(early.rec),rec.score,type = “l”)
rec.score = rec.score%*%diag(sign(coef(lm(rec.score ~ time(early.rec)))[2,]))
#flip the PCs for positive trend #does not materially affect the procedures
reg.early =lm(early.rec[,1:M]~rec.score[,1]+rec.score[,2]+rec.score[,3])
reccoef=coef(reg.early) #4 63
pred.recon = predict(reg.early)
#max(abs(predict(reg.early)-early.rec)) #[1] 5.170858e-07
#3 PCs reverse engineer all 63 reconstructions in early portion to 100% accuracy

#B. scrape late portion in two segments
dat_aws = window(target[,colnames(early.rec)],end=c(2006,12))
#observed data in same columns as recon
not_na = !is.na(dat_aws)
#available values
apply(!not_na,2,sum) #na values by series
trange = t(apply(not_na,2,function(x) range(time(dat_aws)[x])));trange
#figure out first and last times for data for each aws

#B1. construct 1980 to 1995
names(sort(trange[,1])[61:63])
#identify late starters
# 89828 21363 89332
#1996.000 1996.083 1996.167
names1 = names(sort(trange[,1])[61:63])
reg1 = lm(rec.score ~ early.rec[,names1[1]]+ early.rec[,names1[2]]+ early.rec[,names1[3]])
#sapply(summary(reg1), function(A) A$r.squared)
#Response Y1 Response Y2 Response Y3
# 1 1 1
# the three PCs are linear functions of the three recons with no early data
rec.pc1 = cbind(1,window(regem[,names1],start=1980,end =c(1995,12)))%*%coef(reg1)
#matplot(rec.pc1,type = “l”)
#this reconstructs the three PCs from 1980 to 1995

#B2. from 1996 to 2006
#identify earlier finishers
sort(trange[,2])[1:3]
# 89284 89834 89836
#1992.167 1994.750 1995.917
names2 = names(sort(trange[,2])[1:3])

reg2 = lm(rec.score ~ early.rec[,names2[1]]+ early.rec[,names2[2]]+ early.rec[,names2[3]])
rec.pc2 = cbind(1,window(regem[,names2],start=1996))%*%coef(reg2)
#matplot(rec.pc2,type = “l”)
#this reconstructs 3 PCs from 1996 to 2006

#COMBINE to make PCs from 1957 to 2006
working=ts(rbind(rec.score,rec.pc1,rec.pc2),start=c(1957,1),freq=12) ;tsp(working)
#1957.000 2006.917 12.000
estim.rec = ts( cbind(1, working )%*%reccoef,start=c(1957,1),freq=12)
romanm=list(estim=estim.rec,pcs=working,middle3=names1,late3=names2,regem=regem,data=dat_aws,coef=reccoef,svd.early=svd.early)
romanm
}

Now here’s the base case (uncorrected). Download info.

download.file(“http://data.climateaudit.org/data/steig/Data.tab”,”Data.tab”,mode=”wb”)
load(“Data.tab”)
download.file(“http://data.climateaudit.org/data/steig/Info.tab”,”Info.tab”,mode=”wb”)
load(“Info.tab”)
download.file(“http://data.climateaudit.org/data/steig/recon_aws.tab”,”recon_aws.tab”,mode=”wb”)
load(“recon_aws.tab”)
dimnames(recon_aws)[[2]]=Info$aws_grid$id
Data.old=Data

Now do the base case calc:

roman0=romanm(regem=recon_aws,target=Data.old$aws)

Now downloaded new Steig recon and new Harry, R Rock versions:

download.file(“http://faculty.washington.edu/steig/nature09data/data/awsreconcorrected.txt”,”temp.dat”,mode=”wb”)
grid=scan(“temp.dat”,n= -1) # 37800
#600*63 #[1] 3305400
recon_corr=ts(t(array(grid,dim=c(63,600))),start=c(1957,1),freq=12)
dimnames(recon_corr)[[2]]=Info$aws_grid$id

index=c(28,54) #fix Harry,Racer Rock
for( j in 1:2) {
i =index[j]
loc=paste(“http://www.antarctica.ac.uk/met/READER/”,”aws”,”/”,Info$aws$name[i],”.All.temperature.txt”,sep=””)
test=read.table(loc,skip=1,na.strings=”-“)
test=ts( c(t(test[,2:13])),start=test[1,1],freq=12)
temp=(time(Data$aws)>=tsp(test)[1]) & (time(Data$aws)< =tsp(test)[2]); sum(temp)
Data$aws[,i]=NA; Data$aws[temp,i]=test
}
Data.new=Data

Now calc new results:

roman1=romanm(regem=recon_corr,target=Data.new$aws)

If you want to say compare the pre- and post- PC1, then do this (the lists have kept track of the different things in a consistent nomenclature).

ts.plot(roman0$pcs[,1]-roman1$pcs[,1],type=”h”)

Gary A.

Clean Air 1994 Service History

Site:
Clean Air

Performance:
Repaired and moved on 24 January, pressure jumps erratically. expecially during the colder months.

Steve McIntyre

I did a reverse engineering regression of the “PC1” in the pre-1980 period. An r^2 of 0.97.

Here’s a derivation:
A=anom(Data$surface,index=1970:2000)
A=window(A,start=1957,end=c(1979,12))
count=apply(!is.na(A),2,sum)
A=A[,count>=200] #276 24
K=ncol(A)

B=window(roman1$pcs,start=1957,end=c(1979,12));dim(B)
fm=lm(B~A)
name0=Info$surface$name[match(dimnames(A)[[2]],Info$surface$id)];name0

j=1;
GDD(file=paste(“d:/climate/images/2009/steig/rev.eng.PC”,j,”_coef.gif”,sep=””),type=”gif”,w=380,h=300)
par(mar=c(7,3,2,1))
barplot(fm$coef[2:(K+1),j],col=”grey80″,las=3,ylim=c(-3,3),names.arg=name0)
title(paste(“Reverse Engineer PC”,j,” (Roman)Coefficients”,sep=””))
dev.off()

RomanM

I doubt that it has anything to do with seasonal because the data has been deseasonalized by subtracting some sort of monthly mean from each month. In the AWS reconstruction, presumably manned surface stations were used so that if there is some sort of interpretability to the PCs, it is much more likely that it would be geographical in nature.

Steve McIntyre

Here are reverse engineering plots of the coefficients of the “PC1” and “PC3” against surface stations 1957-1979. What does it mean, if anything? Dunno. Right now, and this is still stumbling in the dark a bit, it doesn’t appear that the coefficients of the PC1 in this period have any intelligible meaning. One really wishes the RegEM jockeys did mundane things like report coefficients and residuals, before they held press conferences.

Kenneth Fritsch

With all of the 5509 complete regressions. etc, I ended up saving a 250 MB RData workspace! I’d like to see Ken do that in Excel.

With all of the 5509 complete regressions. etc, I ended up saving a 250 MB RData workspace! I’d like to see Ken do that in Excel.

I couldn’t and that is why I am using R, though is not 250MB the R limit.

Steve McIntyre

#127. Jeff, I don’t think so, tho I know where you’re coming from. However, the pre- and post-Harry AWS recons are different for every station – some more than others. So HArry is used in the AWS reconstruction somehow.

The Harry error is pretty neat for reverse engineering since it gives two recons with only Harry and RR varied. It’s worthwhile comparing them. There are a lot of differences between pre- and post- Harry recons in detail that should help the reverse engineering.

• RomanM
  Posted Feb 8, 2009 at 9:17 PM | Permalink

  Re: Steve McIntyre (#132),

  The Harry error is pretty neat for reverse engineering since it gives two recons with only Harry and RR varied. It’s worthwhile comparing them.

  LOL. I can’t help but see the similarities with the codebreakers trying to decipher messages from the Enigma electronic encrypting machine during World War 2. They were always looking for two encrypted messages with similar content.

• Jeff C.
  Posted Feb 9, 2009 at 11:03 AM | Permalink

  Re: Steve McIntyre (#132),

  Steve said:

  However, the pre- and post-Harry AWS recons are different for every station – some more than others. So HArry is used in the AWS reconstruction somehow.

  The Harry error is pretty neat for reverse engineering since it gives two recons with only Harry and RR varied. It’s worthwhile comparing them. There are a lot of differences between pre- and post- Harry recons in detail that should help the reverse engineering.

  Nots sure if anyone has done it yet, but it may be interesting to use Roman’s recons of the original data and the the post-Harry/RR revised data (i.e. recons without the real, measured data spliced in) and plot the change in trend for each AWS location on the Antarctica map. Since there would be no measured data affecting the results, the change would be entirely based on the presumed “feedback” into the creation of the 3 PCs. Magnitude of changes and proximity to Harry and RR may give some insight into how the data is used.

  It might also be interesting to use the same method to plot the change in the 3 PC coefficients between the two recons for all the same sights. Presumably the coefficients are related to some physical attribute of the AWS locations (lat, long, altitude, distance to coast, etc.) but who knows.

  I’ll try to do this when I get to work (on my break, of course!)

  • Jeff C.
    Posted Feb 9, 2009 at 11:20 AM | Permalink

    Re: Jeff C. (#151),

    On second thought, scratch that second paragraph. If the coefficient are related to physical attributes, they are fixed, it’s the PCs that change. I need coffee.

Steve McIntyre

#133. Roman, while you derived them as PCs, you’re working in reverse. Presumably any 3 linear combinations of the ‘PCs’ would “work” as well. At this point, any interpretation is likely to be post hoc anyway.

Correlations in Mann’s proxy work tend to be entirely opportunistic/data mining and I would assume that this is happening here somewhere – not necessarily intentionally, but certainly without much introspection.

• RomanM
  Posted Feb 8, 2009 at 9:47 PM | Permalink

  Re: Steve McIntyre (#135),

  You’re very right about that. When I was looking at the satellite data, initially, rather than do svd on the 5509 x 600 data matrix, I took smaller subsets of the data (anywhere from 10 to 50 columns) and did princomp on the data. All of the analyses gave me three components, but the three eigenvalues corresponding to one set were not in the same ratios with each other as in the other sets. All set of three PCs could reconstruct the entire 5509 column data set virtually exactly. Each set of PCs was a transformation of any other set which preserved orthogonality.

  I finally did svd on the whole set on the premise that the RegEM would have calculated their PCs from the same matrix.

Steve McIntyre

Jason Smerdon has a new article on RegEM here . He states:

Equation (4) makes it clear that a calibration interval reconstruction is straightforward to calculate using the final m, S, and B outputs of RegEM in this specific case of paleoreconstructions, and therefore calibration statistics should be reported for such results in the future.

This is the sort of thing that reviewers can attend to. I guess the Nature reviewers either didn’t know or didn’t care.

• John Baltutis
  Posted Feb 9, 2009 at 2:55 AM | Permalink

  Re: Steve McIntyre (#140),

  Well, I guess my query WRT to random missing data is answered in Smerdon, et al:

  It should be noted further that this information, especially in terms of the mean, happens to be “additional” only because of a special property of the dataset to which RegEM is applied herein: missing climate data occur during a period with an average temperature that is significantly colder than the calibration period. This property clearly violates an assumption that missing values are missing at random, which is a standard assumption of EM (Schneider 2006). If the missing data within the climate field were truly missing at random, there presumably would not be a significant systematic difference between the M05 and R05 standardizations, and hence corresponding reconstructions. The violation of the randomness assumption, however, is currently unavoidable for all practical problems of climate field reconstructions (CFRs) during the past millennium and thus its role needs to be evaluated for available reconstruction techniques.

  With the ultimate conclusion:

  …further tests of the RegEM algorithm are highly warranted before great confidence can be placed in RegEM-derived climate field reconstructions (CFRs).

  • bender
    Posted Feb 9, 2009 at 11:05 AM | Permalink

    Re: John Baltutis (#142),
    In other words RegEM is an interpolator, not an extrapolator. To use it for extrapolation is either foolish or devious.

    • UC
      Posted Feb 9, 2009 at 1:24 PM | Permalink

      Re: bender (#152),

      In other words RegEM is an interpolator, not an extrapolator.

      That’s my understanding as well. Just like ICE in the calibration exercises we’ve worked with in here. CCE can extrapolate, but this property comes with a price tag, very wide CI.

Jens Andersen

I am one of the very many bystanders, watching with fascination the nimble footwork of the main actors. This is extemely valuable work in my opinion. I am no expert on climate but my instinct is to be very sceptical of complex statistics based on uncertain data! I am an old physics professor and have refereed numerous physics papers over the years. As a referee you try to apply your experience and judgement to evaluate the soundness of a paper and you criticise obscure procedures(!) but it is extremely rare that a referee goes into the kind of detail displayed on this thread. It seems safe to assume that it has not been done in this case either.

I have a small suggestion: In order to get to an interpretation of the three principal components it might be an idea to make an orthogonal transformation with the requirement that the negative weights become as small as possible. If you had just two stations, it would be more transparent to have the independent measurements as a basis rather than their sum and difference.

• RomanM
  Posted Feb 9, 2009 at 7:51 AM | Permalink

  Re: Jens Andersen (#141),

  I have a small suggestion: In order to get to an interpretation of the three principal components it might be an idea to make an orthogonal transformation with the requirement that the negative weights become as small as possible.

  This is an excellent suggestion. Criteria for the rotation of factors (both, orthogonal and oblique) have been developed for exactly this sort of purpose in Factor Analysis.

  For those not familiar with that type of analysis, a dataset is expressed as a linear combination of unseen “factors” (which are sort of like infilled variables with no prior observations). Sometimes, the factors (which are not unique) can be “rotated (i.e. linearly transformed)so that the coefficients for adding the factors together exhibit “interpretable” characteristics. There are a number of methods developed, but factor rotation is basically an ad hoc procedure open to “cherry picking” when determining the type and amount of rotation.

  If I were an author of the original paper, this might be one way to get some insight into understanding the results of the original RegEM analysis. However, here, while we are still in the process of deconstructing the data, I think it would be of more limited value.

  • bernie
    Posted Feb 9, 2009 at 10:41 AM | Permalink

    Re: RomanM (#145), But wouldn’t it also potentially highlight stations that have been spliced or somehow contaminated? After all isn’t that part of the case against the BCPs?

bernie

RomanM:
I wasn’t, honest! I actually have fond memories of it – though it must have been a main frame version of some kind.

Kenneth Fritsch

I looked at the trends for the TIR reconstruction with the monthly data and compared the regression summaries in the table below to those of the AWS reconstruction. I selected the 4 time periods used previously in an attempt to show the effect of start time on the trend and the uncertainty of those later starting trends. The R script for the TIR regressions is listed below.

The table is headed by the adjusted R^2, the trend slope (TS), the standard error of the trend slope (StdErrTS), the probability of the trend slope being different than zero by chance, the lag auto correlation of regression residuals (CorrLag1) and the confidence intervals that were adjusted using Nychka procedure as applied in Santer et al. (2008).

Both reconstructions show that once the starting point moves forward from 1957, the trends become more uncertain and generally more negative. The TIR reconstruction trends differ from AWS in the higher AR1 values and showing a more positive trend in the later trends – although those trends are not significantly different than 0.

Those very small R^2 values bother me as it apparently does Robert (Rabbie?) Burns in the post above. The statistics in the case of 1957-2006 have a low R^2, but the p value says that the chances of the trend being zero are very low. In the other time periods, the trend cannot be shown to be different than zero. Recall that when I did the trend statistics on the annual data for AWS 1957-2006, the p value was greater than 0.05.

Regress TIR Recon Monthly 1957-2006:

>download.file(“http://faculty.washington.edu/steig/nature09data/data/ant_recon.txt”,mode=”wb”,”TempReconTIR”)
> tir=read.table(“TempReconTIR”)
> mattir=as.matrix(tir)
> dim(mattir)
[1] 600 5509
> y=rowMeans(mattir)
> x=1:600
> lmtir=lm(y~x)
> summary(lmtir)
> plot(lmtir)
> acf(residuals(lmtir))$acf[2]

Regress TIR Recon Monthly 1965-2006:

> mattir1965=mattir[c(96:600),]
> y=rowMeans(mattir1965)
> x=96:600
> lmtir1965=lm(y~x)
> summary(lmtir1965)
> plot(lmtir1965)
> corr=acf(residuals(lmtir1965))$acf[2]

Regress TIR Recon Monthly 1970-2006:

> mattir1970=mattir[c(156:600),]
> y=rowMeans(mattir1970)
> x=156:600
> lmtir1970=lm(y~x)
> summary(lmtir1970)
> plot(lmtir1970)
>corr=acf(residuals(lmtir1970))

Regress TIR Recon Monthly 1980-2006:

>mattir1980=mattir[c(276:600),]
> y=rowMeans(mattir1980)
> x=276:600
> lmtir1980=lm(y~x)
> summary(lmtir1980)
>plot(lmtir1980)
> corr=acf(residuals(lmtir1980))$acf[2]

jeff Id

I was able to get the AWS RegEm redonstruction to run last night.

Deconstructing a Reconstruction of Temperature Trend

It’s interesting what happens to the upslope when you use higher order TTLS truncation paramaters. I could use the help of some of you stats experts in moving forward.

• RomanM
  Posted Feb 10, 2009 at 7:27 AM | Permalink

  Re: jeff Id (#166),

  Good work jeff. Keep ploughing ahead, but don’t neglect your day job 😉 .

  One question, are you working with the original AWS data or the corrected? Is there a version of Data$aws floating around with all of the updates?

  • jeff Id
    Posted Feb 10, 2009 at 7:47 AM | Permalink

    Re: RomanM (#167),

    I’m using the one off CA directory with two AWS series removed. – from 65 to 63, I can’t remember the numbers but Matlab won’t converge with them in Data.

• Kenneth Fritsch
  Posted Feb 10, 2009 at 10:17 AM | Permalink

  Re: jeff Id (#166),

  Jeff, your analysis was the cat’s pajamas/neat/cool. The case for the number of PCs to include in the PCA has always been here at CA a matter of relating the selected PCs to some natural occurrence or phenomenon. What your analysis shows is that the conclusions from the Antarctica reconstructions, i.e. the 1957-2006 temperature trends could be susceptible to cherry picking of the PCs. The excerpt listed below from Steig et al. (2009) gives a glimpse of the rationale for the first 2 of the 3 PCs selected, but without much detail.

  I have enjoyed this thread and could not help but learn from participating/observing here.

  The first three principal components are statistically separable and can be meaningfully related to important dynamical features of high-latitude Southern Hemisphere atmospheric circulation, as defined independently by extrapolar instrumental data. The first principal component is significantly correlated with the SAM index (the first principal component of sea-level-pressure or 500-hPa geopotential heights for 20u S–90u S), and the second principal component reflects the zonal wave-3 pattern, which contributes to the Antarctic dipole pattern of sea-ice anomalies in the Ross Sea and Weddell Sea sectors4,8. The first two principal components of TIR alone explain .50% of the monthly and annual temperature variabilities4. Monthly anomalies from microwave data (not affected by clouds) yield virtually identical results4.

Jeff Id

You know the reconstruction was very sensitive, I think it needed to be done on anomaly data to be more reasonable. There isn’t a lot of data in the matrix though.

I tried to leave a comment at RC but it doesn’t even say awaiting moderation, I think I may have experienced a service upgrade.

Does anyone here know what correct values for the truncation parameter should be?

• RomanM
  Posted Feb 10, 2009 at 12:27 PM | Permalink

  Re: Jeff Id (#170),

  The anomaly data as used by Steig is available from his reconstruction by replacing the places where there is no AWS with NAs. Using the variables from my script in the lead post:

  anoms =window(recon_aws, start = 1980)
  anoms[is.na(dat_aws)] =NA

  Try the RegEM using these anomalies.

  Re: Nic L (#168),

  It’s tough sledding trying to figure out the details of what is claimed to be an obvious turnkey operation! I thought it might be a little easier to do if we had a script to check the consistency of their anomaly adjustments (continuing from the script above):

  #check monthly removal:
  diff.recon = dat_aws-anoms

  #function to check range of monthly “means” for a single aws
  #gives two variables:
  #matrix with min and max of monthly adjustments for each month
  #monthly means when range less than .0001 tolerance otherwise NA

  check.range = function(aws) {
  mons =matrix(aws,12,length(aws)/12)
  maxmin =sapply(1:12, function(x) range(mons[x,],na.rm=T))
  rownames(maxmin) = c(“Min”, “Max”)
  colnames(maxmin) = month.abb
  equal = abs(maxmin[1,]-maxmin[2,])<.0001
  monmean = maxmin[1,]
  monmean[!equal] =NA
  list(maxmin = maxmin, monmean = monmean) }

  monthlies = t(sapply(1:63, function(y) check.range(diff.recon[,y])$monmean))

  monthlies is a 63 x 12 matrix of the monthly adjustments made by Steig to calculate anomalies when done consistently for the station. If the value is an NA, then at least two different values were used for that month in different years. In you want to know the ranges for a given aws: e.g. for station 4, use check.range(diff.recon[,4])$maxmin . Some of the NAs could be no months available at all for that aws.

  • Jeff Id
    Posted Feb 10, 2009 at 12:37 PM | Permalink

    Re: RomanM (#174),

    I saw your earlier post and I’ll try that soon. It seemed possible to me that this anomaly was calculated after the fact because of this in the SI:

    and using the READER occupied weather station temperature data to produce a reconstruction of temperature at each AWS site for the period 1957-2006.”

    This was not the case in the Tir data where they clearly state they used anomaly but since so much data is missing in the AWS record I wasn’t sure. I put another question into RC on this but no answer yet. From what I’ve read about RegEm, cyclic errors need to be treated differently.

    • RomanM
      Posted Feb 10, 2009 at 12:59 PM | Permalink

      Re: Jeff Id (#175),

      I would think that you would get convergence faster in the RegEM algorithm be removing the cyclicities (which are of differing amplitudes at the various locations). It can’t hurt to try.

    • Jeff Id
      Posted Feb 10, 2009 at 1:29 PM | Permalink

      Re: RomanM (#176),

      I think you’re right, I’ll give it a try. I’m thinking about pasting the ends on a bunch of Briffa data too. Who knew the implications of imputation could be so implacably interesting? sorry.

  • Nic L
    Posted Feb 11, 2009 at 7:11 AM | Permalink

    Re: RomanM (#174),

    I thought it might be a little easier to do if we had a script to check the consistency of their anomaly adjustments

    Thanks for the neat script, Roman, that is much more elegant than my efforts. Your results are different from mine for some stations, however. This may well be due to my making errors, but in some cases it could be that there are differences in the basic AWS data. I am investigating further.

    • Ryan O
      Posted Feb 11, 2009 at 8:02 AM | Permalink

      Re: Nic L (#199), Nic & Roman, I was confused a bit by how Steig calculated the anomalies as well. He doesn’t use all available data in all cases. So what I did was subtract the recon from READER and find the mode of each month, then used that value to convert READER to anomalies:

      temp2=aws_recon[,x]-ts(temp, start=starty, end=2006+11/12, deltat=1/12, ts.eps=0.03);temp2=matrix(temp2, nrow=12)
      i=1;while(i0, names(sort(-table(temp2[i,1:(length(temp2)/12)])))[1], 0);i=i+1}
      temp=matrix(temp, nrow=12);temp=temp+as.numeric(mode.vector);temp=matrix(temp, ncol=1) #Make anomalies

      .
      “temp” is the READER data.
      .
      After doing this, even Bonaparte Point lies on the zero line . . . except for all the values they didn’t use. 🙂 There were 31 stations where at least 1 READER month was not used – I’m assuming to increase their verification scores – but it’s pretty obvious that they did not use the READER data “as-is”.

jeff Id

I just had a someone call my attention to this.

trunc = min([options.regpar, n-1, p]);
This parameter is then passed as dummy variable r into pttls.m:

% [Xr, Sr, rho, eta] = PTTLS(V, d, colA, colB, r) computes the
% minimum-norm solution Xr of the TLS problem, truncated at rank r

The truncation parameter is equal to the rank of the least squares method used. In my post this means we only see strong positive slopes with rank 1 TTLS.

jeff Id

Last one tonight. For those of us who don’t live with stats, a very cool version of the truncation parameter in ttls.

Click to access Sima.pdf

Ryan O

Sorry. Added comments.

#Make anomalies. Because Steig did not use all available months for the anomalies sometimes,
#making anomalies using monthly means of the actual data results in offsets for each month. To
#prevent, anomalies will be formed by finding the mode of the RECON – ACTUAL
#for each month. This ensures that ACTUAL – RECON for all ACTUAL data used equals
#zero. Any points different than zero indicate ACTUAL data NOT USED for the
#reconstruction.

#Make READER data into a time series and subtract it from the recon data
temp2=aws_recon[,x]-ts(temp, start=starty, end=2006+11/12, deltat=1/12, ts.eps=0.03)
temp2=matrix(temp2, nrow=12) #Turn into matrix, months in rows
i=1;while(i0, names(sort(-table(temp2[i,1:(length(temp2)/12)])))[1], 0) #Steig’s anomaly adjustment is the mode, so let’s find it
i=i+1}
temp=matrix(temp, nrow=12);temp=temp+as.numeric(mode.vector);temp=matrix(temp, ncol=1) #Make anomalies

RomanM

Wow! This has become an R hotbed! It helps to have a lot of eyes looking at the data from many angles. Way to go Ken!

Lost a couple of comments in mid-stream because of two power outages (freezing rain here today – where is GW when you need it? 😉
). I switched to my laptop and wireless and put the internet connection on a back-up power source so will try again.

Jeff C. (#212), what you are doing is somewhat related to the RE and CE statistics that the hockey
team is so fond of. Your code can be simplified – with practice (and a little experience) you can use more of what R was designed
to do well.

For example, you have

for(i in 1:324){
for(j in 1:63){
if(not_na0[i,j] == FALSE){
estim0.rec[i,j]=NA
late0.rec[i,j]=NA
}}}

which can be done as

estim0[is.na(dat_aws)]=NA
late0.rec[is.na(dat_aws)]=NA

although the second line is not necessary for your calculations of diff0.sd.

You can extract the anomalies without a lot of work (and without Data$aws) in a very simple fashion:

#steig = steig reconstruction
#pcrecon = PC reconstruction

anom = steig
anom[abs(steig-pcrecon)<1E-6]=NA

#SD’s of difference between the two
diff = anom-pcrecon
sds = sd(diff,na.rm=T)

This gives the anomalies embedded in the 1957 to 2006 time frame and calcuatesthe SD’s. You can use this to derive the anomalies in the corrected data set as well without knowing what corrections they may have made to the AWS data set. There is a very small chance (about 1 in a million) of an actual anomaly being flagged as a recon value, not enough to worry about.

• Kenneth Fritsch
  Posted Feb 12, 2009 at 11:09 AM | Permalink

  Re: RomanM (#216),

  Roman, my hammeR is seeing nothing but statistical nails. I may soon seek intervention.

• Jeff C.
  Posted Feb 12, 2009 at 1:47 PM | Permalink

  Re: RomanM (#216),

  Roman – thanks for your help in understanding how to simply this. I need to work on breaking my bad Visual Basic habits. Regarding your map plot, I generated something similar and also saw the negative trends on most sites. Since I’m new to R, I figured I must have made a mistake in coding it because it just couldn’t be right.

  • bender
    Posted Feb 12, 2009 at 1:50 PM | Permalink

    Re: Jeff C. (#232),

    it just couldn’t be right

    The other side of reporting bias. Disbelief leading to failure to report.

bender

Do it today!

Kenneth Fritsch

I did a linear detrend of the annual AWS and TIR reconstructions for the period 1957-2006. I have attached the ARn graphs below. I have only listed the R code I used for AWS since TIR follows on it directly.

Does R have detrend and Durbin Watson functions?

AWS detrended:

TIR detrended:

ARn AWS 1957-2006:

load(“AWS1”)
y=rowMeans(mataws)*rep(c(31,28,31,30,31,30,31,31,30,31,30,31)/365,50)
y=colSums(matrix(y,nrow=12,ncol=50,))
x=1957:2006
lmAWS1957=lm(y~x)
yregress=fitted(lmAWS1957)
Meanregress=mean(yregress)
yadj=yregress-Meanregress
ydetrend=y-yadj
lmdetrend=lm(ydetrend~x)
plot(x,ydetrend)
acf(residuals(lmdetrend))$acf[2]
summary(lmdetrend) #Check for detrend with 0 trend slope here.

• RomanM
  Posted Feb 12, 2009 at 2:58 PM | Permalink

  Re: Kenneth Fritsch (#234),

  Does R have detrend and Durbin Watson functions?

  Rather than answer your question directly, let me suggest how you can do wider searches in R. This may not apply to some lower numbered versions but works in R2.8.0 which I am using. You can search R help for general words usning fuzzy matching (simlar to Clicking Help => Search help …) by typing ?? followed without a space by the word you are looking for. If you want to search for a phrase, it should be put in quotes. Examples of such searches are ??Durbin or ??”nonlinear regression” .

  Sometimes the search produces answers which include procedures in other libraries. For instance, one of the results produced by ??durbin looks like car::durbin.watson . This means that the function durbin.watson is in the car library and to get more information (without loading the library), you can type ?car::durbin.watson . I am not sure whether it produces results only from the libraries that you have on your computer or from a broader database.

  • Kenneth Fritsch
    Posted Feb 12, 2009 at 4:45 PM | Permalink

    Re: RomanM (#235),

    Roman, I have R version 2.8.1 and I was aware of the fuzzy searching technique (you had provided me with it and other searching techniques previously and they have been helpful). I attempted your suggested enteries with no success. I have seen the function online using durbin.watson and I was quite sure its was with R. It may, as you suggest, be part of a library I do not have loaded – although I thought at the time of the downloading of R I took the comprehensive options.

    • RomanM
      Posted Feb 12, 2009 at 6:08 PM | Permalink

      Re: Kenneth Fritsch (#245),

      One more piece of information and then the professor hat is off. Old habits are hard to break.

      I haven’t used it much nor do I know how well it works, but when in doubt try

      http://www.rseek.org/

      ..but that’s all… no more advice … I know when to quit …

    • Kenneth Fritsch
      Posted Feb 12, 2009 at 10:37 PM | Permalink

      Re: RomanM (#247),

      Roman, you do remind me of my graduate school professor, but then he was twenty years my senior and I think I am a few years yours. Anyway thanks for the link as it listed all the regression statistics and the packages in which the statistics resides. I downloaded the package car and now have durbin.watson in R. I also loaded the package stats (which I thought was already available) and that should provide me with dwtest.

Jeff C.

Roman – I apologize as it is painfully obvious you know *much* more about this that I do. I didn’t redefine the function each time, but the minor change was changing this line in the function:

recon=window(dat,start=st)

To this:

recon=window(dat,start=st,end=en)

I couldn’t get the end date to change without that modification. I’m probably missing something here.

• RomanM
  Posted Feb 12, 2009 at 4:20 PM | Permalink

  Re: Jeff C. (#239),

  My bad! Good change! All that explanation for nothing. I noticed the change of the default value of st to 1957 in your script and thought that was the change that you were referring to. The change you made was needed. Sorry about that.

• Jeff C.
  Posted Feb 12, 2009 at 4:28 PM | Permalink

  Re: Jeff C. (#239), and Re: RomanM (#238),

  BTW, that wasn’t meant to sound snarky as I am clearly out of my depth here. I appreciate the patient assistance you and the others have shown over the last week as I learn this stuff.

tty

Looking at the AWS maps, it seems that one could make a pretty good case that temperatures in Antarctica (except for the Peninsula) varies in opposition to the rest of the World. Temperatures rose from the fifties up to about 1980 while the rest of the world was cooling, then after 1980 when the rest of the World was warming, temperatures in Antarctica flattened or declined. I wonder, has temperatures in Antarctica started rising again in the last decade, when they have gone flat in the rest of the World?

maksimovich

Just quickly “eyeballing” reader and comparing with Russian data set there seem to be mismatches at Bellingshausen station.

Temporal plots and trends are here http://south.aari.nw.ru/default_en.html

Go to summary tables with yellow icon.

Reader summary here http://www.antarctica.ac.uk/met/gjma/temps.html

Slightly different geographic coordinates as well.

Ryan O

Hm . . . not sure what I’m doing wrong, but when I try to run Steve’s mapping script, I get:

Error in map(“world”, proj = “azequidistant”, orient = c(-90, 0, 0), ylim = c(-90, :
unused argument(s) (proj = “azequidistant”, orient = c(-90, 0, 0), ylim = c(-90, -60), fill = TRUE, col = “white”)

.
Do I need to install something besides just the mapproj library?

Vernon

Jeff C,

the reason you get those results are AVHRR shows warming from 1980 – 2006 and AWS shows no warming (basically flat) from 1980 – 2006. The only way to get the AVHRR and AWS trends to be similar is to use RegEM to backfill the 1856 – 1980 and since the same source is used to back fill both of them, you get similar trends from 1956 – 2006.

• Jeff C.
  Posted Feb 12, 2009 at 11:30 PM | Permalink

  Re: Vernon (#255),

  Yes, I agree and I understand the differences in the trends. My surprise was that for the 1957 to 2006 period, out of 5509 gridsells, 5508 had a positive trend and only one had a negative trend. For the 1980 to 2006 period, 5490 had a positive trend and 19 had a negative trend. This despite the fact the split was roughly 20 pos/43 neg for the AWS recon trends of 1980 to 2006.

  I simplictically viewing the 5509 gridcells as akin to 5509 little thermometers each as a somewhat independent measurement. Since the recon was derived from measured, satellite data I had expected to see at least some regions that have cooled. We know cooling has taken place at some locations based on the READER surface and AWS measurements. Yet from 1980 to 2006, 99.6% of the satellite recon gridcells show warming. I call BS.

  • Phil
    Posted Feb 13, 2009 at 12:24 AM | Permalink

    Re: Jeff C. (#259),

    From Comiso 1999, page 18 (Comiso 1999 is reference 8 in Steig, et al 2009 and Comiso did the satellite reconstructions and was a co-author in Steig 2009):

    Regression analysis of the January 1979–January 1998 dataset yielded a trend of 0.054 ± 0.012°C yr-1 for open ocean, 0.196 ± 0.032°C yr-1 for sea ice, and -0.050 ± 0.043°C yr-1 for the ice sheet. The corresponding values for the July data are 0.035 ± 0.011°C yr-1 for open ocean, 0.078 ± 0.034°C yr-1 for sea ice, and -0.034 ± 0.130°C yr-1 for the ice sheet. The trend analysis of the January and July data combined, with seasonality subtracted by taking the difference of each
    month and the average for the same month over the 20-yr dataset yielded 0.044 ± 0.008°C yr-1 for open ocean, 0.137 ± 0.024°C yr-1 for sea ice, and -0.042 ± 0.067°C yr-1 for the ice sheet.(my bold)

    Comiso goes on to say:

    An entire year of infrared data was processed in 1992, as well as in 1979, to better understand the seasonal characteristics of surface temperature in the region. Annual averages generated from these 2 yr (Figs. 17a and 17b) are shown to be generally in good agreement, especially within the continent. … overall, the general warming in some areas and cooling in other areas is consistent with that observed in station data. The results of a trend analysis on a pixel by pixel basis (Fig. 17d), using data in Figs. 10 and 11, also shows patterns very similar to those shown in Fig. 17c, especially in the continent.(my bold)

    Figures 17c and 17d are on page 21 of Comiso 1999. Figure 17c is the complete yearly average per pixel of 1992 minus 1979. Figure 17d is the trend result per pixel using the Jan and Jul monthly average data from 1979 to 1998. It really is striking to look at how well the 2 figures match up. I suppose this was a way of validating trends computed from just two months per year vs. doing whole years. Figures 17c and 17d show a very different spatial distribution than that shown in your post above. Specifically, compare the third figure in Jeff C. (#251) (decadal trends 1980 to 1990) to Figure 17c in Comiso 1999 (1992 minus 1979), where Comiso 1999 shows strong cooling in the peninsula.

    It would be interesting to hear an explanation of these differences from Comiso himself.

    • Phil
      Posted Feb 13, 2009 at 12:27 AM | Permalink

      Re: Phil (#260), The superscript tags aren’t working. Please subsitute yr^(-1) in the above for yr-1.

    • Jeff C.
      Posted Feb 13, 2009 at 2:57 PM | Permalink

      Re: Phil (#260),

      On page 1 of the SI they seem to dismiss the earlier Comiso data set due to inadequacies in the cloud masking algorithm he used.

      In Comiso (ref. 8), cloud masking was done by a combination of channel differencing and daily differencing, based on the change in observed radiances due to the movement of clouds. Here, we use an additional masking technique in which daily data that differ from the climatological mean by more than some threshold value are assumed to be cloud contaminated and are removed.

      They go on to say:

      The very strong trends in the earlier data set (ref. 8) are not supported by the ground-based observations. In particular, the very strong cooling trends (>1°C/decade) and the strong warming at ~30° to ~ 60°E in East Antarctica results in poor reconstruction statistics.

      It would be interesting to have an explanation. In the main text acknowledgments it says “J.C.C. (Josefino C. Comiso) performed the cloud-masking calculations and provided the updated satellite data set”. Apparently, he agrees with the new conclusions and methodology.

Ryan O

Argh. It’s late. This sentence:

And the reconstructed information at the gridpoints corresponding to the manned stations does NOT appear to be actual data – it appears to be infilled.

should read:

And the reconstructION information at the gridpoints corresponding to the manned stations does NOT appear to be the ground station data – it appears to be infilled.

RomanM

According to Ryan O (#257),

Gavin says at RC:

Only the 42 occupied stations are used to provide the data back to 1957. The AVHRR data, or the AWS data are used only for calculating the co-variance matrices used in the different reconstructions. Thus the reconstruction can either use the covariance with AVHRR to fill the whole interior back to 1957 (the standard reconstruction), or the covariance with the AWS data to fill the AWS locations back to 1957.

Ryan O says:

To me, this means that the surface station data should remain untouched by the infilling and the satellite data provides the spatial map to help RegEM fill in the areas without data. Someone correct me if I’m already off the rails . . .

So if I were to do something totally crazy – like find the gridpoints in the TIR recon that match the station locations – and subtract the ACTUAL data from the TIR recon, I should get a flat line with zero slope and a y-intercept corresponding to the difference between my base period and Steig’s for anomalies, right?

Gavin would be more or less correct in what he says. According to the Steig paper, one reconstruction used the AWS stations and the manned surface stations in one recon and the second used the TIR satellites and the same surface stations. In the first case , the data from the AWS was “spliced” into the recon. In the second, the results have NO actual measurements from anywhere in any of the series. The fit from the 3 pieces is just too good. One would expect a relationship between the manned surface stations and the satellite recon, more so in the pre 1982 time period, because that is the sole source of information used in extending the recon back to 1957. The relationship will extend to other grid points other than just the one in which the station is located. But no “flat line” because it ain’t actual station data!

• Ryan O
  Posted Feb 13, 2009 at 2:02 PM | Permalink

  Re: RomanM (#270),
  Thanks. That’s part of what I was trying to find out. As mentioned in the post above, PCA is unfamiliar territory for me, so how many PCs something has does not have intrinsic meaning to me.
  .
  As far as the second part of what I was trying to find out – how different the recon is from the actual data – I’ll have more tonight. For the moment, they look to be in a similar boat as the AWS data – warming to about 1970 or so and flat thereafter. This contrasts sharply with the satellite recon, which warms throughout. In other words . . . no match.

Kenneth Fritsch

When I said:

Perhaps the TIR data for this period is simply the measured data and then we are looking at difference between TIR and manned station raw measurements.

Re: RomanM (#270),

Roman’s, post above puts this piece of the puzzle to rest for me. The TIR reconstruction data contains no actual measurements over the entire 1957-2006 period, while the AWS does for the 1982-2006 time period and it is spliced to the reconstructed values for the 1957-1981 period.

So can we agree on Roman’s accounting of the reconstructions and on that basis determine whether Ryan’s analysis is compatible with Steigs calibrations. Those colored dots appear to indicate big differences but we need to put numbers on them and do some statistical analyses as Ryan suggests.

RomanM

Well, because of the AWS reconstruction, we all knew it was getting warmer in Antaractica, but until I looked at the satellite reconstruction, I didn’t know that it is getting worse. The trends of the warming are increasing!

How did I find this out? I looked at the trends in the 5509 recons and calculated the trends pre and post 1982, the year in which the satellite observations began. Vitually all the trends were positive. That this was true for the post 1982 throughout the entire continent was certainly a surprise. But a larger one was to come. In a bit more than 45% of the grid points, the decadal trends in the warming after 1982 were larger than they were before 1982 .

The script for calculating this is

#load sat data
#download.file(“http://faculty.washington.edu/steig/nature09data/data/ant_recon.txt”,”ant_recon.txt”)
ant_recon = read.table(“ant_recon.txt”)
ant_recon = ts(ant_recon, start = 1957,frequency=12)

#calculate trend for pre 1982 period
ant_recon.pre = window(ant_recon, en = c(1981,12))
time.pre = time(ant_recon.pre)
pre.trends=sapply(1:5509, function(x) coef(lm(ant_recon.pre[,x] ~ time.pre))[2])

#calculate trend for post 1982 period
ant_recon.post = window(ant_recon, st=1982)
time.post = time(ant_recon.post)
post.trends = sapply(1:5509, function(x) coef(lm(ant_recon.post[,x] ~ time.post))[2])

postgr = post.trends >pre.trends
sum(postgr)/5509 # [1] 0.4523507

co = c(“green”,”red”)
plot(10*pre.trends,10*post.trends,cex =.4,col=co[1+postgr],xlab=”Pre 1982 Trend”,ylab=”Post 1982 Trend”,main = “Relative Trend Post 1982 vs Pre 1982”)
abline(0,1)
legend(“topleft”,c(“Post Pre”), pch =19, col=co, cex =1.2)

The line in the plot is not a regression line but simply a reference line where the two trends are equal.

You can see from the graph that there are some interesting patterns in the plots and I suspect that it is related to the 3PC reconstruction. It also has a relationship to the relative variabilities (Pre and Post) of each of the 5509 series and I will likely put something up on that in a day or so.

However, we aren’t finished here yet. So how does this relate to the location of the grid points? Taking the plotting done by Jeff C. (#251), and adapting it to my purposes:

#download.file(“http://faculty.washington.edu/steig/nature09data/data/Tir_lons.txt”,”tir_lons.txt”)
lon=scan(“tir_lons.txt”)

#download.file(“http://faculty.washington.edu/steig/nature09data/data/Tir_lats.txt”,”tir_lats.txt”)
lat=scan(“tir_lats.txt”)

sat_coord = cbind(lon,lat)

library(mapproj)

ant.plot = function(dat, xy , size =1, head =NULL ) {
# Function to plot anartctic data dat
# xy = (lon, lat)
map(‘world’,proj=’azequidistant’,orient=c(-90,0,0),ylim=c(-90,-60),fill=TRUE,col=”white” )
a=seq(-180,150,30); for(i in 1:length(a)) lines(mapproject(list(x=rep(a[i],181),y= c(-90:90) )),col=4,lty=2)
a=c(-180,0); for(i in 1:length(a)) lines(mapproject(list(x=rep(a[i],181),y= c(-90:90) )),col=4,lty=1)
a=seq(-50,-80,-10);
for(i in 1:4) lines(mapproject(list(y= rep(a[i],72),x= seq(-177.5,177.5,5) )),col=4,lty=2 )
title( main= head)
co = c(“green”,”red”)
xy = mapproject(xy[,1],xy[,2])
points(xy$x,xy$y,pch=19,cex = size,col=co[1+(dat)])
legend(“bottomleft”,c(“Post Pre”), pch =19, col=co, cex =1.2)
#legend(“bottomleft”,c(“.1″,”.2″,”.3″,”.4″, “.5”), pch =19, col=1, pt.cex = sqrt((1:5)/10)*size )
}

head = “Plot of Relative Size of Trends from TIR Recon”
ant.plot(postgr,sat_coord,.4,head)

You can see that it is mainly in West Antarctica where my plot of the AWS station reconstruction in #222, showed cooling in parts of that area after 1980. Combining that with Jeff’s posts in #236, that would indicate considerable disagreement between the two reconstructions. This seems to contradict the statement by Steig at al.: “Results from our AWS-based reconstruction agree well with those from the TIR data (Fig. 2).” I guess it just depends on how you look at the data…

• Kenneth Fritsch
  Posted Feb 13, 2009 at 4:08 PM | Permalink

  Re: RomanM (#275),

  Well, because of the AWS reconstruction, we all knew it was getting warmer in Antaractica, but until I looked at the satellite reconstruction, I didn’t know that it is getting worse. The trends of the warming are increasing!

  Roman, your anlysis is more descriptive in pinpointing differences, but do not we know that there are (visual) differences in trends between AWS and TIR and over the split in time periods by inspection of the trend lines for both. Using the Steig results for AWS and TIR, we know it is getting warmer in the Anarctica starting when? I think the answer is sometime before 1965 as the latter time period trends are not statistically significant and that is the case for both AWS and TIR reconstructions.

  Could one do differences between AWS and TIR reconstruction data and determine if that difference is statistically significant?

  • RomanM
    Posted Feb 13, 2009 at 4:42 PM | Permalink

    Re: Kenneth Fritsch (#276),

    Ken, the two graphs that I have put up are for the TIR reconstruction only, not the AWS. What they show is that about 45% of the grid points supposedly had a LARGER amount of warming in the last 25 years (from 1982-2006) than they did in the first 25 years (1957-1981). This is the opposite of the pattern evident in most of the AWS reconstruction. This would not only indicate that it is continuing to warm, but at an increasing rate. Did we know this already? For the individual grid points?

• Mike B
  Posted Feb 13, 2009 at 4:10 PM | Permalink

  Re: RomanM (#275), Roman and others…thanks very much for all the work you’ve been putting in on this. I confess it’s a bit difficult to follow, mainly becuase it’s spread out over a couple of threads and several hundred posts.

  Any chance you could give us the Reader’s Digest condensed version? Sort of what we know so far, what still needs to be checked, etc.?

  Thanks.

• Kenneth Fritsch
  Posted Feb 13, 2009 at 5:20 PM | Permalink

  Re: RomanM (#275),

  Ken, the two graphs that I have put up are for the TIR reconstruction only, not the AWS.

  I agree that your graphs detail some surprising differences. I guess I was just attempting to segway in some of my own observations.

  This would not only indicate that it is continuing to warm, but at an increasing rate. Did we know this already? For the individual grid points?

  In a word, no.

Kenneth Fritsch

I looked at the Durbin-Watson statistic for lag 1 through lag 25 autocorrelations of the regression residuals to determine the statistically significant lags for the regression residuals for the AWS and TIR reconstructions over the period 1957-2006 for the monthly and annual data. The results are listed below for those lags with p equal to or less than 0.05 along with the R script for the AWS annual calculation.

I plan next to look at the autocorrelation of the residuals for other temperature anomaly time series in order to determine whether these reconstructions have unusual residual autocorrelations.

Monthly AWS: lag1 =0.00; lag2= 0.00; lag4 = 0.01; lag5 = 0.02; lag 11 = 0.00
Monthly TIR: lag1 = 0.00; lag2 = 0.00; lag4 = 0.04; lag10 =0.02; lag11 = 0.00, lag23=0.02; lag24 = 0.00; lag25 =0.01

Annual AWS: lag1 =0.00; lag12 =0.01
Annual TIR: lag3 =0.03

Durbin-Watson AWS Annual 1957-2006:

> load(“AWS1”)
> y=rowMeans(mataws)*rep(c(31,28,31,30,31,30,31,31,30,31,30,31)/365,50)
> y=colSums(matrix(y,nrow=12,ncol=50,))
> x=1957:2006
> lmAWS1957=lm(y~x)
> durbin.watson(lm(y~x),max.lag=25)
Error: could not find function “durbin.watson”
> local({pkg durbin.watson(lm(y~x),max.lag=25)
lag Autocorrelation D-W Statistic p-value
1 0.274450133 1.4182454 0.014
2 -0.065364008 2.0652344 0.868
3 -0.040907103 1.9847145 0.926
4 -0.077302317 2.0113594 0.686
5 -0.256243368 2.3585350 0.060
6 -0.124549529 2.0518936 0.426
7 0.067970411 1.6634797 0.626
8 0.016061778 1.7644956 0.952
9 -0.032691579 1.8565384 0.584
10 -0.047099316 1.8766900 0.478
11 -0.106726377 1.9536365 0.212
12 -0.243363981 2.2221569 0.014
13 -0.079734441 1.8527688 0.256
14 -0.022667051 1.6985708 0.430
15 -0.049721350 1.6781716 0.350
16 0.217915227 0.9934830 0.092
17 0.264636804 0.8763211 0.042
18 0.086155110 1.1854396 0.612
19 0.052095012 1.2011787 0.792
20 -0.033384531 1.2816059 0.770
21 -0.077851762 1.3534064 0.458
22 -0.076377146 1.3064718 0.424
23 -0.007090531 1.1336581 0.804
24 0.033325637 0.9753253 0.714
25 -0.019368345 1.0387367 0.898
Alternative hypothesis: rho[lag] != 0

Jeff C.

That should say Comiso (ref 8). Not sure why that smiley face showed up.

Not Jeff Id

My other stuff made it through, I didn’t see it. thanks.

Re: Jeff C. (#283),

I’m really surprised the sat data didn’t fit better than that.

The RegEm data doesn’t fit either. I’m pretty sure I’ve done everything right but the uptrend I found is greater than the reconstruction. Visually the basic shape of the data is correct though.

If the basic sat trends don’t match, it makes me wonder if the RegEm result I found isn’t also correct. Is it possible that there was some kind of GISS corrected anomaly used? Anyway, after some time looking at the data it looks to me as simple as the pre-1980 surface station data pasted on AWS.

Ryan O

Pics didn’t work. Carbon error. Steve, please feel free to snip out or delete the above post. Here it is again:
.
A bit belated, but here’s the point-by-point comparison with uncertainties for recon gridpoints to actual station data. The method is pretty simple: I find the mean of the common points between 1982-1990 to form the baseline. The reason is that there are the most stations with data available during that time (27), so it allows me to keep the same baseline for a good number of stations. Here’s that:

Red are peninsula stations, blue are west Antarctica, and black are east Antarctica.
.
So now let’s look at 1957-1981:

To me, that looks pretty poor. The recon does a poor job matching stations backwards in time. In general, just like for the AWS recon, the peninsula warming is OVERSTATED in the early period, east looks okay, and there just isn’t much west data, unfortunately.
.
Now 1982-2006:

Looks a little better than the early period. With the exception of one station, the west stations don’t show any statistically significant difference (though they’re all above zero). Of the 27 stations, only 6 fall below the zero line, and all of the ones outside the confidence intervals are above zero.
.
The picture gets even uglier when you break everything up into 8-year blocks. I’ll put these one right after each other so you can see the progression:
.

.
Qualitatively, the actual data from the manned stations is the same as AWS: Except for the peninsula, no statistically significant warming after 1969. The satellite recon, however, significantly understates the peninsula warming trend (by more than a degree C over the analysis period) and overstates the east/west Antarctica warming.
.
In the 1957-1981 period, more of the recon data falls OUTSIDE the 95% confidence level than within. After 1982, about a third of them are out, over half would be out using a 90% confidence level, and most are out on the high side.

Ryan O

1957:2006: AWS TS =0.139; TIR TS = 0.118; GISS TS = 0.114
AWS p = 0.02; TIR p = 0.02; GISS p = 0.00

1965:2006: AWS TS =0.053; TIR TS = 0.086; GISS TS = 0.048
AWS p = 0.46; TIR p = 0.20; GISS p = 0.29

1970:2006: AWS TS =-0.009; TIR TS = 0.064; GISS TS = 0.020
AWS p = 0.91; TIR p = 0.45; GISS p = 0.73

1980:2006: AWS TS =0.-051; TIR TS = 0.112; GISS TS = 0.023
AWS p = 0.60; TIR p = 0.44; GISS p = 0.81

Beautiful.

Kenneth Fritsch

Below are the results of regressions of the reconstruction anomalies (annual) for TIR minus AWS over several time periods. Typical R script is shown below. The trend slopes of the differences (TS) and probability (p) of these slopes being different than zero by chance are listed.

As the trends are moved forward in time from the 1957 start, the TIR-AWS differences become more significant, i.e. p gets smaller while difference trend slope becomes larger. One can see that the 1982-2006 period shows a difference trend that has a p = 0.07. That time period was selected as 1982 was the year the TIR measurements were initiated (AWS started in 1979). Actually coming forward one more year gives a p = 0.01 for the trend slope for the TIR-AWS difference over the period 1983-2006.

Regressing TIR-AWS difference over 1957-2006: TS = -0.02: p = 0.58
Regressing TIR-AWS difference over 1970-2006: TS = 0.074: p = 0.26
Regressing TIR-AWS difference over 1980-2006: TS = 0.169: p = 0.11
Regressing TIR-AWS difference over 1982-2006: TS = 0.223: p = 0.07
Regressing TIR-AWS difference over 1983-2006: TS = 0.318: p = 0.01

load(“AWS1”)
yAWS=rowMeans(mataws)*rep(c(31,28,31,30,31,30,31,31,30,31,30,31)/365,50)
yAWS=colSums(matrix(yAWS,nrow=12,ncol=50,))
load(“TempReconTIR”)
yTIR=rowMeans(mattir)*rep(c(31,28,31,30,31,30,31,31,30,31,30,31)/365,50)
yTIR=colSums(matrix(yTIR,nrow=12,ncol=50,))
Diff=yTIR-yAWS
x=1957:2006
plot(x,Diff)
lmDiff=lm(Diff~x)
summary(lmDiff)

#Standardize AWS anomaly
load(“AWS1”)
yAWS=rowMeans(mataws)*rep(c(31,28,31,30,31,30,31,31,30,31,30,31)/365,50)
yAWS=colSums(matrix(yAWS,nrow=12,ncol=50,))
yAWS=yAWS[c(24:50)]
AnomAWS=yAWS-mean(yAWS)
#Standardize TIR anomaly
load(“TempReconTIR”)
yTIR=rowMeans(mattir)*rep(c(31,28,31,30,31,30,31,31,30,31,30,31)/365,50)
yTIR=colSums(matrix(yTIR,nrow=12,ncol=50,))
yTIR=yTIR[c(24:50)]
AnomTIR=yTIR-mean(yTIR)
#Obtain difference AWS and TIR anomalies and regress over 1980-2006
Diff=AnomTIR-AnomAWS
x=1980:2006
lmDiff=lm(Diff~x)
summary(lmDiff)

Ryan O

Here are graphs of the temporal evolution of the means of the AVHRR reconstruction vs. manned station data. A little explanation here:
.
The way Steig calculates his verification statistics excludes time as a factor. It’s a simple point-by-point comparison. So as long as the series he’s comparing sort of start at the same point and end at the same point, the verification statistics look plausible. This is fine if you’re determining if there is a difference in mean temperature over a time period, but it is not fine if you are looking for temperature trends. They allow you to pass verification even though the trends within the series could be rejected by an appropriate test. This is critical to Steig’s conclusions. Because his verification statistics are insensitive to trends within the data, by choosing the correct starting and ending dates, the satellite reconstruction shows strong warming. Others have already noted the apparent cherry-picking of dates.
.
To show that his trends can be statistically rejected, a test is needed that is sensitive to trends between arbitrary times. The simplest way I can think of to do this is to compare how the difference in means between the ACTUAL data and the RECON data at a common grid point evolves temporally. From there, I can generate a confidence interval for each point and plot that with respect to time.
.
The plots below are generated in the following manner:
.
1. Paired Wilcoxon test (so a non-normal distribution is okay) between the ACTUAL and RECON data for +/- 48 months (time frame doesn’t really matter – the results are pretty much the same whether you use 24-96 months . . . just changes how smooth it is).
2. Normalize the difference in means to the 95% confidence interval (allows multiple series on the same plot).
3. Repeat for every point in the series.
.
Now you have a time series showing where the rolling +/-x months means are as a ratio of the 95% confidence interval. If any points exceed +/- 1.0, the null hypothesis that the +/-x month means are the same at that point in time is rejected at a 95% confidence level.
.
Just to give Steig the best possible chance of passing this test, I added a step to center the ACTUAL and RECON means for the whole series using only common points. I think a more appropriate way would have been to set the means equal at the beginning or end of the series, but even giving Steig this benefit of the doubt doesn’t let him pass. I then grouped the series by location and ran the test for all stations with about 30+ years of data. The results are not good for Mr. Steig.
.

.
Probably the most striking is the understatement of the peninsula warming – the range on most of those series is 2.5 times the 95% confidence interval . . . which means there is an infinitesimal chance of them being correct. You can also see the reason the 1979-1982 start dates are important, as the AVHRR means tend to be outside the lower 95% confidence interval at those times.
.
I am in the process of grouping and plotting the results from the AWS recon vs. the AVHRR recon. They are similar. The null hypothesis that the means temporally evolve in the same way (and, hence, that the derived trends are valid) is rejected at a 95% confidence level.
.
Apologies in advance to Mr. Steig for being blunt, but his AVHRR reconstruction is utter bunk.
.
.
.
Script (be forewarned – running this takes a LONG time. If you do all the stations and the AWS recon at once like I did, you’re doing 600,000 paired Wilcoxon tests on +/-48 months of data apiece):

satrecon=cbind(manned[,1:37], aws_recon[,127:189]);compare=cbind(manned[,38:74], aws_recon[,64:126])
c=matrix(1:600, ncol=1); a=matrix(1:100, ncol=1)

diffmeans=function(x, y, dat1=satrecon, dat2=compare, mo=48) {
#Center the sat recon on the compare mean – this provides the
#best possible case for sat recon – using common points over
#the entire analysis period.
dat1[,y]=dat1[,y]-mean(na.omit(as.vector(dat1[,y]-dat2[,y])))

#Set the period for the Wilcoxon test
hi=ifelse((x+mo)>600, 600, x+mo)
lo=ifelse((x-mo)(hi-lo)/2){
temp=wilcox.test(dat1[lo:hi, y], dat2[lo:hi, y], paired=TRUE, conf.int=TRUE, conf.level=.95)

#Row 1: Diff. in means normalized to the 95% confidence interval
#Row 2: P-value
#Row 3: x (corresponds to row number in recon)
yvalue=c(as.numeric(temp$estimate/(temp$conf.int[2]-temp$estimate)),temp$p.value,x)
} else yvalue=c(NA,NA,x)
yvalue }

getall=function(a) {
temp=cbind(t(apply(c, 1, diffmeans, y=a)))
temp}

#Repeat for all stations
meantest=apply(a, 1, getall)

#Extract normalized values (1:600) and p-values (601:1200)
temp1=meantest[1:600,];temp2=meantest[601:1200,]
ev.means=cbind(temp1,temp2)
colnames(ev.means)=c(mannednames[[1]], awsnames[[1]], paste(“p-value “, mannednames[[1]], sep=””), paste(“p-value “, awsnames[[1]], sep=””)) #Add column names
ev.means=ts(ev.means, start=1957, end=c(2006,12), deltat=1/12) #Time series for plotting

Ryan O

Oh good God . ..
.

Kenneth Fritsch

Ryan, from my layperson’s perspective I judge that your analysis does an excellent job of informing and detailing where differences between TIR and manned stations and AWS and TIR occur. I am curious whether you determined that the paired differences did not have a normal distribution.

I am sure no disrespect was intended on your part, but I would think your reference to Eric Steig should either be Steig or Dr. Steig or Eric Steig.

I think that while your analysis, which shows differences at individual stations and areas in the Antarctica, provides a more complete picture than that reported in Steig et al. (2009), the authors of that paper are going to focus on the overall trends from the accumulated reconstructed stations. They apparently judge that it is important to be able to state that 1957-2006 trends show a warming Antarctica and they will point to those trends (and differences) as a measure of the accumulation of stations and not individual stations.

What I have been repeating ad infinitum here is that given the authors’ perceived intent, I would have to point to the sensitivity of the starting date to trends and view their conclusions in context as either adding nothing new to the science (outside the exploitation and application of the RegEM methodology to these data along with their own unique cloud masking algorithm for adjusting the TIR measurements) or being rather disingenuous about the selection of a starting date. It turns out, of course, that the 1957 start coincides with a number of manned stations temperature measurement starts and thus the start date can be posed as not so arbitrary. But then I would think the authors are obligated even more to talk about the sensitivity of trends to start dates.

• Ryan O
  Posted Feb 17, 2009 at 12:19 PM | Permalink

  Re: Kenneth Fritsch (#302),

  I am sure no disrespect was intended on your part, but I would think your reference to Eric Steig should either be Steig or Dr. Steig or Eric Steig.

  Yes . . . you are correct. I think I will use “Dr.” going forward. 🙂
  .
  Normal distribution: I ran this first using t-tests instead of Wilcoxon tests. The results are very similar except during periods where one time series moves suddenly relative to another. In those times, the differences do not have a normal distribution. When the series do not move suddenly relative to each other, the differences are fairly normally distributed. To avoid any potential problems, though, it made sense just to abandon the assumption of normality and use Wilcoxon. (Oh – I have an extra “0” above – it’s 60,000 paired tests, not 600,000. Haha!)
  .
  As far as the conclusions in the paper go, nothing here seems to debate the 1957-present trend, and I think the 1957 start date is a reasonable choice. It’s the starting point for decently reliable data that is somewhat continuous. If we had such data back to the 1930s then I might question the choice of 1957, but since we do not, I can see no reason why an analysis should not start there.
  .
  The dissection of the reconstructions, I think, highlights a gross conceptual error in the paper that is more important and fundamental than drawing a straight line from 1957 to present and glossing over the fact that the most recent 30+ years shows a flat trend. The error the authors make is switching the dependent and independent variables. Because of this, I think their conclusions do not add anything to the science; rather, they subtract from the science via confusion.
  .
  1. To make the AVHRR reconstruction, they process satellite information in a manner that has not been validated to produce physical results. This is okay – you have to start somewhere. No issue here.
  .
  2. They apply the processed data to the continent along with the manned station data and perform a calibration. Again, okay.
  .
  3. Using the combination of the two data sets, they calculate 3 PCs for analysis. So far so good.
  .
  4. Using the 3 PCs, they infill everything over the entire analysis period and call this the “reconstruction.”
  .
  In my opinion, #4 is where they went off the rails. That is where the dependent variable (the PCs) suddenly becomes the independent variable and the independent variables (actual temperature information) become dependent. In a physical sense, what they are saying is PC causes temperature. They have reversed the cause-and-effect.
  .
  The error does not end there, however. They are fully aware that the method of parsing the satellite data is unproven, so they generate separate reconstructions using the AWS stations to attempt to show that the results are equivalent. They generate verification statistics that indicate no significant difference between the other reconstructions and the AVHRR reconstruction. Unfortunately:
  .
  1. The reconstructions do not constitute a valid reality check because all of the reconstructions are generated using RegEM or PCA with unproven assumptions. In both cases, the unstated (and unproveable in the absence of actual data) assumptions are a) the spatial variance does not change with time; b) the distribution of the missing values with respect to time does not change the results; and c) the distribution of the missing values with respect to location does not change the results. However, in order to show a), you must first show b) and c). They cannot show b) and c) because that information does not exist. In other words, even if the reconstructions were consistent with each other, the only thing you have shown is that they are consistent with each other. That fact confers zero ability to show that they are consistent with reality unless a), b), and c) can be shown – and nothing the authors have done in any way indicates that a), b), and c) hold.
  .
  2. Not only does showing agreement between reconstructions NOT demonstrate agreement with reality, but the reconstructions do not even agree. The mismatch between the AWS and AVHRR trends post-1980 is significant – especially since the trends in the peninsula, west Antarctica, and the coastline east of the bay go in opposite directions.
  .
  3. The verification statistics used are not appropriate tests if subsets of the data are later going to be used for trend analysis. The RE/CE numbers they show are insensitive to how the data changes with time. I feel any test to validate the use of subsets of data must also evaluate the temporal evolution of the data. Time must be explicitly included.
  .
  4. The only real data present – the ground station data – is not consistent with the reconstruction. The one piece of information they have that is actual data has a statistically significant difference in trend.
  .
  To me, this paper is not good science at all. Its conclusions are not supported by the data provided, and the assumptions required to assume that the AVHRR reconstruction reflects actual conditions are both unstated and unproveable. The manner in which the AVHRR reconstruction is presented and discussed confuses cause-and-effect. The method of verification is not sufficient to allow the data to be subsetted for trend analysis.
  .
  Bender had made a comment earlier that this paper deserves a response in the literature, and I find myself agreeing with him. Whether Antarctica is actually warming is not the issue. The issue is that the method used to determine the trend is decidedly erroneous.

  • Kenneth Fritsch
    Posted Feb 17, 2009 at 1:59 PM | Permalink

    Re: Ryan O (#303),

    Ryan, thanks for your comprehensive reply. I am in essential agreement with what you have stated and do not see where it conflicts with the sermon that I continue to preach.

    • Ryan O
      Posted Feb 17, 2009 at 3:37 PM | Permalink

      Re: Kenneth Fritsch (#305), And I wholeheartedly agree with said sermon. 🙂
      .
      Re: MC (#304), I agree. In the effort to get to #5, the authors glazed over 2-4. #5 isn’t where the science is, really. The real science is in 2-4. That’s where the real learning would take place. But this paper seems to have short-circuited that.

  • Jeff C.
    Posted Feb 17, 2009 at 5:32 PM | Permalink

    Re: Ryan O (#303),

    4. Using the 3 PCs, they infill everything over the entire analysis period and call this the “reconstruction.”
    .
    In my opinion, #4 is where they went off the rails. That is where the dependent variable (the PCs) suddenly becomes the independent variable and the independent variables (actual temperature information) become dependent. In a physical sense, what they are saying is PC causes temperature. They have reversed the cause-and-effect.

    Well put. Take a look at my comment in (#283). They showed just a peak at their hand with figures C and D with respect to the satellite data. Setting aside whether or not the cloud masking was done correctly, this is supposedly the satellite data pre-RegEM/3 PC reconstruction. It looks nothing like the reconstruction for the same period. They didn’t infill, they replaced, and then used the dubious AWS stats to bolster it.

    Regarding a response, go for it. Who knows, maybe you can join Steve in the select club who have had their methods called “bizarre” by the Team. It’s a badge of honor. Gavin called my post on distance weighting the input to RegEM, “not worthy of reading on the slowest of days”. He may have been right, but it still felt good.

    • Ryan O
      Posted Feb 17, 2009 at 7:34 PM | Permalink

      Re: Jeff C. (#307), I had forgotten about that one! Yah, that’s a good one, too . . . I was mystified as to how the recon looked so different from that figure until you and Roman, et. al., explained that the recon had no actual data – a possibility that up to that point had never even occurred to me. 🙂

MC

Ryan I’ll go one further. In my mind PCs are good at spotting characteristics but like you say they are the dependent variables. What I would have done is this:
1) Process AWS data and AVHRR and characterise each station in the 1982 to 2006 period (it is 2006 isn’t it?)
2) Describe the characteristics and propose a possible model/simulation with inputs from other measured parameters and physics. End paper.
3) Someone then runs the model, matches the calibration period then extrapolates into the past using the available data as reference tests
4) We see what works and what doesn’t and maybe we all learn something.
5) Then we do the trend analysis

Kenneth Fritsch

Ryan, I know what Steig authors have said/implied, but one needs to look deeper to see what they are really saying, i.e. what they would say they said/concluded under tough questioning at a later date.

Here is the link to the Steve M comparison for GISS gridded versus the Steig recon (appearing as a graph in the Steig paper) for 1979-2003 in the thread titled Steig versus Hansen.

http://www.climateaudit.org/?p=5027#comments

• Ryan O
  Posted Feb 18, 2009 at 1:33 PM | Permalink

  Re: Kenneth Fritsch (#320), One last thing . . . based on the Steig vs. Hansen figure . . . they claim that 1) they replicate that cooling for that period, but 2) that the overall trend is positive. The satellite recon clearly shows a positive trend from 1970-2006 (keeping with the idea that the warming is continuing in recent decades) – albeit the trend for the whole continent has a poor p-value.
  .
  I also don’t agree that what is important is what they say under later questioning. Most people won’t hear that or care. What is important is what they both claim and imply in the paper.

Ryan O

Accidentally hit submit.
.
Note that they start out by saying “in recent decades”. They then go on to explain what they find is different. The implication in their abstract is that the warming is continent wide – not just isolated to the peninsula – and that there is enhanced warming in West Antarctica – and this is relative to recent decades, not just from 1957.
.
That means the trends need to be broken down by region to address their claims. They achieve the enhanced warming in West Antarctica because their satellite recon is outside the confidence intervals of the actual data. I have not done so yet, but I bet if I just regress the West and Bay regions in the TIR recon, I will find statistically significant trends since 1980.
.
They also claim the continent-wide trend is positive. According to the limited surface station data, it is not. And the context for all of this (for most readers, anyway) is the first sentence: in recent decades.
.
Their Figure 2 in the text clearly implies the trend is continuing:
.
http://www.nature.com/nature/journal/v457/n7228/fig_tab/nature07669_F2.html#figure-title
.
As does the sea ice loss rate shown in Figure 4:
.
http://www.nature.com/nature/journal/v457/n7228/fig_tab/nature07669_F4.html#figure-title
.
You do not show a figure correlating sea-ice loss to temperature trends from 1979 – 2003 unless you are implying that the temperature trends are the cause.

Kenneth Fritsch

I also don’t agree that what is important is what they say under later questioning. Most people won’t hear that or care. What is important is what they both claim and imply in the paper.

Ryan, I think we are in reasonably good agreement on what the authors imply and what the astute analyzer can obtain from the authors evidence and the two are the same.

I have read enough papers like Steig et al. (2009) to immediately look at what the authors present as evidence and what they conclude. I am not looking for actual inconsistencies between the two but rather how the casual reader and media types might interpret it and how the authors, under close scrutiny as scientists, and not policy advocates, would be required to answer. I look for deniability and see lots of it in this paper. I am particularly wary when I do not see the authors give confidence intervals and in this case CIs for trends. Their keeping to the mantra of warming over the last 50 years is also a red flag for me.

Meanwhile I will watch with interest as the analyses of the details of the reconstructions continue here at CA

• Ryan O
  Posted Feb 18, 2009 at 5:42 PM | Permalink

  Re: Kenneth Fritsch (#324), Cool. I think maybe we were just looking at it from different angles. I was looking at if from the perspective that there eventually should be something in the literature to show not only the errors in what is actually explicitly claimed but also to show errors in what was implied – without necessarily attributing the implication specifically to Dr. Steig.

Ryan O

Here’s the satellite recon. Average trend is 0.16 deg C/decade with a standard error of 0.05.

Statistics are slope (deg C/decade), standard error, p-value

Byrd:___.2551, .0452, 0 1957-2006
_______.2231, .0703, .0016 1969-2006
.
Doug:___.2217, .0377, 0 1957-2006
________.1953, .0593, .0011 1969-2006
.
Elizabeth:___.2433, .0436, 0 1957-2006
__________.2100, .0679, .0021 1969-2006
.
Erin:___.1582, .0291, 0 1957-2006
_______.1382, .0468, .0033 1969-2006
.
Harry:___.2151, .0369, 0 1957-2006
________.1878, .0582, .0013 1969-2006
.
Mount Siple:___.1683, .0298, 0 1957-2006
_____________.1479, .0463, .0015 1969-2006
.
Siple:___.1900, .0312, 0 1957-2006
_______.1751, .0493, .0004 1969-2006
.
Theresa:___.1898, .0330, 0 1957-2006
__________.1661, .0526, .0017 1969-2006
.
.
Mario Zuchelli:___.0522, .0196, .0080 1957-2006
______________.0287, .0311, .3563 1969-2006
.
McMurdo:___.1574, .0299, 0 1957-2006
__________.1294, .0470, .0062 1969-2006
.
Russkaya:___.1807, .0324, 0 1957-2006
___________.1537, .0582, .0026 1969-2006
.
Avg:___.1906, .0961, 0 1957-2006
______.1649, .0522, .0017 1969-2006

.
Yet the authors claim good correspondence between the 2 recons.

• bender
  Posted Feb 19, 2009 at 6:47 AM | Permalink

  Re: Ryan O (#328),
  Thank you for the stats. Should have been in Steig’s SI.

  the authors claim good correspondence between the 2 recons

  “Good correspondence”?! More like NO correspondence – at least in the low frequency domain.

Robinedwards

Whilst it is very easy to compute linear trends in climate data I still worry about such calculations when simple plotting of the data shows that such a model is often completely inadequate as a compact description if applied over the sorts of time periods that are routinely discussed. Of course, as a start, one must work with “deseasonalised” data. If raw data are used any putative trend will be completely disguised and lost amongst the relatively huge month to month changes that happen everywhere. Thus I submit that the statistics displayed by Ryan O (#327 and #328), whilst illustrating very nicely the points he wants to make may not be based on the most appropriate model. Unless there is a technically plausible hypothesis that a linear trend underlies the observations it may not be the most fruitful way of describing the observations.

It has been my experience over the last 15 years that very abrupt steps in climate can be detected for individual stations over wide geographic ranges and on many time scales. Thus I think that our investigations might be more fruitful if this step behaviour were to be taken into account when examining hard-won observational data. How such step changes are to be detected is another matter (I have a suggestion!) but the techniques of time series analysis might be a good starting point. It is important to note that we are dealing with historical data, and that we have (presumably) no interest in forecasting until a viable model has been identified. My opinion, for what it’s worth – not much perhaps- is that forecasting climate is almost certainly doomed to failure unless we attempt to explain and reproduce some of the features that have occurred in the past. As far as I’m aware current climate models are not good at this.

Robin

• Ryan O
  Posted Feb 19, 2009 at 5:13 PM | Permalink

  Re: Robinedwards (#333), You should look at these, then. I think you will find this more to your liking:
  Ryan O (#300)
  Ryan O (#301)
  The linear trends were to answer bender’s question. 😉

Robinedwards

Re #334, Ryan O.

Well, yes I have looked at those, and others, and some of them are exactly what makes me worry about fitting simple linear trends over longish time periods! Without going into detail it is evident that some of those plots (most?) are of the saw-tooth variety with fairly consistent grand scale slopes, often of about zero. I don’t have to hand the data in the ASCII format that I can work with to check the details so can’t quote the precise statistics. Some do indeed look like long-term linear sequences. And so on a smaller scale do the saw teeth themselves. They are remarkably linear in many cases for periods of some years and then terminate in a step change in both senses.

This is what I find so remarkable. What underlies these linear periods, and what brings about the abrupt step?

This kind of thing happens with a great deal of “climate” data – real temperature observations, tree ring widths, lake deposits, ice core proxies and others.

Does no-one else find these abrupt changes interesting or worthy of technical pursuit?

• Ryan O
  Posted Feb 20, 2009 at 8:30 AM | Permalink

  Re: Robinedwards (#336), Those aren’t plots of temperature. They’re plots of the difference in means normalized to the 95% confidence level – in other words, they show how the satellite reconstruction differs from the ground data. The strange sharpness of the plots is due to the normalization process.
  .
  For example, let’s say that the satellite mean and ground mean differ by 0.5 degrees for a point in 1970, and the 95% confidence interval happens to be +/- 0.6 degrees. This would plot as a “0.833”, because 0.5 degrees does not exceed the 95% confidence interval.
  .
  Now let’s say the means differ by 0.4 degrees in 1980, but the confidence interval is +/- 0.2 degrees. This would plot as a “2.000” because 0.4 degrees is outside the confidence interval.
  .
  This is what gives those graphs their peculiar shapes.

• Kenneth Fritsch
  Posted Feb 20, 2009 at 9:57 AM | Permalink

  Re: Robinedwards (#336),

  This is what I find so remarkable. What underlies these linear periods, and what brings about the abrupt step?

  We have discussed change points in time series previously here at CA. I have looked at some for global zonal temperature series with methods to statistically determine the change points. There have been papers discussed where change points were determined for the global anomalies going back into the late 1800s. I recall that 4 change points can be distinguished. The problem then becomes, to what does one attribute those change points. I would think some a priori rationale would be preferable, like, for example, with the Antarctica trends measured from the Steig data and looking at the 1957-1980 periods and the 1981-2006 period.

  I plan on eventually looking at change points in the Antarctica data (by way of Steig), but my upfront guess would be that that data is too noisy to find change points. In lieu of that analysis I think looking at the trends going forward as I have attempted here at least gives a reasonable picture of the sensitivity of the trend to start dates.

  The reasoning, as I have heard it, for calculating long term trends in temperature series is that a priori the observation that CO2 levels in the atmosphere have been increasing over those time periods. As a layperson I am not at all sure how change points and/or the use of change points are viewed by professional statisticians.

gdn

Just for comparison, I think it would be helpful to see either another map with all of the relevant station types posted in close proximity, or perhaps a thin circle/box around each station so that we can see which stations match as well as the ones that don’t.

Kenneth Fritsch

Below are the results of the AWS reconstruction data broken down into the areas and analyzed in the exact same manner as was done and reported in my previous post on this thread using the TIR reconstruction data.

AWS, like the TIR reconstruction analysis of trends, shows that going forward from 1957 the uncertainty of the data does not allow determining either cooling or warming trends in East and West Antarctica areas. The nominal trends for East and West Antarctica show more cooling for the AWS versus the TIR reconstruction. The Antarctica Peninsula shows warming that is statistically significant in all time periods. The warming trends are nominally higher in the AWS versus the TIR reconstructions – opposite to the AWS versus TIR relationship shown for East and West Antarctica comparisons.

A major difference between the AWS and TIR area relationships are that the TIR reconstruction might allow the implication of the Steig authors that the West Antarctica temperature trends have increased since 1980 in step with increases in the Peninsula. Using the AWS breakdown of those areas would indicate something of the opposite, i.e. that West Antarctica has nominally cooled since 1980 while the Peninsula has continued on nearly the same warming trend value over all the periods 1957-2006, 1965-2006,1970-2006 and 1980-2006.

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%.

download.file(“http://faculty.washington.edu/steig/nature09data/data/ant_recon_aws.txt”,”TempAWS”)
aws=read.table(“TempAWS”)
download.file(“http://faculty.washington.edu/steig/nature09data/data/aws_lats.txt”,”TempAWSlat”)
download.file
(“http://faculty.washington.edu/steig/nature09data/data/aws_lons.txt”,”TempAWSlon”)
lataws=read.table(“TempAWSlat”)
lonaws=read.table(“TempAWSlon”)
Taws=t(aws)
AWS=cbind(lataws,lonaws,Taws)
EWaws=ifelse(test=AWS[,2] dim(Eastaws)
[1] 42 603
Westaws= TAWS[TAWS$EWaws==”Waws”,]
dim(Westaws)
[1] 21 603
Pen1=ifelse(test=Westaws[,3]>285,yes=ifelse(test=Westaws[,3]-76,yes =”Paws”, no=”NoPaws”),no=”NoPaws”),no=”NoPaws”)
PWaws=cbind(Pen1,Westaws)
Penaws=PWaws[PWaws$Pen1==”Paws”,]
WestMinusPen=PWaws[PWaws$Pen1==”NoPaws”,]
dim(Penaws)
[1] 7 604
dim(WestMinusPen)
[1] 14 604
TEaws1957=colMeans(Eastaws[4:603])
DEaws1957=1:600
lmEaws1957=lm(TEaws1957~DEaws1957)
summary(lmEaws1957)
acf(residuals(lmEaws1957))$acf[2]
[1] 0.1943059
plot(DEaws1957,TEaws1957)

Robinedwards

Re 347, for “still” please read “stiff” !

Robin