########################################################################### #Functions to perform MCMC output analysis using CODA. #Author: Ricardo Ehlers Nov 2005 # #Read in MCMC output in the CODA format (produced by WinBUGS/OpenBUGS/JAGS) # #library(coda) # #x = as.mcmc(read.coda(output,index,quiet=T)), where #output: one column file containing the monitored MCMC output #index: file showing which rows of the output file correspond to which variable # #If your output file contains a matrix with one column per variable simply read #it as an MCMC object # #x = as.mcmc(read.table(file=outputfile,header=T)) # ########################################################################### plot.mcmc <- function(x,saveplot=F,smooth=F,file=NULL,type='trace',ask=T, mfcol=NULL,add=F,select=NULL,...){ ########################################################################### #Traces, acf, histograms and density estimates of MCMC output #Set saveplot=T and provide file name for redirecting to a Postscript file. #Use select to specify a subset of variables to monitor. #Ricardo Ehlers Nov 2005 ########################################################################### if (!is.mcmc(x)) x=as.mcmc(x) # stop("argument not of class mcmc\n") if (length(select)>0) x=x[,select] nvar=nvar(x) names=colnames(x) if (!add) { if (saveplot) postscript(file=file,horizontal=F) if (length(mfcol)==0){ if (type=='all') { mfcol=c(3,nvar) } else { mfcol=coda:::set.mfrow(Nchains = nchain(x),Nparms=nvar,nplots=1) } } oldpar = par(mfcol = mfcol) oldpar = c(oldpar, par(ask = ask)) } else { dev.set(dev.cur()) } #if (length(type)==1) { if (type=='trace'){ traceplot(x,smooth=smooth) } else if (type=='autocorr'){ autocorr.plot(x) } else if (type=='hist') { for (j in 1:nvar) hist(x[,j],prob=T,xlab='',main=names[j]) } else if (type=='dens') { for (j in 1:nvar) densplot(x[,j]) } else if (type=='all') { for (j in 1:nvar) { traceplot(x[,j],smooth=smooth) densplot(x[,j]) autocorr.plot(x[,j],auto.layout=F) } } #} if (saveplot) dev.off() } diag.mcmc <- function(x,select=NULL){ ################################################################### #MCMC Diagnostics. #Use select to specify a subset of variables to monitor. #Effective sample size: sample size adjusted for autocorrelation by #estimating the spectral density at frequency 0 using function #spectrum0.ar #Ricardo Ehlers Nov 2005 ################################################################### if (!is.mcmc(x)) stop("argument not of class mcmc\n") if (length(select)>0) x=x[,select] cat('\nEffective sample size for estimating the mean\n') print(effectiveSize(x)) cat('\nGeweke convergence diagnostic\n') print(geweke.diag(x)) cat('\nAutocorrelations\n') print(round(autocorr.diag(x),5)) cat('\nHeidelberger and Welch convergence diagnostic\n') print(heidel.diag(x)) cat('\nRaftery and Lewis diagnostic\n') print(raftery.diag(x)) if (nchain(x)>1){ cat('\nGelman and Rubin convergence diagnostic\n') print(gelman.diag(x)) } } print.summary <- function(x,digits=2,HPD=FALSE,summ=NULL){ ########################################################## #Summary of model parameters using sample x #Ricardo Ehlers, Jun 2006 ########################################################## if (!is.mcmc(x)) stop("argument not of class mcmc\n") options(width=85) if (is.null(summ)) { b=summary(x) } else { b=summ } colnames(b$statistics)=c("Mean","SD","Naive SE","TS SE") print(round(cbind(b$statistics,b$quantiles[,c(1,3,5)]),digits)) if (HPD) { cat('\nHighest Posterior Density Intervals\n') print(round(HPDinterval(x),digits)) } }