Module 2 · Joint modeling and missing data

Lectures 10.1–10.3. Reading: Lecture Notes Part 2, §2.4–2.5. Code: R script — Lecture Video 39 (joint model), building on the aids class activity at the end of R script — Lecture Video 37.

← Back to notes index

Two closing topics. The first brings Part 1 and Part 2 together: a study with both a longitudinal biomarker and a time-to-event outcome, where fitting the Cox model and the mixed model separately is not quite right. The second is the problem that was quietly present in every real longitudinal dataset so far — missing data — and which the aids activity in Notes 6 already showed can change estimates.

Lecture 10.1 — Joint modeling

The setting

The aids study: 467 HIV patients who had failed or were intolerant of AZT, randomized to two antiretrovirals, ddI (didanosine) and ddC (zalcitabine). Two outcomes, both of interest in their own right:

  1. time to death — Figure 2.17, the Kaplan–Meier curves by drug;
  2. CD4 cell count over time — Figure 2.18, the spaghetti plot.

Two research questions:

  • How strong is the association between CD4 count and the risk of death?
  • Is CD4 count a good biomarker? — if treatment improves CD4, does it also improve survival? (A harder question than it looks; see the end of the R section.)

Why not two separate models?

The straightforward approach is one you already know how to do:

  • a Cox model for time to death with CD4 as a time-dependent covariate (the aids data even arrive in (start, stop, event) form, ready for it);
  • a linear mixed model (or GEE) for the CD4 trajectory.

You will see this done all the time. Three subtleties make it not quite right:

  1. The longitudinal measurements exist only while the person is alive. After death, no more CD4. For the Cox model that is fine — it is the nature of a time-varying covariate. For the mixed model it is a problem: dropout here is not someone moving away, it is (probably) death, one of the outcomes of interest. Dropout is not at all random. Contrast a longitudinal measurement external to the person, such as air pollution, where this issue does not arise.
  2. The questions imply a joint process. We talk about the two outcomes as related, then fit them as if they were not. Separate models do not match the conceptual framework.
  3. Biomarkers are measured with error. Biologically, it is a person’s true CD4 level at time $t$ that drives their risk, but what we observe is the true level plus measurement error. The survival model should depend on the true value, not the mis-measured one. Much of the joint-modeling literature grew out of exactly this.

The joint model in three steps

Step 1 — survival, driven by the true biomarker. Let $M_i(t)$ be the true, unobserved longitudinal value at time $t$ and $\mathcal{M}_i(t) = {M_i(s),\ 0 \le s < t}$ the whole true history up to $t$. Specify

\[\lambda(t \mid Z, \mathcal{M}_i(t)) = \lambda_0(t)\exp\{\beta^T Z + \alpha\, M_i(t)\}\]
Symbol Meaning
$Z$ baseline covariates (drug, …), with coefficients $\beta$ as in any Cox model
$M_i(t)$ the true biomarker — a time-varying covariate, except that it is not observed
$\alpha$ the strength of the association between the biomarker and the risk of the event — the parameter research question 1 is about

Structurally this is the time-varying Cox model from §1.9; the difference is that the covariate is latent.

Step 2 — the observed biomarker, with error. Let $Y_i(t)$ be the observed value:

\[Y_i(t) = M_i(t) + \epsilon_i^*(t), \qquad M_i(t) = \theta_0 + X_i(t)\theta_1 + b_{0i} + b_{1i}X_i(t)\]

A linear mixed model with random intercept $b_{0i}$ and random slope $b_{1i}$ — but for the true $M_i(t)$, with $\epsilon_i^*(t)$ now explicitly measurement error. (The coefficients are $\theta$ rather than $\beta$ only because $\beta$ is taken by step 1.) The covariates $X_i(t)$ may be the same as $Z$, overlap with it, or differ: something can be associated with death but not CD4, or vice versa.

Step 3 — tie them together. The hard part. The joint model has the general form

\[p(Y_i(t), T_i, \delta_i) = \int p(Y_i(t) \mid b_i)\;\big\{\lambda(T_i \mid b_i)^{\delta_i}\, S(T_i \mid b_i)\big\}\; p(b_i)\, db_i\]

