##################### ###THIS IS SCRIPT TO SHOW DIFFERENCE BETWEEN PC CALCULATION METHODS #the network used here is simulated to look like NOAMER tree ring network and is in file make.proxy.txt #this is a sample using methods from GRL simulations ###FUNCTIONS USED ##mannomatic #this function re-standardizes a data matrix B on calibration period mean and standard deviation #options are for detrended and undetrended; default is detrended mannomatic<-function(B,M=(79-1),sdmethod="detrended") { N<-nrow(B) m2<-apply(B[(N-M):N,],2,mean) s2<-apply(B[(N-M):N,],2,sd) if (sdmethod=="undetrended") D<-scale(B,center= m2, scale= apply(B[(N-M):N,],2,sd) ) if (sdmethod=="detrended") { D<-scale(B,center= m2, scale= apply(B[(N-M):N,],2,sd) ); sdprox<-apply(D[(N-M):N,],2,sd.detrend) D<-scale(D,center= FALSE,scale= sdprox)} mannomatic<-D mannomatic } ##sd.detrend #this function calculates the ***detrended*** standard deviation of a series #the impact of using *detrended* standard deviations is to icnrease the weight of series with trends #this is secondary to the impact of the decentered mean, but has an additional effect sd.detrend<-function(x) { t<-c(1:length(x)) fm<-lm(x~t) sd.detrend<-sd(fm$residuals) sd.detrend } ## hockeystat #this calculates the "hockey stick index" i.e. (calibration period mean - entire period mean)/std deviation hockeystat<-function(test,N=581,M=79) { hockeystat<- ( mean(test[(N-M+1):N],na.rm=TRUE)- mean(test[1:N],na.rm=TRUE) )/sd(test[1:N],na.rm=TRUE); hockeystat } ssq<-function(x) {ssq<-sum(x^2);ssq} ####################3 ###LOAD PSEUDOPROXY NETWORK url<-"http://www.climateaudit.org/data/zorita" url<-"c:/climate/data/zorita" network.sim<-read.table(file.path(url,"network.sim.txt"),sep="\t",skip=1) dim(network.sim) #581 70 #this has properties of AD1400 NOAMER tree ring nework n<-ncol(network.sim) index<-503:581 tree.mannomatic<-mannomatic(network.sim,sdmethod="detrended") #this does Mannomatic transformation of data max (abs(apply(tree.mannomatic[index,],2,mean))) #2.872614e-16 max (abs(apply(tree.mannomatic[index,],2,sd)-1)) # 0.472871 #this is network after Mann transformation with detrended standard deviation tree.shortnormalize<-mannomatic(network.sim,sdmethod="undetrended") max (abs(apply(tree.shortnormalize[index,],2,mean))) # 2.831991e-16 max (abs(apply(tree.shortnormalize[index,],2,sd)-1)) # 1.110223e-16 #this is network re-standardized on calibration period #the main impact is mean decentering but I'll show the extra effect of detrended standard deviation tree.correlation<-scale(network.sim) #this sets up PC using correlation PCs tree.covariance<-scale(network.sim,scale=FALSE) ### thus we have 4 different scalings and centerings above ###PC ANALYSES ###MBH METHOD svd.mannomatic<-svd(tree.mannomatic) PC1.mannomatic<-svd.mannomatic$u[,1] hockeystat(PC1.mannomatic) #1.687248 #AR1 network -1.437305 pca.mannomatic<-prcomp(tree.mannomatic,center=FALSE) hockeystat(pca.mannomatic$x[,1]) #1.687248 #AR1 network -1.437305 ##ZORITA #there are two methods used in programs standard.f and reconstandard.f #zorita eof calculations use covariance matrix following prior standardizations #thus default covariance in prcomp used ##ZORITA "STANDARD" METHOD #this follows program standard.f in which each series was standardized over full network pca.correlation<-prcomp(tree.correlation) PC1.correlation<-pca.correlation$x[,1] hockeystat(PC1.correlation) # -0.05691923 #AR1 network 0.4888639 ##ZORITA "RECONSTANDARD" METHOD pca.reconstandard<-prcomp(tree.shortnormalize) PC1.reconstandard<-pca.reconstandard$x[,1] hockeystat(PC1.reconstandard) # 0.4265588 # AR1 network 0.1581417 ##PLOT PCs nf<-layout(array(1:2,dim=c(2,1)),heights=c(1.1,1.4)) par(mar=c(0,3,1,1)) x<-PC1.mannomatic plot(1400:1980,sign(hockeystat(x)) * x,type="l",axes=FALSE,xlab="",ylab="") axis(side=1,labels=FALSE);axis(side=2,las=1);box() text(1400,.9*max(x*sign(hockeystat(x))),"MBH",pos=4,font=2) par(mar=c(3,3,0,1)) x<-PC1.correlation;x<-x/ssq(x) plot(1400:1980,x,type="l",axes=FALSE,xlab="",ylab="") axis(side=1);axis(side=2,las=1);box() x<-PC1.reconstandard;x<-x/ssq(x) lines(1400:1980,x,col="red") legend( 1600,.0035,legend=c("Correlation","Zorita"),fill=c("black","red")) text(1400,.9*max(x),"Zorita",pos=4,font=2) combine<-cbind(PC1.mannomatic,PC1.correlation,PC1.reconstandard) #write.table(combine,file="c:/climate/data/zorita/pcs.AR1.txt",sep="\t",quote=FALSE,row.names=FALSE) write.table(combine,file="c:/climate/data/zorita/pcs.sim.txt",sep="\t",quote=FALSE,row.names=FALSE) ##ARIMA.SIM NE x<-arima.sim(n=581*70+100,model=list(ar=0.7)) network.sim<-array(x[101:(581*70+100)],dim=c(581,70)) write.table(network.sim,file="c:/climate/data/zorita/network.sim.AR1.txt",sep="\t",quote=FALSE,row.names=FALSE)