Why Mixed-Effects?

Longitudinal data: same person measured repeatedly leads to correlated outcomes.

Mixed-effects = fixed effects (population-average) + random effects (person-specific deviations).

What Mixed-Effects Captures

Subject-Specific vs Marginal Views

“Marginal” vs “subject-specific” is a contrast of estimands, not methods. The LME delivers both:

Averaging the subject-specific mean over \(b_i \sim N(0, G)\) gives the marginal mean \(X_i\boldsymbol\beta\), so the LME also has a marginal interpretation. GEE is one way to target that marginal mean directly (with robust SEs); the LME targets it while also predicting subject trajectories and leveraging all available data under MAR.

Model Building Blocks

Let \(Y_{ij}\) be the response for subject \(i\) at time \(t_{ij}\).

Random-Intercept (RI) Model

\[ Y_{ij} = \beta_0 + X_{ij}\boldsymbol\beta + b_{0i} + \varepsilon_{ij} \]

where:

• \(b_{0i} \sim N(0, \sigma_b^2)\) is the subject-specific intercept
• \(\varepsilon_{ij} \sim N(0, \sigma^2)\) is the residual error

Random-Slope (RS) Model

\[ Y_{ij} = \beta_0 + \beta_1 t_{ij} + X_{ij}\boldsymbol\beta + b_{0i} + b_{1i} t_{ij} + \varepsilon_{ij} \]

where the random effects follow:

\[ \begin{pmatrix} b_{0i} \\ b_{1i} \end{pmatrix} \sim N\!\left(\mathbf{0}, G\right), \quad G = \begin{pmatrix} \sigma_{b0}^2 & \sigma_{b01} \\ \sigma_{b01} & \sigma_{b1}^2 \end{pmatrix} \]

Implied Covariance (RI Model)

For the random-intercept model:

\[ \mathrm{Var}(Y_{ij}) = \sigma_b^2 + \sigma^2 \]

The shared random intercept \(b_{0i}\) is the only thing two observations on subject \(i\) have in common, so the off-diagonal covariance is

\[ \mathrm{Cov}(Y_{ij}, Y_{ik}) = \sigma_b^2 \quad (j \ne k) \]

Dividing by the variance gives the correlation, which is the same for every pair of times:

\[ \mathrm{Corr}(Y_{ij}, Y_{ik}) = \frac{\sigma_b^2}{\sigma_b^2 + \sigma^2} \quad \text{(Compound Symmetry)} \]

Implied Covariance (RS Model)

For the random-slope model the implied variance changes with time (nonstationary):

\[ \mathrm{Var}(Y_{ij}) = \sigma_{b0}^2 + 2\, t_{ij}\, \sigma_{b01} + t_{ij}^2\, \sigma_{b1}^2 + \sigma^2 \]

This is a quadratic in \(t_{ij}\): trajectories fan out (or in) over time, unlike the constant variance of compound symmetry. This is the nonstationary, “linear-in-time” covariance shown in the RI+RS row of the Link to Chapter 7 table below (it is introduced here, not in Ch. 7).

RI vs RS: When and Why

If slopes vary across persons (common in biology), RS is often essential.

Two-Stage Random Effects Formulation (FLW 8.4)

FLW motivate the LME model in two stages, which makes clear where the random effects come from:

• Stage 1 (within-subject): \(\mathbf{Y}_i = \mathbf{Z}_i \boldsymbol\beta_i + \boldsymbol\varepsilon_i\), with \(\boldsymbol\varepsilon_i \sim N(0, R_i)\). Each subject follows their own regression with coefficients \(\boldsymbol\beta_i\).
• Stage 2 (between-subject): \(\boldsymbol\beta_i = \mathbf{A}_i \boldsymbol\beta + \mathbf{b}_i\), with \(\mathbf{b}_i \sim N(0, G)\). Here \(\mathbf{A}_i\) is a known between-subject design matrix linking each subject’s coefficients to the population means \(\boldsymbol\beta\) (often an identity, or a map of subject-level covariates). The subject-specific coefficients vary around \(\boldsymbol\beta\).

