# Nicolas Nierenberg 4/6/2009 # Updated 4/13/2009 comments to reflect changes to use new SSA # Also added new analysis to the end. # Analyzing Rahmstorf 2007 (R07), and taking a look at an alternative # model of sea level rise. # In his response to comments on his original paper Dr. Rahmstorf # included the data file that he used as well as the matlab script. # Unfortunately the matlab script uses a function called ssatrend which # I can't find anywhere on the internet, nor does there appear to be an # equivalent in R. # Based on a comment from Dr. Rahmstorf on the internet I wrote Aslak Grinsted # and he sent me an unpublished version of ssatrend.m which I have now translated # into R # In his response to comments Dr. Rahmstorf also supplied the Church and White # data set he used in his paper. I didn't find this exact data set at http://www.cmar.csiro.au/sealevel/sl_data_cmar.html # which appears to be their primary data source. It is an annual data set, but # it appears to be very similar to the monthly data set which appears on that site. # In the following code I try to stick as closely as possible to the Rahmstorf # matlab script "sealevel.m" # Install ssa package install.packages("simsalabim",repos="http://R-Forge.R-project.org") library("simsalabim") install.packages("matlab") library("matlab") library("nnet") # This function is an R rewrite of ssatrend by Aslak Grinsted, but not as complete ssatrend = function(x,window) { x.decomp = decompSSA(x,window) # Now we search the eigenvectors for the one representing the trend. e = x.decomp$U indx = which.is.max(abs(sum(e))/sum(abs(e))) # Just need the one vector representing the trend e = e[,indx] # Now create a filter using this vector efilt = convolve(e,flipud(e),type="o") # Normalize efilt = efilt/sum(efilt) # Now we have to pad the original series so that we can get a full trend. # We do this by adding values to the start and end based on the linear trend # of the first window of values. n = length(x) df = data.frame(x=x[1:window],time=seq(1,window)) res = lm(x~time,data=df) left.df = data.frame(time=seq(1-window,0)) x.pad = c(predict.lm(res,left.df),x) df = data.frame(x=x[(n-window+1):n],time=seq(1,window)) res = lm(x~time,data=df) right.df = data.frame(time=seq(window+1,2*window)) x.pad = c(x.pad,predict.lm(res,right.df)) # Now just run it through the filter result = filter(x.pad,efilt) result = result[(window+1):(length(result)-window)] result } # Let's begin by loading the Rahmstorf supplied data church_13221 = read.table("Sea Level/1141283s_codes/church_13221.txt") seayear = church_13221$V1 sealevel = church_13221$V2/10 giss_landocean = read.table("Sea Level/1141283s_codes/giss_landocean.txt") gissyear=giss_landocean$V1 gisstemp=giss_landocean$V2 * .01 # Use ssatrend to do the smoothing to as in the original paper window = 15 smoothsealevel = ssatrend(sealevel,window) smoothgisstemp = ssatrend(gisstemp,window) # For some reason he rescales to 1990 sea level. g1990=(smoothsealevel[120]+smoothsealevel[121])/120 sealevel = sealevel-g1990 smoothsealevel = smoothsealevel-g1990 rateofrise = diff(smoothsealevel) # Now he bins the data rather than just running a regression against the smoothed data # He uses the indexes to sort of align the years between Sea level which begins in 1870 # and temperature which begins in 1880 mtemp = numeric() mrate = numeric() myear = numeric() for(i in 1:24) { mtemp[i] = mean(smoothgisstemp[(5*i-2):(5*i+2)]) mrate[i] = mean(rateofrise[(5*i+7):(5*i+11)]) myear[i] = mean(gissyear[(5*i-2):(5*i+2)]) # I guess this is for plotting } # Rahmstorf reports the correlation, which is .88, note that with a seven year embedding this # drops to .68, at 17 years correlation is at .9. Above 17 years the correlation stays about the same # but the rate starts to drop. cor(mrate,mtemp) # Now we are ready to compute the regression, note that the p-value is very similar to R07 # at 1.3 x 10^-8 res = lm(mrate~mtemp) summary(res) # Plot equivalent to Fig. 2 in Rahmstorf 2007 plot(NULL,xlim=c(-.5,.5),ylim=c(0,3.5),type="n",col="black",xlab="Warming above 1951-1980 mean",ylab="Rate of Sea Level Change (mm/year)") points(mtemp,mrate*10,col="red") # I cheated and just put in the values from the regression rescaled for the graph # to draw the line. lines(c(-.5,.5),c(-.001721,.5*3.39748+1.68153),col="black") # Now just regress the full delta sea level with full smooth temp, the results are basically the # same, but it sets up the plot. res = lm(rateofrise[12:131]~smoothgisstemp[3:122]) summary(res) plot(NULL,xlim=c(1870,2001),ylim=c(0,4),type="n",col="black",xlab="Year",ylab="Rate of Change (mm/yr)") # Note that I use the giss year for the sea level. I didn't spend a ton of time on this but this is # is how R07 lined them up. lines(gissyear[3:122],10*fitted(res)) lines(gissyear[3:122],10*ssatrend(rateofrise[12:131],window),col="red") # Naively compare a model based on the first 12 bins with a model based on the second 12 bins. res = lm(mrate[1:12]~mtemp[1:12]) # coefficient is .44 versus .34 in the full model summary(res) res = lm(mrate[13:24]~mtemp[13:24]) summary(res) # The issue with the simple test above is that due to smoothing the first 12 bins are based # partially on data in the second period. (Likewise for the second 12 bins). This was # noted in a technical correction posted by Dr. Rahmstorf in Science in October 2008, more than a year # after his original response. # So to do this correctly we need to break the data sets up in to two pieces, and then process each # piece separately. This is where the date ranges get interesting. It turns out it matters # that you start with 1882 instead of 1880. This range is different was noted in the technical correction # but it is different than the 1880 to 1940 stated in the original response. I note however that the # way that the binning was done, I believe that both the original paper and the response actually start # in 1882. That made no difference to the full period results. # Another important point is that the reported result in the technical correction can only be achieved # by reducing the window size to 10 years. It makes sense to reduce the window with the smaller data # set, but this reduction was not mentioned in the tehcnical correction. # The technical correction also doesn't report on the result for the second half of the data, using the # split data sets, the result is very different. # These choices yield the same result as the techincal correction window=10 gissrange.per1=3:62 gissrange.per2=63:122 sealevelrange.per1=12:72 sealevelrange.per2=72:132 # Moving the first period back one year changes the coefficient to .26 #gissrange.per1=2:61 #sealevelrange.per1=11:71 # Moving the first period forward one year changes the coefficient to .41 #gissrange.per1=4:62 #sealevelrange.per1=13:73 smoothgisstemp.per1=ssatrend(gisstemp[3:62],10) smoothgisstemp.per2=ssatrend(gisstemp[63:122],10) smoothsealevel.per1=ssatrend(sealevel[12:72],10) smoothsealevel.per2=ssatrend(sealevel[72:132],10) rateofrise.per1=diff(smoothsealevel.per1) rateofrise.per2=diff(smoothsealevel.per2) res = lm(rateofrise.per1~smoothgisstemp.per1) summary(res) # The first period coefficient is .35 As reported res = lm(rateofrise.per2~smoothgisstemp.per2[2:61]) summary(res) # The second period coefficient is .25 As not reported # Just to show that binning makes no difference once again. for(i in 1:12) { mtemp[i] = mean(smoothgisstemp.per1[(5*i-4):(5*i)]) mrate[i] = mean(rateofrise.per1[(5*i-4):(5*i)]) } res = lm(mrate[1:12]~mtemp[1:12]) summary(res) # coefficient is .35, as reported. # Let's try that for the second period for(i in 1:12) { mtemp[i] = mean(smoothgisstemp.per2[(5*i-4):(5*i)]) mrate[i] = mean(rateofrise.per2[(5*i-4):(5*i)]) } res = lm(mrate[1:12]~mtemp[1:12]) summary(res) # coefficient is .25, not reported. # Make a quick estimate of sea level rise given a 3C temperature increase # define f as the linear equation from the model (times ten to make it millimeters/year) f = function(x)(x*.339748+.168153)*10 # The beginning rate is defined by the temperature in the last five year period, which is close enough # to the satellite observed rate at 3.1. The ending temperature is defined as the starting temperature # plus 3. I average the rate and multiply by 100 for the number of years (just multiply by 50) (f(3.31)+f(mrate[24]))*50 # This comes out to 78 cm # If you want to use 4C then (f(4.31)+f(mrate[24]))*50 # It is 95cm