6  Lab 2

6.1 Exercise

The data in hip.txt are from Crowder & Hand (1990), and we can put them in an R data.frame by

hip_url <- url("https://ikosmidis.com/files/APTS-SM-Notes/resources/hip.txt")
hip <- read.table(hip_url,
                  col.names = c("y", "age", "sex", "subj", "time"),
                  colClasses = c(rep("numeric", 4), "factor"))
str(hip)
'data.frame':   88 obs. of  5 variables:
 $ y   : num  47.1 31.1 32.8 44.1 31.5 ...
 $ age : num  66 66 66 70 70 70 44 44 44 70 ...
 $ sex : num  0 0 0 0 0 0 0 0 0 0 ...
 $ subj: num  1 1 1 2 2 2 3 3 3 4 ...
 $ time: Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...

The variable y contains measurements of haematocrit on 30 patients (subj) on up to three occasions (time), one before a hip-replacement operation (coded as 1), and two afterwards (coded as 2 and 3). time is treated here as a categorical variable. The age and sex (0 for male, 1 for female) of the patients are also recorded.

The primary interest in this study is in possible differences in the evolution of haematocrit between males and females and whether there is an age effect.

  1. For each value of sex, plot the time profiles of the response variable for each subject.

    Do you think you think it is necessary to include a random intercept for the subject? What about a random slope for time?

  2. We will analyse these data using linear mixed models (LMMs) of the form \[ \begin{aligned} Y_{ij} \mid x_{ij}, z_{ij} & \stackrel{\text{ind}}{\sim}\mathop{\mathrm{N}}(\mu_{ij}, \sigma^2) \,,\\ \mu_{ij} & = x_{ij}^\top \beta + z_{ij}^\top b_i \,, \\ b_i & \stackrel{\text{ind}}{\sim}\mathop{\mathrm{N}}(0, \Sigma_b^*) \,, \end{aligned} \tag{6.1}\] where \(Y_{ij}\) is the random variable corresponding to the haematocrit measurement for subject \(i\) at time \(j\), and \(x_{ij}\) and \(z_{ij}\) are fixed-effects and mixed-effects covariates, respectively.

    Consider building a set of candidate LMMs for explaining haematocrit, using the fixed-effects of age, sex and time (and possibly interactions of sex to age and time) and random-effects of subj and time. If you want to have models with interaction effects in your candidate set, make sure that you respect the “marginality constraints”, that is the model should have main effects for all terms from which interactions are formed.

    What is your chosen model?

  3. For the model you chose, plot the predicted haematocrit for each patient against time.

    Use your chosen model to predict the full haematocrit profiles of any patients that do not have haematocrit measurements at all three time points.

Hints for plotting

For plotting, you may adapt the code used for producing Figure 2.2 for the rat growth data, or you may use the ggplot2 R package.

Hints for modelling

LMMs for clustered data can be fitted in R using the lmer() function from the lme4 package. For example,

library(lme4)
hip_lmm1 <- lmer(y ~ age + sex + time + (1 | subj), data = hip)

fits the model with \(x_{ij} = (1, \text{age}_i, \text{sex}_i, I(\text{time}_{ij} = 2), I(\text{time}_{ij} = 3))^\top\), and \(z_{ij} = 1\).

The default estimation method in lmer() is REML. If you want to obtain maximum likelihood estimates (for example, for use in model comparison), they can be obtained using the additional argument REML = FALSE.

You might find useful some of the following methods for the object that lmer() returns: summary, fitted, residuals, fixef (fixed effects estimates), ranef (random effects estimates), VarCorr (variance estimates) coef (coefficient estimates at cluster level, incorporating fixed and random effects), AIC, BIC and predict.

For more information, see ?lmer.

Quantities in the general LMM definiton in (2.5)

In terms of the general definition of LMMs in (2.5), in order to get a better understanding of what the formula interface of lmer() does, for hip_lmm1, \(Y\) is

mf1 <- model.frame(hip_lmm1)
model.response(mf1)
    1     2     3     4     5     6     7     8     9    10    11    12    13 
47.10 31.05 32.80 44.10 31.50 37.00 39.70 33.70 24.50 43.30 18.35 36.60 37.40 
   14    15    16    17    18    19    20    21    22    23    24    25    26 
32.25 29.05 45.70 35.50 39.80 44.90 34.10 32.05 42.90 32.05 46.05 28.80 37.80 
   27    28    29    30    31    32    33    34    35    36    37    38    39 
42.10 34.40 36.05 38.25 29.40 30.50 43.00 33.70 36.65 37.80 26.60 30.60 37.25 
   40    41    42    43    44    45    46    47    48    49    50    51    52 
26.50 38.45 27.95 33.95 27.00 32.50 31.95 38.35 32.30 37.90 38.80 32.55 26.85 
   53    54    55    56    57    58    59    60    61    62    63    64    65 
44.65 32.25 34.20 38.00 27.10 37.85 34.00 23.20 25.95 44.80 37.20 29.70 45.95 
   66    67    68    69    70    71    72    73    74    75    76    77    78 
