Hockey Stick CPS Revisited – Part 1

Updated to include Dr. Mann’s words.

—————————-

The last time I did this, R was a brand new language to me. After 6 months messing around in my free time I speak rudimentary R with a C accent. R is a totally free language that anyone can download and learn. This post is a demonstration of the methods behind the Mann08 CPS hockey stick reconstruction. The difference between this and the numerous posts I did before is better R programming and a lot more comments in the code. If you’re serious about understanding vs advocating, you can figure this out. There is not one person I can think of who has ever commented here, incapable of figuring this out.

CPS is composite plus scale, which is an invented method for calibrating proxies to measured temperatures in paleoclimatology reconstructions. In paleoclimatology methods are too often invented to find the signal in the noise – this is not a new problem and it stems from the large signal to noise ration of paleo-data. If you happen to be a paleoclimatologist who does temperature reconstructions, please try your methods on ARIMA data with a known signal before employing it on whatever your proxy is.

I have hundreds of new readers, who didn’t get the day by day experience of my discovery of climatology math. Well some of my early work was a little rough, however it was correct and the specifics still stand uncriticized.

This post was prompted by some people in blogland (despite the complete lack of rational criticism) claiming that my demonstrations of the CPS hockey stick math is faulty rather than the actual hockey stick itself. An oddly reversed situation which could only exist in the new progressive anti-world. Actually, I don’t recall any real criticism of the method or the result other than statements around the AGW crowd that – ‘it’s been proven wrong’.

I’ve cleaned up the code and seriously over-commented it in the hopes that honest people will be able to understand what is going on here. I’ll do a short explanation for it but the primary explanation is in the programming code. To understand it fully you should to read it step by step and run it.

Michael Mann knows full well that this result exists, his explanation is as follows from an RC thread.

Actually, this line of attack is even more disingenuous and/or ill-informed than that. Obviously, if one generates enough red noise surrogate time series (especially when the “redness” is inappropriately inflated, as is often done by the charlatans who make this argument), one can eventually match any target arbitrarily closely.

You can note that this post uses his Actual data rather than red noise data. Dr. Mann’s continued explanation is here.

What this specious line of attack neglects (intentionally, one can safely conclude) is that a screening regression requires independent cross-validation to guard against the selection of false predictors. If a close statistical relationship when training a statistical model arises by chance, as would be the case in such a scenario, then the resulting statistical model will fail when used to make out-of-sample predictions over an independent test period not used in the training of the model. That’s precisely what any serious researchers in this field test for when evaluating the skillfulness of a statistical reconstruction based on any sort of screening regression approach.

So a potentially overstated statistical check is what determines if CPS works. This is in fact false, CPS is incapable of recovering an accurate signal from data. A fact which I will demonstrate in these next few posts. This post however, does not address the statistical validation issue, it does however demonstrate that Dr. Mann is correct that any signal at all can be made using Mann08 math and data. A truth for many of his hockey stick creation methods.

An explanation of Composite Plus Scale (CPS).

Proxy Data

The proxy data is noisy, the noise in proxies comes from other factors than temperature. Tree ring widths are a primary proxy in the Mann08 reconstruction and don’t only grow wider (and narrower) from temperature they also grow according to humidity, bugs, lightning, soil condition, sunlight CO2 and numerous other factors. This means that tree ring widths are naturally noisy data. The same is true for mollusk shells, interpretations of historic documents and a pile of other proxies.

Correlation and deletion.

In engineering, when you have noisy data you might just take an average or an area weighted average and hope you have enough data to make it work. This makes sense but you may choose to try other techniques to extract a signal based on foreknowledge of what created the noise, filtering, PCA, ICA and such. In this case with tree rings and mixed odd proxies, there is very little that can be done which can improve on averaging and filtering but whatever you do you don’t want to use a biased method for signal hunting.

If you are using a known biased signal extracting method as an engineer, and you want your product to actually work, you had better test it. I demonstrate here that Mann 08’s CPS method can pull any signal it wants from noisy proxy data. The CPS method which is likely only used in climatology, naturally distorts an artificial signal in the data but besides that it extracts the best match for whatever shape you wish to find.

Mann starts the process prior to CPS and measures correlation (likeness) of the proxy data to instrumental temperature data and deletes the proxies that don’t follow the measured temperature curve. This is very important, only a few of proxies (Luterbacher) used are demonstrated to be in any way related to temperature and they weren’t really proxies, but rather contain actual temperature data. Really a hard situation to imagine for an engineer.

