Monday, April 13, 2009

Published Comments on Rahmstorf 2007

After Rahmstorf 2007 (R07) was published in science there were two comments published which were each critical of the results. One comment was by Torben Schmith, Søren Johansen, and Peter Thejll (Schmith et al). The other was by Simon Holgate, Svetlana Jevrejeva, Philip Woodworth, and Simon Brewer. (Holgate et al). Dr. Rahmstorf wrote a reply (RR07) to both of these comments. A year later he followed up with a technical correction to the reply. I will touch briefly on Schmith et al. But focus on Holgate et al. In this post I will show that the response to Holgate et al. was not nearly as robust as presented. My conclusion is the same as Holgate et al. "...we do not agree that simplistic projections of the nature presented in [R07] substantially contribute to our understanding of the uncertainties in the nonlinear relationships of the climate system."

Schmith essentially said that due to the fact that there is a trend in both the temperature and sea level series that it violated the "basic assumptions of the statistical methods used." Remarkably they didn't comment on the smoothing. Even more importantly they ignored the fact that the residuals show a very high level of autocorrelation further showing that the model is mis-specified. Even so they concluded that the likelihood of the model was overstated. In RR07 Dr. Rahmstorf attempted to show that even with the trends removed a good fit remained, but this hardly addresses the main points which stand. The rest of the this post will show that the point is largely moot anyway.

The point that Holgate et al. made was simple. R07 didn't test to see if it could predict withheld data by modeling on the rest of the data. There are two problems with their presentation that confused the point when Dr. Rahmstorf responded. First they never duplicated his original result, instead using an SSA algorithm of their own devising. Second they used a different method to build the models. Using their own methods they were unable to predict the second half of the data set using the first half, and of course likewise in reverse. It is true that Dr. Rahmstorf didn't publish his code until he wrote RR07, but it seems like the should have sent him an email so they could start from the same place.

(I also note that at the beginning of his response to comments Dr. Rahmstorf said he was making "the computer code used in the analysis available for use by other researchers." What he neglects to mention is that the computer code he supplied was essentially useless without ssatrend.m. And nothing in his source code indicates where you could find that. The rest of his code is just simple regression and plotting, not much use to other researchers.)

Of course in RR07 Dr. Rahmstorf went back to his original method and reported that he could predict the second half using the first half of the data, and likewise the first half using the second half. In making his response he failed to report certain important points, and more importantly he made a significant error. A year later in October 2008 Science published his "technical correction" for that error, but as I will show that technical correction also fails.

The question is does the first half of the data predict the second half. Or, from my point of view, put more simply would a model built from the first half of the data be similar to a model built from the second half.

To evaluate this I regressed a model using the first twelve of the twenty-four five year bins. In this case the intercept is the same as the full model, but the coefficient is .44 versus .34 in the model built from the full data set. (RR07 reports .42, I'm not sure why the difference but it doesn't really matter. Also it incorrectly reports this as .42mm/year/degree but it is actually .42 cm/year/degree.) Using the second twelve of the five year bins the coefficient is only .24. (RR07 doesn't report this.) So simply put a model based on the first part of the period would have predicted nearly twice the sensitivity to temperature change that was experienced in the second set. RR07 shows graphs that it claims shows that it is making useful predictions, but based on this analysis it seems like a pretty poor match to me.

I've thrown some R code into the bottom of the source file so that the reader can project 2100 sea level with the different coefficients. But suffice it to say it makes a huge difference whether the coefficient is .42 or .24.

In any event I wonder whether the reader of these two posts has noticed the error in this methodology? The problem is that prior to being put into the bins the data had already been smoothed by the ssatrend algorithm with a 15 year window. This means that some of the original data from the second period is actually influencing the smoothed data in the first period. I don't know how Dr. Rahmstorf discovered this but as I said in October 2008 he published a technical correction. In the technical correction he noted the problem;

"This is correct, but it was illustrated by an incorrect figure (Fig. 1),in which the first half of the smoothed sea-level curve (1882 to 1941) was used to predict the sea level for 1942 to 2001. Because the smoothing procedure used a 15-year time window, the smoothed sea-level curve up to 1941 effectively contains sea-level information up to 1948. When this error is corrected and only annual sea-level measurements from 1882 to 1941 are used, the obtained fitgives a sea-level slope of 0.35 mm/year per °C"

(.35 mm/year/°C should, of course, be .35cm/year/°C)

This, he helpfully pointed out, was almost equal to the full trend of the original model even more strongly demonstrating how well a model from the first period predicted the rest of the data.

He failed to note three things. First in the response he said that he trained the data using the period 1880-1940. In the technical correction this was deftly changed to 1882-1941. It turns out the result is highly sensitive to this choice. Second based on my work I have determined that he changed the window size to 10 years from 15. Third using this exact technique the coefficient in the second period is only .22 this is unreported in the technical correction and is still much different than the trend from the first period.

In the data sets supplied by RR07 the sea level data begins in 1870, but the temperature data begins in 1880. As I've said before it isn't clear why he chose to bin the data at all, since it doesn't make any difference to the result, but the way that he did the binning in R07 starts the analysis in 1882. You can only find this out from the code, as R07 says it is using the period 1880-2001. I note that it doesn't change the results of the initial analysis.

In RR07 he says; "...but using only the first half of the data set (1880 to 1940) for deriving the statistical fit." RR07 supplied the code for R07 but not the code for the analysis in the response. This doesn't really matter because as I have duplicated he got the results in RR07 by using the approach of just regressing the first twelve bins. (This effectively started in 1882)

I have duplicated the .35 figure reported in the technical correction. To do this you have to use exactly 1882-1941. You also have to reduce the window to 10 years from 15, this is unreported in the technical correction. The result is non-robust to changes in either the exact year range or the window. If you move the window one year back this lowers the first period coefficient to .26. If you move it one year forward it raises the coefficient to .41. Combining the fact that the first and second periods don't match, to the fact that trivial changes in date ranges and window selections make large changes in the results shows that this model is not well specified. It is certainly not useful for updating sea level predictions beyond the results of AR4.

As a closing note on this analysis I want to say that it may not have been intentional, but R07, RR07, and the technical note were not nearly transparent enough. Instead it seems that they were written to make a point. It also points out, once again, the need for complete code disclosure.

Code for this analysis can be found here.

No comments:

Post a Comment