# R_code.txt : R code to perform ordinal regression for Ceriodaphnia
# epibiont burden (example corresponding to top panel in Table 2 of the
# manuscript)
# 
# Requires JAGS software (http://mcmc-jags.sourceforge.net) and rjags R
# package 
# The rjags interface calls JAGS code from file JAGSmodel.txt 
# Data is read from file supplement4.csv: 1638 observations of 10 
# variables 
# 
# MAR - 18/02/2015 
#################################################################

data <- read.csv("supplement4.csv",head=TRUE)	

# define linear model object to obtain initial values for betas 
lmSol <- lm(y ~ cllen + wk2 + wk3 + depth + fish + depth_fish +
		wk2_fish + wk3_fish,x=TRUE,y=TRUE, data=data) 

# generate input for JAGS model 
jagsInp <- list(y=lmSol$y,x=lmSol$x[,-1], encId=data$encId,
 		N=length(lmSol$y), Mu=rep(0,8),Sig=diag(1.e-08,8)) 

# assign initial parameter values for the two chains 
inits <- list(list(beta=coef(lmSol)[-1]+rnorm(8,0,0.5),
 		sigma=runif(1,0,5), tau0=sort(rnorm(4,0,10))), 
		list(beta=coef(lmSol)[-1]+rnorm(8,0,0.5),
		sigma=runif(1,0,5), tau0=sort(rnorm(4,0,10))))

# load rjags interface to JAGS; CODA package is loaded by default 
require(rjags) 
load.module("dic") 
load.module("glm") 

# adaptive MCMC phase, 1000 iterations 
mod <- jags.model(file="JAGSmodel.txt",data=jagsInp,
		n.adapt=1000,inits=inits, n.chains=2) 

# 100000 MCMC iterations; retain 1 in 50 (thinning) 
n.iter=100000; thin=50 
# results are stored as MCMC list out 
out <- coda.samples(mod, variable.names=c("beta","tau","sigma",
								"psi","deviance", 
								"yrep"),
								n.iter=n.iter,thin=thin) 

# discard first half of iterations (burn in) 
start <- 1000 + thin + n.iter/2; end <- 1000 + n.iter 
halfOut <- window(out,start=start,end=end)

# summary of main parameter estimates 
pars <- 26 
summary(halfOut[,c(1:pars)]) 
plot(halfOut[,1:pars],ask=TRUE) 

# some diagnostics 
acfplot(halfOut[,c(1:pars)]) 
gelman.diag(halfOut[,c(1:pars)])