From Conceptual to Estimation

Until now: models for longitudinal data were conceptual (mean + covariance)

Starting here: focus on estimation of unknown parameters

The Framework

Conditional multivariate normal model:

\[ \mathbb{E}(\mathbf{Y}_i \mid \mathbf{X}_i) = \mathbf{X}_i \beta \]

\[ \text{Cov}(\mathbf{Y}_i \mid \mathbf{X}_i) = \Sigma_i = \Sigma_i(\boldsymbol{\theta}) \]

Here \(\Sigma_i\) is the within-subject covariance for subject \(i\) (FLW’s \(\Sigma_i(\boldsymbol\theta)\)), and \(\boldsymbol{\theta}\) collects its covariance parameters.

Interpreting the Covariance Matrix

\(\Sigma_i\) encodes within-subject covariance of repeated measures

Parameterized by a vector \(\boldsymbol{\theta}\) of covariance parameters

Covariance Structure Options

Why Covariance Structure Matters

Choice of covariance structure directly impacts:

Looking Ahead: Chapters 5-8

Chapters 5-8 will:

• Focus on models for the mean and covariance
• Explore estimation under different covariance assumptions
• Develop tools for inference when \(\Sigma_i\) is known vs estimated

Next up, Lecture 5 (Ch. 5, response profiles) applies this estimation and inference machinery (GLS/ML, REML, the LRT) to the mean response over time.

Discussion Question

What might happen if we assume compound symmetry when the true structure is autoregressive (AR(1))?

Example: Simulating Covariance Structures

set.seed(667)

# Parameters
n_subj <- 25   # enough subjects to SEE constant vs decaying correlation
n_time <- 5
time <- 1:n_time

# Compound symmetry covariance (rho constant)
rho_cs <- 0.8
Sigma_cs <- matrix(rho_cs, n_time, n_time)
diag(Sigma_cs) <- 1

# AR(1) covariance (rho decays with lag)
# row()/col() return the row and column index of each matrix entry, so
# abs(row - col) is the time-lag |j - k| at each cell; rho^lag gives the
# AR(1) correlation that decays as the lag grows.
rho_ar1 <- 0.8
Sigma_ar1 <- rho_ar1 ^ abs(row(Sigma_cs) - col(Sigma_cs))

# Independence covariance (rho is 0)
rho_in <- 0
Sigma_in <- matrix(rho_in, n_time, n_time)
diag(Sigma_in) <- 1

# Simulate data
library(MASS)
sim_cs  <- MASS::mvrnorm(n_subj, mu = rep(0, n_time), Sigma = Sigma_cs)
sim_ar1 <- MASS::mvrnorm(n_subj, mu = rep(0, n_time), Sigma = Sigma_ar1)
sim_in  <- MASS::mvrnorm(n_subj, mu = rep(0, n_time), Sigma = Sigma_in)

# Put in long format
to_long <- function(mat, type){
  data.frame(id = rep(1:nrow(mat), each = n_time),
             time = rep(time, nrow(mat)),
             y = as.vector(t(mat)),
             cov_type = type)
}
df <- rbind(to_long(sim_cs, "CS"), to_long(sim_ar1, "AR1"), to_long(sim_in, "Indep"))

Visualization: CS vs AR(1) vs Independence

library(ggplot2)
ggplot(df, aes(x = time, y = y, group = id, color = cov_type)) +
  geom_line(alpha = 0.75) +
  facet_wrap(~cov_type) +
  labs(title = "Simulated Trajectories under CS vs AR(1) vs Identity Covariance",
       x = "Time", y = "Response") +
  theme_minimal()

Key Takeaways from Simulation

If mis-specified, standard errors are incorrect and inference is invalid (point estimates stay consistent)

ML Estimation Goal

Goal: Jointly estimate both:

• Regression coefficients \(\beta\)
• Covariance parameters \(\boldsymbol{\theta}\)

The Log-Likelihood

Under the multivariate normal model:

\[ \ell(\beta, \boldsymbol{\theta}) = -\tfrac{1}{2} \sum_{i=1}^N \Big[ \log |\Sigma_i| + (\mathbf{Y}_i - \mathbf{X}_i \beta)^\top \Sigma_i^{-1} (\mathbf{Y}_i - \mathbf{X}_i \beta) \Big] \]

Case 1: Independent Observations

When responses are independent (e.g., cross-sectional data):

\[\Sigma_i = \sigma^2 I\]

Log-likelihood simplifies to OLS framework:

\[ \hat{\beta} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{Y} \]

Visualization: Independent Case

set.seed(667)
n <- 50
x <- rnorm(n)
y <- 1 + 2*x + rnorm(n, sd = 1)
fit <- lm(y ~ x)
plot(x, y, pch = 19, col = "steelblue")
abline(fit, col = "darkred", lwd = 2)

Independence: Key Point

With independence:

• OLS line = ML estimate
• All observations contribute equally

Case 2: Non-Independent Observations

For longitudinal data: responses within subject \(i\) are correlated

• \(\Sigma_i\) not proportional to identity matrix
• Need to model covariance explicitly (e.g., CS, AR(1), UN)
• ML must jointly estimate \(\beta\) and \(\boldsymbol{\theta}\)

Key Result: GLS/ML Estimator of \(\beta\)

For a known (or estimated) covariance \(\Sigma_i\), the ML estimator of \(\beta\) is the generalized least squares (GLS) form:

\[ \hat{\beta} = \left( \sum_{i=1}^N \mathbf{X}_i^\top \Sigma_i^{-1} \mathbf{X}_i \right)^{-1} \left( \sum_{i=1}^N \mathbf{X}_i^\top \Sigma_i^{-1} \mathbf{Y}_i \right) \]

Setting \(\Sigma_i = \sigma^2 I\) (independence) drops \(\Sigma_i^{-1}\) to a constant and recovers the OLS form \(\hat{\beta} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{Y}\).

Implications of Non-Independence

If correlation is ignored:

Proper ML uses \(\Sigma_i^{-1}\) weighting to account for correlation

Example: CS vs AR(1) Model Comparison

library(nlme)
data(Orthodont, package = "nlme")
Orthodont <- data.frame(Orthodont)

# Fit ML with compound symmetry vs AR1
fit_cs  <- gls(distance ~ age,
               correlation = corCompSymm(value = 0.5, form = ~1|Subject),
               data = Orthodont, method = "ML")

fit_ar1 <- gls(distance ~ age,
               correlation = corAR1(value = 0.5, form = ~age|Subject),
               data = Orthodont, method = "ML")

AIC(fit_cs, fit_ar1)

Model Comparison: Interpretation

• Compare model fits using AIC
• AR(1) often better when correlation decays with time
• Demonstrates sensitivity of ML to covariance choice

In plain language: When you see an AIC comparison, the model with the lower AIC provides a better balance of fit and parsimony. If AR(1) has lower AIC than CS, this suggests that correlations between measurements do indeed decay as time between measurements increases - a common pattern in longitudinal data.

SAS Equivalent: PROC MIXED (Reference Only)

The same CS vs AR(1) comparison in SAS uses PROC MIXED with a REPEATED statement:

/* Compound symmetry */
proc mixed data=orthodont method=reml;
  class subject;
  model distance = age / solution;
  repeated / type=cs subject=subject;
run;

/* First-order autoregressive */
proc mixed data=orthodont method=reml;
  class subject;
  model distance = age / solution;
  repeated / type=ar(1) subject=subject;
run;

Why ML Matters

ML provides a unified framework for:

Caveat: Can be biased for variance components in small samples - motivates REML

Discussion Prompt

Suppose you analyze repeated measures with independence assumed, but the true covariance is AR(1).

• What happens to \(\hat{\beta}\)?
• What happens to its standard errors?
• How might conclusions differ?

Day 2 Roadmap

Day 1 covered covariance structure and ML estimation. Today:

• Missing data mechanisms (MCAR, MAR, MNAR)
• Statistical inference (Wald, LRT, score tests)
• Restricted maximum likelihood (REML)

The Missing Data Problem

Missing data is ubiquitous in longitudinal studies

The mechanism driving missingness is critical to:

• Validity of parameter estimates
• Interpretation of results

Rubin’s Classification

Core question: Can we ignore the missingness mechanism?

MCAR: Missing Completely at Random

Probability of missingness does not depend on data (neither observed nor unobserved)

Let \(R_i\) be the missingness indicator for subject \(i\) (which responses are observed vs missing). Note this \(R_i\) is the missing-data indicator and is distinct from the residual-covariance matrix \(R_i\) in the master notation.

\[ P(R_i \mid Y_i, X_i) = P(R_i) \]

Example: Lab sample lost due to shipping error

Implication: Complete case analyses remain unbiased, but less efficient

Rare in practice, but a useful baseline

MAR: Missing at Random

Probability of missingness depends only on observed data