Substituting Stage 2 into Stage 1 gives the familiar single-equation form

\[ \mathbf{Y}_i = \mathbf{X}_i \boldsymbol\beta + \mathbf{Z}_i \mathbf{b}_i + \boldsymbol\varepsilon_i, \qquad \Sigma_i = \mathbf{Z}_i G \mathbf{Z}_i^\top + R_i . \]

The random-intercept and random-slope models are special cases: the two-stage view simply makes explicit that \(G\) is the between-subject covariance of the \(\mathbf{b}_i\).

The Link to Chapter 7

Mixed-effects give parsimonious, interpretable \(\Sigma_i\) via random effects rather than directly parameterizing \(\Sigma_i\).

Estimation Overview

Mixed-model equations (Henderson) yield BLUPs of random effects; with unknown variance components they become EBLUPs (empirical BLUPs, the FLW term for predicted random effects: BLUPs computed using the estimated variance components rather than known ones).

Shrinkage Intuition (EBLUP)

Borrowing strength: subject-specific intercepts/slopes are pulled toward population means.

Shrinkage is stronger when:

• Subject has few measurements
• Measurement error \(\sigma^2\) is large
• Random-effect variance is small

This prevents overfitting noisy within-subject patterns.

TLC Data Setup

library(dplyr); library(tidyr); library(nlme); library(ggplot2); library(emmeans)

path <- "../../data/tlc.csv"
tlc_wide <- read.csv(path, check.names = FALSE, stringsAsFactors = FALSE)
names(tlc_wide) <- trimws(names(tlc_wide))

tlc <- tlc_wide %>%
  rename(
    id  = ID,
    trt = `Treatment Group`,
    w0  = `Lead Level Week 0`,
    w1  = `Lead Level Week 1`,
    w4  = `Lead Level Week 4`,
    w6  = `Lead Level Week 6`
  ) %>%
  mutate(trt = factor(trt, levels = c("P","A"), labels = c("Placebo","Succimer"))) %>%
  pivot_longer(c(w0,w1,w4,w6), names_to = "wk", values_to = "y") %>%
  mutate(
    time     = recode(wk, w0="0", w1="1", w4="4", w6="6"),
    time_f   = factor(time, levels = c("0","1","4","6")),
    time_num = as.numeric(as.character(time)),
    id       = factor(id),
    group    = trt
  ) %>%
  select(id, group, time_f, time_num, y) %>%
  arrange(id, time_num)

stopifnot(all(!is.na(tlc$y)))
dplyr::count(tlc, id)

TLC Visuals: Spaghetti Plot

ggplot(tlc, aes(time_num, y, group = id, color = group)) +
  geom_line(alpha = 0.25) +
  stat_summary(aes(group = group), fun = mean, geom = "line", linewidth = 1.1) +
  labs(title = "Lead levels over time (TLC)", x = "Week", y = "Lead") +
  theme_minimal() + guides(color = guide_legend(title = "Group"))

Random-Intercept Model (RI)

lme_ri <- lme(y ~ group * time_f, random = ~ 1 | id, data = tlc, method = "REML")
summary(lme_ri)$tTable[1:8, , drop = FALSE]

Reading VarCorr Output: Random-Intercept Model

RI Model Specification

• Fixed: group x time (factor) = maximal mean for covariance selection
• Random: intercept only (\(b_{0i}\))

Random-Slope Model (RS)

A random slope needs numeric time (time_num) so each subject has a single slope coefficient. (We keep the factor time_f for the population mean in the next model; the consequence of mixing factor and numeric time for model comparison shows up on the AIC slide.)

ctrl <- lmeControl(
  niterEM     = 200,
  msMaxIter   = 200,
  msMaxEval   = 400,
  tolerance   = 1e-6,
  returnObject= TRUE
)

lme_rs <- lme(y ~ group * time_num, random = ~ time_num | id, data = tlc, method = "REML", control = ctrl)
summary(lme_rs)$tTable

Random coefficients allow individual-specific baseline and slope.

