ASR014. LMM for a split-plot design with covariates - Oats varieties

The complete script for this example can be downloaded here:

Dataset

The model that we will fit here is based on the D011 dataset, and the first few rows are presented below:

block	nitrogen	subplot	variety	wplot	yield	column	row	nrate
1	0.6_cwt	1	Marvellous	1	156	1	1	0.6
1	0.4_cwt	2	Marvellous	1	118	2	1	0.4
1	0.2_cwt	3	Marvellous	1	140	1	2	0.2
1	0_cwt	4	Marvellous	1	105	2	2	0.0
1	0_cwt	1	Victory	2	111	1	3	0.0
1	0.2_cwt	2	Victory	2	130	2	3	0.2

Model

In this example we will fit a LMM with the full structure of a split-plot design. The specification of the model is: \[ y = \mu + variety + nitrogen + variety:nitrogen + block + block:wplot +e\ \]

where,

\(y\) is the yield of oats,

\(\mu\) is the overall mean,

\(variety\) is the fixed effect of variety,

\(nitrogen\) is the fixed effect of the level of nitrogen,

\(variety:nitrogen\) is the fixed effect of the interaction between variaty and nitrogen,

\(block\) is the random effect of block, with \(block \sim \mathcal{N}(0,\,\sigma^{2}_{b})\),

\(block:wplot\) is the nested random effect of plot within blocks, with \(block:wplot \sim \mathcal{N}(0,\,\sigma^{2}_{w})\),

\(e\) is the random residual effect, with \(e \sim \mathcal{N}(0,\,\sigma^{2}_{e})\).

Now, let’s take a look at how to write the model with ASReml-R. Note that before fitting the model, variety, nitrogen, block, and wplot need to be set as factors.

asr014 <- asreml(
  fixed = yield ~ variety + nitrogen + variety:nitrogen, 
  random = ~block + block:wplot, 
  residual = ~units,
  data = d011
)

This model can also be written in abbreviated form as:

asr014 <- asreml(
  fixed = yield ~ variety*nitrogen, 
  random = ~block/wplot,
  residual = ~units,
  data = d011
)

Exploring output

Let’s take a look at the ANOVA table:

wald(asr014, denDF = 'numeric')$Wald

                 Df denDF    F.inc           Pr
(Intercept)       1     5 245.1000 1.931823e-05
variety           2    10   1.4850 2.723869e-01
nitrogen          3    45  37.6900 2.457701e-12
variety:nitrogen  6    45   0.3028 9.321988e-01

Here we can see that only the facotr that is significant for yield is nitrogen. Note also the degrees of freedom that describe the hierarchical structure of the experiment as describes a split-plot design.

The variance components estimated from this model are:

summary(asr014)$varcomp

            component std.error  z.ratio bound %ch
block        214.4906 168.66037 1.271731     P 0.1
block:wplot  106.0686  67.88201 1.562543     P 0.0
units!R      177.0946  37.33601 4.743266     P 0.0

Lets obtain the predicted means, based on the above model for the combination of variety:nitrogen:

predict(asr014, classify = 'variety:nitrogen')$pvals

       variety nitrogen predicted.value std.error    status
1  Golden_rain    0_cwt        80.00000  9.106977 Estimable
2  Golden_rain  0.2_cwt        98.50000  9.106977 Estimable
3  Golden_rain  0.4_cwt       114.66667  9.106977 Estimable
4  Golden_rain  0.6_cwt       124.83333  9.106977 Estimable
5   Marvellous    0_cwt        86.66667  9.106977 Estimable
6   Marvellous  0.2_cwt       108.50000  9.106977 Estimable
7   Marvellous  0.4_cwt       117.16667  9.106977 Estimable
8   Marvellous  0.6_cwt       126.83333  9.106977 Estimable
9      Victory    0_cwt        71.50000  9.106977 Estimable
10     Victory  0.2_cwt        89.66667  9.106977 Estimable
11     Victory  0.4_cwt       110.83333  9.106977 Estimable
12     Victory  0.6_cwt       118.50000  9.106977 Estimable

But more relevant are the predictions for the only significant factor:

preds <- predict(asr014, classify = 'nitrogen')$pvals

  nitrogen predicted.value std.error    status
1    0_cwt        79.38889   7.17471 Estimable
2  0.2_cwt        98.88889   7.17471 Estimable
3  0.4_cwt       114.22222   7.17471 Estimable
4  0.6_cwt       123.38889   7.17471 Estimable