## Building a ‘Hokey’ stick

Posted by Jeff Id on October 12, 2008

If this is your first time here recently see this article Will the Real Hockey Stick Please Stand Up?

The first graph of this post is the final result after the 1209 series were magnified and averaged together. Vertical scale is SD units calibrated from 1850-1995. This is my first post which is based on R software.

The point of this post is just to see the types of data going into the latest hokey hockey stick and what shape the groups provide. If you look at all the graphs after the first one, you might see components which add up to the shape of the final graph. Also, you will notice shapes which you can’t find. The fact that you cannot find certain valleys and peaks means very clearly that individual proxies had substantial influence in creating historic trends.

Luterbacher instrument proxies with measured temp data.

Tree ring composite data possibly already HS sorted. Not too many series.

3001 – Conglomeration of marine sediments

4000 -Lake sediments

4001 Rainfall and isotope

5000 – Moisture and grape ripening records – for god’s sakes I bet grape ripening is linear with temperature.

5001-Historic documents!!

6000 – Speleothum – Cave stalactite isotope data.

More cave speleothum data

7000 Isotope ratios in coral

7001 Coral growth (I think) two proxies only

Schweingruber latewood MXD data. I bet $30 this is our little ice age. What tweaks me about this one is the fact that the last 38 years are painted in upslope. (Fake data!) Without it these series don’t get accepted and I bet no little ice age.

Ice core isotope and precipitation data

Ice core isotope records

9000 – Tree ring width data. 76.7% of the 1209 series

Did you catch the proxy which was created from ripening grapes? Imagine converting that plot to temperature.

## Andy said

Erm….if they all measure “temperature” why are they all so different?

As they are all different, which are truthful regarding “temperature” and which are false and how do you tell?

Or is it just those that are Hocky stick shaped measure temperature?

## Jeff Id said

Andy,

You are exactly where I was not that long ago. I have a link which answers some of your questions.

http://noconsensus.wordpress.com/2008/09/11/ten-things-everyone-should-know-about-the-global-warming-hockey-stick/

I hope you like it.

## Demesure said

Hi Jeff