Reading VarCorr Output: Random-Slope Model

RI vs RS: AIC/REML Comparison

data.frame(
  model = c("RI (time_f)", "RS (time_num)"),
  AIC   = c(AIC(lme_ri), AIC(lme_rs)),
  logLik= c(logLik(lme_ri), logLik(lme_rs))
)

Not strictly comparable if fixed effects differ (time factor vs numeric): REML AIC compares random structures only when the fixed effects match (see Estimation Overview).

Factor-Time Fixed Effects with RS Random

Strategy: keep time as factor in fixed effects (rich mean), use time_num in random part to permit RS.

lme_rifs <- lme(y ~ group * time_f, random = ~ time_num | id, data = tlc, method = "REML", control = ctrl)
data.frame(
  model = c("RI", "RI+RS"),
  AIC   = c(AIC(lme_ri), AIC(lme_rifs)),
  logLik= c(logLik(lme_ri), logLik(lme_rifs))
)

Interpreting Fixed Effects with RS

EMMs for Contrasts at Planned Times

emm_ri <- emmeans(lme_ri, ~ group | time_f)
pairs(emm_ri)

Same use of emmeans as in Ch. 7; SEs now reflect random-effects structure.

BLUPs: Subject-Specific Intercepts/Slopes

re <- ranef(lme_rifs)
head(re)

So what? A wide spread of these intercept BLUPs \(b_{0i}\) means substantial between-subject variation in baseline, consistent with a large \(\sigma_{b0}\) in VarCorr: the histogram is the picture of that variance component. (The QQ-plot slide recomputes re from the same fit.)

Heterogeneous Residual Variances + RI/RS

Mixed models can combine random effects + a within-subject residual covariance \(R_i\):

• Heterogeneous residual SDs via weights = varIdent(~ 1 | time_f) (a separate residual SD per visit, instead of one constant \(\sigma\))
• Additional within-subject correlation via corAR1(~ time_num | id) (residual AR(1): correlation that decays with the time gap) if needed

lme_rifs_het <- lme(y ~ group * time_f, random = ~ time_num | id,
                    weights = varIdent(~ 1 | time_f),
                    data = tlc, method = "REML", control = ctrl)
AIC(lme_rifs, lme_rifs_het)

Adding Residual AR(1) to RI/RS

lme_rifs_ar1 <- lme(y ~ group * time_f, random = ~ time_num | id,
                    correlation = corAR1(form = ~ time_num | id),
                    data = tlc, method = "REML", control = ctrl)
data.frame(model = c("RI+RS", "RI+RS + AR(1)"),
           AIC   = c(AIC(lme_rifs), AIC(lme_rifs_ar1)))

Only add AR(1) if random effects do not remove residual serial correlation.

Diagnostics: Residuals vs Fitted

res <- resid(lme_rifs, type = "normalized")
fit <- fitted(lme_rifs)

ggplot(data.frame(fit, res), aes(fit, res)) +
  geom_point(alpha = 0.5) + geom_hline(yintercept = 0, linewidth = 0.3) +
  labs(title = "Residuals vs fitted (normalized)", x = "Fitted", y = "Residual") +
  theme_minimal()

Normalized residuals are Cholesky-standardized: they strip out the modeled within-subject correlation, which is why we check these (not raw residuals) for a correlated-data model.

Diagnostics: ACF of Residuals

acf_out <- ACF(lme_rifs, resType = "normalized")
ggplot(acf_out, aes(lag, ACF)) + geom_col(width = 0.1) +
  geom_hline(yintercept = 0, linewidth = 0.3) +
  labs(title = "ACF of normalized residuals", x = "Lag", y = "ACF",
       subtitle = "Reading rule: bars within the band (flat) = good; slow decay = leftover serial correlation") +
  theme_minimal()

Diagnostics: Semivariogram

