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.
the Air Vent 
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?
Andy,
You are exactly where I was not that long ago. I have a link which answers some of your questions.
https://noconsensus.wordpress.com/2008/09/11/ten-things-everyone-should-know-about-the-global-warming-hockey-stick/
I hope you like it.
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" ...
@3 I meant “don’t forget to enclose it with the tag CODE “.
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.
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.
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.
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.
When I pasted it the tabs disappeared. Makes it tough to read.
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.
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")
}
#
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.