29.10 26.70 41.85 31.95 37.60 38.00 31.65 35.70 42.20 34.00 33.25 39.70 33.45 
   79    80    81    82    83    84    85    86    87    88 
32.65 37.50 28.20 30.30 34.55 30.95 28.75 35.50 24.70 29.75

\(X\) is

model.matrix(hip_lmm1, type = "fixed")
   (Intercept) age sex time2 time3
1            1  66   0     0     0
2            1  66   0     1     0
3            1  66   0     0     1
4            1  70   0     0     0
5            1  70   0     1     0
6            1  70   0     0     1
7            1  44   0     0     0
8            1  44   0     1     0
9            1  44   0     0     1
10           1  70   0     0     0
11           1  70   0     1     0
12           1  70   0     0     1
13           1  74   0     0     0
14           1  74   0     1     0
15           1  74   0     0     1
16           1  65   0     0     0
17           1  65   0     1     0
18           1  65   0     0     1
19           1  54   0     0     0
20           1  54   0     1     0
21           1  54   0     0     1
22           1  63   0     0     0
23           1  63   0     1     0
24           1  71   0     0     0
25           1  71   0     1     0
26           1  71   0     0     1
27           1  68   0     0     0
28           1  68   0     1     0
29           1  68   0     0     1
30           1  69   0     0     0
31           1  69   0     1     0
32           1  69   0     0     1
33           1  64   0     0     0
34           1  64   0     1     0
35           1  64   0     0     1
36           1  70   0     0     0
37           1  70   0     1     0
38           1  70   0     0     1
39           1  60   1     0     0
40           1  60   1     1     0
41           1  60   1     0     1
42           1  52   1     1     0
43           1  52   1     0     1
44           1  52   1     0     0
45           1  52   1     1     0
46           1  52   1     0     1
47           1  75   1     0     0
48           1  75   1     1     0
49           1  75   1     0     1
50           1  72   1     0     0
51           1  72   1     1     0
52           1  72   1     0     1
53           1  54   1     0     0
54           1  54   1     1     0
55           1  54   1     0     1
56           1  71   1     0     0
57           1  71   1     1     0
58           1  71   1     0     1
59           1  58   1     0     0
60           1  58   1     1     0
61           1  58   1     0     1
62           1  77   1     0     0
63           1  77   1     1     0
64           1  77   1     0     1
65           1  66   1     0     0
66           1  66   1     1     0
67           1  66   1     0     1
68           1  53   1     0     0
69           1  53   1     1     0
70           1  53   1     0     1
71           1  74   1     0     0
72           1  74   1     1     0
73           1  74   1     0     1
74           1  78   1     0     0
75           1  78   1     1     0
76           1  78   1     0     1
77           1  74   1     0     0
78           1  74   1     1     0
79           1  74   1     0     1
80           1  79   1     0     0
81           1  79   1     1     0
82           1  79   1     0     1
83           1  71   1     0     0
84           1  71   1     1     0
85           1  71   1     0     1
86           1  68   1     0     0
87           1  68   1     1     0
88           1  68   1     0     1
attr(,"assign")
[1] 0 1 2 3 3
attr(,"contrasts")
attr(,"contrasts")$time
[1] "contr.treatment"

attr(,"msgScaleX")
character(0)

and \(Z\) is

model.matrix(hip_lmm1, type = "random")
88 x 30 sparse Matrix of class "dgCMatrix"
  [[ suppressing 30 column names '1', '2', '3' ... ]]
                                                              
