wardClones.R

library("BRugs")
library("R2WinBUGS")

# CODE WRITTEN BY ERIC WARD, NOV. 2007
# !! make sure to download OpenBUGS, install to default c:\ directory
# data file is in current wd, named "practice.txt"

###################################
# STEP 1: SET PARAMETERS
###################################
numClones<-25;  # must be at least 2
chainLength <- 5000;
burn <- 5000; # should be at least 5000

###################################
# STEP 2: SET UP MCMC PARAMETERS
###################################
dat <- read.table("practice.txt",header=FALSE);
dat[which(dat[,1] < 0),1] <- NA;
dat <- log(dat);

N <- dim(dat)[1];  # number of years
Y <- matrix(0,N,numClones);
for(i in 1:numClones) Y[,i] <- dat[,1];


###################################
# THIS SECTION WRITES THE BUGS CODE
###################################
modelString <- "
model {
   # first year
   for(j in 1:numClones) {
      init[j] ~ dunif(0,10);
   }
   # parameter B
   B ~ dunif(-3,3);
   # Q is process error
   logsigma2.Q ~ dunif(-13.81551,2);
   sigma2.Q <- exp(logsigma2.Q);
   sigma.Q <- sqrt(sigma2.Q);
   invQ <- 1/sigma2.Q;
   # process error model
   for(j in 1:numClones) {
      states[1,j] <- init[j];
      for(i in 2:N) {
         mu[i,j] <- (states[i-1,j] + B);
         states[i,j] ~ dnorm(mu[i,j],invQ);
      }
   }

   # R is observation error
   logsigma2.R ~ dunif(-13.81551,2);
   sigma2.R <- exp(logsigma2.R);
   sigma.R <- sqrt(sigma2.R);
   invR <- 1/sigma2.R;

   for(j in 1:numClones) {
      for(i in 1:N) {
         Y[i,j] ~ dnorm(states[i,j], invR)
      }
   }

}
"

# Write the bugs model to a file
write(modelString,file="cloneModel.txt",append=FALSE)

data <- list("Y","numClones","N");
parameters <- c("B","sigma.Q","sigma.R")
simple <- 0;
try(simple <- bugs (data,inits=NULL,parameters,DIC=TRUE,"cloneModel.txt",n.chains=1,n.burnin=burn, 
   n.thin=1,n.iter=(burn+chainLength),program="openbugs",debug=TRUE),silent=FALSE)
try(attach.bugs (simple, overwrite=TRUE), silent = TRUE)

# calculate the deviance at the mean posterior
pD <- simple$pD;
Dbar <- mean(deviance);
# eval deviance at posterior means
Dhat <- (Dbar - pD)/numClones;