heart | time | A | B | ATP |
---|---|---|---|---|

1 | 0 | 1 | 1 | 85.51 |

1 | 1 | 1 | 1 | 74.56 |

1 | 2 | 1 | 1 | 84.25 |

1 | 3 | 1 | 1 | 84.05 |

1 | 4 | 1 | 1 | 78.41 |

1 | 5 | 1 | 1 | 79.93 |

1 | 6 | 1 | 1 | 71.60 |

1 | 8 | 1 | 1 | 69.90 |

# ASR022. Mixed linear model for repeated measures using an antedependence residual structure - Dog hearts’ enzyme content

The complete script for this example can be downloaded here:

### Dataset

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

### Model

In this example we will fit a LMM for repeated measures using an antedependence residual structure of order \(1\). The specification of the model is:

\[ y = \mu + A + B + A:B + A:time + B:time + heart + e\ \]

where,

\(y\) is the percentage of enzyme content,

\(\mu\) is the overall mean,

\(A\) is the fixed effect of treatment A,

\(B\) is the fixed effect of treatment B,

\(A:B\) is the fixed effect of the interaction between treatment A and B,

\(A:time\) is the fixed effect of the interaction between treatment A and time,

\(B:time\) is the fixed effect of the interaction between treatment B and time,

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

\(e\) is the random residual effect, with \(e \sim \mathcal{N}(0,\mathbf{R})\).

In the above model, \(\mathbf{R}\) is the residual variance-covariance matrix specified as an antedependence model of order \(k = 1\). The structure can be modeled as \(\mathbf{R^{-1} = U D U'}\) where \(\mathbf{U}\) is a unit upper triangular matrix with elements {\(u_{ij}\)}

where,

{\(u_{ii}\)} \(= 1\)

{\(u_{ij}\)} \(= 0, i \gt j\)

{\(u_{ij}\)} \(= 0, j - i \gt k\)

{\(u_{ij}\)} \(=\) {\(u_{ij}\)} \(, 1 \leq j - i \leq k\)

and which can be represented such as: \[ \left[\begin{array}{cc} 1 & u_{ij} & 0 & & 0 \\ 0 & 1 & u_{ij} & & \\ & & \ddots &\ddots & 0 \\ & & & & u_{i,j}\\ 0 & & & 0 & 1 \\ \end{array}\right] \]

\(\mathbf{D}\) is a diagonal matrix such as \(\mathbf{D} = diag(d_1 , ..., d_t)\), with measurement times \(1\) to \(t\), and which can be represented such as:

\[ \left[\begin{array}{cc} d_{1} &0 & & & 0 \\ 0 & d_{2} & & & \\ & & \ddots & & \\ & & & & 0\\ 0 & & & 0 & d_{t} \\ \end{array}\right] \]

Now, let’s take a look at how to write the model with ASReml-R. Note that before fitting the model, `A`

, `B`

, `time`

and `heart`

need to be set as factors.

```
<- asreml(
asr022 fixed = ATP ~ A + B + A:B + A:time + B:time,
random = ~heart,
residual = ~heart:ante(time,1),
data = d018
)<- update(asr022) asr022
```

The command `update()`

was used to request a few extra iterations from the previous fit.

### Exploring output

Let’s first take a look at the estimated variance components from this model.

```
<- summary(asr022)$varcomp
varcomp head(varcomp)
```

```
component std.error z.ratio bound %ch
heart 10.59613794 4.49692944 2.35630514 P 0.0
heart:time!R 1.00000000 NA NA F 0.0
heart:time!time_0:0 0.09596961 0.04325412 2.21873918 U 0.0
heart:time!time_1:0 0.22850098 0.42496932 0.53768818 U 0.3
heart:time!time_1:1 0.06224730 0.02951166 2.10924451 U 0.0
heart:time!time_2:1 -0.02068389 0.35696226 -0.05794418 U 0.7
```

In the first few estimates presented here, V1 is the variance of `heart`

, V3 is the variance of the first level of `time`

in the matrix \(\mathbf{D}\) (i.e., \(d_1\)), and V4 is the covariance between the first and second level of `time`

in the matrix \(\mathbf{U}\) (i.e., \(u_{1,2}\))).

To obtain the \(\mathbf{R}\) matrix of variance-covariance, we can to solve the equation \(\mathbf{R^{-1} = U D U'}\), and then obtain the inverse of \(\mathbf{R^{-1}}\). In this example, the \(\mathbf{R}\) matrix can be visualized such as:

\[ \left[\begin{array}{cc} 10.420 & -2.381 & -0.049 & -0.003 & 0.002 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000\\ -2.381 & 16.609 & 0.344 & 0.023 & -0.012 & 0.003 & 0.002 & 0.001 & 0.000 & 0.000\\ -0.049 & 0.344 & 27.751 & 1.861 & -1.000 & 0.248 & 0.135 & 0.081 & 0.035 & 0.021\\ -0.003 & 0.023 & 1.861& 31.280& -16.818 & 4.175 & 2.265 & 1.365 & 0.583 & 0.355\\ 0.002 &-0.012& -1.000 &-16.818 & 74.010 &-18.372 & -9.968 &-6.006 & -2.566 & -1.564\\ 0.000 & 0.003 &0.248 & 4.175 &-18.372 & 38.848 & 21.078 & 12.700 & 5.426 & 3.308\\ 0.000& 0.002 & 0.135 & 2.265 &-9.968 & 21.078& 108.188 & 65.190 & 27.853 & 16.978\\ 0.000 & 0.001 & 0.081 & 1.365 & -6.006 & 12.700& 65.190& 120.082 & 51.306 & 31.274\\ 0.000 & 0.000 & 0.035 & 0.583 & -2.566 & 5.426 & 27.853& 51.306 &115.451 & 70.375\\ 0.000 & 0.000 & 0.021 & 0.355 & -1.564 & 3.308 & 16.978 & 31.274 & 70.375& 145.111\\ \end{array}\right] \] As expected, the covariances between measurements are decreasing as the measurements were farther apart in time.

Let’s now take a look at the ANOVA table:

`wald(asr022, denDF = 'numeric', ssType = 'incremental')$Wald`

```
Df denDF F.inc Pr
(Intercept) 1 18.6 1.021e+04 0.000000e+00
A 1 18.6 3.879e-01 5.409356e-01
B 1 18.6 2.770e-01 6.048542e-01
A:B 1 18.6 5.738e+00 2.727105e-02
A:time 18 51.2 1.485e+01 2.020606e-14
B:time 9 36.3 2.241e+00 4.146219e-02
```

This is an `incremental`

F-test where all interaction terms are significant.

Finally, of interest will be to have the means for each combination of the factors `A`

and `B`

on each of the time points. This is obtained using:

```
<- predict(asr022, classify = 'A:B:time')$pvals
pp head(pp)
tail(pp)
```

```
A B time predicted.value std.error status
1 1 1 0 77.27397 1.790247 Estimable
2 1 1 1 73.75802 1.997676 Estimable
3 1 1 2 79.00357 2.324947 Estimable
4 1 1 3 74.34428 2.419458 Estimable
5 1 1 4 75.50952 3.358588 Estimable
6 1 1 5 73.32333 2.610551 Estimable
```

```
A B time predicted.value std.error status
35 2 2 4 68.82119 3.358588 Estimable
36 2 2 5 62.54000 2.610551 Estimable
37 2 2 6 58.51675 3.952388 Estimable
38 2 2 8 55.43897 4.138841 Estimable
39 2 2 10 53.20309 4.067176 Estimable
40 2 2 12 50.04651 4.506482 Estimable
```