##BOXPLOT OF TRENDS
    #previously collated Model info (57 A1B runs)

	library(nlme)
	Info=read.table("http://www.climateaudit.info/data/models/santer/2010/info.runs.csv",sep="\t",header=TRUE)   #57 9
	Info=groupedData(T2LT~1|model/run,data=Info)

  ##random effects model
   	fm.lme= lme(T2LT~1,Info)

  #boxplot
 	par(mar=c(10,3,2,1))
	boxplot(T2LT~model, Info[Info$count>=4,],las=3,col="violet",ylim=c(0.1,0.55),cex.axis=.8)
	title("T2LT Trends by Model")
 
 #summary
     summary(fm.lme)

	#Fixed effects: T2LT ~ 1 
	#            Value Std.Error DF t-value p-value
	#(Intercept) 0.341    0.0138 33    24.7       0
	   #this is the ensemble trend

	# Formula: ~1 | run %in% model
	#        (Intercept) Residual
	#StdDev:     0.00443  0.00459
	   #this is the intra-run std deviation

	#Random effects:
	# Formula: ~1 | model
	#        (Intercept)
	#StdDev:      0.0672
           #this is the between-model std deviation

###############
##COMPARE TO OBSERVATIONS
########################################

	source("http://www.climateaudit.info/scripts/utilities.txt")
	download.file("http://www.climateaudit.info/data/models/santer/2010/Sat_201007.tab","temp",mode="wb")
		load("temp") ; names(Sat) #  "T2LT" "T2" 
	(trend.uah=trend(Sat$T2LT$All[,"uah"],method="santer"))
	(trend.rss=trend(Sat$T2LT$All[,"rss"],method="santer"))
	Trend=rbind(trend.uah,trend.rss)

 #GDD(file="d:/climate/images/2010/santer/boxplot_T2LT.png",type="png",h=560,w=480) 
    temp=Info$count>=4
	par(mar=c(10,3,2,1))
	boxplot(T2LT~model, Info[temp,],las=3,col="violet",ylim=c(0,.5),cex.axis=.8,xlim=c(0,9.5))
	title("T2LT Trends by Model")
	points(8,trend.uah["trend"],pch=19,col=2)
	arrows(8, trend.uah["trend"]-2* trend.uah["se"],8,trend.uah["trend"]+2* trend.uah["se"],
	  col=2, 	length=.12,lwd=2,code=3) 

	points(9,trend.rss["trend"],pch=19,col=2)
	arrows(9, trend.rss["trend"]-2* trend.rss["se"],9,trend.rss["trend"]+2* trend.rss["se"], 
	     col=2, 	length=.12,lwd=2,code=3) 

	arrows(7,  fixef(fm.lme)-2* fm.lme$sigma,7, fixef(fm.lme)+2* fm.lme$sigma, 
	     col=2, 	length=.12,lwd=2,angle=90,code=3) 
	axis(side=1,at=c(7,8,9),labels=c("Ensemble","UAH","RSS"),las=3,col=2,cex.axis=.9)
  #dev.off()

#############
## t-tests
##########################

	tropo="T2LT"; K=nrow(Model)
	Model= data.frame(model=unique(Info$model))
	Model$count=tapply(!is.na(Info$model),Info$model,sum)
	Model$tropo_mean= tapply(Info[,tropo],Info$model,mean)
	Model$tropo_sd= tapply(Info[,tropo],Info$model,sd)
	temp=Model$count>=4
	Model$t_uah= NA
	Model$t_rss= NA
	for (j in (1:K)[temp]) {
	   for(i in 1:2) {
	   Model[j,4+i]= (Model$tropo_mean[j]- Trend[i,"trend"])/ sqrt( Model$tropo_sd[j]^2+Trend[i,"se"]^2)  
	}}
	Model[temp,c(1,2,5,6)]
	  #	                model count tropo_mean tropo_sd t_uah t_rss
	  #2  cccma_cgcm_3.1_T47     5      0.347  0.00696  4.01 2.885
	  #7      echam_5_mpi_om     4      0.446  0.01312  5.43 4.376
	  #13            giss_er     5      0.308  0.00183  3.44 2.290
	  #20    mri_cgcm_2.3.2a     5      0.296  0.00756  3.23 2.083
	  #21      ncar_ccsm_3.0     7      0.299  0.00274  3.30 2.150
	  #22         ncar_pcm_1     4      0.209  0.00211  1.92 0.727



   ##Ensemble t-test

	ensemble.stat=rep(NA,2)	 
 	for (i in 1:2) ensemble.stat[i]= (fixef(fm.lme)- Trend[i,"trend"])/ 
		sqrt( fm.lme$sigma^2+Trend[i,"se"]^2) 
	ensemble.stat
	  # 3.93 2.80