4.8 Repeated measures and longitudinal data

When the same units are measured again and again, each unit has its own level and its own path. Models for such data describe the typical path and how units vary around it, and they must respect that measurements of one unit are correlated.

Code: data and style used in this chapter

import warnings

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from matplotlib.axes import Axes

from statly import style

style.use()
rng = np.random.default_rng(1)
warnings.filterwarnings("ignore")

# Fertility (children per woman) by country and year, 1960 to 2011; 40 countries with complete data.
wide = sm.datasets.fertility.load_pandas().data
years = [str(y) for y in range(1960, 2012)]
complete = wide.dropna(subset=years)
chosen = complete.sample(40, random_state=3)
fertility = chosen.melt(
    id_vars="Country Name", value_vars=years, var_name="year", value_name="births"
)
fertility["year"] = fertility.year.astype(int)
fertility["t"] = (fertility.year - 1990) / 10  # decades from 1990
fertility = fertility.rename(columns={"Country Name": "country"})

4.8.1 Paths over time

Longitudinal data follow the same units over time: patients at each visit, pupils each school year, countries each year. The first picture is always the spaghetti plot, one line per unit:

Code

fig, ax = plt.subplots(figsize=(8.5, 3.6))
for _, country in fertility.groupby("country"):
    ax.plot(country.year, country.births, color=style.GREY, lw=0.8, alpha=0.7)
ax.plot(fertility.groupby("year").births.mean(), color=style.CORAL, lw=3)
ax.set(xlabel="year", ylabel="children per woman")
plt.show()

Figure 4.8.1: Fertility (children per woman) in 40 countries, 1960 to 2011, one line per country, with the average across countries (red). Almost every country declines, but they start at different levels, decline at different speeds and at different times.

Two kinds of question arise. How does the typical path look, and how does a predictor change it? And how much do units differ in their level and in their speed? Ordinary regression treats the 2,080 rows as independent, which they are not: a country’s value this year predicts next year’s almost perfectly. Chapter 4.7 introduced mixed models for clustered data; repeated measures are clusters in time.

4.8.2 Growth curves: random intercepts and slopes

A growth curve model is a mixed model in which each unit has its own intercept and its own slope on time, drawn from a common distribution:

$y_{it} = (\beta_0 + b_{0i}) + (\beta_1 + b_{1i})\, t + \varepsilon_{it}, \qquad \begin{pmatrix} b_{0i} \\ b_{1i} \end{pmatrix} \sim N(0, \Sigma)$

$\beta_0$ and $\beta_1$ describe the average path; $\Sigma$ describes how much units vary around it, and whether high starters fall faster (a correlation between $b_0$ and $b_1$ , though a scatter plot of the two shows more than one number can).

Code

growth = smf.mixedlm("births ~ t + I(t**2)", fertility, groups="country", re_formula="~t").fit()
effects = pd.DataFrame(growth.random_effects).T
levels = growth.fe_params["Intercept"] + effects.iloc[:, 0]
speeds = growth.fe_params["t"] + effects.iloc[:, 1]
grid = np.linspace(-3, 2.3, 100)

fig, axes = style.subplots(2, height=3.4, width=10.5)
for name, (b0, b1) in effects.iterrows():
    path = (
        growth.fe_params["Intercept"]
        + b0
        + (growth.fe_params["t"] + b1) * grid
        + growth.fe_params["I(t ** 2)"] * grid**2
    )
    axes[0].plot(1990 + 10 * grid, path, color=style.GREY, lw=0.8, alpha=0.7)
average = (
    growth.fe_params["Intercept"]
    + growth.fe_params["t"] * grid
    + growth.fe_params["I(t ** 2)"] * grid**2
)
axes[0].plot(1990 + 10 * grid, average, color=style.CORAL, lw=3)
axes[0].set(xlabel="year", ylabel="children per woman", title="Fitted paths")
axes[1].scatter(levels, speeds, color=style.BLUE, s=20)
r = np.corrcoef(levels, speeds)[0, 1]
style.note(axes[1], f"correlation {r:.2f}", x=0.55, y=0.95, color=style.INK)
axes[1].set(
    xlabel="level in 1990", ylabel="change per decade", title="Each country's level and speed"
)
fig.tight_layout()
plt.show()

Figure 4.8.2: A growth curve model for fertility with a quadratic average path in time and a random intercept and slope (on time) for each country. Left: the fitted paths of all 40 countries (grey) and the average path (red). Right: each country’s estimated level in 1990 against its rate of change per decade. Countries in the middle, with 3 to 5.5 children per woman in 1990, were falling fastest; those already low and those still high changed slowly. This arch, the demographic transition, is invisible in the correlation of the two random effects, which is close to 0.

The model has a limit worth seeing. It gives every country the same curved average path, shifted and tilted by its own intercept and slope, so each country’s own path is close to a straight line. Real fertility declines are S-shaped: flat, then a fast fall, then flat again, at a different time in each country.

Code

