if (FALSE) {
library(agridat)
data(besag.triticale)
dat <- besag.triticale
dat <- transform(dat, rate=factor(rate), nitro=factor(nitro))
dat <- transform(dat, xf=factor(col), yf=factor(row))
libs(desplot)
desplot(dat, yield ~ col*row,
# aspect unknown
main="besag.triticale")
# Besag & Kempton are not perfectly clear on the model, but
# indicate that there was no evidence of any two-way interactions.
# A reduced, main-effect model had genotype effects that were
# "close to significant" at the five percent level.
# The model below has p-value of gen at .04, so must be slightly
# different than their model.
m2 <- lm(yield ~ gen + rate + nitro + regulator + yf, data=dat)
anova(m2)
# Similar, but not exact, to Besag figure 5
dat$res <- resid(m2)
libs(lattice)
xyplot(res ~ col|as.character(row), data=dat,
as.table=TRUE, type="s", layout=c(1,3),
main="besag.triticale")
if(require("asreml", quietly=TRUE)) {
libs(asreml)
# Besag uses an adjustment based on neighboring plots.
# This analysis fits the standard AR1xAR1 residual model
dat <- dat[order(dat$xf, dat$yf), ]
m3 <- asreml(yield ~ gen + rate + nitro + regulator +
gen:rate + gen:nitro + gen:regulator +
rate:nitro + rate:regulator +
nitro:regulator + yf, data=dat,
resid = ~ ar1(xf):ar1(yf))
wald(m3) # Strongly significant gen, rate, regulator
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 1288255 103.971 < 2.2e-16 ***
## gen 2 903262 72.899 < 2.2e-16 ***
## rate 1 104774 8.456 0.003638 **
## nitro 1 282 0.023 0.880139
## regulator 2 231403 18.676 8.802e-05 ***
## yf 2 3788 0.306 0.858263
## gen:rate 2 1364 0.110 0.946461
## gen:nitro 2 30822 2.488 0.288289
## gen:regulator 4 37269 3.008 0.556507
## rate:nitro 1 1488 0.120 0.728954
## rate:regulator 2 49296 3.979 0.136795
## nitro:regulator 2 41019 3.311 0.191042
## residual (MS) 12391
}
}
Run the code above in your browser using DataLab