ASR010. Random coefficient LMM for longitudinal data (random intercept and slope) - Dog scans

The complete script for this example can be downloaded here:

Dataset

In this example we will use the D006 dataset, and the first few rows are presented below:

dog side day pixel
1 R 0 1045.8
1 R 1 1044.5
1 R 2 1042.9
1 R 4 1050.4
1 R 6 1045.2
1 R 10 1038.9


Model

In this example, we will fit a random coefficient LMM for longitudinal data. The model specification is as follows;

\[ y = \mu + side + day + daysq + dog + dog:day + e\ \] where,

    \(y\) is the pixel intensity,

    \(\mu\) is the overall mean,

    \(side\) is the fixed covariate of side (left and right lymph nodes),

    \(day\) is the fixed covariate of day,

    \(daysq\) is the fixed covariate of the square of day,

    \(dog\) is the random factor of dog, with \(dog \sim \mathcal{N}(0,\,\sigma^{2}_{di})\),

    \(dog:day\) is the random slope of day for each dog, with \(dog:day \sim \mathcal{N}(0,\,\sigma^{2}_{ds})\),

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


Now, let’s take a look at how to write the model with ASReml-R. Here we will use the temporal covariate dayin its linear and quadratic form because our response variable increases and then decreases over time. Note that before fitting the model, dog and side need to be set as factors.

d006$daysq <- d006$day^2
asr010 <- asreml(
  fixed = pixel ~ side + day + daysq,
  random = ~str(~dog + dog:day, ~corgh(2):id(dog)),
  data = d006
)

The str() function is used to join two terms, and therefore it allows to establish a covariance between model terms in this case. This model will fit a regression for each dog and estimate an unique random intercept and a random slope deviations from the overall regression that is defined by the fixed effects. The str() function requires a model formula (e.g., dog + dog:day) and a variance formula (e.g., corgh(2):id(dog)). Notice that the corgh() function takes \(2\) as the argument value referring to a \(2 \times 2\) covariance matrix for the variance of the intercepts and slopes. And in id(dog), the term dog is used to count the number of levels of this factor, but its number can also be provided instead.

Exploring output


The statistical significance of fixed effects can be tested as:

wald(asr010, denDF = 'numeric')$Wald
            Df denDF     F.inc           Pr
(Intercept)  1   9.0 1.689e+04 4.440892e-16
side         1  81.7 3.897e+00 5.173904e-02
day          1   8.3 6.754e-02 8.012880e-01
daysq        1  88.5 4.839e+01 5.783225e-10


The fixed effects (BLUEs) for the intercept, the effects of side, and the slopes of day and daysq are: 

summary(asr010, coef = TRUE)$coef.fixed
                solution   std error    z.ratio
(Intercept) 1075.6488400 10.51615697 102.285354
side_L         0.0000000          NA         NA
side_R        -5.4019608  2.73650557  -1.974036
day            6.0916044  1.13395062   5.372019
daysq         -0.3573287  0.05137289  -6.955588


The variance components estimated from this model are:

summary(asr010)$varcomp
                                         component   std.error   z.ratio bound %ch
dog+dog:day!corgh(2)!2:!corgh(2)!1.cor  -0.4905052   0.3208553 -1.528743     U 0.0
dog+dog:day!corgh(2)_1                 896.4680272 473.2733209  1.894187     P 0.2
dog+dog:day!corgh(2)_2                   3.0978743   2.0526227  1.509227     P 0.1
units!R                                190.9557993  29.8846212  6.389768     P 0.0

Notice here that the correlation between slope and intercept was estimated to be -0.49, followed by their variances.


A portion of the random effects (BLUPs) are presented below:

  • An intercept deviation for each dog (dog):
BLUP <- summary(asr010, coef = TRUE)$coef.random
head(BLUP)
       solution std.error    z.ratio
dog_1 -33.47888  11.10226 -3.0155007
dog_2 -30.04312  11.10226 -2.7060360
dog_3 -26.91620  11.35313 -2.3708180
dog_4 -11.70962  11.35313 -1.0314001
dog_5  29.74348  13.12739  2.2657569
dog_6  10.12405  13.12739  0.7712161
  • An slope deviation for each dog (dog:day):
tail(BLUP)
             solution std.error    z.ratio
dog_5:day  -0.5330376  1.124856 -0.4738720
dog_6:day  -0.9128701  1.124856 -0.8115442
dog_7:day  -1.8781271  1.125645 -1.6684893
dog_8:day   0.6877326  1.125645  0.6109674
dog_9:day  -0.3121488  1.562941 -0.1997188
dog_10:day -0.1718214  1.562683 -0.1099529