Of the 1,209 proxy records in the full dataset, 484 (40%) pass the temperature-screening process over the full (1850–1995) calibration interval – Mann08

So we have noisy data which may or may not be related to temperature and in the case of 71 Luterbacher proxies is temperature data and to sort it we’re going to apply Mannian CPS math to find the hoped for signal hidden in the noise.

Step 1

Measure correlation between the most recent years of proxy data and the most recent years of temperature data. Throw away everyting that doesn’t at least weakly correlate.

– this step in Mann08 deleted 725 of 1209 proxys.

Step 2

All proxies are measured in their own units (i.e. millimeters or parts per million) so you need to recenter (change the mean of) each proxy in the calibration range to match the mean of the measured data.

Step 3

Rescale each proxy so that the standard deviation of the accepted proxies from step 1 match the standard deviation of the calibration signal. Standard deviation is basically how much a graph wiggles up and down around the average. So this matches the amplitude of the wiggles of something which may or may not be temperature to the amplitude of the wiggles of the measured temperature.

Step 4

Take the average.

———

Sounds simple, when I read about step 1 I actually didn’t sleep that night. Throw away the data that doesn’t fit your pre-determined conclusion of a match to temperature. Really, I was so mad I didn’t sleep, remember we’re paying big $$ for this stuff and we must not pretend it’s somehow ok. NOT good.

After the rescaling operation which makes the best possible fit of each proxy to the temperature signal and averaging, the calibration range signal is always found.

Here are some of the signals I was able to extract from Mann08 actual proxy data using their own methods.

Unprecidented Warming

Unprecidented Cooling

And just to really drive the point home.

Unprecidented Sine ----Waving

Inverted Sine Wave

Now I’ve shown that before, but the real point of all this is the code. I’ve cleaned it up and commented the code very thoroughly. It is fully turnkey and anyone interested can put it into R copy it into the script window and run it. Thanks to Steve McIntyre for writing the download (a long time ago) and hosting the data.

Now this is part 1 of the demonstration. Part 2 will uses synthetic ARMA data to demonstrate the distortion of any signal in the data by this method. That post is now available at this link: Historic Hockey Stick – Pt 2 Shock and Recovery

===============================================================

DOWNLOAD – R – HERE — http://www.r-project.org/

DOWNLOAD CODE mann-cps-demo THIS DOWNLOAD AVOIDS FORMATTING PROBLEMS FROM WORDPRESS FOR R

All you need to do is copy and paste from the word document above.

Code is also here below reformatted by wordpress:

##TEMPERATURE CORRELATION
#this calculates statlist with 6 items in list

