# 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