vg <- Variogram(lme_rifs, form = ~ time_num | id, resType = "normalized")
ggplot(vg, aes(dist, variog)) + geom_point() +
  geom_smooth(se = FALSE, span = 0.9) +
  geom_hline(yintercept = 1, linetype = 3) +
  labs(title = "Semivariogram of normalized residuals", x = "Lag (weeks apart)", y = expression(hat(gamma)),
       subtitle = "Reading rule: flat near 1 = correlation captured; rising from a low value = consider AR(1)") +
  theme_minimal()

Diagnostics: Random-Effects QQ Plot

re <- ranef(lme_rifs)
qq <- as.data.frame(re)
ggplot(qq, aes(sample = `(Intercept)`)) +
  stat_qq() + stat_qq_line() +
  labs(title = "QQ plot: predicted random intercepts (BLUPs)", x = "Theoretical", y = "Sample") +
  theme_minimal()

Testing Fixed Effects

FLW Ch. 8 uses large-sample inference: Wald z-statistics and LRTs without a small-sample df adjustment (Tables 8.5 and 8.12 report z-statistics).

lme_rifs_ML <- update(lme_rifs, method = "ML")
lme_rifs_noInt_ML <- lme(y ~ group + time_f, random = ~ time_num | id, data = tlc, method = "ML", control = ctrl)
anova(lme_rifs_noInt_ML, lme_rifs_ML)

Testing Variance Components (Boundary)

Variance-component tests are boundary tests (a variance cannot be negative), so the null distribution is a 50:50 mixture, not a plain \(\chi^2\). Which mixture depends on how many parameters you drop:

• One variance (e.g. \(\sigma_{b0}^2 = 0\): is a random intercept needed at all?): \(\tfrac12\chi^2_0 + \tfrac12\chi^2_1\). Only here is the corrected p the naive \(\chi^2_1\) p-value halved.
• A random slope, dropped from a model that already has a random intercept (our example below: removes \(\sigma_{b1}^2\) and the intercept-slope covariance \(\sigma_{b01}\), two parameters): \(\tfrac12\chi^2_1 + \tfrac12\chi^2_2\) (Stram & Lee). Not a simple halving.

lrt <- anova(lme_ri, lme_rifs)   # REML; same fixed effects, drop the random slope
lrt

The naive anova() \(\chi^2_2\) p-value is conservative (too large). The mixture above is the correct adjustment for this two-parameter test.

Centering Time for Interpretation

With RS, centering time at a meaningful point (e.g., baseline or mid-study) stabilizes covariance and aids interpretation:

• Intercept = mean response at the centering time
• Random-intercept/slope correlation becomes more interpretable

tlc <- tlc %>% mutate(time_c = time_num - mean(time_num))
lme_rifs_c <- lme(y ~ group * time_c, random = ~ time_c | id, data = tlc, method = "REML", control = ctrl)
VarCorr(lme_rifs_c)

Beyond Linear Time: Splines

library(splines)
lme_spl <- lme(y ~ group * ns(time_num, df = 3),
               random = ~ time_num | id, data = tlc, method = "REML", control = ctrl)
AIC(lme_rifs, lme_spl)

Flexible mean over time; keep random slopes (or random spline coefficients if enough data).

Time-Varying Covariates (TVCs)

Mixed models handle TVCs naturally.

Watch out for endogeneity (covariate depends on past \(Y_{ij}\)).

Interactions with time (fixed and/or random slopes) often required.

Missing Data (MAR) and LMMs

Under MAR, ML/REML estimates remain valid with an appropriate mean + random-effects structure.

LMMs use all available data without ad hoc imputation.

Simulation: Shrinkage and EBLUPs

set.seed(8)
N <- 200; m <- 3
ids <- factor(seq_len(N))
times <- 0:2
G <- matrix(c(1.0, 0.3, 0.3, 0.4), 2, 2)  # random-effects covariance (FLW's G)
sigma2 <- 1.2^2                            # residual variance

sim <- do.call(rbind, lapply(ids, function(i){
  b <- MASS::mvrnorm(1, mu = c(0,0), Sigma = G)
  tibble::tibble(id = i, time = times,
                 y = 5 + 0.6*time + b[1] + b[2]*time + rnorm(length(times), sd = sqrt(sigma2)))
}))

