| 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.
asr022 <- asreml(
fixed = ATP ~ A + B + A:B + A:time + B:time,
random = ~heart,
residual = ~heart:ante(time,1),
data = d018
)
asr022 <- update(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.
varcomp <- summary(asr022)$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.7In 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-02This 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:
pp <- predict(asr022, classify = 'A:B:time')$pvals
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