### SANTER FUNCTIONS ### ### df.d1star=function(sd.model,M,se.obs,neff) ### Santer equation 13 made into a function to calculate df for t test ### ### santerf = function(x,end0=2050,start0=1979,M=19, ### trend.target=trend.model.T2LT, sd.target=sd.model.T2LT ) ### calculates trend significance using Santer AR1 method ### ### do.santer = function( X, target, method="lapse") ### wrapper function to do Santer-style analysis against target ### returns 3 endpoints 1999, 2007, 2009 ### $`2009` ### id target trend se sd r1 neff adj_df d1star d1zero pt ptzero # 1 had 0.146 0.1172 0.0560 0.2131 0.9186 15.7413 13.9781 0.4947 2.0931 0.6858 0.9725 # 2 giss 0.146 0.1066 0.0560 0.2099 0.9174 15.9879 14.2152 0.6780 1.9052 0.7457 0.9614 # 3 noaa 0.146 0.1121 0.0519 0.2011 0.9137 16.5533 14.7869 0.6264 2.1604 0.7297 0.9762 ### table1 ### # made from utility function make.table1 to make Santer 2008 Table 1 ### ### table4 ### #from Santer 2008 Table 4 ### ### santer.model.info=model.info ### collation of trend information from Tables 1 and 2 with columns: ### tropo firma trend sd ### ### rss ### collation of RSS series tropical (20-20 series) ### ### ### ## 1. Santer 2008 Table 1 make.table1=function(loc="http://www.climateaudit.info/data/models/santer/santer_2008_table1.dat") { santer=read.table(loc,skip=1) names(santer)=c("item","layer","trend","se","sd","r1","neff") row.names(santer)=paste(santer[,1],santer[,2],sep="_") santer=santer[,3:ncol(santer)] #insert Table III info santer$d1star=NA;santer$dstar=NA santer["UAH_T2LT",c("d1star","dstar")]= c(1.11,7.16) #from Table III for UAH santer["RSS_T2LT",c("d1star","dstar")]= c(0.37,2.25) #from Table III for UAH santer["UAH_T2",c("d1star","dstar")]= c(1.19,6.78) #from Table III for UAH santer["RSS_T2",c("d1star","dstar")]= c(0.44,2.48) #from Table III for UAH make.table1= santer make.table1 } table1=santer=make.table1() # "d:/climate/data/models/santer/santer_2008_table1.dat") # trend se sd r1 neff d1star dstar #HadCRUT3v_TL+O 0.119 0.117 0.197 0.934 8.6 NA NA #Multi-model_mean_TL+O 0.146 0.214 0.274 0.915 11.7 NA NA #... ##2. Santer 2008 Table 4 make.table4=function(loc="http://www.climateaudit.info/data/models/santer/santer_2008_table4.csv") { # "d:/climate/data/models/santer/santer_2008_table4.csv") { make.table4=read.csv("d:/climate/data/models/santer/santer_2008_table4.csv",sep="\t",header=TRUE) #this also has Table VI info make.table4 } table4=make.table4() # Sat Surf Trend SE R1 Neff Dstar D1star #1 HadCRUT3v UAH_T2LT 0.061 0.036 0.642 55.0 17.05 3.50 #2 HadCRUT3v RSS_T2LT -0.046 0.034 0.608 61.5 3.05 0.67 ##3. Santer Model Info tropo=c("T0","T2LT","T2"); firma=c("All","Ocean","Land") make.model.info=function(p=length(tropo),q=length(firma) ) { model.info=data.frame(tropo=gl(p+1,q,labels=c(tropo,"Lapse")),firma=rep(firma,p+1),trend=rep(NA,(p+1)*q),sd=rep(NA,(p+1)*q)) model.info[1,3:4]=table1$trend[2:3] #T0 All from S08 T1 model.info[2,3:4]=table1$trend[7:8] #T0 Ocean from S08 T1 model.info[4,3:4]=table1$trend[11:12] #T2LT All from S08 T1 model.info[7,3:4]=table1$trend[15:16] #T2 All from S08 T1 model.info[10,3:4]=table4$Trend[3:4] #T2LT All Lapse from S08 T4 model.info[11,3:4]=table4$Trend[11:12] #T2LT Ocean Lapse from S08 T4 model.info$trend[5]=model.info$trend[2]-model.info$trend[11] #T2LT Ocean= T0 Ocean minus Lapse Ocean-T2LT k= 0.7822222 #proportion TRP ocean from hadISST mask model.info$trend[3]=(model.info$trend[1]- k*model.info$trend[2])/(1-k) #T0 Land: by difference from T0 All and T0 Ocean model.info$trend[6]=(model.info$trend[4]- k*model.info$trend[5])/(1-k) #T2LT by difference from T2LT All and T2LT Ocean: no difference model.info$trend[12]=(model.info$trend[3]- model.info$trend[6]) #T2LT Land Lapse: T0 Land minus T2LT Land model.info$sd[5:6]=model.info$sd[4] ###interim estimate since not reported by Santer model.info$sd[3]=.07 #interim estimate higher than T0 Ocean model.info$sd[12]=.04 #estimate interim higher than Lapse Ocean return(model.info) } santer.model.info=model.info=make.model.info() # tropo firma trend sd #1 T0 All 0.14600000 0.066 #2 T0 Ocean 0.13000000 0.062 #3 T0 Land 0.20346938 0.070 #4 T2LT All 0.21500000 0.092 #5 T2LT Ocean 0.21500000 0.092 #6 T2LT Land 0.21500000 0.092 #7 T2 All 0.19900000 0.098 #8 T2 Ocean NA NA #9 T2 Land NA NA #10 Lapse All -0.06900000 0.032 #11 Lapse Ocean -0.08500000 0.038 #12 Lapse Land -0.01153062 0.040 ################### ## trend calculation for 1979-2009 period by run ############################### collate.model.information=function(Runs, Info= Info.chad,end0=2009.99 ) { tropo=names(Runs)=c("T0","T2LT","T2") TrendRuns= rep(list(NA),3); names(TrendRuns)=tropo for ( i in 1:3) TrendRuns[[i]]= t (apply(window(Runs[[i]],1979,end0),2, function(x) trend(x,method="santer") )) for(i in 1:3) { TrendRuns[[i]]=data.frame(TrendRuns[[i]]); names(TrendRuns[[i]])=c("trend","neff","se")} #trends for lapse rates TrendRuns[["T0_T2LT"]]=NA; TrendRuns[["T0_T2"]]=NA for ( i in 4:5) { TrendRuns[[i]]= t (apply(window(Runs[[1]]-Runs[[i-2]] ,1979,end0),2, function(x) trend(x,method="santer") )) TrendRuns[[i]]=data.frame(TrendRuns[[i]]); names(TrendRuns[[i]])=c("trend","neff","se") row.names(TrendRuns[[i]])=row.names(TrendRuns[[1]]) } #into Info Info$T0=unlist(TrendRuns[["T0"]]$trend) Info$T2LT=unlist(TrendRuns[["T2LT"]]$trend) Info$T2=unlist(TrendRuns[["T2"]]$trend) Info$T2LT=unlist(TrendRuns[["T2LT"]]$trend) Info$T2=unlist(TrendRuns[["T2"]]$trend) Info$lapse_T2LT=unlist(TrendRuns[["T0_T2LT"]]$trend) Info$lapse_T2=unlist(TrendRuns[["T0_T2"]]$trend) Info=groupedData(T2LT~1|model/run,data=Info) fm.lmlis=rep(list(NA),5);names(fm.lmlis)=names(TrendRuns) for (i in 1:5) { X=Info[,c(1:3,i+3)];names(X)[4]="trend" fm.lmlis[[i]]=lmList(trend~1|model,data=X)} count= tapply(!is.na(Info$model),Info$model,sum) Models=data.frame( array(NA,dim=c(24,7))); names(Models) =c("id","count","T0","T2LT","T2","T0_T2LT","T0_T2") Models$id=unique(Info$model); Models$count=count for (i in 1:5) Models[,i+2]= round(fixef(fm.lmlis[[i]])+ranef(fm.lmlis[[i]]),4) # EXTRACT TARGET MODEL INFORMATION #these are target info that update information in Santer tables of model summaries target.info= cbind( trend =round(apply(Models[,3:7],2,mean ),3), sd= round( apply(Models[,3:7],2,sd ) ,3), count=rep( nrow(ranef(fm.lmlis[["T0"]])),5 ) ) row.names(target.info)=names(TrendRuns); collate.model.information=list(Runs=Runs, Models=Models,Info=Info, TrendRuns=TrendRuns, target.info=target.info,fm=fm.lmlis) return(collate.model.information) } #end function make.trend.information=function(Data) { Trend=list() f=santerf santerf(x=Sat[["T2LT"]][["All"]][,"rss"], target= Data[["loti"]]$target.info["T2LT",]) for(i in 1:3) Surf[[i]]=window(Surf[[i]],start=1979) (Trend[["T2LT"]][["All"]]= cbind( t(apply( Sat[["T2LT"]][["All"]],2, function(x) santerf(x, target= Data[["loti"]]$target.info["T2LT",]) ) ), santer= table1[10:9,"d1star"])) (Trend[["T2"]][["All"]]= cbind( t(apply( Sat[["T2"]][["All"]],2, function(x) santerf(x, target= Data[["loti"]]$target.info["T2",]) ) ), santer= table1[14:13,"d1star"]) ) (Trend[["T2LT"]][["Ocean"]]= cbind( t(apply( Sat[["T2LT"]][["Ocean"]],2, function(x) santerf(x, target= Data[["ocean"]]$target.info["T2LT",]) ) ), santer= table1[c("RSS_T2LT","UAH_T2LT"),"d1star"]) ) (Trend[["T2"]][["Ocean"]]= cbind( t(apply( Sat[["T2"]][["Ocean"]],2, function(x) santerf(x, target= Data[["ocean"]]$target.info["T2",]) ) ), santer= table1[c("RSS_T2","UAH_T2"),"d1star"]) ) (Trend[["T2LT"]][["Land"]]= cbind( t(apply( Sat[["T2LT"]][["Land"]],2, function(x) santerf(x, target= Data[["land"]]$target.info["T2LT",]) ) ), santer= rep(NA,2) ) ) (Trend[["T2"]][["Land"]]= cbind( t(apply( Sat[["T2"]][["Land"]],2, function(x) santerf(x, target= Data[["land"]]$target.info["T2",]) ) ), santer= rep(NA,2) ) ) (Trend[["T0"]][["All"]]= cbind( t(apply( Surf$All,2, function(x) santerf(x, target= Data[["loti"]]$target.info["T0",]) ) ), santer= rep( table1[1,"d1star"],3) )) (Trend[["T0"]][["Ocean"]]= cbind( t(apply( Surf$Ocean,2, function(x) santerf(x, target= Data[["ocean"]]$target.info["T0",]) ) ), santer= table1[4,"d1star"])) (Trend[["T0"]][["Land"]]= cbind( t(apply( Surf$Land,2, function(x) santerf(x, target= Data[["land"]]$target.info["T0",]) ) ), santer= NA)) ### LAPSE RATE TRENDS: Results second section #T2LT All X= window(Surf$All[,"had"],start=1979) - window(Sat[["T2LT"]]$All,end=tsp(Surf$All)[2]) dimnames(X)[[2]]=c("rss","uah") ( Trend[["Lapse_T2LT"]][["All"]]= cbind( t(apply(X,2, function(x) santerf(x, target= Data[["loti"]]$target.info["T0_T2LT",]) ) ), santer= table4[2:1,"D1star"]) ) #Hadcrut #T2 All X= window(Surf$All[,"had"],start=1979) - window(Sat[["T2"]]$All,end=tsp(Surf$All)[2]) dimnames(X)[[2]]=c("rss","uah") ( Trend[["Lapse_T2"]][["All"]]= cbind( t(apply(X,2, function(x) santerf(x, target= Data[["loti"]]$target.info["T0_T2",]) ) ), santer= rep(NA,2) ) ) #T2LT Ocean HadCRU X= window(Surf$Ocean[,"had"],start=1979) - window(Sat[["T2LT"]]$Ocean,end=tsp(Surf$Ocean)[2]) dimnames(X)[[2]]=c("rss","uah") ( Trend[["Lapse_T2LT"]][["Ocean"]]= cbind( t(apply(X,2, function(x) santerf(x, target= Data[["ocean"]]$target.info["T0_T2LT",]) ) ), santer= table4[c(8,5),"D1star"] ) ) #T2LT Ocean HadCRU NOAA X= window(Surf$Ocean[,"noaa"],start=1979) - window(Sat[["T2LT"]]$Ocean,end=tsp(Surf$Ocean)[2]) dimnames(X)[[2]]=c("rss","uah") cbind( t(apply(X,2, function(x) santerf(x, target= Data[["ocean"]]$target.info["T0_T2LT",]) ) ), santer= table4[c(10,7),"D1star"]) #NOAA #T2LT Land HadCRU X= window(Surf$Land[,"had"],start=1979) - window(Sat[["T2LT"]]$Land,end=tsp(Surf$Land)[2]) dimnames(X)[[2]]=c("rss","uah") ( Trend[["Lapse_T2LT"]][["Land"]]= cbind( t(apply(X,2, function(x) santerf(x, target= Data[["land"]]$target.info["T0_T2LT",]) ) ), santer= NA ) ) return(Trend) } # function santerf= function(x,target,period="monthly" ) { #monthly trend.model=target["trend"]; sd.model=target["sd"];M=target["count"] N=length(x); sd.obs=sd(x,na.rm=T) fm= lm (x~ I( (1:length(x))/120 )) trend.obs= fm$coef[2] r= arima (fm$residuals,order=c(1,0,0))$coef[1]; # neff= N * (1-r)/(1+r) ; se.obs= sqrt((N-2)/(neff-2))* summary(fm)$coef[2,"Std. Error"] d1star = (trend.model-trend.obs)/ sqrt( (sd.model/sqrt(M-1))^2 + se.obs^2) adj_df= ( (sd.model/sqrt(M-1))^2 + se.obs^2) / ( (sd.model/sqrt(M-1))^2 /(M-1) + se.obs^2/(neff-2) ) # d1zero =trend.obs/se.obs pt.obs= pt(d1star,adj_df) # ptzero=pt(d1zero,adj_df) santerf= round(c(trend.obs,trend.model,d1star,pt=pt.obs),3) names(santerf)=c("obs","target","d1star","pt") santerf } ## df.d1star df.d1star=function(sd.model,M,se.obs,neff) { #df calculation from Santer equation 13 made into function C= (sd.model/sqrt(M-1))^2 + se.obs^2;C # 0.0195 D=(sd.model/sqrt(M-1))^2 /(M-1) + se.obs^2/(neff-2) ;D # 0.00155 df.d1star= C/D # 12.6 df.d1star } df_effective= function(se.obs,se.target,df.obs,df.target) { C = se.target^2 + se.obs^2;C # 0.0195 D= se.target^2 /df.target + se.obs^2/df.obs ;D # 0.00155 df_effective = C/D # 12.6 df_effective } ### # function to compare trend against a target with standard error and do t-test allowing for AR1 autocorrelation # x is a time series # (se.model.T2LT= sd.model.T2LT/sqrt (nrow(info.model)-1) ) # 0.02168461 adjusted_ttest =function(fm.obs , fm.target ) { f=rep(NA,4); names(f)=c("trend","target", "t_adj","pt_adj") year=c(time(x))/10; N=length(x); year=year-mean(year) fm= lm (x~year) trend.obs=f["trend"]= fm$coef[2] r= arima (fm$residuals,order=c(1,0,0))$coef[1]; # neff= N * (1-r)/(1+r) ; se.obs= sqrt((N-2)/(neff-2))* summary(fm)$coef[2,"Std. Error"] se.target= summary(fm.target)$coef [1,"Std. Error"] f["target"]= fm.target$coef f["t_adj"]= (target-trend.obs)/ sqrt (se.target ^2 + se.obs^2) f["df_adj"]= df_effective(se.obs, se.target,df.obs= neff-2,df.target ) f["pt_adj"]=pt(f["t_adj"],f["df_adj"]) return(f) } # fm.obs= lm(Y[,"rich"]~X[,"had"]) anom= function(x, start0= 1979, end0=1999.99) { test= ts (t( array(x, dim=c (12, length(x)/12 ))), start=tsp(x)[1]) m0= apply( window (test,start=start0,end=end0),2,mean,na.rm=T) anom= round( scale( test,center=m0,scale=FALSE) ,3) anom= ts(c(t(anom)),start=tsp(x)[1],freq=12) return(anom) } trend =function(x,method="short" ) { if (method=="short") return(round(lm(x~ I((1:length(x))/120))$coef[2],3)) if (method=="santer") { fm= lm(x~ I((1:length(x))/120)) se.ols=summary(fm)$coef[2,2] N=length(x) (trend.obs= fm$coef[2])# 1.116659 fm.arima=arima(fm$residuals,order=c(1,0,0));fm.arima r=fm.arima$coef[1] neff= N * (1-r)/(1+r) ;neff # 102.2358 se.obs= sqrt((N-2)/(neff-2))* se.ols; se.obs # 0.071758 trend= round(c(trend=unlist(trend.obs), neff=neff, se=se.obs),4) trend }} do.santer=function( X, target, method="lapse", endpoints = c(1999,2007,2009) ) { if ( method=="level") { K=length(endpoints) Obs=rep(list(NA),K);names(Obs)=paste(endpoints); q=ncol(X) for (k in 1:K) Obs[[k]]=data.frame( id=dimnames(X)[[2]],target=rep( target["trend"],q),trend=NA,se=NA,sd=NA,r1=NA,neff=NA,adj_df=NA,d1star=NA,d1zero=NA,pt=NA,ptzero=NA) for (j in 1:q) { temp1=!is.na(X[,paste(Obs[[1]]$id[j])] ); end0=min( max(time(X)[temp1]) ) x=window( X[,paste(Obs[[1]]$id[j])],start=1979,end=end0 ) Obs[[1]][j,3:12]=round(f(x,end0=1999,trend.model=target["trend"] ,sd.model=target["sd"]),4) Obs[[2]][j,3:12]=round(f(x,end0=2007,trend.model=target["trend"] ,sd.model=target["sd"]),4) Obs[[3]][j,3:12]=round(f(x,end0=2009,trend.model=target["trend"] ,sd.model=target["sd"]),4) } for(k in 1:3) row.names(Obs[[k]])=paste("#",row.names(Obs[[k]])) return(Obs) } else { Surf=X$surf; Sat=X$sat; Lapse=rep(list(NA),3);names(Lapse)=c(1999,2007,2009) p=ncol(Surf); q=ncol(Sat) for (k in 1:3) { Lapse[[k]]=data.frame( surf=gl(p,q,labels=dimnames(Surf)[[2]]), sat=rep(dimnames(Sat)[[2]],p),target=rep( target["trend"],p*q), trend=NA,se=NA,sd=NA,r1=NA,neff=NA,adj_df=NA,d1star=NA,d1zero=NA,pt=NA,ptzero=NA) } #santerf= function(x,end0=2050,start0=1979,M=19, #trend.model=trend.model.T2LT, sd.model=sd.model.T2LT ) { for (j in 1:(p*q)) { temp1=!is.na(Surf[,paste(Lapse[[1]]$surf[j])] ); temp2=!is.na(Sat[,paste(Lapse[[1]]$sat[j])] ) ; end0=min( max(time(Surf)[temp1]), max(time(Sat)[temp2]) ) x=window( Surf[,paste(Lapse[[1]]$surf[j])],start=1979,end=end0 ) y=window( Sat[,paste(Lapse[[1]]$sat[j])],start=1979,end=end0 ) lapse=x-y Lapse[[1]][j,4:13]=round(f(lapse,end0=1999,trend.model=target[1],sd.model=target[2]),4) Lapse[[2]][j,4:13]=round(f(lapse,end0=2007,trend.model=target[1],sd.model=target[2]),4) Lapse[[3]][j,4:13]=round(f(lapse,end0=2009,trend.model=target[1],sd.model=target[2]),4) } for(k in 1:3) row.names(Lapse[[k]])=paste("#",row.names(Lapse[[k]])) return(Lapse) } }