if (FALSE) {
library(agridat)
data(holshouser.splitstrip)
dat <- holshouser.splitstrip
dat$spacing <- factor(dat$spacing)
dat$pop <- factor(dat$pop)
# Experiment layout and field trends
libs(desplot)
desplot(dat, yield ~ col*row,
out1=block, # unknown aspect
main="holshouser.splitstrip")
desplot(dat, spacing ~ col*row,
out1=block, out2=cultivar, # unknown aspect
col=cultivar, text=pop, cex=.8, shorten='none', col.regions=c('wheat','white'),
main="holshouser.splitstrip experiment design")
# Overall main effects and interactions
libs(HH)
interaction2wt(yield~cultivar*spacing*pop, dat,
x.between=0, y.between=0,
main="holshouser.splitstrip")
## Schabenberger's SAS model, page 497
## proc mixed data=splitstripplot;
## class block cultivar pop spacing;
## model yield = cultivar spacing spacing*cultivar pop pop*cultivar
## spacing*pop spacing*pop*cultivar / ddfm=satterth;
## random block block*cultivar block*cultivar*spacing block*cultivar*pop;
## run;
## Now lme4. This design has five error terms--four are explicitly given.
libs(lme4)
libs(lucid)
m1 <- lmer(yield ~ cultivar * spacing * pop +
(1|block) + (1|block:cultivar) + (1|block:cultivar:spacing) +
(1|block:cultivar:pop), data=dat)
vc(m1) ## Variances match Schabenberger, page 498.
## grp var1 var2 vcov sdcor
## block:cultivar:pop (Intercept) 2.421 1.556
## block:cultivar:spacing (Intercept) 1.244 1.116
## block:cultivar (Intercept) 0.4523 0.6725
## block (Intercept) 3.037 1.743
## Residual 3.928 1.982
}
Run the code above in your browser using DataLab