Can you post your R script please ? (don’t forget to enclose it with

`)`

I would like to change the graphs' title to "Icecores : x / 1029 proxies", "Tree ring width: y / 1029 proxies" ...

## Demesure said

@3 I meant "don’t forget to enclose it with the tag CODE ".

## Jeff Id said

Demesure,

I tried to upload it today and apparently I mad a mistake and didn’t save the final version of my code. I have something close though so I will fix it and post tonight after work. That’s the main reason I am learning R, so people can review my work easily. C++ would be huge to go through.

Perhaps you can help me, I needed to make an indexed array of time series so I could join them in a ts.union. The function becomes very slow when the union grows large. I gave up after a couple of hours but what I needed was to address a bunch of time series like this x[i]. I ended up just typing in 4 series like xa, xb, xc, xd. It was a pain but it got the job done.

You’ll see what I mean later.

## Gary said

Jeff, You are doing amazing stuff. Keep it up. Is there some reason that we can’t just be shown a graph of the actual proxies results without CPS, EIV or any other selection process? I did read your original blog on how the data gets chosen. But I would really like to just see all the data in it’s original form as a graph. Thanks for your efforts.

## Jeff Id said

The original data has values as high as 3000 (tree ring) and others around 1 (pre-calibrated isotope). The only way to look at it together is to standardize. With all 1209 proxy series on a single graph it is just a mess.

There is an answer though. On Climate Audit, (Link on the right) Steve McIntyre made gif based movies of the proxies which play one after another on the screen in standardized units, he even had a way to download the movie file.

The posts were around the beginning of September08 and there were several of them. It was a bit jaw dropping for me.

I hope that helps.

## Jeff Id said

R script used to generate above graph. If anyone knows the correct way to replace chronA,B,C,D into an array like chron[i], it would be appreciated.

———————

library(gplots)

source(“http://www.climateaudit.org/scripts/utilities.txt”)

download.file(“http://www.climateaudit.org/data/mann.2008/mann.tab”,”temp.dat”,mode=”wb”);load(“temp.dat”)

download.file(“http://www.climateaudit.org/data/mann.2008/details.tab”,”temp.dat”,mode=”wb”);load(“temp.dat”)

codelist=NULL

for(j in 1:10000)

{

temp=(details$code==j)

K=sum(temp)

if(K>0)

{

codelist[length(codelist)+1]=j

}

}

for(j in 1:length(codelist))

{

temp=header$Data.type.code==codelist[j]

K=sum (temp)

index=(1:1209)[temp]

chronA=NULL

chronB=NULL

chronC=NULL

chronD=NULL

for ( i in 1:K)

{

X=g(mann[[ index[i] ]]) ##convert to time series

temp2=(time(X)>=1850)&(time(X)<=1995) ##collect mann proxy calibration into temp

m0=mean(X[temp2],na.rm=T)

sd0=sd(X[temp2],na.rm=T) ##Calc mean and standard deviation

Xscaled=(scale(X[],center=m0,scale=sd0)) ##Linear scaling of temp

Xscaled=ts(Xscaled,start=time(X)[1]) ##Cram it back into time series

if(i<250)

{

chronA=ts.union(chronA,Xscaled ) ##stick everything into time series array

}

else

{

if(i<500)

{

chronB=ts.union(chronB,Xscaled ) ##stick everything into time series array

}

else

{

if(i<750)

{

chronC=ts.union(chronC,Xscaled ) ##stick everything into time series array

}

else

{

chronD=ts.union(chronD,Xscaled ) ##stick everything into time series array

}

}

}

}

if(length(chronB)!=0)

{

chronA=ts.union(chronA,chronB)

}

if(length(chronC)!=0)

{

chronA=ts.union(chronA,chronC)

}

if(length(chronD)!=0)

{

chronA=ts.union(chronA,chronD)

}

mannavg=ts(apply(chronA,1,mean,na.rm=T),start=tsp(chronA)[1]) ##average time series

Y=mannavg;##//ts.union(annual,mannavg) ##cram time series into Y

f=function(x) filter.combine.pad(x,truncated.gauss.weights(11) )[,2] ##filter functionh

plot(time(Y),f(Y),col=2,type=”l”,lwd=2,las=1,ylab=””,xlab=”Year”,xlim=c(0,2000), ylim = c(-5,5))

title(paste(“Average of Proxy Type =”,codelist[j]))

##savePlot(paste(“c:/agw/r/Average of Proxy Type =”,codelist[j]),type=”jpg”)

}

It should run as is.

## Jeff Id said

When I pasted it the tabs disappeared. Makes it tough to read.

## Andy said

Thanks Jeff, if only I could be where you are, I’m afraid my non-existant computer/maths skills are dwarfed by those here and on other blogs.

Keep up the good work.

## Demesure said

Jeff,

I can’t get your code running (“header$Data.type.code”, “g(mann[[ index[i] ]])” and “filter.combine.pad” don’t run) so I tweaked it my way. I don’t know why your chronA..D so don’t know where I can help.

Now things are much clearer to me : (see my graphs here

http://pichuile.free.fr/ecad/mann2008 )

You see strange things when data are not filtered in code 5000 and 8001. And of course, that 927 proxies of 1209 are tree rings.

library(gplots)

source("http://www.climateaudit.org/scripts/utilities.txt")

download.file("http://www.climateaudit.org/data/mann.2008/mann.tab","temp.dat",mode="wb");load("temp.dat")

download.file("http://www.climateaudit.org/data/mann.2008/details.tab","temp.dat",mode="wb");load("temp.dat")

proxiname=c("Miscellaneous","Composite","Marine sediments","Miscellaneous","Miscellaneous",

"Miscellaneous","Miscellaneous","Miscellaneous" ,"Miscellaneous" ,"Miscellaneous" ,"Miscellaneous" ,

"Tree latewood density","Oxygen isotopes","Miscellaneous","Tree width")

codelist=as.numeric(levels(details$code))

for(j in 1:length(codelist))

{

y=mann[details$code== codelist[j]]

chronA=NULL

for ( dat in y)

{

X=ts(dat[,2], start=dat[1,1]) ##convert to time series

temp2=(time(X)>=1850)&(time(X)<=1995) ##collect mann proxy calibration into temp

m0=mean(X[temp2],na.rm=T)

sd0=sd(X[temp2],na.rm=T) ##Calc mean and standard deviation

Xscaled=(scale(X[],center=m0,scale=sd0)) ##Linear scaling of temp

Xscaled=ts(Xscaled,start=time(X)[1]) ##Cram it back into time series

chronA=ts.union(chronA,Xscaled ) ##stick everything into time series array

}

`Y=ts(apply(chronA,1,mean,na.rm=T),start=tsp(chronA)[1]) ##average time series`

plot(Y,col=2,type="l",las=1,xlab="Year",ylab="values scaled to 1 sigma mean",xlim=c(0,2000), ylim = c(-3,3))

title(paste(proxiname[j]," (code=",codelist[j], ") / ",length(y)," proxies",sep=""))

abline(h=seq(-3,3,1), v=seq(0,2000,500),col="gray", lty="dotted", add=T)

lines(Y,col=2)

savePlot(paste("c:/agw/Average of Proxy Type =",codelist[j]),type="jpg")

}

#

## Jeff Id said

I finally got a chance to play around with your code. Much cleaner. The chronA..D is a method for increasing the speed of the union for the tree ring proxies. It cuts the time down approximately 2^4 or 16 times. I just figured out how to dim the array beforehand. Works great.

Thanks for the example.