BIOS 667 - Lecture 10: Residual Analyses and Diagnostics (Ch. 10)

Lecture Objectives

By the end of this lecture, you will be able to:

• Compute raw and transformed (Cholesky-standardized) residuals and explain why the transform decorrelates them
• Use the Mahalanobis distance \(d_i = r_i^{*\top} r_i^{*}\) to flag outlying individuals
• Use the semi-variogram of transformed residuals to diagnose covariance misspecification
• Use random-effect (empirical BLUP) plots to spot individuals with unusual profiles (NOT to test RE normality)
• Read each diagnostic plot: what a good vs problematic pattern looks like

Prerequisites / Recall

This lecture builds directly on two earlier ones.

Roadmap

Notation reference

Notation (this lecture)

Part 1: Raw, Transformed, and Outlier Diagnostics

Why Diagnostics for LMMs?

Linear mixed-effects models rest on assumptions more complex than OLS.

Raw Residuals (FLW eq 10.1)

The raw (marginal) residual for subject \(i\) is

\[ r_i \;=\; \mathbf{Y}_i - \mathbf{X}_i\hat{\boldsymbol\beta}. \tag{10.1}\]

• One residual vector per subject, length \(n_i\).
• If the mean model \(\mathbf{X}_i\boldsymbol\beta\) is correct, these are centered at zero.
• Problem: raw residuals are correlated and have unequal variances (they inherit \(\Sigma_i\)), so ordinary regression diagnostics do not apply directly.

Transformed Residuals (FLW eq 10.2)

Take the Cholesky factorization of the estimated marginal covariance, \(\hat\Sigma_i = L_i L_i^\top\) (\(L_i\) lower-triangular). Define

\[ r_i^{*} \;=\; L_i^{-1} r_i \;=\; L_i^{-1}\!\left(\mathbf{Y}_i - \mathbf{X}_i\hat{\boldsymbol\beta}\right). \tag{10.2}\]

• First element = standardized baseline residual.
• \(k\)-th element estimates \(\dfrac{Y_{ik}-E(Y_{ik}\mid Y_{i1},\dots,Y_{i,k-1})}{\sqrt{\operatorname{Var}(Y_{ik}\mid Y_{i1},\dots,Y_{i,k-1})}}\).

Mahalanobis Distance for Outlying Individuals (FLW eq 10.3)

A single summary of how far subject \(i\)’s whole profile is from its fitted values:

\[ d_i \;=\; r_i^{*\top}\, r_i^{*}. \tag{10.3}\]

• p-values give a common metric when \(n_i\) varies across subjects.
• Caution (FLW p.268): because we focus on the most extreme \(d_i\), the true tail behaves worse than a plain \(\chi^2_{n_i}\); a Bonferroni correction (\(p \times N\)) is available but very conservative. Expect about \(0.05N\) “significant” \(d_i\) by chance.

Part 2: The Semi-Variogram

Variogram: FLW Definition (eq 10.4)

The semi-variogram is one-half the expected squared difference between two residuals on the same individual, at time-lag \(h_{ijk}\):

\[ \gamma(h_{ijk}) \;=\; \tfrac12\,\mathbb{E}\big[(r_{ij}-r_{ik})^2\big]. \tag{10.4}\]

Since residuals have mean zero,

\[ \gamma(h_{ijk}) \;=\; \tfrac12\operatorname{Var}(r_{ij}) + \tfrac12\operatorname{Var}(r_{ik}) - \operatorname{Cov}(r_{ij},r_{ik}). \]

Variogram: the Geostatistics Motivation (flagged)

How to Read a Variogram

Part 3: Case Study (Real FLW Data)

Case Study: Dental Growth Data

• 27 children (11 girls, 16 boys), distance (mm) from pituitary to pteryomaxillary fissure at ages 8, 10, 12, 14.
• Fit the Ch. 8 LME, then: transformed residuals -> Mahalanobis outlier scan -> semi-variogram.