fertility["fitted"] = growth.fittedvalues
fertility["residual"] = fertility.births - fertility.fitted
fertility["size"] = fertility.residual.abs()
examples = (
    fertility.groupby("country")["size"].max().nlargest(4).index
)  # the worst-fitted countries

fig, axes = style.subplots(2, height=3.3, width=10.5)
for name, color in zip(examples, (style.BLUE, style.CORAL, style.GREEN, style.GOLD)):
    country = fertility[fertility.country == name]
    axes[0].scatter(country.year, country.births, s=6, color=color, alpha=0.7)
    axes[0].plot(country.year, country.fitted, color=color, lw=2, label=name)
axes[0].legend(fontsize=7.5)
axes[0].set(xlabel="year", ylabel="children per woman", title="Four countries")
for _, country in fertility.groupby("country"):
    axes[1].plot(country.year, country.residual, color=style.GREY, lw=0.7, alpha=0.6)
axes[1].axhline(0, color=style.INK, lw=0.8)
axes[1].set(xlabel="year", ylabel="observed − fitted", title="Residuals of all countries")
fig.tight_layout()
plt.show()

Figure 4.8.3: Where the growth curve model falls short. Left: four countries’ observed fertility (dots) and their fitted paths (lines); the straight-ish paths cut across the S-shaped declines, too low early and late and too high in the middle of each country’s transition, or the other way round. Right: the residuals of all 40 countries over time; they run in long waves instead of scattering around zero. Country-specific curves, for instance a spline in time (Chapter 6.5) with random coefficients, would follow the paths better.

4.8.3 Why the correlation matters

Ignoring the repeated structure does not usually bias the average path, but it badly understates the uncertainty. With 52 yearly values per country, ordinary regression acts as if it had 2,080 independent observations; in truth there are 40 countries, each a cluster of nearly redundant values.

Code

formula = "births ~ t + I(t**2)"
ols = smf.ols(formula, fertility).fit()
robust = smf.ols(formula, fertility).fit(
    cov_type="cluster", cov_kwds={"groups": pd.factorize(fertility.country)[0]}
)
gee = smf.gee(formula, "country", fertility, cov_struct=sm.cov_struct.Exchangeable()).fit()
ANALYSES = {
    "ordinary least squares": (ols.params["t"], ols.bse["t"]),
    "cluster-robust": (robust.params["t"], robust.bse["t"]),
    "GEE": (gee.params["t"], gee.bse["t"]),
    "growth curve (mixed)": (growth.fe_params["t"], growth.bse_fe["t"]),
}

fig, ax = plt.subplots(figsize=(8, 2.9))
for i, (name, (est, se)) in enumerate(ANALYSES.items()):
    ax.errorbar(
        est,
        i,
        xerr=1.96 * se,
        fmt="o",
        color=style.GREY if i == 0 else style.BLUE,
        capsize=4,
        lw=2,
    )
    ax.text(est + 1.96 * se + 0.02, i, f"± {1.96 * se:.2f}", va="center", fontsize=8.5)
ax.set_yticks(range(len(ANALYSES)), list(ANALYSES))
ax.invert_yaxis()
ax.set(xlabel="change in children per woman per decade (at 1990)")
plt.show()

Figure 4.8.4: The average change in fertility per decade (at 1990), with 95% intervals, from four analyses of the same data. Ordinary least squares treats every country-year as independent and gives an interval half as wide as it should be. Cluster-robust standard errors, generalised estimating equations (GEE, exchangeable correlation) and the growth curve model all account for the 40 countries and agree on a much wider interval.

4.8.4 Two kinds of answer: subject-specific and population-averaged

Mixed models and generalised estimating equations (GEE) answer slightly different questions.

A mixed model describes each unit: “for a given country, fertility falls by so much per decade.” The random effects are part of the model, and individual paths can be predicted.
GEE describes the population average: “across countries, average fertility falls by so much per decade.” It fits only the average path, with a working correlation structure (independent, exchangeable, autoregressive) for the repeated values, and uses robust standard errors, which stay valid even if that structure is wrong.

For a straight-line outcome the two give the same coefficients. For a yes/no or count outcome they differ: a subject-specific odds ratio is further from 1 than the population-averaged one, because averaging S-shaped curves over units flattens them.

4.8.5 The classic: repeated measures ANOVA

In an experiment where every subject is measured under each condition, or at a fixed set of times, the traditional analysis is the repeated measures ANOVA. It splits the variation into between-subject and within-subject parts, and tests the time (or condition) effect against the within-subject noise only. It assumes sphericity: equal variances of all differences between time points. The mixed model with a random intercept gives the same test under that assumption, and handles unequal spacing, missing visits and other correlation structures as well.

Code

patients, visits = 20, 4
level = rng.normal(50, 10, patients)
scores = level[:, None] - 2.0 * np.arange(visits) + rng.normal(0, 3, (patients, visits))
long = pd.DataFrame(
    {
        "patient": np.repeat(np.arange(patients), visits),
        "visit": np.tile(np.arange(visits), patients),
        "score": scores.ravel(),
    }
)
oneway = smf.ols("score ~ C(visit)", long).fit()
repeated = sm.stats.AnovaRM(long, "score", "patient", within=["visit"]).fit().anova_table

