# R code for Exam 2

# Setup
library(arm)
options(digits=7)
regenerate = read.table("regenerate.txt",header=T)
regenerate$treatment = with(regenerate, reorder(regenerate$treatment,stems) )
regenerate$site = with(regenerate, reorder(regenerate$site,stems) )

# Problem 1a
print( xyplot(stems ~ treatment | site,data=regenerate) )

# Problem 1b
# Other models are acceptable
fit.1b = glm(formula = stems ~ treatment*site + spacing, family = "poisson",data = regenerate)
coef.1b = coef(fit.1b)[1:18]

# Simulation
set.seed(23532)
#sim.1b = sim(fit.1b,10000)

# Note that sim.1b$beta has 19 columns, but the 19th is ignored in the model
# Code for point estimates using predict and direct multiplication for a check
Cb.LS = exp( sum(coef.1b * c(1,0,0,0,1,0,0,1,8,0,0,0,0,0,0,0,0,1) ) )
#Cb.LSb = exp(predict(fit.1b,data.frame(site="LaSelva",treatment="Cb",spacing=8)))
Ta.LS = exp( sum(coef.1b * c(1,1,0,0,0,0,0,1,8,0,0,0,0,0,1,0,0,0) ) )
#Ta.LSb = exp(predict(fit.1b,data.frame(site="LaSelva",treatment="Ta",spacing=8)))
Cb.P4 = exp( sum(coef.1b * c(1,0,0,0,1,0,0,0,4,0,0,0,0,0,0,0,0,0) ) )
#Cb.P4b = exp(predict(fit.1b,data.frame(site="Paniagua",treatment="Cb",spacing=4)))
Ta.P4 = exp( sum(coef.1b * c(1,1,0,0,0,0,0,0,4,0,0,0,0,0,0,0,0,0) ) )
#Ta.P4b = exp(predict(fit.1b,data.frame(site="Paniagua",treatment="Ta",spacing=4)))
Cb.P32 = exp( sum(coef.1b * c(1,0,0,0,1,0,0,0,32,0,0,0,0,0,0,0,0,0) ) )
#Cb.P32b = exp(predict(fit.1b,data.frame(site="Paniagua",treatment="Cb",spacing=32)))
Ta.P32 = exp( sum(coef.1b * c(1,1,0,0,0,0,0,0,32,0,0,0,0,0,0,0,0,0) ) )
#Ta.P32b = exp(predict(fit.1b,data.frame(site="Paniagua",treatment="Ta",spacing=32)))
Cb.Q = exp( sum(coef.1b * c(1,0,0,0,1,0,1,0,16,0,0,0,1,0,0,0,0,0) ) )
#Cb.Qb = exp(predict(fit.1b,data.frame(site="Quesada",treatment="Cb",spacing=16)))
Ta.Q = exp( sum(coef.1b * c(1,1,0,0,0,0,1,0,16,1,0,0,0,0,0,0,0,0) ) )
#Ta.Qb = exp(predict(fit.1b,data.frame(site="Quesada",treatment="Ta",spacing=16)))

# Code to get the 95% confidence intervals
# La Selva, spacing = 8

print(round(quantile(exp(sim.1b$beta[,1:18] %*% c(1,0,0,0,1,0,0,1,8,0,0,0,0,0,0,0,0,1))
               - exp(sim.1b$beta[,1:18] %*% c(1,1,0,0,0,0,0,1,8,0,0,0,0,0,1,0,0,0)),c(0.025,0.5,0.975)),1))

# Paniagua, spacing = 4
print(round(quantile(exp(sim.1b$beta[,1:18] %*% c(1,0,0,0,1,0,0,0,4,0,0,0,0,0,0,0,0,0))
               - exp(sim.1b$beta[,1:18] %*% c(1,1,0,0,0,0,0,0,4,0,0,0,0,0,0,0,0,0)),c(0.025,0.5,0.975)),1))

# Paniagua, spacing = 32
print(round(quantile(exp(sim.1b$beta[,1:18] %*% c(1,0,0,0,1,0,0,0,32,0,0,0,0,0,0,0,0,0))
               - exp(sim.1b$beta[,1:18] %*% c(1,1,0,0,0,0,0,0,32,0,0,0,0,0,0,0,0,0)),c(0.025,0.5,0.975)),1))

# Quesada, spacing = 16
print(round(quantile(exp(sim.1b$beta[,1:18] %*% c(1,0,0,0,1,0,1,0,16,0,0,0,1,0,0,0,0,0))
               - exp(sim.1b$beta[,1:18] %*% c(1,1,0,0,0,0,1,0,16,1,0,0,0,0,0,0,0,0)),c(0.025,0.5,0.975)),1))


#
# Problem 2
#

# Read in the data

hamster = read.table("hamster.txt",header=T)
hamster$id = with(hamster, reorder(hamster$id,10*(day=="Long")+conc))

# Problem 2a
dotplot(id ~ conc,groups=day,data=hamster,auto.key=list(space="right"))

# Problem 2b
fit.2b = lmer(conc ~ day + (1 | id), data=hamster)
display(fit.2b,digits=3)
fit.2bi = lmer(conc ~ day - 1 + (1 | id), data=hamster)
display(fit.2bi,digits=3)

# Problem 2c
set.seed(34323)
sim.2b = mcmcsamp(fit.2b,10000)
sim.2b[1,]
sigma = sqrt(exp(sim.2b[,3]))
sigmaA = sqrt(exp(sim.2b[,4]))
hamEff1 = rnorm(10000,0,sigmaA)
meas1 = sim.2b[,1] + hamEff1 + rnorm(10000,0,sigma)
meas2 = sim.2b[,1] + hamEff1 + rnorm(10000,0,sigma)
quantile(meas1-meas2,c(0.025,0.975))

# Problem 2d
muLong = rnorm(10000,sim.2b[,1],sigmaA)
muShort = rnorm(10000,sim.2b[,1] + sim.2b[,2],sigmaA)
xbarLong = rnorm(10000,muLong,sigma/2)
xbarShort = rnorm(10000,muShort,sigma/2)
round(quantile(xbarLong - xbarShort,c(0.025,0.975)),3)