fit <- lme(y ~ time, random = ~ time | id, data = sim, method = "REML", control = ctrl)
blup <- ranef(fit)
head(blup)

Note BLUPs near 0 when subject has little information (shrinkage).

Visual: Shrinkage with Few vs Many Points

set.seed(81)
ids_few <- sample(levels(sim$id), 6)
pred <- expand.grid(id = ids_few, time = seq(0,2,by=.05))
pred$pop <- predict(fit, newdata = pred, level = 0)
pred$ss  <- predict(fit, newdata = pred, level = 1)

ggplot(sim %>% dplyr::filter(id %in% ids_few),
       aes(time, y, group = id)) +
  geom_point() +
  geom_line(aes(y = pop, group = id), data = pred, linetype = 2) +
  geom_line(aes(y = ss,  group = id), data = pred) +
  facet_wrap(~ id, ncol = 3) +
  labs(title = "Shrinkage stronger with sparse/noisy subjects", x = "Time", y = "Y") +
  theme_minimal()

Prediction and Intervals

Prediction intervals require accounting for both residual and random-effect uncertainty, so there is no one-line predict() interval for a subject-specific prediction. The standard route (stated as reference, not run here):

Influence and Outliers

• Check standardized residuals, random-effect outliers, leverage
• Consider robust mixed models or transformations if heavy tails/outliers

Random-Effects Covariance Interpretation

When to Add Residual Correlation

If RI/RS + heterogeneous residual variances still leave serial correlation in residuals (ACF/variogram), add corAR1 or similar.

Do not over-layer: random slopes often absorb a lot of serial correlation.

Practical Modeling Sequence

Code Template for Exercise

# 1) RI
fit_ri <- lme(y ~ group * time_f, random = ~ 1 | id, data = tlc, method = "REML")

# 2) RI+RS (time_num in random part)
fit_rs <- lme(y ~ group * time_f, random = ~ time_num | id, data = tlc, method = "REML", control = ctrl)

# 3) Add hetero residual SDs by visit
fit_rs_het <- update(fit_rs, weights = varIdent(~ 1 | time_f))

# Compare REML AIC
data.frame(model=c("RI","RI+RS","RI+RS+het"),
           AIC=c(AIC(fit_ri), AIC(fit_rs), AIC(fit_rs_het)))

Common Pitfalls

Common Mistakes to Avoid

RI/RS vs Pure Covariance Models (Ch. 7)

• Many covariance patterns can be mimicked by random effects
• Mixed models interpret \(\Sigma_i\) through heterogeneity in baseline/slopes
• For prediction and patient-specific effects, mixed models are natural

Reporting Checklist (LMM)

Quick Reference: nlme::lme Syntax

# Random-intercept
lme(y ~ fixed_terms, random = ~ 1 | id, data = dat)

# Random-slope (linear time)
lme(y ~ fixed_terms, random = ~ time | id, data = dat)

# Add hetero residual SDs by visit
lme(y ~ fixed_terms, random = ~ time | id,
    weights = varIdent(~ 1 | time_factor), data = dat)

# Add residual AR(1)
lme(y ~ fixed_terms, random = ~ time | id,
    correlation = corAR1(~ time | id), data = dat)

The Same Model in SAS (PROC MIXED)

R is our runnable language (used in the slides above and in the homework). For students who work in SAS, here is the equivalent, shown for reference only and not executed. This mirrors FLW’s Computing section (8.9).

/* Not run here. SAS counterpart of the random-slope fit above, on the TLC data. */
PROC MIXED DATA=tlc METHOD=REML;
  CLASS id trt;
  MODEL y = time trt time*trt / SOLUTION;
  RANDOM intercept time / SUBJECT=id TYPE=UN G;
RUN;

• RANDOM ... / SUBJECT=id is the SAS counterpart of R’s random = ~ time | id (or (time | id) in lmer).
• METHOD=REML matches the nlme::lme/lmer default; use METHOD=ML for LRTs of fixed effects.
• TYPE=UN is an unstructured \(G\) (correlated intercept and slope); TYPE=VC makes them independent.