Split \(Y_i\) into its observed and missing parts: \(Y_{i,\text{obs}}\) is the responses actually recorded for subject \(i\), and \(Y_{i,\text{mis}}\) is the responses that are missing.

\[ P(R_i \mid Y_i, X_i) = P(R_i \mid Y_{i,\text{obs}}, X_i) \]

Example: Dropout depends on age or baseline health status (observed)

Implication: Valid inference possible if model conditions on observed predictors

MAR: Methods for Valid Inference

MAR: A Critical Caveat

MAR does not mean automatically safe

Valid inference requires that the variables driving missingness are part of your analysis (or imputation) model

Otherwise: bias creeps in, even if the true mechanism is MAR

MAR Examples: Correct vs Incorrect

Takeaway: Always include predictors of missingness in your model if they are observed

MNAR: Missing Not at Random

Probability of missingness depends on unobserved data

\[ P(R_i \mid Y_i, X_i) = P(R_i \mid Y_{i,\text{mis}}, Y_{i,\text{obs}}, X_i) \]

Example: Patients with worsening symptoms more likely to drop out

MNAR: Implications

Ignorability fails; specialized models needed:

Assumptions must be explicit and untestable

Visual Illustration of Mechanisms

library(ggplot2)
set.seed(667)

n <- 200
x <- rnorm(n)
y <- 0.5 + 1.2*x + rnorm(n)

# Create missingness
mar  <- ifelse(x > 0.5, NA, y)           # MAR: depends on observed X
mcar <- ifelse(runif(n) < 0.3, NA, y)    # MCAR: random
mnar <- ifelse(y > 1, NA, y)             # MNAR: depends on unobserved Y

df <- data.frame(x, y, mcar, mar, mnar)
df_long <- tidyr::pivot_longer(df, cols = c(y, mcar, mar, mnar),
                               names_to = "Type", values_to = "Value")

ggplot(df_long, aes(x, Value, color = Type)) +
  geom_point(alpha = 0.6) +
  facet_wrap(~Type) +
  theme_minimal() +
  labs(title = "Illustrating MCAR, MAR, and MNAR")

Consequences of Ignoring Mechanisms

Incorrect assumption can change clinical conclusions

Discussion Prompt

Suppose a trial of a new cancer therapy has 20% dropout:

• Healthier patients drop due to travel burden
• Sicker patients drop due to toxicity

• Which mechanism is most plausible?
• How would this affect analysis?

Practical Takeaways for Missing Data

Goals of Inference

After estimation, we want to make inferences about:

• Regression coefficients \(\beta\)
• Covariance parameters \(\boldsymbol{\theta}\)

Why Longitudinal Inference Differs

Inference in longitudinal models differs from OLS because:

Standard OLS Inference (Review)

In cross-sectional linear regression:

\[\hat{\beta} \sim N(\beta, \sigma^2 (\mathbf{X}^\top \mathbf{X})^{-1})\]

Assumes:

• Independent, identically distributed errors
• Homoscedasticity

Inference with Correlated Data

With correlated data:

\[\text{Var}(\hat{\beta}) = \left( \sum_i \mathbf{X}_i^\top \Sigma_i^{-1} \mathbf{X}_i \right)^{-1}\]

Key difference: Dependence on \(\Sigma_i\)

If \(\Sigma_i\) mis-specified: incorrect standard errors, invalid inference (the point estimate \(\hat{\beta}\) stays consistent)

Three Approaches to Hypothesis Testing

Likelihood Ratio Tests (LRT)

For testing hypotheses about parameters in nested models:

\[\Lambda = -2(\ell_{\text{restricted}} - \ell_{\text{full}})\]

Asymptotically \(\chi^2\)-distributed

Used for testing fixed effects and comparing covariance structures

LRT: Pros and Cons

Wald Tests

Based on asymptotic normality of \(\hat{\beta}\):

\[W = (\hat{\beta} - \beta_0)^\top \text{Var}(\hat{\beta})^{-1} (\hat{\beta} - \beta_0)\]

Score Tests

Based on derivative of likelihood at restricted model

Comparison of Inference Approaches

The two tests of FLW Ch. 4 are the LRT and the Wald test. The score test (last row) is an instructor extension, not Ch. 4 content.

When to Use Which Test?

LRT and Wald are the FLW Ch. 4 tests; the score-test row is the instructor extension.

Satterthwaite Degrees of Freedom

In finite samples, \(\chi^2\) and \(t\) approximations may be too liberal

