Module 2 · Longitudinal data — notation & visualization

Lectures 7.1–7.2. Reading: Lecture Notes Part 2, §2.1. Code: R script — Lecture Video 32.

← Back to notes index

Part 2 leaves survival behind. The outcome is no longer a time-to-event but a measurement repeated on the same people over time — weight, a blood biomarker, a cognitive score. Some examples will feel like the time-varying covariates of §1.9, but nothing here is about hazards.

Lecture 7.1 — What longitudinal data is, and why it is hard

Longitudinal data are repeated observations taken on individuals. The analysis usually has two goals:

  1. Characterize the change in the response over time — not just whether it changes but the shape of the change.
  2. Identify factors that influence that change.

The complication is the one sentence the whole module turns on:

Observations that belong to the same subject tend to be more similar than observations from different subjects.

Your weight tomorrow looks like your weight today far more than it looks like anyone else’s weight. That within-subject correlation must be accounted for to get valid inference — ignoring it treats $N \times n$ correlated numbers as if they were $N \times n$ independent ones.

Longitudinal data are a special case of hierarchical data — observations nested in levels: pupils within schools (within regions), patients within hospitals. Here the nesting is measurements within individuals.

Notation

Symbol Meaning
$Y_{ij}$ response for individual $i = 1, \dots, N$ at time $t_{ij}$, $j = 1, \dots, n$
$\mu_j = E(Y_{ij})$ mean at occasion $j$ — the location of the distribution
$\sigma^2j = \text{Var}(Y{ij})$ variance at occasion $j$ — the spread
$\sigma_{jk} = \text{Cov}(Y_{ij}, Y_{ik})$ covariance of two measurements on the same person at occasions $j$ and $k$; $\sigma_{jk} = 0$ means no linear dependence
$r_{jk} = \sigma_{jk} / (\sigma_j \sigma_k)$ correlation — scale-free, in $[-1, 1]$
$\mathbf{Y}i = (Y{i1}, \dots, Y_{in})^T$ the vector of person $i$’s repeated measures

So $N$ people, $n$ time points. The covariance of the vector collects all of the above into a symmetric $n \times n$ matrix,

\[\text{Cov}(\mathbf{Y}_i) = \Sigma_i = \begin{pmatrix} \sigma^2_1 & \sigma_{12} & \cdots & \sigma_{1n} \\ \sigma_{12} & \sigma^2_2 & \cdots & \cdot \\ \cdot & \cdot & \ddots & \cdot \\ \sigma_{1n} & \cdot & \cdots & \sigma^2_n \end{pmatrix}\]

— variances on the diagonal, covariances off it — and the correlation matrix is the same thing standardized, with $1$ on the diagonal:

\[\mathbf{r}_i = \begin{pmatrix} 1 & r_{12} & \cdots & r_{1n} \\ r_{12} & 1 & \cdots & \cdot \\ \cdot & \cdot & \ddots & \cdot \\ r_{1n} & \cdot & \cdots & 1 \end{pmatrix}\]

If the off-diagonals were all zero we would have longitudinal data in form but not in substance — ordinary methods would do. On the same person they are essentially never zero.

Four empirical facts about $r_{jk}$

With longitudinal data it is usually the case that the correlations

  1. are positive;
  2. decrease with increasing time separation — weight today vs. tomorrow is more correlated than weight today vs. six months out;
  3. rarely approach zero, however far apart the occasions;
  4. rarely approach one, however close together.

Why this matters — facts 3 and 4 are what make the problem hard. If the correlation died out after a week, measurements far enough apart could be treated as independent and everything would simplify. If nearby measurements were perfectly correlated, they would be redundant and everything would simplify the other way. Neither extreme happens, so the correlation structure has to be modelled rather than assumed away. Sections 2.2 and 2.3 are two different ways of doing that.

The dental growth data

The running example (Potthoff & Roy 1964): 16 boys and 11 girls, each measured at ages 8, 10, 12 and 14 — the distance (mm) from the centre of the pituitary gland to the pteryomaxillary fissure, a quantity that matters in orthodontic therapy.

Variable Content
id unique child ID, 1–27
sex Girl / Boy (coded 0 / 1 in the notes)
y8, y10, y12, y14 the distance at each age

So $N = 27$, $n = 4$, everyone measured at the same four times.

The three research questions

