## Not run: 
  
  library(agridat)
  data(besag.bayesian)
  dat <- besag.bayesian

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

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

  libs(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)
  m1 <- update(m1)
  m1 <- update(m1)
    
  # Visualize trends, similar to Besag figure 2.
  # Need 'as.vector' because asreml4 uses a named vector
  dat$res <- unname(m1$resid)
  dat$geneff <- coef(m1)$fixed[as.numeric(dat$gen)]
  dat <- transform(dat, fert=yield-geneff-res)
  libs(lattice)
  xyplot(geneff ~ rrow|col, dat, layout=c(1,3), type='s',
         main="besag.bayesian - Variety effects", ylim=c(5,15 ))
  xyplot(fert ~ rrow|col, dat, layout=c(1,3), type='s',
         main="besag.bayesian - Fertility", ylim=c(-2,2))
  xyplot(res ~ rrow|col, dat, layout=c(1,3), type='s',
         main="besag.bayesian - Residuals", ylim=c(-4,4))


## End(Not run)