if (FALSE) {
library(agridat)
data(hughes.grapes)
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
libs(lattice)
foo <- bwplot(rate ~ vine|block*trt, dat, main="hughes.grapes",
xlab="vine")
libs(latticeExtra)
useOuterStrips(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
libs(lme4)
m1 <- glmer(cbind(diseased, total-diseased) ~ trt + block + (1|plot/vine),
data=dat, family=binomial)
m1
# Switch from binomial counts to bernoulli data
libs(aod)
bdat <- splitbin(cbind(diseased, total-diseased) ~ block+trt+plot+vine+shoot,
data=dat)$tab
names(bdat)[2] <- 'y'
# Using lme4
m2 <- glmer(y ~ trt + block + (1|plot/vine), data=bdat, family=binomial)
m2
# Now using MASS:::glmmPQL
libs(MASS)
m3 <- glmmPQL(y ~ trt + block, data=bdat,
random=~1|plot/vine, family=binomial)
m3
}
Run the code above in your browser using DataLab