##This is a file of support functions #These functions are used to read jags files and write jags files #The write jags files is useful for running native command line jags with the data generated in R. read.jagsdata <- function(file){ e <- new.env(); eval(parse(file), e); return(as.list(e)) } write.jagsdata <- function(x,file){jagsdata<-lapply(x,as.double); attach(jagsdata) ; try(dump(names(jagsdata),file=file)) ; detach(jagsdata)} #This funciton is used to generate column names from lists i.e. #list(a=c("tom","bill"),b="george") #gives a.tom,a.bill,b.george in this order unlist.names<-function(x) { #x is a list with names paste(rep(names(x),sapply(x,length)),unlist(x),sep=".") } #These plots are generated using tplot jplots<-function(plotfun,data,title,mfrow=c(4,3),...) { x11() par(mfrow=mfrow) par(oma=c(0,0,1,0)) plotfun(data,auto.layout=F,ask=F,...) title(outer=T,main=title) } #This allows us to plot our coefficients nicely it works with plots `tplot` <- function (x, trace = TRUE, density = TRUE, smooth = TRUE, bwf, auto.layout = TRUE, ask = par("ask"), mfrow=NULL,...) { oldpar <- NULL on.exit(par(oldpar)) if (auto.layout & is.null(mfrow)) { mfrow <- set.mfrow(Nchains = nchain(x), Nparms = nvar(x), nplots = trace + density) oldpar <- par(mfrow = mfrow) } else oldpar<-par(mfrow=mfrow) for (i in 1:nvar(x)) { if (trace) traceplot(x[, i, drop = FALSE], smooth = smooth, ...) if (density) { if (missing(bwf)) densplot(x[, i, drop = FALSE], ...) else densplot(x[, i, drop = FALSE], bwf = bwf, ...) } if (i == 1) oldpar <- c(oldpar, par(ask = ask)) } } environment(tplot)<-environment(coda::traceplot) #This function converts an mcmc list correctly into an mcarray (It may not be needed). `as.mcmc.list.mcarray` <- function (x, ...) { if (is.null(dim(x)) || is.null(names(dim(x)))) { NextMethod() } xdim <- dim(x) ndim <- length(xdim) dn <- names(xdim) which.iter <- which(dn == "iteration") if (length(which.iter) != 1) { stop("Bad iteration dimension in mcarray") } which.chain <- which(dn == "chain") if (length(which.chain) > 1) { stop("Bad chain dimension in mcarray") } niter <- xdim[which.iter] if (length(which.chain) == 0) { perm <- c(which.iter,(1:ndim)[-which.iter]) x <- matrix(aperm(x, perm), nrow = niter) ans <- mcmc.list(mcmc(x)) } else { nchain <- xdim[which.chain] ans <- vector("list", nchain) len <- prod(xdim[-which.chain]) perm <- c(which.iter, (1:ndim)[-c(which.iter, which.chain)], which.chain) x <- aperm(x, perm) for (i in 1:nchain) { ans[[i]] <- mcmc(matrix(x[1:len + (i - 1) * len], nrow = niter)) } ans <- mcmc.list(ans) } return(ans) } environment(as.mcmc.list.mcarray)<-environment(rjags:::as.mcmc.list.mcarray) #coda ###The following functions are used for the posterior sampling and plotting: samplergpd<-function(x=jagsgpdout, newcov=NULL,N=NULL,u=threshold.gpd,mask=list(LogRate=c("INT","soi"),Sigma=c("INT","nao","sst"),Xi=c("INT","soi")),yp=c(5,10,20,50,100,200,500,1000),etafun=c(exp,exp,identity)) { if(length(dim(x))!=2) stop("need matrix for input") um<-setdiff(unlist(mask),c("INT","Yr")) #remove INT and YR if(is.null(newcov)) { stop("Need Covariate Matrix, use loss dataframe input") } if(is.matrix(newcov)) newcov<-data.frame(newcov) missingcov<-setdiff(um,names(newcov)) if(length(missingcov)) stop(paste("Missing covariates:",paste(missingcov,collapse=","),)) newcov<-data.frame(newcov[um]) ###Generate samples based on Normal distribution of covariates. This is for adding prediction uncertainty into the covariates. if(!is.null(N)) { if(verbose) print(paste("Sampling new covariates independently using first row as covariate and second row as standard deviations with",N,"samples")) if(N<=1) newcov<-newcov[1,,drop=F] else newcov<-data.frame(lapply(newcov,function(x,N=N) rnorm(N,mean=x[1],sd=x[2]))) } #For fixed quantities we just do the following: newcov$INT<-1 # parmask<-rep(names(mask),sapply(mask,length)) parmask<- split(unlist.names(mask),rep(names(mask),sapply(mask,length))) Ssplit<-lapply(names(mask),function(nn,p=parmask,pp=mask,xx=x) { y<-xx[,parmask[[nn]]]; colnames(y)<-pp[[nn]]; return(y)}) ;names(Ssplit)<-names(mask) Nsplit<-lapply(mask,function(nn,xx=newcov) { xx[,nn,drop=F] }) rates<-NULL for(i in 1:length(Ssplit)) { rates[[i]]<-as.matrix(Ssplit[[i]])%*%t(as.matrix(Nsplit[[i]])) } names(rates)<-names(mask) dimnames(rates[[1]])<-list(Sample=paste(1:nrow(x)),Predictor=make.add.text(cbind(newcov[,um]))) ratesa<-as.array.list(rates,name="Parameter") ###For each sample Year (thats nsamples from WinBugs by M covariates) returnlevels<-apply(ratesa,c(1,2), returnlevels.sample,u=u,yp=yp,etafun=etafun) #Function will always return a 3 dimensional array , each sample generates only 1 return level. #We could also generate seasonal samples, but this does not give us enouogh seasons to examine return levels. dn<-dimnames(returnlevels)[2:3] dimnames(returnlevels)<-c(list(yp=yp),dn) return(aperm(returnlevels,c(2,3,1))) } returnlevels.sample<-function(x,u,yp,etafun=c(exp,exp,identity)) { #x[1]=tc, x[2]=ls, x[3]=xi rate<-etafun[[1]](x[1]) sigma<-etafun[[2]](x[2]) xi<-etafun[[3]](x[3]) if(xi==0) return(u+sigma*log(rate*yp)) else return(u+sigma/xi*((rate*yp)^xi-1)) } make.add.text<-function (x, sep = "=", collapse = ",", digits = 3) { sapply(1:nrow(x),function(i,xx=x,ss=sep,cc=collapse,dd=digits,cn=colnames(xx)) paste(cn, round(xx[i,], dd), sep = ss, collapse = cc)) } jagsgpdposterior<-function(newcov=NULL,x=jagsgpdoutput,N=NULL,thin=1,u=threshold.gpd,yp=c(5,10,20,50,100,200,500,1000),plot=c("Boxplot","Dotplot"),mask=parmaskgpd,etafun=c(exp,exp,identity),mfrow=c(2,3),probs=c(.01,.025,.05,.25,.5,.75,.9,.95,.975,.99),yrange=c(8,14),yl=14,raw=T,type=field,newplot=F) { #Used to generate the samples from the posterior samples. if(thin>1) x<-x[seq(1,nrow(x),round(thin)),] deviance<-try(x[,"deviance"],silent=T) if(class(deviance)=="try-error") stop("deviance is missing from dataset") mindeviance<-which(min(deviance)==deviance) returnlevels<-samplergpd(x=x, newcov=newcov,N=N,u=u,mask=mask,yp=yp,etafun=etafun) n<-nrow(newcov) murl<-apply(returnlevels,c(2,3),mean) rlsmle<-apply(returnlevels,c(2,3),function(x,index=mindeviance) mean(x[index])) stats<-apply(returnlevels,c(2,3),quantile,probs=probs) if(newplot) x11() par(mfrow=mfrow) coltext<-dimnames(stats)[[2]] plot<-match.arg(plot) if(plot=="Boxplot") { for(i in 1:n) { add.text=paste("\n for",coltext[i]) rls<-returnlevels[,i,] tomboxes(rls,yr=yrange,title=paste("Return level quantiles with MLE for",type,"log loss from POT model",add.text=add.text),ylab= substitute(paste(Log[10]," 2005 Dollars")),col="blue2",mufun=mean,yl=yl) points(colnames(rls),rlsmle[i,,drop=T],col="red2",pch=16) abline(h=u,col="green") legend(x=1.2,y=yrange[2]+.2,legend=c("Posterior Mean","MLE","Threshold"),pch=16,col=c("blue2","red2","green"),bg="white") } } if(raw) out<-list(quantiles=stats,mle=rlsmle,postmean=murl,cov=newcov,returnlevels=returnlevels) else out<-list(quantiles=stats,mle=rlsmle,postmean=murl,cov=newcov) if(plot=="Dotplot") plotextreme(y=out,text=T,ylim=yrange,mfrow=mfrow) out } tomboxes<-function(data, qs = c(0.975, 0.75, 0.5, 0.25, 0.025), xlab = "Years", ylab = "kt", title = "Median Climate Return Level for Gulf Coast", xr = c(1, 1000), yr = c(0, 400), cap.shrink = 1., wd = 0.2, yl = 16, col = 7,mufun=mean) { scol <- par("col") ##xr,yr x and y ranges x>=1, xlab, ylab, labels, title, main label ##data, array with numberic column ids, from the call to get.samp.full, this is the med.yrl output ##cap.shrink controls size of cap to size of mid box, wd is width, in log plots, relative to log_10(x) x <- as.numeric(colnames(data)) namex <- colnames(data) stats <- apply(data, 2, quantile, qs) statmean <- apply(data,2,mufun) plot(type = "n", xr, yr, axes = F, lab = c(9, 17, 8), xlab = xlab, ylab = ylab, log = "x") abline(v = c(1, x), lty = 2, col=scol) abline(h = seq(yr[1], yr[2], length = yl+1), lty = 2,col=scol) axis(side = 1, at = c(1, x), tick = T, pos = yr[1], lwd=1,tcl=-.8) axis(side = 2, at = seq(yr[1], yr[2], length = yl/2 + 1), labels = T, tick = F) axis(side = 2, at = seq(yr[1], yr[2], length = yl + 1), labels = F, tick = T, pos = 1 ,lwd=1,tcl=-.8) title(title) par(col = col) width <- x * wd n <- length(x) del <- matrix(width, 10, n, byrow = T) small <- c(1, 5, 6, 10) del[small, ] <- cap.shrink * del[small, ] del[1:5, ] <- - del[1:5, ] xx <- rep(x, rep(10, n)) + del #whiskers, ticks at max & min, UQ and LQ segments(xx[1:5, ], stats, xx[6:10, ], stats) #sides of boxes segments(xx[c(2, 7), ], stats[c(2, 2), ], xx[c(2, 7), ], stats[c(4, 4), ]) xx <- rep(x, rep(2, n)) segments(xx, stats[c(5, 2), ], xx, stats[c(4, 1), ]) points(x, statmean, pch = 16) par(col=scol) invisible(NULL) } unlist.names<-function(x) { #x is a list with names paste(rep(names(x),sapply(x,length)),unlist(x),sep=".") } as.array.list<-function(x,name=NULL) { dim1<-dim(x[[1]]) if(is.null(dim1)) { dim1<-length(x[[1]]) dn1<-list(names(x[[1]])) } else dn1<-dimnames(x[[1]]) lx<-length(x) dimx<-c(dim1,lx) if(!is.null(dn1)) { dimnx<-c(dn1,list(1:lx)) } else { dimnx<-lapply(dimx,function(x) 1:x) } ldim<-length(dimnx) if(!is.null(nx<-names(x))) dimnx[[ldim]]<-nx ndx<-paste("Dim",c(1:ldim),sep="") if(!is.null(name)) ndx[ldim]<-name if(!is.null(dn1) && !is.null(nxx<-names(dn1))) ndx<-c(nxx,ndx[ldim]) names(dimnx)<-ndx xx<-unlist(lapply(x, function(x) if(is.data.frame(x)) as.matrix(x) else x) ,recursive=F) dim(xx)<-dimx dimnames(xx)<-dimnx xx } plotextreme<-function(x=c(5,10,20,50,100,200,500,1000) , y=jagsPosteriorSamples,qtls=c("2.5%","5%","25%","50%","75%","95%","97.5%"),clrs=c("red","orange","green","black","green","orange","red"), text=F ,ylim=c(9,14),mfrow=c(1,1)) { par(mfrow=mfrow,las=1,ps=13) newq<-y$quantiles ncov<-y$cov lq<-length(qtls) xi<-as.character(x) for(i in 1:nrow(ncov)) { plot(x,rep(9,8),xlim=c(3,1000),ylim=ylim,type="n",bty="n",cex=1.4,log="x",ylab="Log Damage Loss",xlab="Return Period (yr)") grid(col="grey",nx=NULL,ny=NULL) abline(col="grey",v=c(5,20,50,200,500)) for(j in seq(along=x)) { points(rep(x[j],lq),newq[qtls,i,xi[j]],col=clrs,pch=19,cex=1.0) } if(text) mtext(maketext(ncov[i,]),side=3,padj=+1,line=1,cex=1.0) } xpds<-par("xpd") par(xpd=NA) if(mfrow[1]==2 && mfrow[2]==2) { x=.255 y= ylim[2]+.87*(ylim[2]-ylim[1]) } else { x= 3 y= ylim[2] } legend(list(x=x,y=y),legend=rev(qtls),pch=20,col=clrs,inset=c(-.5,+.5),plot=T,bg="white") par(xpd=xpds) } maketext<-function(x,sig=3) { x1<-names(x) x1[x1=="AMO"]<-"SST" x2<-as.character(signif(unlist(x),sig)) paste(x1,"=",x2," ",sep="",collapse=",") }