with $b_i = (b_{0i}, b_{1i})$. Read it left to right: the density of the observed longitudinal data given the random effects; the survival contribution $\lambda^{\delta} S$ from §1.3’s censored likelihood, also given the random effects; the distribution of the random effects; integrated over them. The random effects are what explain the interdependence between the two processes — that is the entire modeling assumption.

Practicalities:

  • For $S(T_i \mid b_i)$ the advice is a parametric but flexible model, e.g. splines for the baseline hazard. A fully unspecified Cox baseline sometimes fails.
  • Estimation is maximum likelihood (the pieces are parametric now) via an iterative procedure, usually the EM algorithm; or Bayesian.
  • The integral is the difficulty. It must actually be computed numerically, and joint fits sometimes fail because the integration does not converge.
  • R: package JM (and JMbayes for the Bayesian route).

References: Rizopoulos, Joint Modeling of Longitudinal and Survival Data in R; Tsiatis & Davidian (2004), Statistica Sinica.

Lecture 10.2 — Joint modeling in R

From R script — Lecture Video 39.

Two versions of the data ship with JM:

Object Rows Use
aids 1405 — one per CD4 measurement, with start/stop/event the longitudinal model
aids.id 467 — one per patient, with Time and death the survival model

The two component models, fit first on their own:

lmeFit <- lme(CD4 ~ obstime*drug, random = ~ obstime | patient, data = aids)
coxFit <- coxph(Surv(Time, death) ~ drug, data = aids.id, x = TRUE)

lme (package nlme) is a different mixed-model function from lmer — the one JM wants. The fixed-effects formula looks the same; the random effects go in a separate argument, random = ~ obstime | patient for intercepts and slopes (~ 1 | patient for intercepts only). In the Cox model only drug appears — CD4 is deliberately left out, because the joint model will supply it as $M_i(t)$ — and x = TRUE tells coxph to keep the design matrix in the object, which jointModel needs. Alone, the Cox model shows no drug effect (HR $= 1.23$, $p = 0.15$).

The joint model:

jointFit <- jointModel(lmeFit, coxFit, timeVar = "obstime", method = "spline-PH-GH")
summary(jointFit)

The two fits go in, timeVar names the time variable it should find in both, and method has three parts:

Part Options Meaning
baseline piecewise, spline, weibull, Cox how $\lambda_0(t)$ is handled: piecewise-constant, B-spline approximation of the log baseline hazard, Weibull, or unspecified
model PH, AFT proportional hazards or accelerated failure time
integration GH Gauss–Hermite quadrature for the integral over $b_i$

So "spline-PH-GH" is a PH model with a spline-approximated baseline hazard — parametric but flexible, step 3’s advice. "Cox-PH-GH" feels most familiar but the integration sometimes fails; "weibull-PH-GH" is more restrictive. The fit can take a while (the script pastes its output in case). Reading the summary top to bottom:

  • Data descriptives — 1405 observations in 467 groups; 188 events (40.3%).
  • Variance components — SDs 4.52 (intercept), 0.17 (slope), corr $-0.05$, residual 1.88. The random-effects block, as in any mixed model.

  • Longitudinal process — the CD4 model within the joint framework:

      Value SE $p$
    (Intercept) 7.071 0.173 $< 0.0001$
    obstime $-0.181$ 0.022 $< 0.0001$
    drugddI 0.302 0.270 0.26
    obstime:drugddI 0.004 0.031 0.90

    CD4 falls over time; no evidence of a treatment effect on CD4.

  • Event process — the survival model within the joint framework:

      Value SE $z$ $p$
    drugddI 0.356 0.158 2.25 0.024
    Assoct $-0.300$ 0.039 $-7.79$ $< 0.0001$
    bs1 … bs9       spline basis for the baseline hazard

    Assoct is $\hat\alpha$ — the answer to research question 1. Each unit higher true CD4 multiplies the hazard of death by $e^{-0.30} = 0.74$; there is a strong association between CD4 and survival. And there is evidence of a treatment effect on survival (ddI worse, HR $= 1.43$), which the Cox model alone did not show.

  • Integration / optimization — Gauss–Hermite with 15 quadrature points; Convergence: 0 means it worked.

