if (FALSE) {
library(agridat)
data(mcconway.turnip)
dat <- mcconway.turnip
dat$densf <- factor(dat$density)
# Table 2 of Lee et al.
m0 <- aov( yield ~ gen * densf * date + block, dat )
summary(m0)
## Df Sum Sq Mean Sq F value Pr(>F)
## gen 1 84.0 83.95 8.753 0.00491 **
## densf 3 470.4 156.79 16.347 2.51e-07 ***
## date 1 233.7 233.71 24.367 1.14e-05 ***
## block 3 163.7 54.58 5.690 0.00216 **
## gen:densf 3 8.6 2.88 0.301 0.82485
## gen:date 1 36.5 36.45 3.800 0.05749 .
## densf:date 3 154.8 51.60 5.380 0.00299 **
## gen:densf:date 3 18.0 6.00 0.626 0.60224
## Residuals 45 431.6 9.59
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Boxplots suggest heteroskedasticity for date, density
libs("HH")
interaction2wt(yield ~ gen + date + densf +block, dat,
x.between=0, y.between=0,
main="mcconway.turnip - yield")
libs(nlme)
# Random block model
m1 <- lme(yield ~ gen * date * densf, random= ~1|block, data=dat)
summary(m1)
anova(m1)
# Multiplicative variance model over densities and dates
m2 <- update(m1,
weights=varComb(varIdent(form=~1|densf),
varIdent(form=~1|date)))
summary(m2)
anova(m2)
# Unstructured variance model over densities and dates
m3 <- update(m1, weights=varIdent(form=~1|densf*date))
summary(m3)
anova(m3)
# Table 3 of Piepho, using transformation
m4 <- aov( yield^.235 ~ gen * date * densf + block, dat )
summary(m4)
}
Run the code above in your browser using DataLab