# Analysis of kiwi data set

# Read in the data

kiwi = read.table("kiwi.txt",header=T)

# Examine structure

str(kiwi)

# Plot yield versus shade, group by block.

library(lattice)
print( xyplot(yield ~ shade | block, data=kiwi, aspect=1, layout = c(3,1), groups=shade) )

# Reorder factors by yield to make trends easier to see

kiwi$shade = with(kiwi, reorder(shade,yield))
kiwi$block = with(kiwi, reorder(block,yield))

# Replot
print( xyplot(yield ~ shade | block, data=kiwi, aspect=1, layout = c(3,1), groups=shade) )

# Load arm library and fix the digits option
library(arm)
options(digits=7)

# For comparison, here we fit a standard model
fit.0 = lm(yield ~ shade + block + plot,data=kiwi)
cat("fit.0: standard model\n")
display(fit.0)

# Fit a multilevel model
# Block level variation
# Shade level variation
# Plot level variation
# Vine level variation

fit.1 = lmer(yield ~ (1 | shade) + (1 | block) + (1 | plot),data=kiwi)
cat("fit.1: multi-level model\n")
display(fit.1)

# Alternative model with shade as fixed
# Here we remove the intercept so that the shade effects are easy to see
fit.2 = lmer(yield ~ shade - 1 + (1 | block) + (1 | plot), data=kiwi)
cat("fit.2: model with shade fixed \n")
display(fit.2)

# Prediction interval for difference between two measurements at two vines
# in the no shade treatment in the north block

# See notes from the next lecture for more details about prediction intervals.