fig, axes = style.subplots(2, height=3.2, width=10)
for p in range(patients):
    axes[0].plot(range(visits), scores[p], color=style.GREY, lw=0.8)
axes[0].plot(range(visits), scores.mean(axis=0), "o-", color=style.CORAL, lw=3)
axes[0].set(
    xticks=range(visits), xticklabels=[f"visit {v + 1}" for v in range(visits)], ylabel="score"
)
tests = {
    "ANOVA ignoring patients": (oneway.fvalue, oneway.f_pvalue),
    "repeated measures": (repeated["F Value"].iloc[0], repeated["Pr > F"].iloc[0]),
}
axes[1].barh(
    list(tests), [f for f, _ in tests.values()], color=[style.GREY, style.BLUE], alpha=0.7
)
for i, (f, p) in enumerate(tests.values()):
    axes[1].text(f, i, f"  F = {f:.1f}, p = {p:.2g}", va="center", fontsize=8.5)
axes[1].invert_yaxis()
axes[1].set(xlabel="F statistic for visits", xlim=(0, max(f for f, _ in tests.values()) * 1.6))
fig.tight_layout()
plt.show()

Figure 4.8.5: Twenty simulated patients measured at four visits during a treatment. Left: each patient’s path (grey) and the visit means (red). Patients differ a lot in level, but each improves. Right: the F-test for visits from a one-way ANOVA that ignores who is who, and from a repeated measures analysis (a mixed model with a random intercept per patient). Removing the between-patient differences makes the same change clearly visible.

4.8.6 Dropout

Long studies lose people, and often not at random: patients who do badly leave. The average of those still present then drifts upwards and suggests an improvement that did not happen. A mixed model, fitted to all available visits, uses each patient’s early values and stays valid if dropout depends only on what was observed (missing at random, Chapter 5.3). Carrying each patient’s last value forward is a common but biased shortcut.

Code

N_PAT, N_VISIT = 300, 6
person = rng.normal(0, 6, N_PAT)
true = 50 + person[:, None] + rng.normal(0, 3, (N_PAT, N_VISIT))
present = np.ones((N_PAT, N_VISIT), bool)
for v in range(1, N_VISIT):
    leave = rng.uniform(size=N_PAT) < 1 / (1 + np.exp((true[:, v - 1] - 42) / 3))
    present[:, v] = present[:, v - 1] & ~leave
observed = pd.DataFrame(
    {
        "patient": np.repeat(np.arange(N_PAT), N_VISIT),
        "visit": np.tile(np.arange(N_VISIT), N_PAT),
        "score": true.ravel(),
        "present": present.ravel(),
    }
)
seen = observed[observed.present]
mixed = smf.mixedlm("score ~ C(visit)", seen, groups="patient").fit()
mixed_path = (
    mixed.fe_params["Intercept"]
    + np.r_[0, [mixed.fe_params[f"C(visit)[T.{v}]"] for v in range(1, N_VISIT)]]
)
carried = pd.DataFrame(np.where(present, true, np.nan)).ffill(axis=1).to_numpy()

fig, ax = plt.subplots(figsize=(8, 3.2))
ax.plot(
    range(1, N_VISIT + 1),
    seen.groupby("visit").score.mean(),
    "o-",
    color=style.GREY,
    lw=2,
    label="mean of those still present",
)
ax.plot(
    range(1, N_VISIT + 1),
    carried.mean(axis=0),
    "o-",
    color=style.GOLD,
    lw=2,
    label="last value carried forward",
)
ax.plot(range(1, N_VISIT + 1), mixed_path, "o-", color=style.BLUE, lw=2, label="mixed model")
ax.axhline(true.mean(), color=style.CORAL, ls="--", lw=1.2)
style.note(
    ax,
    f"still present at visit 6: {present[:, -1].mean():.0%}",
    x=0.03,
    y=0.95,
    color=style.INK,
)
ax.legend(fontsize=8, loc="upper left", bbox_to_anchor=(0, 0.88))
ax.set(xlabel="visit", ylabel="average score")
plt.show()

Figure 4.8.6: Simulated patients whose true average score is flat over six visits. After each visit, patients with a low score are more likely to drop out. Grey: the mean of the patients still in the study rises, a false improvement. Gold: carrying the last observation forward. Blue: the mixed model’s average path, which uses all observed visits and recovers the flat line (dashed).

4.8.7 In one picture

Start with the spaghetti plot: one line per unit.
Growth curve (mixed) models give each unit its own intercept and slope around an average path.
Ignoring the repeated structure gives far too narrow intervals; use mixed models, GEE or cluster-robust errors.
Mixed models are subject-specific, GEE population-averaged; they differ for yes/no and count outcomes.
Repeated measures ANOVA is the classic special case; it needs sphericity.
Dropout related to the outcome biases simple averages and last-value-carried-forward; mixed models on all visits are valid under missing at random.

Next: Part 5, how the data are collected, and when they can show cause and effect.