Notes 6 — Linear Mixed Models & GEE
Module 2 · Linear mixed models and generalized estimating equations
Lectures 8.1–8.4. Reading: Lecture Notes Part 2, §2.2–2.3. Code: R script — Lecture Video 35 (mixed models), R script — Lecture Video 37 (GEE).
Picks up from Notes 5: $N$ people, each with a vector $\mathbf{Y}i = (Y{i1}, \dots, Y_{in})^T$ of repeated measures and a within-person covariance $\Sigma_i$ that cannot be ignored. Two ways to handle it — model it through random effects (this page’s first half) or model only the mean and treat the correlation as a nuisance (the second half) — and the three research questions (time effect, group effect, time × group) answered with each.
Lecture 8.1 — Linear mixed models
Three assumptions, and three ways to handle $\Sigma_i$
To do inference we assume:
- A model for the mean, $E(\mathbf{Y}_i) = X_i\beta$.
- Multivariate normality, $\mathbf{Y}_i \sim N(X_i\beta,\ \Sigma_i(\theta))$, with $\theta$ a vector of covariance parameters.
- Some structure for $\Sigma_i(\theta)$. Three families:
| Structure | What it assumes | Parameters |
|---|---|---|
| Unstructured | nothing — every variance and covariance is free | $n(n+1)/2$ |
| Compound symmetry | constant variance $\sigma^2$; every pair of occasions has the same correlation $\rho$ | 2 |
| Autoregressive (AR-1) | constant variance; correlation $\rho^{\lvert j-k \rvert}$ decays with time separation | 2 |
With these, $\hat\beta$ is the maximum-likelihood (generalized least squares) estimate. In practice REML — restricted maximum likelihood — is used instead of plain ML, because ML’s variance-component estimates are biased downward, especially with few observations per person. R does this automatically.
Three approaches to the covariance:
- assume unstructured — many parameters;
- assume a pattern like CS or AR-1;
- use a random-effects covariance structure.
Approaches 1 and 2 are covariance pattern models: they soak up every source of within-subject correlation into one matrix without separating between-subject from within-subject variability, and most of them need equally spaced, everyone-measured- at-the-same-time data. Approach 3 — the linear mixed model — is the most common, and the course’s focus.
The mixed-model idea
Each individual has their own mean trajectory. The mean response is a combination of fixed effects $\beta$ shared by everyone (population characteristics) and random effects unique to each person — hence mixed. It is a hierarchical model.
Random intercepts model
\[Y_{ij} = X_{ij}\beta + b_i + \epsilon_{ij}, \qquad b_i \sim N(0, \sigma_b^2), \quad \epsilon_{ij} \sim N(0, \sigma^2), \quad b_i \perp \epsilon_{ij} \tag{2.1}\]| Term | Meaning |
|---|---|
| $X_{ij}\beta$ | the fixed effects — the population mean |
| $b_i$ | the subject effect: person $i$’s deviation from the population intercept, after covariates. A random variable; we estimate its variance, not the individual values |
| $\epsilon_{ij}$ | within-subject error |
Two different “mean responses” come out of this model, and the distinction runs through everything that follows:
\[\underbrace{E(Y_{ij} \mid b_i) = X_{ij}\beta + b_i}_{\text{conditional — one individual's trajectory}} \qquad\qquad \underbrace{E(Y_{ij}) = X_{ij}\beta}_{\text{marginal — the population average}}\]In the pictures: the marginal mean is the one thick line; the conditional means are the thin parallel lines, one per person, each shifted up or down by its own $b_i$ (Figure 2.10, left panel).
The induced covariance. Because $b_i$ is shared by all of person $i$’s measurements,
\[\text{Var}(Y_{ij}) = \sigma_b^2 + \sigma^2, \qquad \text{Cov}(Y_{ij}, Y_{ik}) = \sigma_b^2\]so the random intercepts model induces compound symmetry — same variance at every occasion, same covariance for every pair. $\sigma_b^2$ and $\sigma^2$ are the variance components. Often they are nuisance parameters on the way to $\beta$, but they also give the intraclass correlation,
\[\text{ICC} = \frac{\sigma_b^2}{\sigma_b^2 + \sigma^2}\]— between 0 and 1, the proportion of total variance “explained” by between-subject variability; equivalently the common correlation $\rho$ between any two measurements on the same person.
Estimation notes. With the normality assumptions the likelihood can be written
out, including both fixed and random terms; REML maximizes it. Degrees of freedom
are tricky — treating a random effect per person as a parameter over-penalizes,
so approximations are used; R’s lmerTest reports Satterthwaite df.
Answering the three questions
Everything in the dental analysis is a random intercepts model with different fixed effects. Age is coded with indicators — $X_{ij1}, X_{ij2}, X_{ij3}$ for ages 10, 12, 14 with age 8 as reference — which assumes nothing about the shape of the trajectory (not forced linear), at the cost of only working when everyone is measured at the same times.
| Question | Model | $H_0$ |
|---|---|---|
| Time effect — is the mean response varying with time? | $Y_{ij} = \beta_0 + \beta_1 X_{ij1} + \beta_2 X_{ij2} + \beta_3 X_{ij3} + b_i + \epsilon_{ij}$ | $\beta_1 = \beta_2 = \beta_3 = 0$ |
| Group effect — is the mean response different between groups? | add $\beta_4 X_{ij4}$, $X_{ij4} = 1$ for boys | $\beta_4 = 0$ |
| Time × group — is the change over time different by group? | add $\beta_5 X_{ij1}X_{ij4} + \beta_6 X_{ij2}X_{ij4} + \beta_7 X_{ij3}X_{ij4}$ | $\beta_5 = \beta_6 = \beta_7 = 0$ |
In the group model $E(Y_{ij} \mid \text{boy}) - E(Y_{ij} \mid \text{girl}) = \beta_4$, so with two groups the test is a single coefficient; with more groups it would be a joint (ANOVA-type) test like the time effect.
Linear time
Replace the indicators with age as a number. The interaction model becomes
\[Y_{ij} = \beta_0 + \beta_1 X_{ij1} + \beta_2 X_{ij2} + \beta_3 X_{ij1}X_{ij2} + b_i + \epsilon_{ij}\]with $X_{ij1}$ the sex indicator and $X_{ij2}$ age. Table 2.1 in the notes:
| Estimate | SE | $p$ | |
|---|---|---|---|
| (Intercept) | 17.373 | 1.184 | $< 2 \times 10^{-16}$ |
| sexBoy | $-1.032$ | 1.537 | 0.504 |
| age | 0.480 | 0.093 | $2 \times 10^{-6}$ |
| sexBoy:age | 0.305 | 0.121 | 0.014 |
| $\sigma_b^2$ (id intercept) | 3.299 | ||
| $\sigma^2$ (residual) | 1.922 |
Reading it exactly as ordinary regression with an interaction:
- 17.37 — mean distance for girls (reference) at age 0. Not meaningful (no teeth at age 0), but that is where the intercept lives.
- $-1.03$ — boys minus girls at age 0. Same caveat.
- 0.48 — change in mean distance per year of age for girls.
- 0.30 — the additional per-year change for boys, so the boys’ slope is $0.48 + 0.30 = 0.78$ mm/year.
Conclusion: a time effect and a time × group interaction. Because the model is linear, the interpretation of $\beta$ is identical to ordinary regression — the random effect changes the inference, not the meaning of the coefficients. (That sentence stops being true for non-linear models; see GEE below.)
Predicting the random effects: BLUPs
Usually the three questions are answered by inference on $\beta$ alone and the $b_i$ are never looked at. But sometimes you want subject-specific response profiles, $X_i\hat\beta + \hat{b}_i$. The estimate of $b_i$ is its conditional mean given the data, the best linear unbiased predictor (BLUP), also called the empirical Bayes estimate because the same formula comes out of a fully Bayesian derivation. For the random intercepts model it has the form
\[\hat{b}_i = \frac{n\,\sigma_b^2}{n\,\sigma_b^2 + \sigma^2}\,(\bar{Y}_{i\cdot} - \bar{Y}_{\cdot\cdot})\]Read it in two pieces. $\bar{Y}{i\cdot} - \bar{Y}{\cdot\cdot}$ — person $i$’s mean minus the grand mean — is the naive “how different is this person” estimate. The fraction in front shrinks it toward zero, which is why BLUPs are called shrunken estimates:
| Shrinkage | Intuition | |
|---|---|---|
| more measurements per person ($n$ large) | less | with 100 points on you, believe your own difference |
| $\sigma_b^2$ large relative to $\sigma^2$ | less | people genuinely differ; trust the individual |
| $\sigma^2$ large relative to $\sigma_b^2$ | more | mostly noise within people; pull everyone toward the average |
So a subject-specific line is not the same as fitting a regression to that person alone — it still borrows strength from everyone through $\hat\beta$ and the shrinkage.
Random intercepts and slopes
\[Y_{ij} = \beta_0 + X_{ij}\beta + b_i + \alpha_i X_{ij} + \epsilon_{ij}, \qquad \begin{pmatrix} b_i \\ \alpha_i \end{pmatrix} \sim N\!\left(\begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} \sigma_b^2 & \sigma_{b\alpha} \\ \sigma_{b\alpha} & \sigma_\alpha^2 \end{pmatrix}\right) \tag{2.2}\]Each person now has their own intercept and their own slope, bivariate normal — so there are two new things to estimate, not one: $\sigma_\alpha^2$ and the covariance $\sigma_{b\alpha}$. Regrouping,
\[E(Y_{ij} \mid b_i, \alpha_i) = (\beta_0 + b_i) + (\beta + \alpha_i)X_{ij}, \qquad E(Y_{ij}) = \beta_0 + \beta X_{ij}\]$b_i$ is the deviation from the population intercept, $\alpha_i$ from the population slope (Figure 2.10, right panel). A random-slopes-only model (middle panel) drops $b_i$.
Lecture 8.2 — Mixed models in R
From R script — Lecture Video 35.
Same packages and the same wide-to-long dental_long as Notes 5; the new ones are
lme4 (fits), lmerTest (Satterthwaite tests), emmeans, car::Anova, Epi::ci.lin.
Time effect
lin_age = lmer(distance ~ measurement + (1 | id), data = dental_long)
summary(lin_age)
The formula reads like lm except for + (1 | id): “give each id its own
random intercept.” The output has three parts worth recognizing:
- a line saying REML and Satterthwaite’s method — the two estimation notes above;
- Random effects —
id (Intercept)variance $= \hat\sigma_b^2 = 4.46$ andResidual$= \hat\sigma^2 = 2.08$, so $\text{ICC} = 4.46 / (4.46 + 2.08) = 0.68$; - Fixed effects — the usual estimate / SE / df / $t$ / $p$ table:
| Estimate | SE | $p$ | |
|---|---|---|---|
| (Intercept) | 22.19 | 0.49 | |
| age 10 | 0.98 | 0.39 | 0.015 |
| age 12 | 2.46 | 0.39 | $< 0.001$ |
| age 14 | 3.91 | 0.39 | $< 0.001$ |
Mean distance at age 8 (the reference) is 22.18 mm; ages 10, 12, 14 are 0.98, 2.46, 3.91 mm higher. Three helpers:
ci.lin(lin_age) # 95% CIs for the fixed effects
ranova(lin_age) # is the random intercept needed?
tidy(emmeans(lin_age, "measurement"), conf.int = TRUE) # the mean AT each age, with CI
Anova(lin_age) # H0: all three age betas = 0
ranova tests the random effect: LRT $= 60.5$, $p = 7 \times 10^{-15}$ — the
random intercepts are needed (if this were not significant you could perhaps justify a
fixed-effects-only model and ignore the correlation). emmeans turns the
“difference from age 8” coefficients into what you actually want: the estimated mean
at each age (22.2, 23.2, 24.6, 26.1) with SE and CI. Anova is the Wald test of the
time effect, $\chi^2_3 = 114$, $p < 10^{-16}$: the trajectory is not flat.
Plotting the estimate over the spaghetti plot uses predict with re.form = ~ 0,
meaning “no random effects” — the marginal mean, one blue dot per age:
dental_fit = bind_cols(dental_long, pred_age = predict(lin_age, re.form = ~ 0))
ggplot(dental_fit, aes(age, distance)) + geom_line(aes(group = factor(id))) +
geom_point(aes(y = pred_age), col = "blue", size = 2)
Drop re.form = ~ 0 and you get a dot per person instead.
Group effect
lin_agesex = lmer(distance ~ measurement + sex + (1 | id), data = dental_long)
sexBoy $= 2.32$ (95% CI 0.83–3.81, $p = 0.005$): boys’ mean distance is 2.3 mm
higher than girls’, adjusting for time. Girl is the reference. With two groups this
single coefficient is the test; no Anova needed.
tidy(emmeans(lin_agesex, c("measurement", "sex"))) gives the eight means (girls
20.8 → 24.7, boys 23.1 → 27.0). Plotted with col = sex, the blue dots sit uniformly
above the red — parallel by construction, since there is no interaction yet.
Time × group interaction
lin_agesexinter = lmer(distance ~ measurement*sex + (1 | id), data = dental_long)
Anova(lin_agesexinter, type = 3)
measurement*sex expands to main effects plus interaction. The table is now busy —
three age betas, one sex beta, three interaction betas. Anova(type = 3) tests each
term with its higher-order terms left in:
| $\chi^2$ | df | $p$ | |
|---|---|---|---|
| measurement | 25.6 | 3 | $10^{-5}$ |
| sex | 3.6 | 1 | 0.059 |
| measurement:sex | 7.1 | 3 | 0.069 |
Marginal evidence of an interaction. In the plot the blue and red dots are no longer parallel — the boys’ slope is steeper.
Linear time
lin_agec = lmer(distance ~ age + (1 | id), data = dental_long) # age numeric
lin_agecsexinter = lmer(distance ~ sex*age + (1 | id), data = dental_long)
Using age (a number) instead of measurement (a factor) fits one slope: 0.66
mm/year. The interaction model is Table 2.1 — and now the interaction is clearly
significant ($p = 0.014$ vs. 0.069 before). Two reasons: more power, since three
betas are estimated instead of seven; and it is a different model, one that forces
the trajectory to be linear. Plotting needs a grid of hypothetical people:
pred_agecsexinter = expand.grid(age = seq(8, 14, .5), sex = levels(dental_long$sex))
pred_agecsexinter$pred = predict(lin_agecsexinter, newdata = pred_agecsexinter, re.form = ~ 0)
ggplot(pred_agecsexinter, aes(age, pred, col = sex)) + geom_line()
Subject-specific trajectories
Children 10 (a girl) and 21 (a boy) from lin_agec. The one change: omit
re.form = ~ 0 and predict includes the BLUP $\hat{b}_i$:
pred.data = expand.grid(age = seq(8, 14, .5), id = c(10, 21))
pred.data$indiv_pred = predict(lin_agec, newdata = pred.data) # X beta + b_i
pred.data$marg_pred = predict(lin_agec, newdata = pred.data, re.form = ~ 0) # X beta
| age 8 | age 14 | $\hat{b}_i$ | |
|---|---|---|---|
| marginal | 22.04 | 26.00 | — |
| child 10 | 17.09 | 21.05 | $-4.96$ |
| child 21 | 26.96 | 30.92 | $+4.91$ |
The plot is Figure 2.15: the thick marginal line, and two parallel lines
$\hat{b}_i$ above and below it — the vertical gap is the BLUP. left_join-ing the
children’s actual measurements adds the observed points; the NAs it produces at the
half-year grid ages are harmless.
Random slopes and intercepts
lin_agecr = lmer(distance ~ age + (age | id), data = dental_long)
lin_agecsexinterr = lmer(distance ~ age*sex + (age | id), data = dental_long)
VarCorr(lin_agecr) # just the random-effects block
(age | id) instead of (1 | id) asks for random intercepts and slopes. The
random-effects block grows accordingly — for the interaction model:
| Variance | SD | Corr | |
|---|---|---|---|
| id (Intercept) | 5.77 | 2.40 | |
| id age | 0.032 | 0.18 | $-0.67$ |
| Residual | 1.72 | 1.31 |
Three things estimated where there was one: $\sigma_b^2$, $\sigma_\alpha^2$, and their correlation. The fixed effects are identical to the random-intercepts fit (17.37, $-1.03$, 0.48, 0.30) — only the SEs and $p$-values move (interaction $p = 0.033$ vs. 0.014) — same conclusion. Plotting every child’s line in grey under the two marginal lines (Figure 2.16) shows them crossing: different intercepts and different slopes.
Random slopes only is the clunky (0 + age | id) — the 0 + removes the
intercept; without it R assumes you want both. Only $\sigma_\alpha^2$ and $\sigma^2$
remain; conclusions unchanged.
Lecture 8.3 — Generalized estimating equations
Why marginal models exist
For the linear models so far, the interpretation of $\beta$ did not depend on how the correlation was handled — random intercepts, random slopes, it made no difference to what “0.48 mm per year for girls” meant. Here is why. With a random intercept,
\[E(Y_i) = E\{E(Y_i \mid b_i)\} = E(X_i\beta + b_i) = X_i\beta + E(b_i) = X_i\beta\]because expectation is linear: the $b_i$ term drops out and the marginal mean is the subject-specific mean averaged over people. For a non-linear link — logistic, Poisson —
\[E(Y_i) = E\{h(X_i\beta + b_i)\} \ne h(X_i\beta)\]in general. The random effect no longer averages away, so what you assume about the correlation structure changes what the coefficients mean. This forces a distinction:
| Marginal model | Mixed-effects model | |
|---|---|---|
| also called | population-average | subject-specific |
| $\beta$ describes | the population mean response | an individual’s response, given their $b_i$ |
| estimated by | GEE | (RE)ML |
A marginal model is one whose mean depends only on the covariates — no random effects, no previous responses. And it requires no distributional assumption, only a regression model for the mean; dropping the likelihood is what leads to generalized estimating equations.
The three-part specification
A marginal model for longitudinal data is specified by:
- The mean, through a known link $g$: $g(\mu_{ij}) = \eta_{ij} = X_{ij}\beta$, where $\mu_{ij} = E(Y_{ij} \mid X_{ij})$.
- The variance, as a function of the mean: $\text{Var}(Y_{ij}) = \phi\, v(\mu_{ij})$, with $v$ a known variance function and $\phi$ a scale parameter, known or estimated.
- The within-subject association, as a function of extra parameters $\alpha$.
Nowhere is a distribution named — but three things are assumed. “No distributional assumptions” is not “no assumptions”; nothing is free.
Linear case, side by side with what we already know:
| Part | Specification | Equivalent to |
|---|---|---|
| 1 | identity link, $E(Y_{ij}) = X_{ij}\beta$ | the usual linear model |
| 2 | $v(\mu) = 1$, $\text{Var}(Y_{ij}) = \phi$ | constant variance over time, $\sigma^2$ |
| 3 | $\text{Corr}(Y_{ij}, Y_{ik}) = \alpha$ for all $j \ne k$ | compound symmetry again |
The random intercepts model, the CS pattern model, and this GEE all land on the same matrix. In the linear case the point is almost invisible; in the non-linear case it is everything.
Binary case — $Y_{ij} \in {0,1}$ (were you happy when you woke up, every day for 30 days):
| Part | Specification |
|---|---|
| 1 | $\text{logit}(\mu_{ij}) = X_{ij}\beta$ — looks exactly like logistic regression |
| 2 | $\text{Var}(Y_{ij}) = \mu_{ij}(1 - \mu_{ij})$, $\phi = 1$ — forced, a Bernoulli’s variance is a function of its mean |
| 3 | a common log odds ratio $\log \text{OR}(Y_{ij}, Y_{ik}) = \alpha$ — same flavour as “all off-diagonals equal,” in odds-ratio form because $Y$ is binary |
Estimation
No likelihood, so no MLE. The GEE estimator of $\beta$ can be thought of as minimizing
\[\sum_{i=1}^{N} \{\mathbf{Y}_i - \mu_i(\beta)\}^T V_i^{-1} \{\mathbf{Y}_i - \mu_i(\beta)\}\]where $\mu_{ij}(\beta) = g^{-1}(X_{ij}\beta)$ comes from part 1 and $V_i$, the working covariance matrix, from parts 2 and 3. The minimizer, if it exists, solves the generalized estimating equations
\[\sum_{i=1}^{N} D_i^T V_i^{-1}(\mathbf{Y}_i - \mu_i) = 0, \qquad D_i = \partial\mu_i / \partial\beta\]by an iterative two-stage procedure that software does for you. The resulting $\hat\beta$ describes the effect of covariates on the population mean response.
The difference in interpretation: ECG example
67 subjects, active drug vs. placebo, binary outcome $Y = 1$ if the ECG is abnormal, measured at two time points. $X_{ij1}$ = drug (1 active, 0 placebo), $X_{ij2}$ = time (0 at time 1, 1 at time 2).
Marginal model, estimated by GEE, with a common log odds ratio $\alpha$ for the within-subject association:
\[\text{logit}(\mu_{ij}) = \beta_0 + \beta_1 X_{ij1} + \beta_2 X_{ij2}, \qquad \hat\beta_1 = 0.57\]The odds of an abnormal ECG are $e^{0.57} = 1.77$ times higher when treated with the active drug versus placebo.
Generalized linear mixed model, random intercept:
\[\text{logit}\{E(Y_{ij} \mid b_i)\} = \beta_0^* + \beta_1^* X_{ij1} + \beta_2^* X_{ij2} + b_i, \qquad b_i \sim N(0, \sigma_b^2), \qquad \hat\beta_1^* = 0.89\]Each patient has an underlying propensity $b_i$ for an abnormal ECG; the drug multiplies that patient’s odds by $e^{\beta_1^*}$ regardless of their propensity.
A patient’s odds of an abnormal ECG are $e^{0.89} = 2.4$ times higher when treated with the active drug than when treated with placebo.
Same data, same question, 1.77 vs. 2.4. The marginal coefficient is a statement about the population; the mixed-model coefficient is a statement about an individual. The subject-specific effect is typically larger in magnitude, and neither is wrong — they answer different questions. That is why the distinction between population-average and subject-specific models had to be made, and why it only bites for non-linear links.
Mixed model or GEE?
- Either way you must assume something about the correlation between time points. GEE assumes less (no distribution) but not nothing.
- Want variance components — between vs. within, the ICC? → mixed model. GEE has no between/within decomposition.
- Want subject-specific trajectories (the grey lines)? → mixed model. GEE cannot produce them.
- Want to avoid distributional assumptions? → GEE.
- Only care about covariate effects on the population mean? → GEE.
- For a non-linear model, remember the interpretations differ. For a linear model they are equivalent — demonstrated next.
Lecture 8.4 — GEE in R
From R script — Lecture Video 37.
New package geepack, function geeglm (generalized estimating equation,
generalized linear model). Same dental_long; the script jumps straight to the
interaction model with linear age.
gee_inter = geeglm(distance ~ age + sex + sex*age, data = dental_long,
id = id, corstr = "exchangeable")
summary(gee_inter)
The formula is the familiar part. What is missing is any (1 | id) — there are no
random effects. What is new is id = (which rows belong to the same cluster) and
corstr = (part 3 of the specification). "exchangeable" is geepack’s name for
compound symmetry: same variance on the diagonal, same covariance everywhere off it,
two parameters. Alternatives: "independence", "ar1", "unstructured",
"userdefined" — more flexible, more parameters, the usual trade-off. family
defaults to gaussian (identity link).
| Estimate | SE | Wald | $p$ | |
|---|---|---|---|---|
| (Intercept) | 17.373 | 0.725 | 573.9 | $< 10^{-16}$ |
| age | 0.480 | 0.063 | 57.7 | $3 \times 10^{-14}$ |
| sexBoy | $-1.032$ | 1.378 | 0.56 | 0.45 |
| age:sexBoy | 0.305 | 0.117 | 6.80 | 0.009 |
| scale $\hat\phi$ | 4.905 | |||
| $\hat\alpha$ | 0.618 | 0.131 |
Below the coefficients the output reports exactly the other two parts of the
specification: the scale (part 2, the common variance) and alpha (part 3, the
common correlation, 0.62 — compare the ICC of 0.63 from the random-intercepts model).
tidy(gee_inter, conf.int = TRUE) adds CIs; ggpredict(gee_inter, terms = c("age",
"sex")) builds the prediction grid for a two-line plot, as expand.grid did for lmer.
GEE vs. mixed model on the same data
lin_agecsexinterr = lmer(distance ~ age*sex + (age | id), data = dental_long)
ci.lin(lin_agecsexinterr)
tidy(gee_inter, conf.int = TRUE)
| GEE estimate (SE) | Mixed estimate (SE) | |
|---|---|---|
| (Intercept) | 17.373 (0.725) | 17.373 (1.228) |
| age | 0.480 (0.063) | 0.480 (0.104) |
| sexBoy | $-1.032$ (1.378) | $-1.032$ (1.595) |
| age:sexBoy | 0.305 (0.117) | 0.305 (0.135) |
The coefficients are identical; the standard errors and $p$-values are not. The
point estimates agree because this is a linear model, where the marginal and
subject-specific means coincide. The inference differs because the two methods handle
the repeated-measures correlation differently. With a logistic or Poisson model the
coefficients themselves would differ, as in the ECG example. You can tell which model
someone fit from the output alone: a random-effects block means lmer; scale and
alpha mean geeglm.
Class activity: aids from the JM package
A randomized trial of 467 HIV patients (ddI vs. ddC), with CD4 count measured
repeatedly (obstime), plus survival time, drug, gender, prevOI (previous
opportunistic infection: AIDS vs. noAIDS) and AZT (intolerance vs. failure). The
survival part waits for joint modelling; here only the CD4 trajectories are used. Unlike
the dental data, measurement is unbalanced — 1 to 5 visits per patient, 1405 rows.
- Spaghetti plot — the Notes 5 code with
obstime,CD4,col = factor(patient), legend removed. Not everyone is measured at every time. - Time effect —
lmer(CD4 ~ obstime + (1 | patient)):obstime$= -0.149$, $p < 10^{-16}$. CD4 falls over time. ($\hat\sigma_b^2 = 20.3$, $\hat\sigma^2 = 3.85$, ICC $\approx 0.84$ — most variation is between patients.) - Group effects —
+ drug:drugddI$= 0.54$, $p = 0.21$, no drug effect.+ prevOI:prevOIAIDS$= -4.78$, $p < 10^{-16}$ — patients with a previous infection have much lower CD4. - Interactions —
obstime*drug: $p = 0.67$;obstime*prevOI: $p = 0.94$. CD4 declines over time, and differs byprevOI, but the rate of decline differs by neither. -
GEE version —
geeglm(CD4 ~ obstime*prevOI, id = patient, corstr = "exchangeable"):Mixed (lmer) GEE obstime $-0.1500$ $-0.1518$ prevOIAIDS $-4.7755$ $-4.7818$ obstime:prevOIAIDS $-0.0019$ $-0.0036$ Same conclusions (no interaction, $p = 0.90$). But the coefficients are close, not identical — unlike the dental comparison. The reason is missing data: with unbalanced follow-up the two methods use all the available data but weight it differently when estimating the covariance structure, so with incomplete profiles even a linear model’s $\hat\beta$ diverges slightly. The dental data were complete, so there the match was exact. This is the first appearance of a theme §2.5 takes up properly.