Case Study: Load and Reshape

suppressPackageStartupMessages({
  library(nlme); library(ggplot2); library(dplyr)
})

# Download fallback (FLW site); local copy preferred
data_file <- "../../data/dental.txt"
if (!file.exists(data_file)) {
  dir.create("../../data", showWarnings = FALSE, recursive = TRUE)
  download.file("https://content.sph.harvard.edu/fitzmaur/ala2e/dental.txt",
                data_file)
}

# FLW dental.txt has a comment-header block and is WIDE (one row per child)
raw <- readLines(data_file)
dat_lines <- grep("^[[:space:]]*[0-9]+[[:space:]]+[FM]", raw, value = TRUE)
wide <- read.table(text = dat_lines, header = FALSE,
                   col.names = c("id", "gender", "age8", "age10", "age12", "age14"))

# reshape WIDE -> LONG
long <- reshape(wide, varying = c("age8", "age10", "age12", "age14"),
                v.names = "distance", timevar = "age",
                times = c(8, 10, 12, 14), idvar = "id", direction = "long")
long <- long[order(long$id, long$age), ]
long$id    <- factor(long$id)
long$age_c <- long$age - 8          # center at first occasion
head(long, 6)

Case Study: Fit the Ch. 8 LME

fit <- lme(distance ~ age_c * gender,
           random = ~ 1 + age_c | id,
           data = long, method = "REML")
round(fixef(fit), 3)

VarCorr(fit)

Case Study: Transformed Residuals (eq 10.2)

uids   <- levels(long$id)
r_marg <- resid(fit, level = 0)          # RAW residual Y - X*beta (level=0)
long$r_marg <- r_marg

long$r_star <- NA_real_
d_i <- setNames(numeric(length(uids)), uids)   # Mahalanobis distances

for (u in uids) {
  idx   <- which(long$id == u)
  Sig_i <- getVarCov(fit, individual = u, type = "marginal")[[1]]  # Sigma_i hat
  L_i   <- t(chol(Sig_i))                # lower-triangular: Sigma_i = L_i L_i'
  r_i   <- long$r_marg[idx]
  rstar <- solve(L_i, r_i)               # r* = L_i^{-1} r_i   (eq 10.2)
  long$r_star[idx] <- rstar
  d_i[u] <- sum(rstar^2)                 # d_i = r*' r*        (eq 10.3)
}

c(var_transformed = round(var(long$r_star), 3))   # should be ~1

Case Study: QQ Plot of Transformed Residuals

Case Study: Mahalanobis Outlier Scan (eq 10.3)

Case Study: Semi-Variogram of Transformed Residuals (eq 10.4)

Case Study: Does Dropping the Random Slope Change the Variogram? (dental)

Case Study: Standalone SAS (PROC MIXED)

/* Dental growth LME with residual + influence diagnostics.
   Data already reshaped to LONG: id gender age_c distance.            */
proc mixed data=dental_long method=reml;
    class id gender;
    model distance = age_c gender age_c*gender /
          solution outpm=marg_resid outp=cond_resid
          residual influence(effect=id);
    random intercept age_c / subject=id type=un g;
run;

/* outpm= : marginal residuals (Y - X*beta)         -> variogram input
   outp=  : conditional residuals (Y - X*beta - Z*b)
   influence(effect=id) : leave-one-cluster-out Cook's D / restricted
                          likelihood distance per child (Level-2)       */

Part 4: Software Diagnostics (nlme / lme4)

Two Residual Types in Software

FLW’s split is raw vs transformed. Software (nlme, lme4) instead exposes two levels of residual. Keep the vocabulary distinct:

Conditional Residuals vs Fitted

Conditional Residual ACF (beyond FLW)

Random-Effect (Empirical BLUP) Plots: Correct Use

Part 5: Beyond-FLW Diagnostics (flagged)

