dat <- besag.elbatan
# Besag figure 1
xyplot(yield~row, dat, groups=block, type=c('l'),
main="besag.elbatan")
desplot(yield~block*row, dat, main="besag.elbatan")
# RCB
m1 <- lm(yield ~ 0 + gen + block, dat)
p1 <- coef(m1)[1:50]
# Add smooth trend with GAM
require(gam)
m2 <- gam(yield ~ 0 + gen + block + lo(row), data=dat)
plot(m2, residuals=TRUE, se=TRUE, col=dat$block)
p2 <- coef(m2)[1:50]
# Compare estimates
plot(p1, p2, xlab="RCB", ylab="RCB with smooth trend", type='n')
text(p1, p2, 1:50, cex=.5)
abline(0,1,col="gray")Run the code above in your browser using DataLab