snedecor.asparagus {agridat} R Documentation

## Asparagus yields for different cutting treatments

### Description

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

```## Not run:

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)
# 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)

## # 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

## End(Not run)
```

[Package agridat version 1.18 Index]