Separate vs. joint. The separate time-varying Cox model with observed CD4 gives $-0.19$ for the CD4 coefficient; the joint model’s $\alpha = -0.30$ is larger in magnitude — the usual pattern once measurement error in the biomarker is accounted for, since error attenuates the naive coefficient. The conclusions are the same either way; the estimates are not.

On research question 2. “Is CD4 a good biomarker?” is not really answered by the joint model, and the answer depends on what is meant. If it means prediction — can CD4 trajectories predict death? — this is not a prediction tool. If it means surrogacy — can a treatment effect on CD4 stand in for the treatment effect on survival, so trials need not wait for deaths? — that is the surrogate-marker literature, a different framework. The joint model gives a rigorous measure of association between a trajectory and a time-to-event outcome; no more.

Lecture 10.3 — Missing data

With longitudinal data, missing data are the rule, not the exception: people drop out, withdraw, move, or just miss a visit. In the survival setting censoring is itself a kind of missingness, and covariates can be missing too — the standard Cox software silently does a complete-case analysis, dropping anyone missing any covariate.

Two implications:

  1. Loss of efficiency — less data, less precision, less power. Unavoidable.
  2. Bias — misleading inference about the mean response. This is the serious one, and whether it happens depends on why the data are missing: the missing data mechanism.

Notation: $Y_i$ and $R_i$

Alongside the response vector $\mathbf{Y}i = (Y{i1}, \dots, Y_{in})^T$ define the response indicator vector $\mathbf{R}i = (R{i1}, \dots, R_{in})^T$, with $R_{ij} = 1$ if $Y_{ij}$ is observed and $0$ if missing. Four kinds of person, five time points:

  $\mathbf{Y}_i$ $\mathbf{R}_i$ Pattern
1 50, 51, 55, 40, 42 1 1 1 1 1 completer
2 42, 56, 51, NA, NA 1 1 1 0 0 dropout
3 NA, NA, 40, 92, 70 0 0 1 1 1 late entry
4 39, NA, 54, NA, 49 1 0 1 0 1 intermittent

(First question when you receive data: how is missingness coded? NA, ., 99, 999 … A 99 code is a trap once a later variable can genuinely take the value 99.)

Split each person’s vector into $\mathbf{Y}_i^O$, the observed components, and $\mathbf{Y}_i^M$, the missing ones — the latter hypothetical, since by definition we do not have it, but needed to state the mechanisms.

Three mechanisms

Mechanism Definition Observed data a random sample of the complete data? Which methods are valid
MCAR — missing completely at random $P(\mathbf{R}_i \mid \mathbf{Y}_i^O, \mathbf{Y}_i^M) = P(\mathbf{R}_i)$ yes anything valid for complete data; just less efficient
MAR — missing at random $P(\mathbf{R}_i \mid \mathbf{Y}_i^O, \mathbf{Y}_i^M) = P(\mathbf{R}_i \mid \mathbf{Y}_i^O)$ no mixed models with a correctly specified correlation structure; not GEE (needs a correction); not mixed models with a mis-specified structure
MNAR — missing not at random depends on $\mathbf{Y}_i^M$ (and possibly $\mathbf{Y}_i^O$) — cannot be simplified no essentially nothing, without modeling the joint distribution of $(\mathbf{R}_i, \mathbf{Y}_i^O, \mathbf{Y}_i^M)$

Examples:

  • MCAR — subjects leave after a pre-determined number of measurements; lab samples lost to equipment or staff error. Rare in practice.
  • MAR — protocol removes patients whose response exceeds a threshold (toxicity); physicians give rescue medication to non-responders. Missingness depends on values you have seen, not on the ones you then miss.
  • MNAR — drug-addiction studies where those who relapse stop reporting; the alcohol-use participant who skips Monday’s visit after a weekend bender; quality-of-life questionnaires left blank precisely when quality of life is poor. Missingness depends on the value that is missing.