Quick Reference: ML vs REML

Side Note: lme vs lmer

Choose based on needs; results for fixed effects are often similar.

Worked Example: ML vs REML

fit_full_ML <- lme(y ~ group * time_f, random = ~ time_num | id, data = tlc, method = "ML", control = ctrl)
fit_noint_ML <- lme(y ~ group + time_f, random = ~ time_num | id, data = tlc, method = "ML", control = ctrl)
anova(fit_noint_ML, fit_full_ML)

Subject-Specific Predictions Table

sub_ids <- head(levels(tlc$id), 6)
new6 <- tlc %>% distinct(id, group) %>% filter(id %in% sub_ids) %>%
  mutate(time_num = 6, time_f = factor("6", levels = c("0","1","4","6")))
data.frame(
  id = sub_ids,
  pop = predict(lme_rifs, newdata = new6, level = 0),
  ss  = predict(lme_rifs, newdata = new6, level = 1)
)

Model Criticism: Overdispersion/Heavy Tails

If residuals or random effects are heavy-tailed, consider:

• Transformations, variance functions, robust LMMs (e.g., t-distributed random effects/errors)
• Check sensitivity to outliers

Time as Factor vs Continuous: Guidance

Growth Curves vs Spline LMMs

• Polynomial growth curves are easy but can oscillate
• Natural cubic splines are more stable
• Choose df by AIC/CV and subject-matter considerations

TVCs and Lagged Effects

• Incorporate lagged covariates if plausible causal ordering
• Beware of post-treatment covariates when estimating treatment effects

Beyond Gaussian LMMs

• GLMMs for non-Gaussian outcomes (binary, counts)
• Similar principles; estimation via Laplace/AGQ
• Not covered in this chapter

Cheat-Sheet: Model Choice

Cheat-Sheet: Testing

Do not mix ML and REML across compared models.

What to Include in Your Report

• Data structure and missingness; time coding
• Final mean, random, residual structures; estimation method
• Key estimates (fixed) with CIs; variance components
• Diagnostics; sensitivity; clear interpretations

Recipe Card: Fitting Mixed-Effects Models

# Ensure data is in long format with:
# - Subject ID as factor
# - Time as both factor (for rich mean) and numeric (for slopes)
dat <- dat %>%
  mutate(id = factor(id),
         time_f = factor(time),
         time_num = as.numeric(time),
         time_c = time_num - mean(time_num))  # centered

# Start with random intercept
fit_ri <- lme(y ~ group * time_f, random = ~ 1 | id,
              data = dat, method = "REML")

# Add random slope
fit_rs <- lme(y ~ group * time_f, random = ~ time_num | id,
              data = dat, method = "REML")

# Add heterogeneous residual variances
fit_rs_het <- update(fit_rs, weights = varIdent(~ 1 | time_f))

AIC(fit_ri, fit_rs, fit_rs_het)  # Lower is better

plot(fit_rs)                    # Residuals vs fitted
qqnorm(resid(fit_rs, type = "normalized"))
ACF(fit_rs)                     # Check for leftover correlation
qqnorm(ranef(fit_rs)[,1])       # BLUPs: spot unusual subjects, NOT a normality test (FLW p.273)

fit_full_ML <- update(fit_rs, method = "ML")
fit_reduced_ML <- lme(y ~ group + time_f, random = ~ time_num | id,
                      data = dat, method = "ML")
anova(fit_reduced_ML, fit_full_ML)  # LRT for interaction

summary(fit_rs)$tTable           # Fixed effects
VarCorr(fit_rs)                  # Variance components
emmeans(fit_rs, ~ group | time_f) %>% pairs()  # Contrasts

Further Reading and For Next Time

• FLW Ch. 8: the Case Studies (8.8), the Computing section on PROC MIXED (8.9), and Further Reading (8.10).
• Laird & Ware (1982), the original two-stage random-effects formulation behind 8.4.
• Next: FLW Ch. 9, Fixed Effects versus Random Effects Models and the Hausman test (Lecture 9).