Satterthwaite adjustment: estimates a fractional degrees of freedom from the variance of the estimated variance components, so the \(t\) reference distribution matches the small-sample behavior of \(\hat{\beta} / \text{SE}(\hat{\beta})\). lmerTest computes it automatically; you do not specify the formula by hand.

Satterthwaite: Key Features

Inference Summary

The FLW Ch. 4 toolkit is the Wald test and the LRT; the score test is an instructor extension.

Motivation for REML

Estimating variance components by ML can be biased

Because ML does not account for loss of degrees of freedom from estimating \(\beta\)

The REML Solution

REML modifies the likelihood to use only error contrasts of the data

This eliminates \(\beta\) from the variance estimation

Mathematical Definition: ML

Suppose \(\mathbf{Y} \sim N(\mathbf{X}\beta, \boldsymbol\Sigma)\)

ML log-likelihood:

\[ \ell(\beta, \boldsymbol{\theta}) = -\tfrac{1}{2}\left\{\log|\boldsymbol\Sigma| + (\mathbf{Y}-\mathbf{X}\beta)^\top \boldsymbol\Sigma^{-1}(\mathbf{Y}-\mathbf{X}\beta) \right\} \]

The OLS Intuition: \(n\) vs \(n-p\)

Start with the familiar OLS case (here \(n\) is the generic OLS sample size, the number of data points, not the longitudinal occasion count):

Dividing by \(n\) ignores that we spent \(p\) degrees of freedom estimating \(\beta\). The \(-p\) correction is what makes the estimate unbiased. REML carries this same correction into the general covariance setting.

Mathematical Definition: REML

REML log-likelihood integrates out \(\beta\):

\[ \ell_R(\boldsymbol{\theta}) = -\tfrac{1}{2}\left\{ \log|\boldsymbol\Sigma| + \underbrace{\log|\mathbf{X}^\top \boldsymbol\Sigma^{-1}\mathbf{X}|}_{\text{multivariate analogue of the } -p \text{ correction}} + (\mathbf{Y}-\mathbf{X}\hat{\beta})^\top \boldsymbol\Sigma^{-1}(\mathbf{Y}-\mathbf{X}\hat{\beta}) \right\} \]

The \(\log|\mathbf{X}^\top \boldsymbol\Sigma^{-1}\mathbf{X}|\) term is the matrix generalization of subtracting \(p\) from the denominator: it debits the likelihood for the \(p\) fixed effects that were estimated.

Intuition: Why REML Works

Advantages of REML

Disadvantages of REML

ML vs REML: Decision Guide

Visualization: Bias in Variance Estimates

set.seed(667)
library(nlme)

# Simulate small dataset
dat <- data.frame(id = rep(1:20, each = 3),
                  time = rep(1:3, times = 20))
# rnorm(20) draws one random intercept per subject; [dat$id] then repeats
# each subject's intercept across that subject's 3 rows (a subject-level effect).
dat$y <- 1 + 0.5*dat$time + rnorm(60, sd = 1) + rnorm(20)[dat$id]

# Fit ML and REML
fit_ml   <- lme(y ~ time, random = ~1|id, data = dat, method = "ML")
fit_reml <- lme(y ~ time, random = ~1|id, data = dat, method = "REML")

c("ML Var(Intercept)"   = VarCorr(fit_ml)[1, "Variance"],
  "REML Var(Intercept)" = VarCorr(fit_reml)[1, "Variance"])

ML vs REML: Interpreting the Output

In plain language: Notice that REML typically produces a larger variance estimate than ML. This is because REML corrects for the downward bias inherent in ML estimation. In small samples, this difference can be substantial. The REML estimate is analogous to using \(n-p\) instead of \(n\) in the denominator when estimating residual variance in ordinary regression.

REML Takeaways

Common Mistakes to Avoid

• Using REML for fixed effects LRT: Use ML when comparing models with different fixed effects; REML likelihoods are not comparable
• Ignoring correlation structure entirely: Even though \(\hat{\beta}\) stays consistent, standard errors are incorrect and inference is invalid - potentially leading to false conclusions
• Assuming MAR when dropout depends on unobserved outcomes: If sicker patients drop out due to worsening symptoms (unobserved), treating this as MAR leads to bias
• Not including predictors of missingness in your model: Under MAR, validity requires conditioning on observed variables that predict missingness
• Using overly complex covariance structures with small samples: Unstructured covariance with \(n_i(n_i+1)/2\) parameters can fail to converge or produce unstable estimates

Lecture Summary