Stated for the dental data, then in the general form that applies to any longitudinal study:

  Dental data In general
1. Time effect How does dental growth change over time? What is the shape of the trajectory? What is the shape of the trajectory of the mean response over time?
2. Group effect Is there a difference in dental growth between boys and girls? What is the average difference between groups of individuals?
3. Time × group interaction How does the relationship between dental growth and time vary by sex? How does the relationship between the response and time vary according to group?

“Shape” in question 1 is deliberate: the answer is not “yes, it changes” but increasing, decreasing, increasing then plateauing, and so on. Question 2 says nothing about time; question 3 is where the two meet. Every model in the next sections is built to answer one of these three.

Lecture 7.2 — Looking at the data

Always plot the data, but for longitudinal data especially: the plots illustrate both the goals and the patterns — the changes over time, and whether some factor goes with a differential change. Each plot below corresponds to one piece of the notation.

Means over time → $\mu_j$

Summary statistics of the distance at each age:

Age 8 10 12 14
mean (mm) 22.19 23.17 24.65 26.09
sd 2.43 2.16 2.82 2.77

Before any analysis: the mean increases with age, by about 4 mm from 8 to 14, and the spread grows a little too. A boxplot per age with the individual points jittered on top (Figure 2.1 in the notes) shows the same thing with the interquartile ranges visible.

By sex:

  8 10 12 14
Girls (n = 11) 21.18 22.23 23.09 24.09
Boys (n = 16) 22.88 23.81 25.72 27.47

Both increase; boys are higher throughout and the gap widens — a first hint at all three questions at once. Figure 2.2 groups the boxplots girls-left / boys-right with age in the legend; Figure 2.3 is a bar chart of just the means (less information than the boxplot, no spread).

Correlations → $\mathbf{r}_i$

The estimated correlation matrix of the four measurements:

  y8 y10 y12 y14
y8 1 0.626 0.711 0.600
y10   1 0.635 0.759
y12     1 0.795
y14       1

All positive and far from zero — an 8-year-old’s distance clearly carries information about their distance at 10. Facts 1, 3 and 4 above are visible; fact 2 is only loosely so here (0.600 at six years apart vs. 0.626 at two), which is a reminder that these are estimates from 27 children.

The ggpairs plot (Figure 2.4) packs a lot into one panel: densities on the diagonal, the numeric correlations with significance stars in the upper triangle, and scatter plots with a fitted line and confidence band in the lower triangle. Done by sex (Figure 2.5):

Pair All Girls Boys
y8–y10 0.626 0.830 0.437
y8–y12 0.711 0.862 0.558
y8–y14 0.600 0.841 0.315
y10–y12 0.635 0.895 0.387
y10–y14 0.759 0.879 0.631
y12–y14 0.795 0.948 0.586

Girls’ measurements are far more correlated across time than boys’ — girls’ trajectories are more parallel; boys’ cross more. That is itself a group difference, in the covariance rather than the mean.

Trajectories → the shape of change

Three views of the same 108 points:

  • Mean trajectory by sex with confidence intervals (Figure 2.6) — the mean at each age for each sex, joined by lines, with error bars. The girls’ points are nudged 0.05 years left so the bars at age 8 do not sit on top of each other.
  • One panel per child (Figure 2.7) — feasible with 27, not with hundreds. Shows what an individual trajectory actually looks like: mostly rising, not perfectly smooth, a few dips.
  • Spaghetti plot (Figure 2.8) — every child’s line on one set of axes. By sex with a smooth curve overlaid (Figure 2.9): the boys’ fan is visibly wider than the girls’, matching the correlation table.

The boxplots describe $\mu_j$, the pairs plot describes $\mathbf{r}_i$, and the spaghetti plot shows the raw $\mathbf{Y}_i$ vectors themselves.

In R

From R script — Lecture Video 32. This is where the course switches from base graphics to the tidyverse and ggplot2 — more packages, much better longitudinal plots. The install list is long (tidyverse, labelled, rstatix, ggpubr, GGally, car, Epi, lme4, lmerTest, emmeans, multcomp, geepack, ggeffects, gt) and can take anywhere from minutes to overnight; most are for the modelling sections to come. The script avoids the pipe (%>% / |>) on purpose and does everything step by step. Reference: R for Data Science, r4ds.had.co.nz.

