Asparagus yields for different cutting treatments

Asparagus yields for different cutting treatments, in 4 years.

Format

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

block

block factor, 4 levels

year

year, numeric

trt

treatment factor of final cutting date

yield

yield, ounces

Details

Planted in 1927. Cutting began in 1929. Yield is the weight of asparagus cuttings up to Jun 1 in each plot. Some plots received continued cuttings until Jun 15, Jul 1, and Jul 15.

In the past, repeated-measurement experiments like this were sometimes analyzed as if they were a split-plot experiment. This violates some indpendence assumptions.

Source

Snedecor and Cochran, 1989. Statistical Methods.

References

Mick O'Neill, 2010. A Guide To Linear Mixed Models In An Experimental Design Context. Statistical Advisory & Training Service Pty Ltd.

Examples

# \dontrun{
  
  library(agridat)
  data(snedecor.asparagus)
  dat <- snedecor.asparagus

  dat <- transform(dat, year=factor(year))
  dat$trt <- factor(dat$trt,
                    levels=c("Jun-01", "Jun-15", "Jul-01", "Jul-15"))

  # Continued cutting reduces plant vigor and yield
  libs(lattice)
  dotplot(yield ~ trt|year, data=dat,
          xlab="Cutting treatment", main="snedecor.asparagus")

  # Split-plot
  if(0){
    libs(lme4)
    m1 <- lmer(yield ~ trt + year + trt:year +
                 (1|block) + (1|block:trt), data=dat)
  }

  # ----------

  libs(asreml,lucid) # asreml4

  # Split-plot with asreml
  m2 <- asreml(yield ~ trt + year + trt:year, data=dat,
               random = ~ block + block:trt)
#> Model fitted using the gamma parameterization.
#> ASReml 4.1.0 Mon Jan 11 17:08:50 2021
#>           LogLik        Sigma2     DF     wall    cpu
#>  1      -245.589       735.637     60 17:08:51    0.0
#>  2      -242.013       605.606     60 17:08:51    0.0
#>  3      -239.057       500.049     60 17:08:51    0.0
#>  4      -237.707       441.384     60 17:08:51    0.0
#>  5      -237.333       411.824     60 17:08:51    0.0
#>  6      -237.307       405.253     60 17:08:51    0.0
#>  7      -237.307       404.663     60 17:08:51    0.0
  # vc(m2)
  ##    effect component std.error z.ratio bound 
  ##     block     354.3     405      0.87     P 0.1
  ## block:trt     462.8     256.9    1.8      P 0  
  ##   units!R     404.7      82.6    4.9      P 0  

  ## # Antedependence with asreml.  See O'Neill (2010).
  dat <- dat[order(dat$block, dat$trt), ]
  m3 <- asreml(yield ~ year * trt, data=dat,
               random = ~ block,
               residual = ~ block:trt:ante(year,1),
               max=50)
#> Model fitted using the sigma parameterization.
#> ASReml 4.1.0 Mon Jan 11 17:08:51 2021
#>           LogLik        Sigma2     DF     wall    cpu
#>  1      -239.093           1.0     60 17:08:51    0.0
#>  2      -237.801           1.0     60 17:08:51    0.0
#>  3      -233.259           1.0     60 17:08:51    0.0
#>  4      -223.667           1.0     60 17:08:51    0.0
#>  5      -220.355           1.0     60 17:08:51    0.0
#>  6      -220.139           1.0     60 17:08:51    0.0
#>  7      -220.091           1.0     60 17:08:51    0.0
#>  8      -220.079           1.0     60 17:08:51    0.0
#>  9      -220.076           1.0     60 17:08:51    0.0
#> Warning: Some components changed by more than 1% on the last iteration.

  ## # Extract the covariance matrix for years and convert to correlation
  ## covmat <- diag(4)
  ## covmat[upper.tri(covmat,diag=TRUE)] <- m3$R.param$`block:trt:year`$year$initial
  ## covmat[lower.tri(covmat)] <- t(covmat)[lower.tri(covmat)]
  ## round(cov2cor(covmat),2) # correlation among the 4 years
  ## #      [,1] [,2] [,3] [,4]
  ## # [1,] 1.00 0.45 0.39 0.31
  ## # [2,] 0.45 1.00 0.86 0.69
  ## # [3,] 0.39 0.86 1.00 0.80
  ## # [4,] 0.31 0.69 0.80 1.00

  ## # We can also build the covariance Sigma by hand from the estimated
  ## # variance components via: Sigma^-1 = U D^-1 U'
  ## vv <- vc(m3)
  ## print(vv)
  ## ##            effect component std.error z.ratio constr
  ## ##   block!block.var  86.56    156.9        0.55    pos
  ## ##        R!variance   1              NA      NA    fix
  ## ##  R!year.1930:1930   0.00233   0.00106    2.2   uncon
  ## ##  R!year.1931:1930  -0.7169    0.4528    -1.6   uncon
  ## ##  R!year.1931:1931   0.00116   0.00048    2.4   uncon
  ## ##  R!year.1932:1931  -1.139     0.1962    -5.8   uncon
  ## ##  R!year.1932:1932   0.00208   0.00085    2.4   uncon
  ## ##  R!year.1933:1932  -0.6782    0.1555    -4.4   uncon
  ## ##  R!year.1933:1933   0.00201   0.00083    2.4   uncon
  
  ## U <- diag(4)
  ## U[1,2] <- vv[4,2] ; U[2,3] <- vv[6,2] ; U[3,4] <- vv[8,2]
  ## Dinv <- diag(c(vv[3,2], vv[5,2], vv[7,2], vv[9,2]))
  ## # solve(U 
  ## solve(crossprod(t(U), tcrossprod(Dinv, U)) )
  ## ##          [,1]      [,2]      [,3]      [,4]
  ## ## [1,] 428.4310  307.1478  349.8152  237.2453
  ## ## [2,] 307.1478 1083.9717 1234.5516  837.2751
  ## ## [3,] 349.8152 1234.5516 1886.5150 1279.4378
  ## ## [4,] 237.2453  837.2751 1279.4378 1364.8446

# }