Question 20. Toenail infection GLMM
2024 examination, Question 4 (20 marks).
In this trial, 294 patients were randomised to terbinafine (250 mg/day) or itraconazole (200 mg/day), and assessed at weeks 0,4,8,12,24,36,48. The binary outcome is moderate or severe onycholysis (1) versus none or mild (0). The command and summary output for the first fit are reproduced below.
fit1 = glmer(outcome ~ treatment*weeks + (1|patientID), family=binomial)
Random effects:
Groups Name Variance Std.Dev.
patientID (Intercept) 21.05 4.588
Number of obs: 1908, groups: patientID, 294
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.54669 0.78251 -3.255 0.00114 **
treatmentTerbinafine -0.25799 0.69136 -0.373 0.70902
weeks -0.10343 0.01214 -8.517 <2e-16 ***
treatmentTerbinafine:weeks -0.04059 0.01840 -2.207 0.02735 *(a) [2 marks] Why is a patient-specific random intercept sensible?
(b) [6 marks] Write the GLMM corresponding to the first command.
(c) [5 marks] Interpret the five reported values, on the scale of the model.
(d) [2 marks] Is there a treatment difference in the first fit, and which treatment looks better?
The second fit allows a patient-specific week slope:
fit2 = glmer(outcome ~ treatment*weeks + (weeks|patientID), family=binomial)
Random effects:
Groups Name Variance Std.Dev. Corr
patientID (Intercept) 7380.733 85.911
weeks 9.124 3.021 -0.30
Number of obs: 1908, groups: patientID, 294
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -12.68683 1.21837 -10.413 <2e-16 ***
treatmentterbinafine -0.06899 1.56848 -0.044 0.965
weeks -2.80367 0.31638 -8.862 <2e-16 ***
treatmentterbinafine:weeks -0.27229 0.27409 -0.993 0.321(e) [2 marks] How does the second model differ from the first?
(f) [3 marks] Does the second fit show a treatment difference? How can the two conclusions be reconciled?
Definitions. Let j index patients, i index visits, t_{ji} be week, and T_j=1 for terbinafine or 0 for itraconazole. The response Y_{ji}=1 indicates moderate or severe onycholysis. The conditional probability is p_{ji}(\alpha_j)=\mathbb{P}(Y_{ji}=1\mid T_j,t_{ji},\alpha_j). The logit is \log\{p/(1-p)\} and its inverse is \operatorname{expit}(\eta)=e^\eta/(1+e^\eta). Fixed effects describe the population predictor; random effects describe patient departures from it.
(a) Each patient is assessed at seven visits, so the responses are repeated measures on the same person. A patient intercept is a draw from a population of patients: one variance, rather than a separate mean for each of the 294 named people. The trial is about the treatments in that population. Two visits on the same patient share the intercept, so they are dependent. Conditional on the intercept they are independent. Assigning the treatment at random does not remove that dependence.
For two distinct visits on the same patient, the law of total covariance gives
\begin{aligned} \operatorname{Cov}(Y_{ji},Y_{jk}\mid X) &=\mathbb{E}\{\operatorname{Cov}(Y_{ji},Y_{jk}\mid\alpha_j,X)\mid X\}\\ &\quad+\operatorname{Cov}\{p_{ji}(\alpha_j),p_{jk}(\alpha_j)\mid X\}. \end{aligned}
The first term is zero under conditional independence. The second is generally nonzero because both probabilities change with the same latent intercept. This is the mechanism by which the model accounts for repeated-measure dependence.
(b) Let T_j=1 for terbinafine and T_j=0 for itraconazole, and let visit i of patient j have response Y_{ji}. Naming binomial and not naming a link selects the logit. The model is
\begin{aligned} Y_{ji}\mid \alpha_j &\sim \operatorname{Bernoulli}(p_{ji}),\\ \operatorname{logit}(p_{ji}) &= \beta_0+\beta_1 T_j+\beta_2\,\mathrm{weeks}_{ji} +\beta_3 T_j\,\mathrm{weeks}_{ji}+\alpha_j,\\ \alpha_j &\stackrel{\mathrm{iid}}{\sim} \operatorname{N}(0,\sigma_A^2), \end{aligned}
with visits conditionally independent given \alpha_j, patients independent of one another, and random intercepts independent of the assigned treatment. Equivalently, p_{ji} is the expit of the linear predictor. The normal distribution is a model for a latent intercept on the log-odds scale. It is not a claim that a normal measurement was observed. The likelihood integrates \alpha_j out, and there is no closed form.
For patient j the likelihood contribution is the one-dimensional integral of the product of the visit-specific Bernoulli masses against the normal random-intercept density. Integration is performed once per patient, since all visits share that intercept. Similarly, the population-averaged probability is \mathbb{E}_{\alpha_j}\{p_{ji}(\alpha_j)\}. Because expit is nonlinear, this generally differs from setting \alpha_j=0 in the predictor.
(c) Every coefficient below holds the patient’s intercept fixed. A unit change is a change in log-odds, and \exp(\beta) multiplies the odds. It is not a change in probability.
The standard deviation 4.588 is the standard deviation of \alpha_j on the log-odds scale. Holding fixed covariates constant, a one-standard-deviation increase in the latent intercept multiplies conditional odds by \exp(4.588)\approx98. It is a scale of between-patient heterogeneity, not a standard deviation of observed binary responses or disease probabilities.
The intercept -2.547 is the log-odds for a typical patient, \alpha_j=0, on itraconazole at week 0. Week 0 is a real visit. The probability is
\frac{\exp(-2.547)}{1+\exp(-2.547)}\approx 0.073.
The terbinafine coefficient -0.258 is the week-0 contrast for a given patient. Terbinafine multiplies the odds of moderate or severe onycholysis by \exp(-0.258)\approx 0.77 relative to itraconazole. The z statistic is -0.37 and p=0.71, so this is not a detected difference.
The weeks coefficient -0.103 is the per-week change for a patient on itraconazole. One extra week multiplies the odds by \exp(-0.103)\approx 0.90, at that patient’s intercept. The z statistic is -8.5.
The interaction -0.0406 is the extra weekly change under terbinafine. Each extra week multiplies the terbinafine-versus-itraconazole odds ratio by \exp(-0.0406)\approx 0.96. The weekly multiplier itself is \exp(-0.144)\approx 0.87 on terbinafine and \exp(-0.103)\approx 0.90 on itraconazole. The z statistic is -2.21 and p=0.027, past the usual 5\% line.
One can obtain all these interpretations by writing the two fitted predictors separately:
\eta_I(t,\alpha)=-2.547-0.103t+\alpha,\qquad \eta_T(t,\alpha)=-2.805-0.144t+\alpha.
Subtracting them at a common week and random-intercept value gives \log\operatorname{OR}_{T:I}(t)=\beta_1+\beta_3t. Differentiating with respect to week gives the two log-odds slopes and their difference \beta_3. This avoids confusing a baseline contrast, a time slope and a difference in slopes.
(d) There is evidence of a difference in the rate, with a steeper decline under terbinafine. The week-0 contrast has p=0.71 and describes baseline differences before treatment, not the difference in improvement. The treatment-by-weeks interaction has p=0.027 and addresses that difference. At week 48, which is a visit in the trial, the subject-specific log-odds contrast is
-0.258+48(-0.0406)=-2.21,
so the odds of moderate or severe onycholysis under terbinafine are \exp(-2.21)\approx 0.11 times the odds under itraconazole, for a given patient.
This is a fitted contrast, not a significance test at week 48. Its standard error would require \operatorname{Var}(\hat\beta_1)+48^2\operatorname{Var}(\hat\beta_3)+96\operatorname{Cov}(\hat\beta_1,\hat\beta_3), and the covariance is not printed. The interaction test itself is H_0:\beta_3=0; equality of the treatment predictors at every week would require the joint null \beta_1=\beta_3=0.
(e) (weeks|patientID) adds a random slope on weeks and a covariance between that slope and the intercept. In the notation of (b),
\operatorname{logit}(p_{ji}) = X_{ji}^{\top}\beta+Z_{ji}^{\top}\alpha_j, \qquad \alpha_j\sim \operatorname{N}(0,\Sigma_A),
where Z_{ji}=(1,\mathrm{weeks}_{ji})^{\top} and \Sigma_A holds two variances and one covariance. The first model is the special case with the slope variance equal to zero. Patients may then differ in baseline and in rate. The fixed effects remain the population constants. The \alpha_j are draws.
Writing \alpha_j=(\alpha_{0j},\alpha_{1j})^\top, patient j’s log-odds slope is \beta_2+\beta_3T_j+\alpha_{1j}. If the two random-effect variances are \tau_0^2,\tau_1^2 and their covariance is \tau_{01}, the covariance between the random parts of that patient’s predictors at weeks s,t is \tau_0^2+(s+t)\tau_{01}+st\tau_1^2. The model therefore changes the time-dependent covariance structure as well as allowing different trajectories.
(f) In the second fit the treatment coefficient has p=0.965 and the treatment-by-weeks interaction has p=0.321. It therefore gives no clear evidence that the two drugs differ in either baseline log-odds or rate of change. Both fitted weekly slopes are negative, but this alone does not establish a causal improvement for either drug. The first fit forces all patients in a treatment group to share a time slope; the second allows patient-specific slopes. Some variation attributed to the fixed interaction in the first fit can instead be represented by those individual trajectories. The interaction standard error rises from 0.018 to 0.274, so the second estimate is much less precise.
The fitted intercept–slope correlation is -0.30. On the model’s log-odds scale, patients with a higher random intercept tend to have a more negative random slope, hence a steeper decline; the association is modest and does not itself explain a treatment contrast. Alan’s central point is that allowing different patient trajectories changes the treatment comparison. The very large fitted random-effect standard deviations also warrant a fit check. Neither fit proves that the drugs are equally effective.
For example, an approximate 95\% Wald interval for the second interaction is -0.27229\pm1.96(0.27409), or about (-0.810,0.265). It includes zero and substantial nonzero differences in slopes. The changed conclusion is therefore about the precision of the treatment comparison under a different dependence model. Adding random slopes does not universally increase every standard error, and the change in a p-value alone does not establish which random-effects structure is adequate.