1  1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2  1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3  1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4  . 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5  . 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6  . 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7  . . 1 . . . . . . . . . . . . . . . . . . . . . . . . . . .
8  . . 1 . . . . . . . . . . . . . . . . . . . . . . . . . . .
9  . . 1 . . . . . . . . . . . . . . . . . . . . . . . . . . .
10 . . . 1 . . . . . . . . . . . . . . . . . . . . . . . . . .
11 . . . 1 . . . . . . . . . . . . . . . . . . . . . . . . . .
12 . . . 1 . . . . . . . . . . . . . . . . . . . . . . . . . .
13 . . . . 1 . . . . . . . . . . . . . . . . . . . . . . . . .
14 . . . . 1 . . . . . . . . . . . . . . . . . . . . . . . . .
15 . . . . 1 . . . . . . . . . . . . . . . . . . . . . . . . .
16 . . . . . 1 . . . . . . . . . . . . . . . . . . . . . . . .
17 . . . . . 1 . . . . . . . . . . . . . . . . . . . . . . . .
18 . . . . . 1 . . . . . . . . . . . . . . . . . . . . . . . .
19 . . . . . . 1 . . . . . . . . . . . . . . . . . . . . . . .
20 . . . . . . 1 . . . . . . . . . . . . . . . . . . . . . . .
21 . . . . . . 1 . . . . . . . . . . . . . . . . . . . . . . .
22 . . . . . . . 1 . . . . . . . . . . . . . . . . . . . . . .
23 . . . . . . . 1 . . . . . . . . . . . . . . . . . . . . . .
24 . . . . . . . . 1 . . . . . . . . . . . . . . . . . . . . .
25 . . . . . . . . 1 . . . . . . . . . . . . . . . . . . . . .
26 . . . . . . . . 1 . . . . . . . . . . . . . . . . . . . . .
27 . . . . . . . . . 1 . . . . . . . . . . . . . . . . . . . .
28 . . . . . . . . . 1 . . . . . . . . . . . . . . . . . . . .
29 . . . . . . . . . 1 . . . . . . . . . . . . . . . . . . . .
30 . . . . . . . . . . 1 . . . . . . . . . . . . . . . . . . .
31 . . . . . . . . . . 1 . . . . . . . . . . . . . . . . . . .
32 . . . . . . . . . . 1 . . . . . . . . . . . . . . . . . . .
33 . . . . . . . . . . . 1 . . . . . . . . . . . . . . . . . .
34 . . . . . . . . . . . 1 . . . . . . . . . . . . . . . . . .
35 . . . . . . . . . . . 1 . . . . . . . . . . . . . . . . . .
36 . . . . . . . . . . . . 1 . . . . . . . . . . . . . . . . .
37 . . . . . . . . . . . . 1 . . . . . . . . . . . . . . . . .
38 . . . . . . . . . . . . 1 . . . . . . . . . . . . . . . . .
39 . . . . . . . . . . . . . 1 . . . . . . . . . . . . . . . .
40 . . . . . . . . . . . . . 1 . . . . . . . . . . . . . . . .
41 . . . . . . . . . . . . . 1 . . . . . . . . . . . . . . . .
42 . . . . . . . . . . . . . . 1 . . . . . . . . . . . . . . .
43 . . . . . . . . . . . . . . 1 . . . . . . . . . . . . . . .
44 . . . . . . . . . . . . . . . 1 . . . . . . . . . . . . . .
45 . . . . . . . . . . . . . . . 1 . . . . . . . . . . . . . .
46 . . . . . . . . . . . . . . . 1 . . . . . . . . . . . . . .
47 . . . . . . . . . . . . . . . . 1 . . . . . . . . . . . . .
48 . . . . . . . . . . . . . . . . 1 . . . . . . . . . . . . .
49 . . . . . . . . . . . . . . . . 1 . . . . . . . . . . . . .
50 . . . . . . . . . . . . . . . . . 1 . . . . . . . . . . . .
51 . . . . . . . . . . . . . . . . . 1 . . . . . . . . . . . .
52 . . . . . . . . . . . . . . . . . 1 . . . . . . . . . . . .
53 . . . . . . . . . . . . . . . . . . 1 . . . . . . . . . . .
54 . . . . . . . . . . . . . . . . . . 1 . . . . . . . . . . .
55 . . . . . . . . . . . . . . . . . . 1 . . . . . . . . . . .
56 . . . . . . . . . . . . . . . . . . . 1 . . . . . . . . . .
57 . . . . . . . . . . . . . . . . . . . 1 . . . . . . . . . .
58 . . . . . . . . . . . . . . . . . . . 1 . . . . . . . . . .
59 . . . . . . . . . . . . . . . . . . . . 1 . . . . . . . . .
60 . . . . . . . . . . . . . . . . . . . . 1 . . . . . . . . .
61 . . . . . . . . . . . . . . . . . . . . 1 . . . . . . . . .
62 . . . . . . . . . . . . . . . . . . . . . 1 . . . . . . . .
63 . . . . . . . . . . . . . . . . . . . . . 1 . . . . . . . .
64 . . . . . . . . . . . . . . . . . . . . . 1 . . . . . . . .
65 . . . . . . . . . . . . . . . . . . . . . . 1 . . . . . . .
66 . . . . . . . . . . . . . . . . . . . . . . 1 . . . . . . .
67 . . . . . . . . . . . . . . . . . . . . . . 1 . . . . . . .
68 . . . . . . . . . . . . . . . . . . . . . . . 1 . . . . . .
69 . . . . . . . . . . . . . . . . . . . . . . . 1 . . . . . .
70 . . . . . . . . . . . . . . . . . . . . . . . 1 . . . . . .
71 . . . . . . . . . . . . . . . . . . . . . . . . 1 . . . . .
72 . . . . . . . . . . . . . . . . . . . . . . . . 1 . . . . .
73 . . . . . . . . . . . . . . . . . . . . . . . . 1 . . . . .
74 . . . . . . . . . . . . . . . . . . . . . . . . . 1 . . . .
75 . . . . . . . . . . . . . . . . . . . . . . . . . 1 . . . .
76 . . . . . . . . . . . . . . . . . . . . . . . . . 1 . . . .
77 . . . . . . . . . . . . . . . . . . . . . . . . . . 1 . . .
78 . . . . . . . . . . . . . . . . . . . . . . . . . . 1 . . .
79 . . . . . . . . . . . . . . . . . . . . . . . . . . 1 . . .
80 . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 . .
81 . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 . .
82 . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 . .
83 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 .
84 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 .
85 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 .
86 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
87 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
88 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1