Cluster Influence: Leave-One-Out Refit (beyond FLW)

# REAL leave-one-cluster-out refit (N=27 children -> fast, runs in full)
beta_full <- fixef(fit)
ctrl  <- lmeControl(maxIter = 200, msMaxIter = 200, opt = "optim")
delta <- setNames(rep(NA_real_, length(uids)), uids)
for (u in uids) {
  sub <- long[long$id != u, ]
  rf  <- tryCatch(
    lme(distance ~ age_c * gender, random = ~ 1 + age_c | id,
        data = sub, method = "REML", control = ctrl),
    error = function(e) NULL)
  if (!is.null(rf)) delta[u] <- max(abs(fixef(rf) - beta_full))  # max |change|
}
delta <- delta[!is.na(delta)]   # drop any child whose deletion-refit failed
plot(seq_along(delta), delta, type = "h", lwd = 2,
     xlab = "Child (id) dropped", ylab = "max |delta beta|",
     main = "Leave-one-child-out: max |change| in fixed effects")
text(which.max(delta), max(delta), labels = names(which.max(delta)), pos = 3, col = "red")

CUSUM (recursive): beyond FLW (flagged)

suppressPackageStartupMessages(library(strucchange))
# Recursive CUSUM is an OLS/structural-change construct -> run on an OLS fit
fm <- efp(distance ~ age_c + gender, type = "Rec-CUSUM", data = long)
plot(fm, main = "Recursive CUSUM (OLS mean model) - beyond-FLW demo")

What FLW Actually Does: Aggregated Sums (Section 10.4)

FLW Section 10.4 (“aggregating residuals”) uses the Lin, Wei & Ying (2002) cumulative / moving sums of raw residuals over a covariate \(x_{ijk}\) or fitted value:

\[ W_k(x) = \frac{1}{\sqrt{N}}\sum_{i=1}^{N}\sum_{j=1}^{n_i} I(x_{ijk}\le x)\, r_{ij}. \]

• Under a correct mean model, \(W_k(\cdot)\) behaves like a zero-mean Gaussian process.
• Assess by simulating many null realizations and comparing the observed \(\max|W_k|\) to the simulated maxima (a supremum test, p-value = proportion exceeding).
• Robust: valid regardless of the true covariance among responses.

Part 6: Checklist and Pitfalls

Practical Diagnostic Checklist (FLW Ch. 10)

RI vs RS: Comparison (boundary caveat)

Common Mistakes to Avoid

Diagnostic Decision Tree

Diagnostic Principles Across Models

Summary

• Raw residual (eq 10.1) \(r_i=\mathbf Y_i-\mathbf X_i\hat{\boldsymbol\beta}\) - correlated, unequal variance.
• Transformed residual (eq 10.2) \(r_i^{*}=L_i^{-1}r_i\) - uncorrelated, unit variance; the chapter’s defining tool.
• Mahalanobis (eq 10.3) \(d_i=r_i^{*\top}r_i^{*}\sim\chi^2_{n_i}\) flags outlying individuals.
• Semi-variogram (eq 10.4) of transformed residuals centers at 1 under a correct covariance model.
• Empirical-BLUP plots spot unusual individuals; they do NOT test RE normality (FLW p.273).
• ACF, cluster-influence, CUSUM are useful but beyond FLW Ch. 10 (flagged).

Further Reading

• Fitzmaurice, Laird & Ware (2011), Applied Longitudinal Analysis, 2nd ed., Chapter 10 (Residual Analyses and Diagnostics).
• Lin, Wei & Ying (2002), “Model-checking techniques based on cumulative residuals,” Biometrics 58, 1-12 (FLW’s aggregated-residuals method).
• Pinheiro & Bates (2000), Mixed-Effects Models in S and S-PLUS (nlme residual / Variogram / ACF tooling).
• Verbeke & Molenberghs (2000) and Diggle et al. (2002) for the geostatistics variogram (nugget / sill).