load(url("http://alecri.github.io/downloads/data/dental.RData"))
dental

Wide → long

dental arrives wide: one row per child, the four measurements in four columns. For plotting and modelling you usually want long: one row per measurement, so four rows per child.

Wide (27 rows)   Long (108 rows)
id sex y8 y10 y12 y14 → id sex measurement distance
1 Girl 21.0 20.0 21.5 23.0   1 Girl y8 21.0 / 1 Girl y10 20.0 / …

Three steps, each worth running on its own to see what changed:

dental_long = pivot_longer(dental, cols = starts_with("y"),
                           names_to = "measurement", values_to = "distance")

pivot_longer takes the columns whose names start with y, puts those column names into a new column measurement and their values into distance.

dental_long = mutate(dental_long,
                     age = parse_number(measurement),
                     measurement = fct_inorder(paste("Measure at age", age)))

parse_number("y10") is 10 — it pulls the number out of a string and ignores the rest. The measurement column is then rewritten as a label ("Measure at age 10") and fct_inorder fixes the factor levels in the order they appear, so plots run 8, 10, 12, 14 rather than alphabetically. Redundant with age, but it becomes the axis text.

dental_long = set_variable_labels(dental_long,
    age = "Age of the child at measurement",
    measurement = "Label for time measurement",
    distance = "Measurement")

Variable labels that plotting functions pick up automatically. When a plot later “looks perfect”, this preparation is why.

Summaries and the anatomy of a ggplot

get_summary_stats(group_by(dental_long, age), distance)
get_summary_stats(group_by(dental_long, sex, measurement), distance, show = c("mean", "sd"))

The first gives the four-row table above (min, max, median, quartiles, mean, sd …); the second the eight-row by-sex version, trimmed to mean and sd.

The first plot, built up one layer at a time — the way to understand any ggplot:

ggplot(dental_long, aes(measurement, distance, fill = measurement)) +
  geom_boxplot() +
  geom_jitter(width = 0.2) +
  guides(fill = "none") +
  labs(x = "", y = "Dental growth, mm") +
  theme(text = element_text(size = 18))
Layer What it does Without it
ggplot(..., aes(x, y, fill)) sets the axes and a fill colour per measurement empty axes, nothing drawn
geom_boxplot() the boxes —
geom_jitter(width = 0.2) the raw points, wiggled horizontally so ties don’t overprint geom_point() stacks tied values on top of each other; width = 0 is the same
guides(fill = "none") removes the legend a legend that just repeats the x-axis labels
labs(x = "", y = ...) axis titles lowercase variable names
theme(text = element_text(size = 18)) bigger text default small text

By sex, two ways:

# boxplots, girls on one side / boys on the other, age in the legend (legend now useful)
ggplot(dental_long, aes(sex, distance, fill = measurement)) + geom_boxplot() +
  labs(x = "", y = "Dental growth, mm", fill = "")

# bar chart of means only, with the value printed on each bar
ggplot(summarise(group_by(dental_long, sex, measurement),
                 mean_distance = mean(distance), .groups = "drop"),
       aes(sex, mean_distance, fill = measurement, label = round(mean_distance))) +
  geom_col(position = "dodge") +
  geom_text(position = position_dodge(width = 0.9), vjust = -0.5) +
  coord_flip()

Correlations

cov_obs = cov(dplyr::select(dental, starts_with("y")))   # the 4x4 covariance matrix
cov2cor(cov_obs)                                          # ... as correlations

This uses the wide data: cov() wants one column per occasion. cov2cor converts $\Sigma$ to $\mathbf{r}$.

ggpairs(dplyr::select(dental, starts_with("y")), lower = list(continuous = "smooth"))
ggpairs(dental, mapping = aes(colour = sex), columns = 3:6, lower = list(continuous = "smooth"))

The lower = list(continuous = "smooth") argument adds the fitted line and band to the scatter plots; drop it and you get plain points. The by-sex version takes the whole data frame (it needs the sex column) and columns = 3:6 to say which columns to correlate.

Trajectories

The mean-with-CI plot needs a small summary table first:

plot_data = group_by(dental_long, sex, age)                       # 8 groups
plot_data = summarise(plot_data, mean = list(mean_ci(distance)), .groups = "drop")
plot_data = unnest_wider(plot_data, mean)                         # -> y, ymin, ymax
plot_data = mutate(plot_data, agex = age - .05 + .05*(sex == "Boy"))   # the wiggle

