if (FALSE) {
library(agridat)
data(vaneeuwijk.drymatter)
dat <- vaneeuwijk.drymatter
dat <- transform(dat, year=factor(year))
dat <- transform(dat, env=factor(paste(year,site)))
libs(HH)
HH::interaction2wt(y ~ year+site+variety,dat,rot=c(90,0),
x.between=0, y.between=0,
main="vaneeuwijk.drymatter")
# anova model
m1 <- aov(y ~ variety+env+variety:env, data=dat)
anova(m1) # Similar to VanEeuwijk table 2
m2 <- aov(y ~ year*site*variety, data=dat)
anova(m2) # matches Sahai table 5.5
# variance components model
libs(lme4)
libs(lucid)
m3 <- lmer(y ~ (1|year) + (1|site) + (1|variety) +
(1|year:site) + (1|year:variety) + (1|site:variety),
data=dat)
vc(m3) # matches Sahai page 266
## grp var1 var2 vcov sdcor
## year:variety (Intercept) 0.3187 0.5645
## year:site (Intercept) 7.735 2.781
## site:variety (Intercept) 0.03502 0.1871
## year (Intercept) 6.272 2.504
## variety (Intercept) 0.4867 0.6976
## site (Intercept) 6.504 2.55
## Residual 0.8885 0.9426
}
Run the code above in your browser using DataLab