#MCINTYRE AND MCKITRICK, E&E (2005) #This shows scripts used to produce figures and calculations. #it is hoped that the outline and the script itself will be reasonably self-explanatory. ########################### #DIRECTORIES AND FUNCTIONS #mbhscripts<-"c:/climate/scripts" #"http://www.climateaudit.org/scripts" mbhscripts<-"http://www.climate2003.com/scripts/MM05_EE" #edited for internet access Oct 2006 source(file.path(mbhscripts,"ee2005.functions.txt") ) #mbhdata<- "c:/climate/data/mann" #"http://www.climateaudit.org/data/MBH98" #commented out OCt 2006 mbhscripts<-"http://www.climateaudit.org/scripts" #edited for internet access OCt 2006 source(file.path(mbhscripts, "MBH98config4.txt")) ########################################## #FIGURE 1 ####################### proxy<-proxy.MBH case<-"MBH98" NH0F<-NHbeta(proxy,case,series=series.MBH,directory=directory.MBH,flag="new",method3="scale.RPC", method4="scale.NH")#use new function emulation.config0<-ts(NH0F[[1]],start=1400,end=1980) #CONFIRM THAT THIS IS THE SERIES ORIGINALLY ARCHIVED MM05b<-read.table("http://www.climate2003.com/data/MM05.EE/emulation.txt",sep="\t",header=TRUE) # this was archived in March 2005 # THIS IS VERSION USED IN EE[2005] #FIGURE 1 PANEL 1 max(MM05b[,2]-NH0F[[1]]) ## 6.661338e-16 #thus the same #COMPARE EMULATION TO ARCHIVED MBH (**NOT** TO NH TEMPERATURE) verification.stats(emulation.config0,NH.MBH,c(1902,1980),c(1400,1901)) # RE.cal RE.ver R2.cal R2.ver CE sign.test prod.test #1 0.772 0.829 0.804 0.677 0.359 406 14 #used on p. 71: Relative to the MBH98 reconstruction, it has a Reduction of Error (“RE”) statistic of 0.83 in the 1400–1901 period (R2 – 0.68) – both hockeystat4(NH0F[[1]]) # 1.53 hockeystat4(NH.MBH) # 1.77 # CASE 2 - ELIMINATE GASPE EXTRAPOLATION proxy[1:4,53]<-NA #eliminate Gaspe extrapolation NH2<-NHbeta(proxy,case,series=series.MBH,directory=directory.MBH,flag="new",method3="scale.RPC", method4="scale.NH")#use new function hockeystat4(NH2[[1]]) # 1.26 # CASE 3 - FIGURE 1 BOTTOM PANEL - ELIMINATE GASPE EXTRAPOLATION AND USE CENTERED COVARIANCE TREE RING PCS case<-"WDCP.persist" NH3<-NHbeta(proxy,case,series=series.MBH,directory=directory.MBH,flag="new",method3="scale.RPC", method4="scale.NH")#use new function max(MM05b[,4]-NH3[[1]]) ##4.996004e-16 c(hockeystat4(NH3[[1]]),hockeystat4(NH3[[3]])) # -0.3474396 -0.5402734 verification.stats(ts(NH3[[3]],start=1400),sparse,c(1902,1980),c(1854,1901))[1:5] # RE.cal RE.ver R2.cal R2.ver CE #1 0.05381303 -0.7983807 0.1804436 4.778862e-05 -3.215624 # CASE 4 - FIGURE 1 MIDDLE PANEL - CENTERED NOAMER TREE RING PCS ONLY proxy[1:4,53]<-proxy[5,53] NH2A<-NHbeta(proxy,case,series=series.MBH,directory=directory.MBH,flag="new",method3="scale.RPC", method4="scale.NH")#use new function max(MM05b[,3]-NH2A[[1]]) ## 4.996004e-16 hockeystat4(NH2A[[1]]) # 0.2860645 # DRAW FIGURE 1 ##smooth series emulation.config0<-filter.combine.pad(emulation.config0,truncated.gauss.weights(25)) emulation.config2<-filter.combine.pad(ts(NH2[[1]],start=1400,end=1980),truncated.gauss.weights(25)) emulation.config3<-filter.combine.pad(ts(NH3[[1]],start=1400,end=1980),truncated.gauss.weights(25)) #postscript(file.path("c:/climate/scripts/ee.2004","figure1.eps"), width = 7.5, height = 3.5, # horizontal = FALSE, onefile = FALSE, paper = "special", # family = "Helvetica",bg="white",pointsize=8) nf <- layout(array(1:3,dim=c(3,1)),heights=c(1.1,1,1.1)) layout.show(nf) par(mar=c(0,4,3,2)) plot(1400:1980,emulation.config0[,1],axes=FALSE,type="l",ylab="",font.lab=2,ylim=c(-.6,.5),col="grey70") lines(1400:1980,emulation.config0[,2]) axis(side=2,las=2,cex.axis=0.75,font=2)#at=seq(-.4,.2,.2), box() text(1415,0.45,"MBH98 Emulation",font=2,cex=.75,pos=4) par(mar=c(0,4,0,2)) plot(1400:1980,emulation.config2[,1],axes=FALSE,type="l",ylab="",font.lab=2,ylim=c(-.6,.5),col="grey70") lines(1400:1980,emulation.config2[,2]) axis(side=2,las=2,cex.axis=0.75,font=2)#at=seq(-.4,.2,.2), box() text(1415,0.45,"Archived Gaspe",font=2,cex=.75,pos=4) # par(mar=c(3,4,0,2)) plot(1400:1980,emulation.config3[,1],axes=FALSE,type="l",ylab="",font.lab=2,ylim=c(-.6,.5),col="grey70") lines(1400:1980,emulation.config3[,2]) axis(side=2,las=2,cex.axis=0.75,font=2)#at=seq(-.4,.2,.2), box() axis(side=1, at = seq(1400,2000,100),labels=TRUE,font=2,las=1) text(1415,0.45,"Archived Gaspe and Centered PCs",font=2,cex=.75,pos=4) # dev.off() #SAVE DATA emulation<-cbind(nhmann[1:581,2],emulation.config0[,1],emulation.config2[,1],NH2A[[1]], emulation.config3[,1]) dimnames(emulation)[[2]]<-c("MBH98","emulation","Gaspe.only","NOAMER.only","corrected") #write.table(emulation,file="c:/climate/scripts/ee.2004/emulation.MM05b.txt",sep="\t",quote=FALSE) ######END FIGURE 1 ################################# ################### ###FIGURE 2 - IMPACT OF +0.5 IN 15th CENTURY #1. REPLICATE MBH98 CALCULATIONS #CALCULATE NOAMER PC1 USING MANN DATA AND SVD noamer.file<-"ftp://holocene.evsc.virginia.edu/pub/MBH98/TREE/ITRDB/NOAMER/BACKTO_1400" noamer.file<-"c:/climate/data/new/TREE/ITRDB/NOAMER/BACKTO_1400" MBH.PC1<-read.table(file.path(noamer.file,"pc01.out")) MBH.PC1<-ts(MBH.PC1[,2],start=1400,end=1980) mbhdata<-"c:/climate/data/mann" tree<-read.table(file.path(mbhdata,"UVA/NOAMER.txt"),sep="\t",header=TRUE) tree<-ts(tree[,2:ncol(tree)],start=tree[1,1],end=tree[nrow(tree),1]) tree0<-tree temp<-!is.na(tree[1,]) tree<-tree[,temp] tree<-tree[1:581,] Tree<-extend.persist(tree) N<-nrow(tree) Tree.mannomatic<-apply(Tree,2,mannomatic) #MAKE MATRIX WITH INCREASED VALUES FROM 1400-1450 #IDENTIFY CENSORED SERIES url<-"c:/climate/data/new/TREE/ITRDB/NOAMER" id<-list.files(url) #175 is censored list; 176 is uncensored list roster1<-read.table(file.path(url,id[175]));length(roster1[,2])#192 roster2<-read.table(file.path(url,id[176]) ); length(roster2[,2])#212 length(roster2[,2][!is.na(match(roster2[,2],roster1[,2]))]) #192 series length(roster2[,2][is.na(match(roster2[,2],roster1[,2]))]) #20 series roster2[,2][is.na(match(roster2[,2],roster1[,2]))] #this lists the following 20 series # [1] az510.txt ca528.txt ca529.txt ca530.txt ca533.txt # [6] ca534.txt co511.txt co522.txt co523.txt co524.txt #[11] co525.txt co535.txt co545.txt co547.txt nv510.txt # [16] nv511.txt nv512.txt nv513.txt nv514.txt nv516.txt graybill.list<-as.character(roster2[,2][is.na(match(roster2[,2],roster1[,2]))]) graybill.logical<-!is.na(match(dimnames(Tree)[[2]],graybill.list)) Tree.adj<- Tree #creates mirror object for testing Tree.adj[1:51,!graybill.logical]<- Tree.adj[1:51,!graybill.logical]+0.5 # adds 0.5 to all non-Graybill sites (mean is 1) Tree.adj.mannomatic<-apply(Tree.adj,2,mannomatic) ##Oct 2006 comment - the purpose of this diagram was to illustrate the perverse impact of the MBH method #in the case where there was an aribtrary increment to non-bristlecone chronologies #in all other calculations where we tried to illustrate the effect of conventional chronologies, #we used principal components algorithms princomp or prcomp #here in order to match scales, we illustrated by using svd on the original matrix standardized to a mean of 1000 #an alternate version using prcomp is added in below PCA.mannomatic.svd<-svd(Tree.mannomatic,nu=9,nv=9) #max(PCA.mannomatic.svd$u[,1]-MBH.PC1)## 5.955034e-08 PC1.mannomatic<-ts(PCA.mannomatic.svd$u[,1],start=1400) PCA.adj.mannomatic<- svd(Tree.adj.mannomatic)#svd on data matrix PC1.adj.mannomatic<- ts(PCA.adj.mannomatic$u [,1],start=1400) #selects PC1 from svd model PC1.base<- ts( svd(Tree)$u[,1], start=1400) PC1.adj.base<- ts( svd(Tree.adj)$u[,1], start=1400) f<-function(x){f<-filter.combine.pad(x,truncated.gauss.weights(21))[,2];f} par(mar=c(3,3,1,1)) ylim0<-c(-.065,.065) plot(1400:1980,f(PC1.mannomatic),type="l",xlab="",ylab="",ylim=ylim0) lines(1400:1980,-f(PC1.adj.mannomatic),lty=2) lines(1400:1980,-f(PC1.base),col="grey50") lines(1400:1980,-f(PC1.adj.base),col="grey50",lty=2) mean(PC1.adj.mannomatic)# 0.03574446 mean(PC1.adj.base)#-0.04100004 mean(PC1.base)# -0.0412064 ## OCt 2006: alternate version of Figure 2 using prcomp PCA.mannomatic<-prcomp(Tree.mannomatic,center=FALSE) PC1.mannomatic.rev<-ts(PCA.mannomatic$x[,1],start=1400) PCA.adj.mannomatic<- prcomp(Tree.adj.mannomatic,center=FALSE)#svd on data matrix PC1.adj.mannomatic.rev<- ts(PCA.adj.mannomatic$x [,1],start=1400) #selects PC1 from svd model PC1.base.rev<- ts( prcomp(Tree)$x[,1], start=1400) PC1.adj.base.rev<- ts( prcomp(Tree.adj)$x[,1], start=1400) f<-function(x){f<-filter.combine.pad(x,truncated.gauss.weights(21))[,2];f} par(mar=c(3,3,1,1)) ylim0<-c(-12,4) plot(1400:1980,f(PC1.mannomatic.rev),type="l",xlab="",ylab="",ylim=ylim0) lines(1400:1980,-f(PC1.adj.mannomatic.rev),lty=2) #note that PC1.adj is flipped to match lines(1400:1980,f(PC1.base.rev),col="grey50") lines(1400:1980,-f(PC1.adj.base.rev),col="grey50",lty=2) mean(PC1.adj.mannomatic.rev)#5.7787 mean(PC1.adj.base.rev)#5.26954e-17 mean(PC1.base.rev)# -6.356149e-17 cor(-PC1.base.rev,PC1.base)# 0.8807274 combine<-cbind(-PC1.base,PC1.base.rev) dimnames(combine)[[2]]<-c("Figure 2", "prcomp") plot.ts(combine,main="") #end update ######## END FIGURE 2 ##PERMUTATIONS AND COMBINATIONS #CASE 5- USE CENSORED MBH98 PCs AND ARCHIVED GASPE #this is done by setting NOAMER.CENSORED=TRUE proxy<-proxy.MBH proxy[1:4,53]<-NA NH5<-NHbeta(proxy,"MBH98",NOAMER.CENSORED=TRUE, series=series.MBH,directory=directory.MBH,flag="new",method3="scale.RPC", method4="scale.NH")#use new function hockeystat4(NH5[[1]]) # -0.1965684 #this is similar to the impact of using covariance PCs #rm(NH5) ##this looks very similar to the base case cor(NH0F[[1]],NH5[1])# 0.8886202 verification.stats(ts(NH0F[[3]],start=1400),sparse,c(1902,1980),c(1854,1901)) # RE.cal RE.ver R2.cal R2.ver #1 0.1103655 0.3378769 0.2125288 0.01851249 verification.stats(ts(NH5[[3]],start=1400),sparse,c(1902,1980),c(1854,1901)) # RE.cal RE.ver R2.cal R2.ver #1 0.06507186 -0.8050038 0.1866221 0.002099475 #CASE 6 - USE 5 NOAMER series in AD1400 and AD1450 STEPS #this is done by defining new series schedule retaining 5 PCs while using WDCP.persist directory of covariance PCs proxy<-proxy.MBH proxy[1:4,53]<-NA series.new<-series.MBH series.new[3,10:11]<-5 #retain 5 NOAMER PCs NH6<-NHbeta(proxy,"WDCP.persist", series=series.new,directory=directory.MBH,flag="new",method3="scale.RPC", method4="scale.NH")#use new function c(hockeystat4(NH6[[1]]),hockeystat4(NH6[[3]])) #[1] 0.7790334 0.5165209 #this is similar to the impact of using covariance PCs verification.stats(ts(NH6[[3]],start=1400),sparse,c(1902,1980),c(1854,1901)) #RE.cal RE.ver R2.cal R2.ver #1 0.1307584 0.02364327 0.2247425 0.001648292 ###WA case 5(c) 4 PCs - RE.cal 0.29; RE ver: 0.14 #CASE 7 - SAME AS CASE 6, BUT WITHOUT THE PC4 #this is done by setting bristlecone.tag=TRUE which excludes the PC4 NH7<-NHbeta(proxy,"WDCP.persist",bristlecone.tag=TRUE, series=series.new,directory=directory.MBH,flag="new",method3="scale.RPC", method4="scale.NH")#use new function c(hockeystat4(NH7[[1]]),hockeystat4(NH7[[3]])) #[1] -0.1426421 -0.2920432 #this is almost the same as the censored PCs verification.stats(ts(NH7[[3]],start=1400),sparse,c(1902,1980),c(1854,1901))[1:4] # RE.cal RE.ver R2.cal R2.ver #1 0.07021073 -0.6462009 0.1894767 6.904155e-08 #CASE 8 - SAME AS CASE 7 (WITHOUT THE PC4), BUT WITH 7 PCs in NOAMER EARLY PART #this is done by setting bristlecone.tag=TRUE which excludes the PC4 series.new<-series.MBH series.new[3,9:11]<-7 NH8<-NHbeta(proxy,"WDCP.persist",bristlecone.tag=TRUE, series=series.new,directory=directory.MBH,flag="new",method3="scale.RPC", method4="scale.NH")#use new function c(hockeystat4(NH8[[1]]),hockeystat4(NH8[[3]])) #[1] -0.2366295 -0.3065319 #this is almost the same as the censored PCs verification.stats(ts(NH8[[3]],start=1400),sparse,c(1902,1980),c(1854,1901))[1:4] # RE.cal RE.ver R2.cal R2.ver #1 0.1182854 -0.9948159 0.2172316 0.002296584 #CASE 9 - CORRELATION PCs #this is done by using pre-calculated PCs in subdirectory "WDCP.correlation" proxy<-proxy.MBH proxy[1:4,53]<-NA NH9<-NHbeta(proxy,"WDCP.correlation", series=series.MBH,directory=directory.MBH,flag="new",method3="scale.RPC", method4="scale.NH")#use new function c(hockeystat4(NH9[[1]]),hockeystat4(NH9[[3]])) #[1] 0.6185650 0.4253359 #this is similar to the impact of using covariance PCs verification.stats(ts(NH9[[3]],start=1400),sparse,c(1902,1980),c(1854,1901))[1:4] # RE.cal RE.ver R2.cal R2.ver #1 0.2035161 0.02208175 0.2710988 0.004947316 #this corresponds to WA case 5(b) RE cal 0.32; RE ver 0.18 # MAKE BIG PROXY dataversion.MBH<-rep("UVA",6) #or "MBH98" or "UVA" method<-"persist" #or "remaining" #step is denominated as 1-11 tree.bigproxy<-function(dataversion=dataversion.MBH) { tree.bigproxy<-NULL id<-NULL for (m in 1:6) { loc<-file.path (mbhdata,dataversion[m], paste( regionid[m],"txt",sep=".") ) tree<-read.table(loc,header=TRUE) tree<-ts(tree[,2:ncol(tree)],start=tree[1,1],end=tree[nrow(tree),1]) idtemp<-dimnames(tree)[[2]] tree.bigproxy<-ts.union(tree.bigproxy,tree) id<-c(id,idtemp) } dimnames(tree.bigproxy)[[2]]<-id tree.bigproxy } tree.big<-tree.bigproxy(dataversion.MBH) dim(tree.big) #[1] 985 334 tsp(tree.big) # 1011 1995 1 tree.big<-ts(tree.big[(1400-1010):(1980-1010),],start=1400) tree.big<-extend.persist(tree.big) load("c:/climate/data/mann/proxy.tab") proxy[1:4,53]<-NA proxy<- ts( proxy[1:581,c(1:68,100:112)],start=1400) proxy<-extend.persist(proxy) big.proxy<-ts.union(proxy,tree.big) dim(big.proxy) # 581 415 dimnames(big.proxy)[[2]]<-c(dimnames(proxy)[[2]],dimnames(tree.big)[[2]]) NHno1<-NHnoPC(big.proxy,flag="new",method3="scale.RPC", method4="scale.NH") c(hockeystat4(NHno1[[1]]),hockeystat4(NHno1[[3]])) #[1] 1.946827 1.960265 #this is similar to the impact of using covariance PCs verification.stats(ts(NHno1[[3]],start=1400),sparse,c(1902,1980),c(1854,1901))[1:4] # RE.cal RE.ver R2.cal R2.ver # 0.1320199 0.359942 0.2255092 0.002976207 #DO with 95 temp<-!is.na(big.proxy)[1,];sum(temp) dimnames(big.proxy)[[2]][temp] #this lists 109 #the WA uses only 95 #we use 94 because Gaspe not in test<-(1:415)[temp][95:109];test big.proxy[1:51,test]<-NA #set to NA to reduce to 95 NHno2<-NHnoPC(big.proxy,flag="new",method3="scale.RPC", method4="scale.NH") c(hockeystat4(NHno2[[1]]),hockeystat4(NHno2[[3]])) #[1] 2.166292 2.012479 #this is similar to the impact of using covariance PCs verification.stats(ts(NHno2[[3]],start=1400),sparse,c(1902,1980),c(1854,1901))[1:4] # RE.cal RE.ver R2.cal R2.ver #1 0.1624460 0.4779063 0.2443978 0.009379043 ##ALTERNATE GRAYBILL LIST load("c:/climate/data/tree/northamerica.details.tab") id0<-strip(id[19:length(id)],2) id0<-c(id[1:18],id0) temp<-!is.na(match(details$id,id0)); sum(temp) #69 fred<-substr(as.character(details$author[temp]),1,6) temp1<-(fred=="DONALD"); sum(temp1) #26 test<-as.character(details$id[temp][temp1]) temp1<-!is.na(match(dimnames(big.proxy)[[2]],paste(test,"txt",sep=".")));sum(temp1) #26 NHno3<-NHnoPC(big.proxy[,!temp1],flag="new",method3="scale.RPC", method4="scale.NH") c(hockeystat4(NHno3[[1]]),hockeystat4(NHno3[[3]])) #[1] 0.7378430 0.5727606 #this is similar to the impact of using covariance PCs verification.stats(ts(NHno3[[3]],start=1400),sparse,c(1902,1980),c(1854,1901))[1:4] # RE.cal RE.ver R2.cal R2.ver #1 0.2490273 -0.7025898 0.3023037 0.03716841 temp<-!is.na(big.proxy[1,]) #95 id<-dimnames(big.proxy)[[2]][temp] #one australian data.frame(id,round(fred1$weights[temp,11]*1000,2)) h3<-apply(big.proxy[,temp],2,hockeystat4) plot(h3,fred1$weights[temp,11]) order1<-order(fred1$weights[temp,11],decreasing=TRUE) id[order1] #### [1] "ca529.txt" "co545.txt" "nv512.txt" "ca528.txt" [5] "ca532.txt" "arge013.txt" "ca534.txt" "nv513.txt" [9] "nv511.txt" "az510.txt" "co523.txt" "co525.txt" [13] "ca530.txt" "V14" "V12" "V43" [17] "V67" "ga002.txt" "co524.txt" "nv060.txt" [21] "V47" "nv061.txt" "ar053.txt" "V20" [25] "nv514.txt" "co522.txt" "ca531.txt" "ca087.txt" [29] "ca065.txt" "or009.txt" "az550.txt" "swmxdfew11.dat" [33] "swmxdfew09.dat" "co067.txt" "V64" "nv510.txt" [37] "swmxdflw11.dat" "ca533.txt" "nm025.txt" "nm026.txt" [41] "swmxdflw01.dat" "co511.txt" "swmxdfew01.dat" "V63" [45] "V109" "wy023x.txt" "nc008.txt" "co535.txt" [49] "ar049.txt" "wy006.txt" "V62" "nm572.txt" [53] "nv515.txt" "ut509.txt" "swmxdflw09.dat" "nm560.txt" [57] "ga003.txt" "or012.txt" "nv517.txt" "ga004.txt" [61] "V68" "V66" "nv037.txt" "ca073.txt" [65] "nm559.txt" "ar050.txt" "la001.txt" "or015.txt" [69] "V13" "ca535.txt" "sd017.txt" "va021.txt" [73] "ut508.txt" "V102" "ut023.txt" "co076.txt" [77] "nv056.txt" "co509x.txt" "ar052.txt" "V18" [81] "mt006.txt" "V103" "co509.txt" "nv049.txt" [85] "V108" "nv053.txt" "V15" "wy005b.txt" [89] "az086.txt" "az082.txt" "ca084.txt" "nv516.txt" [93] "ca555.txt" "co547.txt" "wy023.txt" temp<-!is.na(big.proxy[1,]) #95 temp1<-!is.na(match(dimnames(big.proxy)[[2]],graybill.list));sum(temp1) temp2<-temp&!temp1;sum(temp2) #75 fred2<-NHnoPC(big.proxy[,!temp1],flag="new",method3="scale.RPC", method4="scale.NH") hockeystat4(fred2[[1]]) #[1] 2.151194 data.frame(dimnames(big.proxy)[[2]][temp2],round(fred1$weights[temp2,11]*1000,2)) h3<-apply(big.proxy[,temp2],2,hockeystat4) plot(h3,fred1$weights[temp2,11]) order1<-order(fred1$weights[temp2,11],decreasing=TRUE) id[order1] ############# ###ROBUSTNESS CALCULATIONS ########### ##WA UPDATE #Only 14 such chronologies account for over 93% of the variance in the PC1 #see GRL #The PC1 in turn accounts for 38% of the total variance. #see GRL #In a centered calculation on the same data, the influence of the bristlecone pines drops to the PC4 #see GRL #The PC4 in a centered calculation only accounts for only about 8% of the total variance, #see GRL #If a centered PC calculation on the North American network is carried out (as we #advocate), then MM-type results occur if the first 2 NOAMER PCs are used in the #AD1400 network (the number as used in MBH98) #MBH-type results occur if the NOAMER network is expanded to 5 PCs in the AD1400 segment #If de-centered PC calculation is carried out (as in MBH98), then MM-type results #still occur regardless of the presence or absence of the PC4 if the bristlecone pine #sites are excluded, ##C. CENSORED CALCULATION #show impact of new information on PC retention rosters count.directory<-"treepc.byregion.count.May2004.txt" #this is new information loc.directory<-"treepc.byregion.directory.new.txt" #this is new information bristlecone.tag<-FALSE NOAMER.CENSORED<-TRUE #this uses CENSORED directory NH2<-NHnew(proxy,case,count.directory,loc.directory,bristlecone.tag)#use new function par(mfrow=c(3,1)) plot.ts(MBH.NH) plot.ts(NH0[[1]]) plot.ts(NH2[[1]]) #MBH-type results occur if bristlecone pine sites (and PC4) are included #If the data are transformed as in MBH98, but the principal components are #calculated on the covariance matrix, rather than directly on the de-centered data, the #results move about halfway from MBH to MM. #If the data are not transformed #(MM), but the principal components are calculated on the correlation matrix rather #than the covariance matrix, the results move part way from MM to MBH, with #bristlecone pine data moving up from the PC4 to influence the PC2. #If no North American PC1 is used at all in the AD1400 calculations (which occurs #if PC calculations are done over the maximum period in which all sites are #available, as done in MM03), then MM-type results occur under both centered and #decentered PC calculations, with and without bristlecone pines. #If, as is suggested in Mann et al. [2004a, 2004b], no PC calculations are applied to #the North American and Stahle/SWM networks and the sites are instead used as #individual proxies (while otherwise carrying on with MBH98 methods), then #MBH-type results are obtained regardless of whether the Gaspé series is duplicated #or extrapolated. In this case, the MBH temperature reconstruction becomes little #more than an index of bristlecone pine growth. However, if the bristlecone pine #sites are excluded from this network, then MM-type results are obtained. ################# ####### ##FIGURES 3 AND 4 ########################## ###FUNCTIONS corfunction<-function(x) {corfunction<-cor(x[503:581],pc[1:79,1]); corfunction} load("c:/climate/data/mann/proxy.tab") proxy<-ts(proxy[1:581,],start=1400,end=1980) #LOAD MBH98 PC SERIES NOMENCLATURE pcstore<-"c:/climate/data/pcs" count.directory.MBH04<-file.path(pcstore,"treepc.byregion.count.SI.2004.txt") #this is new information loc.directory.MBH04<-file.path(pcstore,"treepc.byregion.directory.SI.2004.txt") #this is new information count.directory<-count.directory.MBH04 loc.directory<-loc.directory.MBH04 count.directory0<-count.directory.MBH04 loc.directory0<-loc.directory.MBH04 series<-read.table ( count.directory0,sep="\t",header=TRUE) #table of retained PCs by period and region directory<-read.table(loc.directory0,sep="\t",header=TRUE) #table of look-up directories for above series<-series[2:ncol(series)] directory<-directory[2:ncol(directory)] case<-"MBH98" #read in principal component series #tags for special studies on NOAM.CENSORED and PC4 effect k<-11 for (m in 1:6) { if (series[m,k]>0) {load(file.path(pcstore,case,paste(regionid[m],directory[m,k],"tab",sep=".")))} #this proides for using CENSORED version if (series[m,k]>0) { proxy[,(mannid[m]-1+(1:series[m,k]))]<-pca[,(1:series[m,k])] } } sum(!is.na(proxy[1,])) #[1] 22 temp<-!is.na(proxy[1,]); proxy[,temp]<-extend.persist(proxy[,temp]) proxy[,temp]<-mannomatic(proxy[,temp],sdmethod="undetrended") #this is identical to below # m1<-apply(proxy[503:581,],2,mean,na.rm=TRUE);## means # s1<-apply(proxy[503:581,],2,sd,na.rm=TRUE);##standard deviations # for (j in 1:ncol(proxy)) { # temp<-!is.na(proxy[,j]); # proxy[temp,j]<-(proxy[temp,j]-m1[j])/s1[j] # } TPC<-ts( scale(pc[1:79,]),start=1902,end=1980) #CALCULATE WEIGHTING FACTORS roster<-temp M<-sum(roster) N<-1 G<-array(c(rep(NA,N*ncol(proxy))),dim=c(N,ncol(proxy))) selection<-c(1) U<-TPC[1:79,selection] Ut<-t(U) B<-proxy[503:581,roster] G[selection,roster]<- solve(Ut %*% U ) %*% Ut %*% B; #this is array of regression coefficients in calibration G[selection,roster] # [1] 0.3136919 -0.0947182 0.3176608 -0.2580482 -0.1708455 0.1609345 0.3034515 0.1812157 0.5876245 0.0006883 #[11] 0.0343040 0.0969501 -0.0629002 0.3016627 -0.0605463 0.0764932 0.4639370 0.2271881 -0.1395499 -0.1966504 #[21] -0.2379293 0.0342403 apply(proxy[,temp],2,corfunction) # V12 V13 V14 V15 V18 V20 V43 V47 V53 # 0.3136919 -0.0947182 0.3176608 -0.2580482 -0.1708455 0.1609345 0.3034515 0.1812157 0.5876245 # V62 V63 V64 V66 V67 V68 V72 V84 V85 # 0.0006883 0.0343040 0.0969501 -0.0629002 0.3016627 -0.0605463 0.0764932 0.4639370 0.2271881 # V102 V103 V108 V109 #-0.1395499 -0.1966504 -0.2379293 0.0342403 P<-diag(weights1[roster]^0.5) #this is used in weighted regression C <- P %*% ( t(G[selection,roster] %*% P) %*% solve ((G[selection,roster] %*% P) %*% t((G[selection,roster] %*%P))) ) c(t(C)) #calculate weighting factors # [1] 0.1754800 -0.0529856 0.1777002 -0.1443528 -0.1911428 0.0900272 0.3395030 0.1520587 0.3287185 0.0002541 #[11] 0.0126652 0.0357945 -0.1055596 0.3375017 -0.0677395 0.0855809 0.5190548 0.2541791 -0.0515226 -0.0726044 #[21] -0.1330982 0.0191541 cor(C,G[selection,roster]) # 0.9291 #plot(C,G[selection,roster],xlab="Calculated weights",ylab="TPC1 Correlations", # main="") #cbind(C,G[selection,roster],dimnames(proxy)[[2]][roster]) #legend(locator(1),legend="Gaspe",bty="n", text.width=10) #legend(locator(1),legend="NOAMER PC1",bty="n") ## lines(G[selection,roster],fitted(fm)) fm<-lm(C~G[selection,roster]-1) coef(fm) # 0.7704 test3<-apply(proxy[,roster],2,hockeystat3) weight.1400<-G[1,roster]*test3 order.1400<-order(weight.1400) weight.1400.ordered<-weight.1400[order.1400] ##FIGURE 3 ##CALCULATE BASE SCATTER PLOT #postscript(file.path("c:/climate/scripts/ee.2004","figure3.eps"), width = 7.5, height = 3.5, # horizontal = FALSE, onefile = FALSE, paper = "special", # family = "Helvetica",bg="white",pointsize=8) par(mfrow=c(1,2),mar=c(4,4,2,2)) test3<-apply(proxy[,roster],2,hockeystat3) plot(G[1,roster],test3,main="",xlab="TPC1 Correlation",font.lab=2, ylab="Standardized Difference in Mean",axes=FALSE) abline(h=0) abline(v=0) gaspe<-!is.na(match(1:22,9)) NOAMER<-!is.na(match(1:22,17:18)) axis(side=2, las=2,cex.axis=.75,font=2) axis(side=1, cex.axis=.75,font=2) text(0.32,1.51,"NOAMER PC1",cex=.6) text(0.5,1.4,"GASPE",cex=.6) box() plot(cumsum(weight.1400.ordered),xlab="Proxy Network",ylab="Cumulative Contribution" ,font.lab=2,axes=FALSE) abline(h=0) axis(side=2, las=2,cex.axis=.75,font=2) axis(side=1, cex.axis=.75,font=2) #text(-1.9,1.15,"Centered",font=2) box() # dev.off() ##FIGURE 4 ##STUDY DIFFERENCES #LOAD SHEENJEK INSTEAD OF GASPE load("c:/climate/data/mann/proxy4.tab") sheenjek<-ts(proxy[1:581,55],start=1400,end=1980) sheenjek<-extend.persist(as.matrix(sheenjek))[,1] sheenjek<-normalize(ts(sheenjek,start=1400,end=1980),1902,1980) load("c:/climate/data/mann/proxy.tab") proxy<-ts(proxy[1:581,],start=1400,end=1980) temp<-!is.na(proxy[1,]); proxy[,temp]<-extend.persist(proxy[,temp]) proxy[,temp]<-mannomatic(proxy[,temp],sdmethod="undetrended") #this is identical to below #CALCULATE CENTERED PC SERIES tree<-read.table(file.path(mbhdata,"UVA/NOAMER.txt"),sep="\t",header=TRUE) tree<-ts(tree[,2:ncol(tree)],start=tree[1,1],end=tree[nrow(tree),1]) temp<-!is.na(tree[1,]) tree<-tree[,temp] tree<-tree[1:581,] tree<-extend.persist(tree) dim(tree) #581 70 tree0<-tree z<-princomp(tree) pca<-ts (z$scores[,1:9],start=1400,end=1980) pca<-mannomatic(pca,sdmethod="undetrended") test7<-apply(pca,2,hockeystat3) # Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9 #-0.23654102 -0.25882634 -0.61488051 1.09831635 0.02862578 -0.13700836 -0.89699010 -0.13989351 0.30520733 test8<-apply(pca,2,corfunction) # Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9 # 0.107522502 -0.040820749 -0.153178226 0.485981264 -0.007063538 -0.285181326 0.093584201 0.081530458 0.223798131 #LOAD AND SELECT NOAMER BACKTO1400 SITES par(mfrow=c(1,2)) plot(G[1,roster][gaspe],test3[gaspe],main="",xlab="TPC1 Correlation", ylab="Standardized Difference in Mean", xlim=c(-.25,.6),ylim=c(-.75,1.5)) abline(h=0) abline(v=0) points( corfunction(sheenjek),hockeystat3(sheenjek),col="blue",pch="+") points(G[1,roster][NOAMER],test3[NOAMER]) points(test8[1:2],test7[1:2], col="red",pch="+") #CALC NEW CUMSUM proxy.v1<-proxy proxy.v1[,55]<-sheenjek proxy.v1[1:4,53]<-NA proxy.v1[,84:85]<-pca[,1:2] temp<-!is.na(proxy.v1[1,]) test3<-apply(proxy.v1[,temp],2,hockeystat3) test4<-apply(proxy.v1[,temp],2,corfunction) weight.1400<-test3*test4 order.1400<-order(weight.1400) weight.1400.ordered<-weight.1400[order.1400] #FIGURE 4 #postscript(file.path("c:/climate/scripts/ee.2004","figure4.eps"), width = 7.5, height = 3.5, # horizontal = FALSE, onefile = FALSE, paper = "special", # family = "Helvetica",bg="white",pointsize=8) par(mfrow=c(1,2),mar=c(4,4,2,2)) plot(test4,test3,main="",xlab="TPC1 Correlation", ylab="Standardized Difference in Mean", xlim=c(-.25,.6),ylim=c(-.75,1.5),axes=FALSE,font.lab=2) abline(h=0) abline(v=0) axis(side=2, las=2,cex.axis=.75,font=2) axis(side=1, cex.axis=.75,font=2) box() plot(cumsum(weight.1400.ordered),xlab="Proxy Network",ylab="Cumulative Contribution",ylim=c(-.7,1),font.lab=2,axes=FALSE) abline(h=0) axis(side=2, las=2,cex.axis=.75,font=2) axis(side=1, cex.axis=.75,font=2) #text(-1.9,1.15,"Centered",font=2) box() # dev.off() ##END FIGURES 3 AND 4 ######################3 ### ##FIGURE 5 ########################### #Figure 5. Solid: Weighted average of 6 MBH98 AD1400 step proxies: 4 Quelccaya #series (averaged to one series), Tornetrask temperature reconstruction and Tasmania #temperature reconstruction; Dashed – average of MBH98 NOAMER PC1 and #Gaspé. All series smoothed with 25-year Gaussian filter. #load("c:/climate/data/mann/prname.tab") #prname[1:112,3] #4 Quelc are 12-15; Tasmania 43; Tornetrask 67; Gaspe 53. f<-function(x){f<-filter.combine.pad(x,truncated.gauss.weights(25))[,2];f} load("c:/climate/data/mann/proxy.tab") temp<-c(1:68,100:112) proxy<-ts(proxy[1:581,],start=1400) proxy[,temp]<-extend.persist(proxy[,temp]) proxy[,temp]<-normalize(proxy[,temp],1902,1980) noamer.file<-"c:/climate/data/new/TREE/ITRDB/NOAMER/BACKTO_1400" MBH.PC1<-read.table(file.path(noamer.file,"pc01.out")) MBH.PC1<-ts(MBH.PC1[,2],start=1400,end=1980) proxy[,84]<-normalize(MBH.PC1,1902,1980) test<- proxy[,c(12:15,43,67)] %*% as.matrix ( c ( rep(.25/3,4),1/3,1/3) ) series1<- normalize (ts( test,start=1400),1902,1980) series2<- ts( apply(cbind(proxy[,c(53,84)]),1,mean) ,start=1400) series2<-normalize(series2,1902,1980) ylim0<-c(-3.8,1.5) plot(1400:1980,f(series1),xlab="",ylab="" ,type="l",ylim=ylim0) lines(1400:1980,f(series2),lty=2) #END FIGURE 5 ######################## ### RE of PC1 against Briffa load("c:/climate/data/mann/proxy.tab") Briffa<-proxy[,49] g<-read.table("c:/climate/data/new/TREE/ITRDB/NOAMER/BACKTO_1400/pc01.out") PC1<-ts(g[,2],start=g[1,1],end=g[nrow(g),1]) hockeystat(PC1) #[1] 1.618660 #scale both to 1902-1980 PC1<-normalize(PC1,1902,1980) Briffa<-normalize(Briffa,1902,1980) par(mfrow=c(1,1),mar=c(4,4,2,2)) plot.ts(PC1,ylab="S.D. Units",xlab="") lines(Briffa,col="red") RE(PC1,Briffa,c(1902,1980),c(1600,1901))# -7.716685 Briffa<-filter.combine.pad(Briffa,truncated.gauss.weights(21)) PC1<-filter.combine.pad(PC1,truncated.gauss.weights(21)) ylim0<-range(c(PC1[,2],Briffa[,2]),na.rm=T) par(mar=c(3,3,1,1)) plot(1400:1980,PC1[,2],type="l",lty=2,xlab="",ylab="",ylim=ylim0,lwd=2) lines(1400:1994,Briffa[,2],col="grey50",lwd=2) ###END FIGURE 5 w.emulation<-svd(tree.mannomatic) MBH.eof1<-read.table("c:/climate/data/new/TREE/ITRDB/NOAMER/BACKTO_1400/eof01.out") plot(MBH.eof1,type="l") lines(w.emulation$v[,1],col="red") temp<-(MBH.eof1[,2]-w.emulation$v[,1])>0.005 dimnames(tree)[[2]][temp] # "nm572.txt" "nv515.txt" test<-read.table("c:/climate/data/new/TREE/ITRDB/NOAMER/nm572.txt") test<-ts(test[,2],start=test[1,1],end=test[nrow(test),1]) combine<-ts.union(test,ts(tree[,"nm572.txt"],start=1400,end=1981)) tree.base<-tree z<-princomp(tree) z$sdev^2/sum(z$sdev^2) # Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 # 0.1818731340 0.0896893585 0.0804297123 0.0698970200 0.0546745916 Data1<-apply( tree,2,arima.data) hist(Data1,main="",xlim=c(0,1)) eof1<-read.table("c:/climate/data/new/TREE/ITRDB/NOAMER/BACKTO_1400/eof01.out") par(mar=c(4,4,2,2)) plot(eof1[,2],Data1,xlab="MBH98 EOF1 Weighting",ylab="AR1 Coefficient") temp<-!is.na(match(dimnames(tree)[[2]],graybill.list)) #20 sites points(eof1[temp,2],Data1[temp],pch="+") acf(tree[ !is.na(tree[,19]),19],lag.max=350) lines(0:350,Data1[19]^(0:350),col="red",lwd=2) MM.PC1<-ts(-z$scores[,1],start=1400,end=1980) #SIGN-ADJUSTED MM.PC1<-norm(MM.PC1,1902,1980) MM.PC1<-as.numeric(MM.PC1) par(mfrow=c(3,1)) m0<-mean(NOAMER.PC1[503:581]) sd0<-sd(NOAMER.PC1[503:581]) MBH<-ts((NOAMER.PC1-m0)/sd0,start=1400,end=1980) ts.plot(MBH,ylab="MBH") abline(h=0) ts.plot(MM.PC1,ylab="MM") abline(h=0) PC4<-ts(z$scores[,4],start=1400,end=1980) #SIGN-ADJUSTED PC4<-norm(PC4,1902,1980) ts.plot(PC4,ylab="PC4") abline(h=0) ############## ##FIGURE 1 load("c:/climate/data/MBH99/hockeysticks.tab") par(mfrow=c(2,1),mar=c(3,4,1,2)) fred<-norm(hockeysticks[,74],1902,1980) fred<-hockeysticks[,74] plot.ts(fred,ylab="(a)", main="",axes=FALSE) #abline(h=0) plot.ts(NH.MBH,ylab="(b)",main="",axes=FALSE) #abline(h=0) axis(side=1, at=seq(1400,2000,by=100), labels=c(1400,"",1600,"",1800,"",2000)) loc.out<-"c:/climate/data/ee/example.png" png(file=loc.out,width=360,height=240,pointsize=5,bg="white"); par(mfrow=c(2,1),mar=c(2,4,2,2)) plot.ts(fred,ylab="(a)", main="",axes=FALSE) #abline(h=0) plot.ts(NH.MBH,ylab="(b)",main="",axes=FALSE) #abline(h=0) axis(side=1, at=seq(1400,2000,by=100), labels=c(1400,"",1600,"",1800,"",2000)) dev.off(); ##FIGURE acf(ca534S) ########################## ############################ #SKIP THIS STEP EXCEP IN SETUP #REDO TREE RING PC CALCULATIONS USING PRINCOMP AND STORE COVARIANCE PC CALCULATIONS # used in reconstruction variation ############################ #A. LOAD PARAMETERS url<-"c:/climate/data/mann/WDCP" dataversion<-c(rep("WDCP",5),"UVA") #or "MBH98" or "UVA" method<-"persist" #or "remaining" #B. DO CALCULATIONS FOR 6 REGIONS for (m in 1:6) { loc<-file.path (mbhdata,dataversion[m], paste( regionid[m],"txt",sep=".") ) tree<-read.table(loc,header=TRUE) tree<-ts(tree,start=tree[1,1],end=tree[nrow(tree),1]) tree<-tree[,2:ncol(tree)] tree0<-tree for ( k in 1:length(mbh98.backto[[m]]) ) { pca<- array ( rep(NA,mannpc[m]*581), dim=c(581,mannpc[m] ) ) temp<-!is.na(tree0[ (mbh98.backto[[m]][k]-tsp(tree0)[1]+1) ,] ) tree<-tree0[,temp] #dim(tree)#[1] 593 70 FOR NOAMER BACKTO_1400 dates<-tsp(tree)# [1] 1400 1992 1 #EXTEND BY PERSISTENCE AND TRUNCATE TO 1980 AND end0<-min(1980,dates[2]) if (method=="persist") tree<-extend.persist(tree) if (method=="remaining") tree<-extend.remaining(tree) tree<-tree[1:(end0-dates[1]+1),] z<-princomp( tree[ ( mbh98.backto[[m]][k]-tsp(tree0)[1]+1):(end0-tsp(tree0)[1]+1), ] ) # pca[(mbh98.backto[[m]][k]-1400+1):(end0-1400+1),(1:min(mannpc[m],ncol(z$scores)))]<-z$scores[,1:min(mannpc[m],ncol(z$scores))] #EXTEND TO 1980 BY PERSISTENCE IF NECESSARY if (end0<1980) pca [(end0+1-1400+1):581,] <- t( array( rep(pca[(end0-1400+1),],(1980-end0)) , dim=c(ncol(pca),(1980-end0) )) ) #save(pca,file=file.path(pcstore,paste(dataversion[1],method,sep="."),paste(regionid[m],mbh98.backto[[m]][k],"tab",sep=".") ) ) } #end of k iteration } # end of m-iteration # end saving of tree ring PC series ###########################################