agridat (version 1.16)

besag.elbatan: RCB experiment of 50 varieties of wheat in 3 blocks with strong spatial trend.

Description

RCB experiment of 50 varieties of wheat in 3 blocks with strong spatial trend.

Format

A data frame with 150 observations on the following 4 variables.

yield

yield of wheat

gen

genotype, factor with 50 levels

col

column/block

row

row

Details

RCB experiment on wheat at El Batan, Mexico. There are three single-column replicates with 50 varieties in each replicate.

Plot dimensions are not given by Besag.

References

Wilkinson 1984.

Besag & Seheult 1989.

Examples

# NOT RUN {
data(besag.elbatan)
dat <- besag.elbatan

if(require(desplot)){
  desplot(yield~col*row, dat,
          num=gen, # aspect unknown
          main="besag.elbatan - wheat yields")
}

# }
# NOT RUN {
  # Besag figure 1
  library(lattice)
  xyplot(yield~row|col, dat, type=c('l'),
         layout=c(1,3), main="besag.elbatan wheat yields")
# }
# NOT RUN {
# RCB
m1 <- lm(yield ~ 0 + gen + factor(col), dat)
p1 <- coef(m1)[1:50]

# }
# NOT RUN {
  # Formerly used gam package, but as of R 3.1, Rcmd check --as-cran
  # is complaining
  # Calls: plot.gam ... model.matrix.gam -> predict -> predict.gam -> array
  # but it works perfectly in interactive mode !!!
  # Remove the FALSE to run the code below
  if(is.element("gam", search())) detach(package:gam)
  library(mgcv)
  m2 <- mgcv::gam(yield ~ -1 + gen + factor(col) + s(row), data=dat)
  plot(m2, residuals=TRUE, main="besag.elbatan")
  pred <- cbind(dat, predict(m2, dat, type="terms"))
  # Need to correct for the average loess effect, which is like
  # an overall intercept term.
  adjlo <-  mean(pred$"s(row)")
  p2 <- coef(m2)[1:50] + adjlo

  # Compare estimates
  lims <- range(c(p1,p2))
  plot(p1, p2, xlab="RCB prediction",
       ylab="RCB with smooth trend (predicted)",
       type='n', xlim=lims, ylim=lims,
       main="besag.elbatan")
  text(p1, p2, 1:50, cex=.5)
  abline(0,1,col="gray")
# }
# NOT RUN {
# }