What you can and cannot check. With the observed data you can explore MCAR vs. MAR — fit a model for $R_{ij}$ on the observed earlier $Y$’s and see whether missingness is associated with them. You cannot empirically distinguish MAR from MNAR, because the predictor that would reveal it, $\mathbf{Y}_i^M$, is the thing you do not have. That distinction rests on subject-matter argument.

Three ways to handle it

  1. Complete-case analysis — keep only completers. Valid only under MCAR. For longitudinal data it can be brutal (100 people, 10 completers) and is what the Cox software does by default with missing covariates.
  2. Available-data analysis — use everything observed. Mixed models and GEE both do this. More efficient than complete-case; valid under MCAR, and under MAR only for methods whose conditional means and covariances are correctly specified (the mixed-model column in the table above).
  3. Imputation — fill in the missing values. Many ways.

Imputation

Single imputation fills in once, then analyses the completed dataset as if the filled values were real. The problem common to all its variants: it ignores the uncertainty introduced by guessing, so reported standard errors are too small.

  • Last observation carried forward (LOCF) — for dropout, repeat the last value (person 2 becomes 42, 56, 51, 51, 51). Extensively criticized; believed conservative, which it is not; useless for late entry; and nonetheless widely used and accepted by regulators such as the FDA. Variants: baseline observation carried forward (42, 42), worst observation carried forward (56, 56 if higher is worse).
  • Mean imputation — the person’s own mean of observed values, or the mean of everyone’s value at that occasion.

Multiple imputation is widely argued to be the method of choice. Replace the missing values with $m$ plausible draws, giving $m$ completed datasets; run the analysis on each; combine. For a regression parameter, with $\hat\beta^{(k)}$ and $\widehat{\text{Cov}}(\hat\beta^{(k)})$ from the $k$-th dataset:

\[\bar\beta = \frac{1}{m}\sum_{k=1}^{m}\hat\beta^{(k)}\] \[\widehat{\text{Cov}}(\bar\beta) = \underbrace{\frac{1}{m}\sum_{k=1}^{m}\widehat{\text{Cov}}(\hat\beta^{(k)})}_{\text{within-imputation}} \;+\; \Big(1 + \frac{1}{m}\Big)\underbrace{\frac{1}{m-1}\sum_{k=1}^{m}(\hat\beta^{(k)} - \bar\beta)(\hat\beta^{(k)} - \bar\beta)^T}_{\text{between-imputation}}\]

The point estimate is just the average — easy. The variance has two pieces: the average of the usual variances, plus a term measuring how much the $m$ estimates disagree with each other. If the 100 imputed $\hat\beta$’s are all nearly the same, the imputation barely mattered and the second term is small; if they scatter widely, the imputation matters a lot and the variance says so. That second term is exactly what single imputation throws away.

How to generate the $m$ draws is the hard part, and it must involve randomness — LOCF has only one way to be done, so it cannot be “multiple”. The idea is to draw $\mathbf{Y}_i^M$ from the conditional distribution $f(\mathbf{Y}_i^M \mid \mathbf{Y}_i^O, X_i)$. If someone handed you that distribution (say, $N(\text{age} + 8 + \text{female},\ 2)$), you would sample from it $m$ times per missing value and be done.

  • Predictive mean matching builds it from regressions: with values missing only at time 5, fit $E(Y_{i5}) = a_0 + a_1 Y_{i1} + \dots + a_4 Y_{i4} + a_5 X_i$ on the complete cases; then, in a Bayesian spirit, draw a random set of $a$’s from their estimated distribution, plug in the person’s observed values to get an imputed $Y_{i5}$, and repeat $m$ times with fresh draws of the $a$’s.
  • EM algorithm — model the joint distribution $f(\mathbf{Y}_i \mid X_i)$ directly, estimate its parameters iteratively, and impute from the fitted conditional distribution.

None of this is specific to longitudinal data: the three mechanisms and multiple imputation are general tools, used whenever data are missing, inference must be unbiased, and the mechanism is plausibly MAR — and even under MCAR, because complete-case analysis wastes so much.

References: Crippa, A review of Longitudinal Data Analysis in R; Potthoff & Roy (1964); Fitzmaurice, Laird & Ware (2012), Applied Longitudinal Analysis.