agridat (version 1.16)

besag.bayesian: Spring barley in United Kingdom

Description

An experiment with 75 varieties of barley, planted in 3 reps.

Format

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

col

column (also blocking factor)

row

row

yield

yield

gen

variety/genotype

Details

RCB design, each column is one rep.

References

Davison, A. C. (2003). Statistical Models. Cambridge University Press. Pages 534-535.

Examples

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

# Yield values were scaled to unit variance
var(dat$yield, na.rm=TRUE)

# Besag Fig 2. Reverse row numbers to match Besag, Davison
dat$rrow <- 76 - dat$row
require("lattice")
xyplot(yield ~ rrow|col, dat, layout=c(1,3), type='s',
       xlab="row", ylab="yield", main="besag.bayesian")

# ----------------------------------------------------------------------------

# }
# NOT RUN {
  # asreml3
  require(asreml)
  
  # Use asreml to fit a model with AR1 gradient in rows
  
  dat <- transform(dat, cf=factor(col), rf=factor(rrow))
  m1 <- asreml(yield ~ -1 + gen, data=dat, random=~ar1v(rf))
  m1 <- update(m1)
  
  # Visualize trends, similar to Besag figure 2.
  dat$res <- m1$residuals
  dat$geneff <- coef(m1)$fixed[as.numeric(dat$gen)]
  dat <- transform(dat, fert=yield-geneff-res)
  xyplot(geneff ~ rrow|col, dat, layout=c(1,3), type='s',
         main="Variety effects", ylim=c(5,15 ))
  xyplot(fert ~ rrow|col, dat, layout=c(1,3), type='s',
         main="Fertility", ylim=c(-2,2))
  xyplot(res ~ rrow|col, dat, layout=c(1,3), type='s',
         main="Residuals", ylim=c(-4,4))
  
# }
# NOT RUN {
# ----------------------------------------------------------------------------

# }
# NOT RUN {
  ## require(asreml4)
  
  ## # Use asreml to fit a model with AR1 gradient in rows
  
  ## dat <- transform(dat, cf=factor(col), rf=factor(rrow))
  ## m1 <- asreml(yield ~ -1 + gen, data=dat, random=~ar1v(rf))
  ## m1 <- update(m1)
  ## m1 <- update(m1)
  ## m1 <- update(m1)
  ## m1 <- update(m1)
  ## m1 <- update(m1)
  ## m1 <- update(m1)
  
  ## # Visualize trends, similar to Besag figure 2.
  ## dat$res <- resid(m1)
  ## dat$geneff <- coef(m1)$fixed[as.numeric(dat$gen)]
  ## dat <- transform(dat, fert=yield-geneff-res)
  ## xyplot(geneff ~ rrow|col, dat, layout=c(1,3), type='s',
  ##        main="Variety effects", ylim=c(5,15 ))
  ## xyplot(fert ~ rrow|col, dat, layout=c(1,3), type='s',
  ##        main="Fertility", ylim=c(-2,2))
  ## xyplot(res ~ rrow|col, dat, layout=c(1,3), type='s',
  ##        main="Residuals", ylim=c(-4,4))

# }
# NOT RUN {
# }