###########################################################################
##################### TURNKEY DATA DOWNLOADER #############################
###########################################################################
#Mann PRoxy Data and Info
method=”online”
if (method==”offline”) {load(“d:/climate/data/mann08/mann.tab”)
load(“d:/climate/data/mann08/mannorig.tab”)
load(“d:/climate/data/mann08/idorig.tab”)
names(mannorig)=idorig$idorig
load(“d:/climate/data/mann08/details.tab”)
source(“d:/climate/scripts/mann.2008/instrumental.txt”)} else {
url=”http://www.climateaudit.info/data/mann.2008″
download.file(file.path(url,”mann.tab”),”temp.dat”,mode=”wb”);load(“temp.dat”)
download.file(file.path(url,”mannorig.tab”),”temp.dat”,mode=”wb”);load(“temp.dat”)
download.file(file.path(url,”details.tab”),”temp.dat”,mode=”wb”);load(“temp.dat”)
download.file(file.path(url,”idorig.tab”),”temp.dat”,mode=”wb”);load(“temp.dat”)
names(mannorig)=idorig$idorig
source(“http://www.climateaudit.info/scripts/mann.2008/instrumental.txt”)}
############################################################################
################ FUNCTIONS ################################################
############################################################################

######## CREATE TIME SERIES ##########
gg = function(X)
{
ts(X[,2],start=X[1,1])
}

######## GAUSSIAN FILTER FUNCTION TESTED AT CA #########
######## set to 11 years ###############################
source(“http://www.climateaudit.info/scripts/utilities.txt”)

ff=function(x)
{
filter.combine.pad(x,truncated.gauss.weights(11) )[,2]
}
##########################
##### END FUNCTIONS ######
##########################

library(gplots)
use0=”pairwise.complete.obs”

#variable names
###########################################################################
########### mannorig is 1357 infilled series,,, ###########################
########### mann is the 1209 series after RegEM infilling #################
###########################################################################
####### THIS SECTION LETS YOU USE INFILLED OR NON-INFILLED DATA #####
#######try either of these for fun. comment out EITHER LINE#######

mano=mannorig #copy non-infilled PROXY series to different memory
#mano=mann #copy REGEM infilled PROXY series to different memory

#####################################################################
#####################################################################

lent=length(mano) ##lent = number of series chosen above

#filter all series using gaussian CA filter

for(i in 1:lent)
{
mano[[i]][,2]=ff(mano[[i]][,2])
}

#apply functions to look for change values in this section
#change values here and run from here down only

#############################################################################################
#############################################################################################
# MANN O8 MATH ALLOWS YOU TO FIND WHATEVER PATTERN YOU WANT IN PROXY DATA
# THIS SECTION GENERATES DATA REPRESENTING WHAT YOU WANT TO FIND
# I’VE PROVIDED SEVERAL POSSIBLE FUNCTIONS YOU CAN CHOOSE FROM BY DELETING THE ‘#’ SIGNS
# I’VE ALSO PROVIDED A SECTION BY WHICH YOU CAN SET THE PARAMETERS FOR EACH FUNCTION
# YOU CAN MAKE YOUR OWN POSITIVE AND NEGATIVE HOCKEY STICKS FROM THE ACTUAL PROXY DATA
# FOLLOW THE INSTRUCTIONS AND YOU CAN DO IT YOURSELF
#############################################################################################

######### CHANGE THESE VALUES TO DETERMINE THE CALIBRATION PERIOD
######### THESE VALUES ARE EQUIVALENT TO GROUND TEMPERATURES IN M08
######### INSTEAD OF TEMPERATURES WE CAN LOOK FOR ANYTHING AND STILL FIND IT
######### “CALIBRATION PERIOD” DETERMINES THE FUNCTION TIME PERIOD YOU ARE LOOKING FOR IN YEARS
######### APPROPRIATE VALUES ARE INCLUDED IN THE REMARKS AS SUGGESTIONS FOR EACH LINE
#########
######### AFTER TOP IS LOADED RUN FROM HERE DOWN
endyr=1995 #end year for function correlation 100-2000 (INTEGER)
leng=200 #lenght of function to correlate in years 50 – 1000 (INTEGER)
cycles=4 #cycles of sine wave (ANY)
##################

################## FUNCTIONS CHOOSE ONE BY COMMENTING OUT OTHER 3 WITH ‘#’
################## TRY THESE AND MAKE YOUR OWN — ‘FAKE CALIBRATION TEMPERATURE DATA’
#m=c(1:leng)/-leng ###### 0 to -1 NEGATIVE linear slope
#m=c(1:leng)/leng ###### 0 to 1 POSITIVE linear slope
#m=sin( (1:leng)/(leng/cycles)*pi) #sin wave with # cycles
m=-sin( (1:leng)/(leng/cycles)*pi) #-sin wave with # cycles
##################

#############################################################################################
# THIS LINE DETERMINES HOW WELL EACH INDIVIDUAL PROXY MUST CORRELATE TO THE FUNCTION
# CORRELATION IS ALWAYS BETWEEN 0 AND 1
# CORRELATION VALUES GREATER THAN THIS NUMBER ARE ACCEPTED
# MANN 08 USED 0.1 HIGHER NUMBERS MEAN LESS PROXIES ACCEPTED AND BETTER FINAL FIT
# LOWER NUMBERS MEAN MORE ACCEPTED AND WORSE FINAL FIT

r=.3 #sorting correlation criteria -1 TO 1 POSITIVE VALUES MAKE MORE SENSE

#############################################################################################
#############################################################################################
#############################################################################################

leng=length(m) ## RECALCULATES leng SO THAT PEOPLE WHO MESS AROUND WITH THE FUNCTIONS ABOVE DON’T HAVE TROUBLE
m=ts(m,endyr-leng) ## CONVERT FAKE CALIBRATION TEMPERATURE DATA TO TIME SERIES STARTING AT endyr MINUS LENT

corv = array(0,lent) #set corv AS A VECTOR OF LENGTH EQUAL TO NUMBERR OF PROXIES

stat = rep(NA,lent) #INITIALIZE STAT VECTOR, THANKS ROMAN M

for(i in 1:lent)## CYCLE BETWEEN { } ONE TIME FOR EACH PROXY
{
y=gg(mano[[i]]) ## assign y as a timeseries version of each individual proxy
y=ts.union(m,y) ## make a union – two column proxy of the fake temperature function and the individual proxy
val = window(y, start(m)[1],end=end(m)[1])
stat[i]=cor(val,use=use0)[1,2]
}

pass= stat > r # BUILD VECTOR OF PASS FAILS, ONE ELEMENT FOR EACH PROXY .. TRUE == GREATER THAN r

#############################################################################################
# NOW WE HAVE SELECTED OUR PROXIES BASED ON HIGHEST CORRELATION TO THE CHOSEN
# ‘FAKE CALIBRATION TEMP DATA’
#############################################################################################
#############################################################################################
#######NOW WE CALCULATE THE MEAN AND STANDARD DEVIATION FOR EACH PROXY AND SUBTRACT AND DIVIDE
#######THIS IS THE BUSINESS PORTION OF MANN 08 CPS – COMPOSITE PLUS SCALE
#######calculate SD for sorting function m — used in CPS to match standard
#######deviation and offset of all proxies to the function sd and offset
#############################################################################################
#############################################################################################

msig=mean (m,na.rm=T) ## CALCULATE THE MEAN OF FAKE CALIBRATION TEMP FUNCTION
ssig=sd(m,na.rm=T) ## CALCULATE THE STANDARD DEVIATON OF FAKE CALIBRATION TEMP FUNCTION

#clear values
y=NULL;
x=NULL;
Z=NULL

for(j in 1:lent)#LOOP THROUGH ALL PROXIES
{
if(pass[j]==TRUE) #ONLY USE PROXIES THAT PASS CORRELATION REMEMBER pass IS A TRUE/FALSE FOR EACH PROXY
{
y=(gg(mano[[j]])) #MAKE Y A TIME SERIES FOR EACH PROXY
val = window(y, start(m)[1],end=end(m)[1]) #TRUNCATE TIME SERIES TO YEARS WHERE FAKE CALIBRATION TEMP FUNCTION EXISTS
#CPS USES ONLY CALIBRATION RANGE TO SCALE AND OFFSET PROXY
#THIS DETAIL IS THE PRIMARY POINT WHERE IT GOES WRONG
m0=mean(val,na.rm=T)# CALCULATE MEAN OF PROXY SERIES IN FAKE CALIBRATION TEMP YEAR RANGE
sd0=sd(val,na.rm=T) # CALCULATE STANDARD DEVIATION OF PROXY SERIES IN FAKE CALIBRATION TEMP YEAR RANGE
################# CPS METHOD IS APPLIED HERE #######################
yscale=y-m0 #CENTER TEMP PROXIES BY SUBTRACTING MEAN OF PROXY IN FK CAL TEMP RANGE
yscale=yscale/sd0 #SCALE INDIVIDUAL PROXY BY DIVIDING SD OF PROXY IN FK CAL TEMP RANGE – more wiggles stronger division
yscale=yscale*ssig #SCALE INDIVIDUAL PROXY BY MULTIPLYING BY SD OF FK CAL TEMP – this scales the proxies to match the wiggles of fk cal temp data
yscale=yscale+msig #CENTER INDIVIDUAL PROXY BY ADDING TO OFFSET MATCH FK CAL TEMP

x=ts.union(x,ts (yscale,time(y)[1])) #ADDS EACH INDIVIDUAL PROXY INTO ONE BIG MATRIX – EACH TIME THROUGH ONE MORE IS ADDED
}
}

adrecon = window(x,start=0,end=1995) #A COUPLE OF PROXIES EXTEND FAR IN THE PAST SO USE YEARS 0 to 1995 TO MATCH MANN 08
Z=rowMeans(adrecon,na.rm=TRUE) #AVERAGE ALL REMAINING SCALED SERIES YEARS TOGETHER PROXIES ARE IN COLUMNS BY YEARS AT THIS POINT SO AVERAGE ROWS
Z=ts(Z,start=0) #MAKE Z INTO A TIME SERIES STARTING AT TIME 0AD

###plot###
tit=paste(“Our Mann 08 CPS Temperature Reconstruction \nUsing Artificial Calibration Data\nPercent Proxies Retained (r >”,r,”) =”,round(sum(pass)/length(pass)*100,2),”%” )
plot(Z,type=”l”,ylim=c(-1.5,1.5),main = tit,xlab=”Year”, ylab=”Degrees C”)
lines(m,col=”red”,lwd=3)
grid()
smartlegend(x=”left”,y=”bottom”,inset=0,fill=1:2,legend=c(“Mann 2008 Actual Proxy Average”,”Fk Calibration Temperature”) ,cex=.7)
#savePlot(“c:/agw/r=0_3 NEG SINE”,type=”jpg”)