# R script used to fit Bayesian changepoint models (via WinBUGS) for POD data. # Jim Thomson. Monash University. August 2008 # requiers WinBUGS 1.4 with "jump" add-in to be installed # requires R package R2WinBUGS setwd("C:/NCEAS") #set working directory library(R2WinBUGS) stdz=function(x){ return((x-mean(x))/sd(x)) } data=read.delim("poddata.txt") data=data[1:40,] species.cols=c(2,3,4,6,8,9,10,12,14,15,17,18,19,21,23,24,25,27,31,33,35,37,39,41,43) haserr=c(4,6,10,12,15,19,21,25,27,31,33,35,37,39,41,43) colnames(data)[species.cols] test.sp=species.cols[c(1,3,5,7,9:11,13)] probs=NULL ppps=NULL fitted=NULL betas=NULL Rhats=NULL for(spcol in test.sp){ #loop thru response variables spname=colnames(data)[spcol] y=data[,spcol] errest=length(which(haserr==spcol)) if(errest>0){ se=ifelse(data[,spcol+1]==0,0.0001,data[,spcol+1]) yvar=1/(se^2) yvar=ifelse(is.na(yvar),mean(yvar,na.rm=T),yvar) } if(errest==0){yvar=rep(10000,length(y))} # time=data[,1]-1967 start=min(which(is.finite(y))) time=time[start:length(y)] yvar=yvar[start:length(y)] y=y[start:length(y)] N=length(y) T=N Xtime=NULL for(t in 1:(T-2)){ Xtime=cbind(Xtime,c(rep(0,t),rep(1,T-1-t))) } covar=solve(t(Xtime)%*%Xtime) scalef=mean(diag(covar))*3.96 kmax=c(5,5) # maximum number of step and slope changes pc=c(.5,.5) # prior proportion of kmax years that have a changepoint order=1 # order of piecewise spline for trend bugsdat=list("N","y","time","kmax","pc","yvar","scalef","order") #list of data sent to WinBUGS parameters=c("k","mu","y.act","alpha","beta","change","ppp")# nodes to monitor in WinBUGS yspan=diff(range(y,na.rm=T)) # for scaling priors prior.names=c("gamma","point","uniform","cauchy") for(prior.type in 1:4){ # loop thru (or specify here) different priors for(scaleff in c(0.5,1,2)){ # loop thru relative scaling factors for priors if(prior.type==1){ # gamma prior inits=function(){list(gmean=1,beta.prec=c(1,1),tau=1)} #initial values scalef=ifelse(scaleff==0.5,0.1,ifelse(scaleff==1,.01,.001))#scaling factor for priors } if(prior.type==2){ # point prior inits=function(){list(gmean=1,tau=1)} scalef=c(1/((scaleff*yspan/1.96)^2),1/(((scaleff*yspan)/(4*1.96))^2)) # prior prob 0.95 magnitude < scalef } if(prior.type==3){ # unif prior inits=function(){list(gmean=1,beta.sd=scalef/2,tau=1)} scalef=c(0.8*yspan*scaleff,0.8*scaleff*yspan/4) } if(prior.type==4){ # half-Caucy prior inits=function(){list(gmean=1,tauZ.eta=c(1,1),tau=1)} scalef=c(100/((scaleff*yspan)^2),100/(((scaleff*yspan)/4)^2)) } model.name=paste("changepoints_spline_",prior.names[prior.type],"_obserr.txt",sep="") # fit model in WinBUGS # fit=bugs(bugsdat,inits,parameters,model.name,n.iter=500000, debug=F, working.directory="C:/NCEAS/WinBUGS model files") # store results in various matrices probs=cbind(probs,rbind(matrix(NA,ncol=2,nrow=40-N),t(fit$mean$change))) ppps=rbind(ppps,c(spname,fit$mean$ppp)) fitted=cbind(fitted,c(rep(NA,40-N),fit$mean$mu),c(rep(NA,40-N),fit$mean$y.act)) betas=cbind(betas,c(rep(NA,40-N),fit$mean$alpha),c(rep(NA,40-N),fit$mean$beta)) Rhats=rbind(Rhats,c(spname,quantile(fit$summary[,"Rhat"]))) colnames(probs)[(ncol(probs)-1):ncol(probs)]=spname colnames(fitted)[(ncol(fitted)-1):ncol(fitted)]=spname colnames(betas)[(ncol(betas)-1):ncol(betas)]=spname } } } # PLOT ALL RESULTS par(mfrow=c(2,2)) time=data[,1] plotcols=seq(1,ncol(fitted)-1,2) probc=1 for(c in plotcols){ y=data[,colnames(betas)[c]] plotdata=cbind(fitted[,c:(c+1)],betas[,c:(c+1)]) plot(y~time,main=colnames(fitted)[c]) lines(plotdata[,1]~time) lines(plotdata[,3]~time,lty=2,col="red") lines(plotdata[,4]+mean(plotdata[,3],na.rm=T)~time,lty=2,col="blue") plot(y~time,main=colnames(fitted)[c]) lines(plotdata[,2]~time) lines(plotdata[,3]~time,lty=2,col="red") lines(plotdata[,4]+mean(plotdata[,3],na.rm=T)~time,lty=2,col="blue") plot(c(0,probs[,c])~time, type="h", main=colnames(probs)[probc],ylim=c(0,1)) plot(c(0,probs[,c+1])~time,type="h", main=colnames(probs)[probc],ylim=c(0,1)) } ##SAVE RESULTS - change NAME probs.NAME=probs fitted.NAME=fitted betas.NAME=betas save(fitted.NAME,file="fitted.NAME.R") save(betas.NAME,file="betas.NAME.R") save(probs.NAME,file="probs.NAME.R")