RCB experiment of spring barley in United Kingdom

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.

Used with permission of David Higdon.

Source

Besag, J. E., Green, P. J., Higdon, D. and Mengersen, K. (1995). Bayesian computation and stochastic systems. Statistical Science, 10, 3-66. https://www.jstor.org/stable/2246224

References

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

Examples

# \dontrun{
  
  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)
#> Loading required package: Matrix
#> Offline License checked out Mon Jan 11 17:07:42 2021

  # 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))
#> Offline License checked out Mon Jan 11 17:07:43 2021
#> Model fitted using the gamma parameterization.
#> ASReml 4.1.0 Mon Jan 11 17:07:43 2021
#>           LogLik        Sigma2     DF     wall    cpu
#>  1      -94.4798      0.689559    149 17:07:43    0.0 (2 restrained)
#>  2      -95.8277      0.763261    149 17:07:43    0.0 (2 restrained)
#>  3      -96.2472      0.772025    149 17:07:43    0.0 (2 restrained)
#>  4      -95.2684      0.758604    149 17:07:43    0.0 (1 restrained)
#>  5      -88.7865      0.645542    149 17:07:43    0.0 (2 restrained)
#>  6      -95.1889      0.757622    149 17:07:43    0.0 (1 restrained)
#>  7      -90.7329      0.677440    149 17:07:43    0.0 (2 restrained)
#> Warning: Oscillating parameter(s) reset to average value (iteration 7).
#>  8      -95.6188      0.762867    149 17:07:43    0.0 (1 restrained)
#>  9      -91.3575      0.705741    149 17:07:43    0.0
#> 10      -87.7909      0.644795    149 17:07:43    0.0 (1 restrained)
#> 11      -83.9513      0.601796    149 17:07:43    0.0
#> 12      -81.8901      0.548584    149 17:07:43    0.0 (1 restrained)
#> 13      -82.1499      0.568466    149 17:07:43    0.0
#> Warning: Log-likelihood not converged
#> Warning: Some components changed by more than 1% on the last iteration.
  m1 <- update(m1)
#> Model fitted using the gamma parameterization.
#> ASReml 4.1.0 Mon Jan 11 17:07:43 2021
#>           LogLik        Sigma2     DF     wall    cpu
#>  1      -82.2785      0.567157    149 17:07:43    0.0
#>  2      -81.8595      0.561558    149 17:07:43    0.0
#>  3      -81.3047      0.543601    149 17:07:43    0.0
#>  4      -81.2374      0.541386    149 17:07:43    0.0
#>  5      -81.1807      0.538231    149 17:07:43    0.0
#>  6      -81.1162      0.532665    149 17:07:43    0.0
#>  7      -81.1079      0.528363    149 17:07:43    0.0
#>  8      -81.1147      0.531189    149 17:07:43    0.0
#>  9      -81.1090      0.530548    149 17:07:43    0.0
#> 10      -81.1089      0.530215    149 17:07:43    0.0
#> 11      -81.1068      0.529576    149 17:07:43    0.0
#> 12      -81.1077      0.529867    149 17:07:43    0.0
#> 13      -81.1065      0.529367    149 17:07:43    0.0
#> Warning: Log-likelihood not converged
#> Warning: Some components changed by more than 1% on the last iteration.
  m1 <- update(m1)
#> Model fitted using the gamma parameterization.
#> ASReml 4.1.0 Mon Jan 11 17:07:43 2021
#>           LogLik        Sigma2     DF     wall    cpu
#>  1      -81.1071      0.529705    149 17:07:43    0.0
#>  2      -81.1062      0.529401    149 17:07:43    0.0
  m1 <- update(m1)
#> Model fitted using the gamma parameterization.
#> ASReml 4.1.0 Mon Jan 11 17:07:43 2021
#>           LogLik        Sigma2     DF     wall    cpu
#>  1      -81.1063      0.529138    149 17:07:43    0.0
#>  2      -81.1062      0.529216    149 17:07:43    0.0
    
  # 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))

# }