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.
Monday, April 13, 2009
Saturday, April 11, 2009
Duplicating Rahmstorf 2007
I had quite an interesting time replicating the calculations from "A Semi-Empirical Approach to Projecting Future Sea-Level Rise" (Science 1/19/2007 R07) by Stefan Rahmstorf. The general concept in the paper is very simple. He estimates a linear equation for the rate of sea level rise based on global temperature from measured data. The difficulty came in because the paper relies on an algorithm that is not generally available. To demonstrate this it seems that although there were two critical comments on his paper published in Science neither replicated the original calculations of R07. It only became possible to replicate the results when Dr. Rahmstorf published his code as part of his response to those comments. Even then you had to do a little hunting. In this post I will show that I have duplicated those results, so that my subsequent comments make sense.
R07 is quite simple in concept. His theory is that the rate of sea level rise has a linear correlation with surface temperature. The idea is that you start with a situation where sea level is stable. Then you raise the temperature causing sea level to rise until is reaches a new equilibrium. I won't explain the whole paper, but this makes perfect sense. There are, of course, multiple inputs to sea level rise, but pretty much all of them should respond to a change in surface temperature, the question is how fast and how much.
To estimate this linear equation he uses sea level data from Church and White (IPCC), as well as the GISS global temperature record. But this is where things get a little tricky. He doesn't do a simple regression on the temperature and sea level because this data is "noisy." Instead runs both the temperature and the sea level data through a process to separate the trend from the noise. This isn't explained in the main text of the paper, but it is referred to in the text below figures 2 and 3. "Both temperature andsea-level curves were smoothed by computing nonlinear trend lines, with an embedding period of 15 years (14)." He uses the word "smoothed" but this is not a typical smoothing algorithm. Use of a typical polynomial fit will not get as "good" a result, as I discovered in early attempts.
I might have discovered this more quickly if I had read all the way to the bottom of Dr. Rahmstorf's reply to comments where a link to the code and data was published. The code and data weren't linked from the original paper which is where I started. Both published comments also didn't have the benefit of this code and data which is clear from reading them, and I guess after that everyone called it a day. But I will get to this stuff in a later post.
Reference 14 from the paper is "J. C. Moore, A. Grinsted, S. Jevrejeva, Eos 86, 226 (2005)." This paper is titled "New Tools for Analyzing Time Series Relationships and Trends." It is actually a short review article of several new mathematical techniques. At first glance it certainly isn't immediately apparent that this tells us how the "smoothing" was done. But there is a section titled "Nonlinear Trends in Sea Level and Temperature." This section refers to the use of Single Spectrum Analysis (SSA) to extract a nonlinear trend. It refers to Ghil et. al. 2002. So this is at best an indirect reference.
At this point I still hadn't discovered the code at the bottom of the reply, but I did have a lead. So I looked into SSA and discovered software at UCLA that was available. To make a long story short SSA is necessary but not sufficient to duplicate the R07 results. In fact further searches, and eventually the code led me to the fact that R07 relied on a matlab function called ssatrend.m. Dr. Rahmstorf on Real Climate indicated that this was available from Aslak Grinsted.
I wrote Dr. Grinsted who wrote me back very promptly, and sent me the source code to ssatrend.m. He also commented that he had no idea how Dr. Rahmstorf had gotten a copy of it, and that he had never meant for it to be distributed. I think that he was just concerned about it being unsupported. My own view is that it is pretty strange to use some random piece of code in a published paper without making the code your own. Especially where, as in this case, the strength of ther result depends on using this specific algorithm. It isn't available from any of the usual Matlab repositories which is the first place I looked.
In any event with a little effort I found an SSA algorithm written in R on source forge. I checked its results agains the UCLA program and they were identical so I knew the foundation was good. Then with some effort, I translated ssatrend.m into R so that I could run it on an open platform.
Finally I translated the code supplied in the reply written by Dr. Rahmstorf into R so that I could at see if I was starting from the same place. This was successful.
Using the supplied input files for sea level and temperature I get a correlation coefficient (R) of .88 exactly as reported in R07. In addition the following graphs are identical to figures two and three (top) of R07. The code is here.
At the end of the day the use of this unpublished algorithm made this much harder than it should have been as the underlying concept is very simple. I have no idea how a reviewer could have evaluated whether the use of this algorithm made sense. Having said that I don't really see any problem with it, and I look forward to trying ssatrend with other data.
But now that I can get the exact results from the paper, I have a couple of comments that go beyond what has already been written in Science.
R07 is quite simple in concept. His theory is that the rate of sea level rise has a linear correlation with surface temperature. The idea is that you start with a situation where sea level is stable. Then you raise the temperature causing sea level to rise until is reaches a new equilibrium. I won't explain the whole paper, but this makes perfect sense. There are, of course, multiple inputs to sea level rise, but pretty much all of them should respond to a change in surface temperature, the question is how fast and how much.
To estimate this linear equation he uses sea level data from Church and White (IPCC), as well as the GISS global temperature record. But this is where things get a little tricky. He doesn't do a simple regression on the temperature and sea level because this data is "noisy." Instead runs both the temperature and the sea level data through a process to separate the trend from the noise. This isn't explained in the main text of the paper, but it is referred to in the text below figures 2 and 3. "Both temperature andsea-level curves were smoothed by computing nonlinear trend lines, with an embedding period of 15 years (14)." He uses the word "smoothed" but this is not a typical smoothing algorithm. Use of a typical polynomial fit will not get as "good" a result, as I discovered in early attempts.
I might have discovered this more quickly if I had read all the way to the bottom of Dr. Rahmstorf's reply to comments where a link to the code and data was published. The code and data weren't linked from the original paper which is where I started. Both published comments also didn't have the benefit of this code and data which is clear from reading them, and I guess after that everyone called it a day. But I will get to this stuff in a later post.
Reference 14 from the paper is "J. C. Moore, A. Grinsted, S. Jevrejeva, Eos 86, 226 (2005)." This paper is titled "New Tools for Analyzing Time Series Relationships and Trends." It is actually a short review article of several new mathematical techniques. At first glance it certainly isn't immediately apparent that this tells us how the "smoothing" was done. But there is a section titled "Nonlinear Trends in Sea Level and Temperature." This section refers to the use of Single Spectrum Analysis (SSA) to extract a nonlinear trend. It refers to Ghil et. al. 2002. So this is at best an indirect reference.
At this point I still hadn't discovered the code at the bottom of the reply, but I did have a lead. So I looked into SSA and discovered software at UCLA that was available. To make a long story short SSA is necessary but not sufficient to duplicate the R07 results. In fact further searches, and eventually the code led me to the fact that R07 relied on a matlab function called ssatrend.m. Dr. Rahmstorf on Real Climate indicated that this was available from Aslak Grinsted.
I wrote Dr. Grinsted who wrote me back very promptly, and sent me the source code to ssatrend.m. He also commented that he had no idea how Dr. Rahmstorf had gotten a copy of it, and that he had never meant for it to be distributed. I think that he was just concerned about it being unsupported. My own view is that it is pretty strange to use some random piece of code in a published paper without making the code your own. Especially where, as in this case, the strength of ther result depends on using this specific algorithm. It isn't available from any of the usual Matlab repositories which is the first place I looked.
In any event with a little effort I found an SSA algorithm written in R on source forge. I checked its results agains the UCLA program and they were identical so I knew the foundation was good. Then with some effort, I translated ssatrend.m into R so that I could run it on an open platform.
Finally I translated the code supplied in the reply written by Dr. Rahmstorf into R so that I could at see if I was starting from the same place. This was successful.
Using the supplied input files for sea level and temperature I get a correlation coefficient (R) of .88 exactly as reported in R07. In addition the following graphs are identical to figures two and three (top) of R07. The code is here.
At the end of the day the use of this unpublished algorithm made this much harder than it should have been as the underlying concept is very simple. I have no idea how a reviewer could have evaluated whether the use of this algorithm made sense. Having said that I don't really see any problem with it, and I look forward to trying ssatrend with other data.
But now that I can get the exact results from the paper, I have a couple of comments that go beyond what has already been written in Science.
New consensus on sea level rise?
Recently because of some posts on stoat I have become interested in how the consensus is moving on projected sea level rise. There are clearly some scientists who now believe that 1M above 1990 by 2100 is likely. These include Stefan Rahmstorf, and Aslak Grinsted who have published papers based on empirical analysis to come to this conclusion. There are a number of issues with the paper by Dr. Rahmstorf some of which were published as comments in Science. I will post some additional comments which follow the publishing of the code and data. The paper by Dr. Grinsted is more interesting to me, but because of the time periods he uses to train and test his model much of the input data and assumptions are necessarily pretty fuzzy.
In AR4 the consensus estimates for sea level rise are on page 821 figure 10.33. Table 10.7 on the previous page breaks down the components. The range for the various scenarios including the error bars is roughly .2 meters to .6 meters. Many people seem to feel that they just threw up their hands at faster ice sheet discharge but table 10.33 includes figures for that at the bottom. They declined to add these into the projections because they couldn't assess the likelihood of this happening. It is important to note that these would have only added .09 to .17 meters to the high end of the range bringing the top to about .8 meters. (It would have slightly lowered the bottom of the range as well.)
So I wonder how far the consensus has actually changed in the last couple of years. Or are there simply some scientists who believe the consensus is low? Should the best estimate now be considered 1M? Or should we stay with the IPCC conclusions? Of course I don't attend the conferences with the types of experts who are called on to determine the consensus view. But I note that it doesn't seem to me that any of the lead authors of Chapter 10 are authors of these recent studies. (I could easily be wrong as cross checking that type of thing is difficult.)
Anyway I am going to write a couple of posts looking at these empirical papers. I am also looking at building an empirical model of my own to see if it will improve on the results of R07.
In AR4 the consensus estimates for sea level rise are on page 821 figure 10.33. Table 10.7 on the previous page breaks down the components. The range for the various scenarios including the error bars is roughly .2 meters to .6 meters. Many people seem to feel that they just threw up their hands at faster ice sheet discharge but table 10.33 includes figures for that at the bottom. They declined to add these into the projections because they couldn't assess the likelihood of this happening. It is important to note that these would have only added .09 to .17 meters to the high end of the range bringing the top to about .8 meters. (It would have slightly lowered the bottom of the range as well.)
So I wonder how far the consensus has actually changed in the last couple of years. Or are there simply some scientists who believe the consensus is low? Should the best estimate now be considered 1M? Or should we stay with the IPCC conclusions? Of course I don't attend the conferences with the types of experts who are called on to determine the consensus view. But I note that it doesn't seem to me that any of the lead authors of Chapter 10 are authors of these recent studies. (I could easily be wrong as cross checking that type of thing is difficult.)
Anyway I am going to write a couple of posts looking at these empirical papers. I am also looking at building an empirical model of my own to see if it will improve on the results of R07.
Subscribe to:
Posts (Atom)