ASR023. LMM specifying a generalized unstructured correlation matrix - Effect of ozone on lung capacity

The complete script for this example can be downloaded here:

Dataset

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

rat	time	resp
1	After	9.4
1	Before	8.7
2	After	9.8
2	Before	7.9
3	After	9.9
3	Before	8.3
4	After	10.3
4	Before	8.4

Model

In this example we will fit a LMM specifying a generalized unstructured correlation matrix (corgh()), and using initial values for the estimation of its variance components. The specification of the model is: \[ y = \mu + time + e\ \]

where,

\(y\) is the lung capacity in milliliters,

\(\mu\) is the overall mean,

\(time\) is the fixed effect of the measurement time (i.e., Before and After),

\(e\)

\(e \sim \mathcal{N}(0,\mathbf{R})\)

In the above model, the residual matrix \(\mathbf{R = U \otimes I_m}\) represents a direct product between a 2 \(\times\) 2 matrix (\(\mathbf{U}\)) describing the structure of each individual over their two measurements. This implies a variance for each time point and a correlation between residuals. In addition, the second matrix (\(\mathbf{I_m}\)) is describing the independence between the \(m\) individuals.

Before fitting our model, we need to obtain the initial values that we will use in the corgh() structure. For this, we will fit two simpler models to obtain the variance components associated with the time levels of “Before” and “After”. We will be using the ASReml-R option subsetto select each measurement separately. For the correlation between the two measurement, we will use a rough estimate of 0.5. Further more, the data have to be re-ordered to match the order specified by the residual formula, which in this case is by time level. Note that before fitting the models, rat, and time need to be set as factors.

asr023_before <- asreml(
  fixed = resp ~ 1, 
  subset = time == 'Before', 
  data = d019
)

summary(asr023_before)$varcomp

        component std.error  z.ratio bound %ch
units!R 0.2663636 0.1135778 2.345208     P   0

asr023_after <- asreml(
  fixed = resp ~ 1, 
  subset = time == 'After', 
  data = d019
)

summary(asr023_after)$varcomp

        component std.error  z.ratio bound %ch
units!R 0.9753788 0.4159029 2.345208     P   0

So, now using the results from the above models, we can set up the initial values for our next model. This corresponds to three terms. The first is a guess of the correlation between the residual effects, followed by our previous estimates of the residual variance by time level.

initR <- c(0.5, 0.266, 0.975)

Now, let’s take a look at how to write our final model with ASReml-R, but note that we are first sorting the data:

d019 <- d019[order(d019$time),]
asr023 <- asreml(
  fixed = resp ~ time, 
  residual = ~corgh(time, init = initR):id(rat),
  data = d019
)

Exploring output

The variance components estimated from this model are:

summary(asr023)$varcomp

                                      component std.error  z.ratio bound %ch
time:rat!R                           1.00000000        NA       NA     F 0.0
time:rat!time!Before:!time!After.cor 0.07936761 0.2994146 0.265076     U 0.1
time:rat!time_After                  0.97537227 0.4142945 2.354297     P 0.3
time:rat!time_Before                 0.26636364 0.1135774 2.345216     P 0.0

Here, V2 corresponds to the correlation of the residuals between “Before” and “After” with a relatively low value of 0.08. The components V3 and V4 correspond to their variance components.

Let’s take a look at the ANOVA table:

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

            Df denDF   F.inc           Pr
(Intercept)  1    11 4070.00 1.776357e-15
time         1    11   15.09 2.541295e-03

This indicates that there are statistically significant differences between time points.

Finally, we can request the predictions for the factor `time’:

predict(asr023, classify = 'time')$pvals

    time predicted.value std.error    status
1  After        9.658333 0.2850992 Estimable
2 Before        8.450000 0.1489865 Estimable

These results indicate that after 30 days of exposure to ozone, the rats have increased, on average, their lung capacity.