ggplot(plot_data, aes(agex, y, col = sex, shape = sex)) +
  geom_point() +
  geom_errorbar(aes(ymin = ymin, ymax = ymax), width = 0.2) +
  geom_line() +
  labs(x = "Age, years", y = "Mean Dental growth, mm", shape = "Sex", col = "Sex")

mean_ci() returns the mean and CI as a list in one cell; unnest_wider splits that into three columns. .groups = "drop" just forgets the grouping once the summary is done. agex shifts girls 0.05 years left and boys stay put, so the error bars at the same age don’t overlap.

Individual trajectories:

# one panel per child
ggplot(dental_long, aes(age, distance, col = factor(id))) +
  geom_point() + geom_line() + facet_wrap(~ id) +
  labs(x = "Age, years", y = "Dental growth, mm", col = "Child id") +
  guides(col = guide_legend(nrow = 3)) + theme(legend.position = "bottom")

# spaghetti plot: same, without facet_wrap
ggplot(dental_long, aes(age, distance, col = factor(id))) + geom_line() + ...

# spaghetti by sex with a smooth
ggplot(dental_long, aes(age, distance)) +
  geom_line(aes(group = factor(id))) +
  geom_smooth() +
  facet_grid(~ sex)

facet_wrap(~ id) makes a panel per child; facet_grid(~ sex) puts the sexes side by side. factor(id) is needed so that colour is categorical (27 colours) rather than a gradient over the numbers 1–27. In the last plot the lines are grouped by child but not coloured by child, which is the trick the class activity asks for. geom_smooth() warns that it has few points per x — expected with four ages.

Class activity: body fat

bodyfat — 162 girls from the MIT Growth and Development Study, pre-menarcheal and non-obese at entry, measured roughly annually until four years after menarche. Response pbf (percent body fat, by bioelectric impedance); occasion 1–10; also agec, agem (age at menarche) and time (years relative to menarche).

Two differences from the dental data: it is already long (one row per measurement), and occasion does not mean age — everyone’s occasion 1 is a different age, so a mean “at each age” makes no sense and the visualizations go by occasion. The design is unbalanced: 162 girls at occasions 1–3, dropping to 9 at occasion 9 and a single girl at occasion 10. (A real analysis would use time relative to menarche, which is the study’s clock.)

  1. Boxplot per occasion — no pivot needed, but occasion is numeric and must become a factor for grouping:

    bodyfat = mutate(bodyfat, occasion_f = as.factor(occasion))
ggplot(bodyfat, aes(occasion_f, pbf, fill = occasion_f)) + geom_boxplot() +
  geom_jitter(width = 0.2) + guides(fill = "none") + labs(x = "", y = "Percent Body Fat, %")

  2. Correlation plot — ggpairs needs wide data, so pivot the other way, then drop the id column (correlating IDs with body fat is meaningless) and occasion 10 (one girl; no correlation can be computed):

    bodyfat_wide = pivot_wider(bodyfat[, c(1, 5, 6)], names_from = "occasion", values_from = "pbf")
ggpairs(bodyfat_wide[, -c(1, 11)], lower = list(continuous = "smooth"))

    Adjacent occasions correlate around 0.8–0.9; all positive — facts 1–4 again.

  3. Spaghetti plot, no legend, black lines — with 162 girls a per-person legend is useless. theme(legend.position = "none") removes it; and colouring is what makes the rainbow, so moving id from col = to group = gives black lines:

    ggplot(bodyfat, aes(occasion, pbf)) + geom_line(aes(group = id)) +
  labs(x = "Occasion", y = "Percent Body Fat, %") + theme(legend.position = "none")

What comes next

With the questions framed and the data visualized, §2.2–2.5 of the notes supply the analysis: linear mixed models (random intercepts and slopes, variance components, the ICC, REML, BLUPs), generalized estimating equations (marginal vs. subject-specific models and why their coefficients differ for binary outcomes), joint modelling of a longitudinal biomarker with a time-to-event outcome, and missing data (MCAR / MAR / MNAR, LOCF, multiple imputation). Those follow in the next notes pages as the lectures and scripts arrive, reusing the dental data and these same plots to display the fitted models.