dat <- hughes.grapes
dat <- transform(dat, rate = diseased/total, plot=trt:block)
# Trt 1 has higher rate, more variable, Trt 3 lower rate, less variable
foo <- bwplot(rate ~ vine|block*trt, dat)
if (require(latticeExtra)) useOuterStrips(foo) else print(foo)
# Table 1 of Piepho 1999
tapply(dat$rate, dat$trt, mean) # trt 1 does not match Piepho
tapply(dat$rate, dat$trt, max)
# Piepho model 3. Binomial data. May not be exactly the same model
# Use the binomial count data with lme4
library(lme4)
m1 <- lmer(cbind(diseased, total-diseased) ~ trt + block + (1|plot/vine),
data=dat, family=binomial)
m1
# Switch from binomial counts to bernoulli data
require(aod)
bdat <- splitbin(cbind(diseased, total-diseased) ~ block+trt+plot+vine+shoot,
data=dat)$tab
names(bdat)[2] <- 'y'
# Using lme4
m2 <- lmer(y ~ trt + block + (1|plot/vine), data=bdat, family=binomial)
m2
# Now using MASS:::glmmPQL
require(MASS)
m3 <- glmmPQL(y ~ trt + block, data=bdat,
random=~1|plot/vine, family=binomial)
m3Run the code above in your browser using DataLab