Response Vector

Each subject has a vector of responses:

\[\mathbf{Y}_i = (Y_{i1}, Y_{i2}, \dots, Y_{in_i})^\top\]

Dimension: \(n_i \times 1\)

Covariate Information

Covariates for subject \(i\) at visit \(j\) form a row vector \(X_{ij}\) (dimension \(1 \times p\)):

\[X_{ij} = (x_{ij1}, x_{ij2}, \dots, x_{ijp})\]

The mean model is written \(E(Y_{ij}) = X_{ij}\boldsymbol\beta\), the row vector times the coefficient vector (\(\boldsymbol\beta\) is the vector of population-average regression coefficients).

Covariates include:

• Time-invariant: sex, treatment arm, genotype
• Time-varying: weight, adherence, current dose

Each row is one occasion \(X_{ij}\); each column is a covariate. So \(X_i\) is \(n_i \times p\) (here \(3 \times 3\)).

Design Matrix

Stack the covariate row vectors for subject \(i\):

\[\mathbf{X}_i = \begin{pmatrix} X_{i1} \\ X_{i2} \\ \vdots \\ X_{in_i} \end{pmatrix}\]

Dimension: \(n_i \times p\). The mean for the whole subject is \(E(\mathbf{Y}_i) = \mathbf{X}_i\boldsymbol\beta\).

Complete Data Structure

The basic longitudinal data structure:

\[(\mathbf{Y}_i, \mathbf{X}_i), \quad i = 1, \dots, N\]

Key features:

• Multiple responses per subject
• Within-subject correlation
• Unequal \(n_i\) possible

Wide vs. Long Format

Wide format (each row = subject):

Wide vs. Long Format

Long format (each row = observation):

Most R modeling functions expect long format.

Reshaping: Wide to Long

df_wide <- tibble(
  id = c(1, 2),
  sex = c(0, 1),
  bp_t1 = c(120, 135),
  bp_t2 = c(122, 137),
  bp_t3 = c(125, 140)
)

df_long <- df_wide |>
  pivot_longer(cols = starts_with("bp"),
               names_to = "visit",
               values_to = "bp") |>
  # build a numeric time index from the visit label:
  #   gsub("bp_t", "", visit) strips the "bp_t" prefix, leaving "1", "2", "3"
  #   as.numeric() coerces that text to a number
  #   - 1 shifts the index so the first visit is baseline = 0
  mutate(time = as.numeric(gsub("bp_t", "", visit)) - 1)

df_long

Interpretation: The output shows the same data restructured so each row represents one observation (one visit for one subject). Subject 1 now has 3 rows (one per visit) instead of 3 columns. The time variable is coded as 0, 1, 2 to represent the three visits. This long format is required by most R modeling functions like lme() and gls().

Covariance and Correlation in Real Data

The chapter is built around one object: the \(n_i \times n_i\) covariance matrix of a subject’s repeated measures (and its standardized form, the correlation matrix). Let us estimate it from the TLC lead-exposure trial, placebo group (50 children measured at weeks 0, 1, 4, 6). These reproduce FLW Tables 2.2-2.3.

# tlc.csv is WIDE: ID, Treatment Group, then lead level at weeks 0, 1, 4, 6
tlc_file <- "../../data/tlc.csv"
if (file.exists(tlc_file)) {
  tlc <- read.csv(tlc_file)
} else {
  # fallback: the FLW site ships a whitespace-delimited, header-free file
  tlc <- read.table("https://content.sph.harvard.edu/fitzmaur/ala2e/tlc-data.txt",
                    header = FALSE,
                    col.names = c("ID", "Treatment.Group",
                                  "Lead.Level.Week.0", "Lead.Level.Week.1",
                                  "Lead.Level.Week.4", "Lead.Level.Week.6"))
}

# placebo group only: one row per subject, one column per occasion
pbo <- tlc[tlc$Treatment.Group == "P", 3:6]
colnames(pbo) <- c("Week 0", "Week 1", "Week 4", "Week 6")
dim(pbo)          # 50 subjects x 4 occasions

Source 1: Between-Individual Heterogeneity

• Some individuals are consistently high responders
• Others are consistently low responders
• Creates positive correlation across all time points

Statistical approach: Random effects (Ch. 8)

Visualizing Heterogeneity

set.seed(667)
df <- tibble(
  id = rep(1:8, each = 6),
  time = rep(1:6, 8),
  # Random intercept per subject
  y = rep(rnorm(8, 10, 3), each = 6) + rnorm(48, 0, 1)
)

ggplot(df, aes(x = time, y = y, group = id, color = factor(id))) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  theme_minimal() +
  theme(legend.position = "none") +
  labs(title = "Between-Subject Heterogeneity",
       subtitle = "Each subject has a different baseline level",
       x = "Time", y = "Outcome")

Interpretation: Notice how each line (representing one subject) stays roughly at its own level throughout the study. Some subjects sit consistently high, others stay low. This separation illustrates between-subject heterogeneity: some individuals are consistently higher or lower responders. This creates positive correlation among all measurements within each person.

Source 2: Within-Individual Biological Variation

• Natural fluctuations around homeostatic set point
• Driven by circadian rhythms, diet, environment
• Serial correlation: Closer measurements -> higher correlation

Serial Correlation

AR(1) (first-order autoregressive) means each measurement equals \(\rho\) times the previous one plus noise, so the correlation at lag \(k\) is \(\rho^k\) (\(\rho\) is the within-subject correlation parameter; for AR(1) the lag-\(k\) correlation is \(\rho^k\)).

set.seed(667)

# Simulate AR(1) process
n <- 50
rho <- 0.8
bp <- arima.sim(model = list(ar = rho), n = n, sd = 5) + 120

tibble(time = 1:n, bp = as.numeric(bp)) |>
  ggplot(aes(x = time, y = bp)) +
  geom_line(linewidth = 0.8, color = "steelblue") +
  geom_point(size = 2, color = "steelblue") +
  theme_minimal() +
  labs(title = "Serial Correlation in Repeated Measures",
       subtitle = expression(paste("AR(1) with ", rho, " = 0.8")),
       x = "Time", y = "Blood Pressure")

Interpretation: This plot shows blood pressure measurements over 50 time points from a single simulated individual. Notice how the trajectory does not jump erratically: when blood pressure is high at one time point, it tends to remain elevated at the next. This “smoothness” reflects serial correlation: measurements close in time are more similar than measurements far apart. The AR(1) model with \(\rho = 0.8\) captures this pattern mathematically.

Autocorrelation Function

acf(bp, main = "Autocorrelation of Repeated Measures",
    col = "steelblue", lwd = 2)

Interpretation: The autocorrelation function (ACF) plot shows how correlation changes with time lag. The x-axis shows the lag (number of time steps apart), and the y-axis shows the correlation. At lag 0, the correlation is 1 (a measurement is perfectly correlated with itself). At lag 1, the correlation is about 0.8, meaning measurements one time step apart are highly correlated. As the lag increases, the bars get shorter, showing that correlation decays with time separation. The dashed blue lines indicate statistical significance; bars extending beyond these lines suggest meaningful correlation.

Source 3: Measurement Error

• Imprecision in instruments/procedures
• Distinct from biological variability
• Attenuates observed correlations

Reliability = Var(True) / Var(Total)

Reliability Examples

Lower reliability -> weaker observed correlations.

Attenuation Effect

Reliability \(=\) Var(true) / Var(total), and the observed correlation \(=\) true correlation \(\times\) reliability, so low reliability attenuates (shrinks) the observed correlation toward 0.

set.seed(667)
n <- 300
true_data <- matrix(rnorm(n * 2, 10, 2), ncol = 2)

add_error <- function(data, reliability) {
  v_true <- var(as.vector(data))
  v_err <- v_true * (1 - reliability) / reliability
  data + matrix(rnorm(length(data), 0, sqrt(v_err)), ncol = 2)
}

rels <- c(0.98, 0.85, 0.50)
df_plot <- do.call(rbind, lapply(rels, function(r) {
  d <- add_error(true_data, r)
  tibble(t1 = d[,1], t2 = d[,2], Reliability = paste("Reliability =", r))
}))

ggplot(df_plot, aes(t1, t2)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  facet_wrap(~Reliability) +
  theme_minimal() +
  labs(title = "Attenuation of Correlation with Lower Reliability",
       x = "Measurement at Time 1", y = "Measurement at Time 2")

Summary: Three Sources

All three combine to determine observed correlation structure.

Common Mistakes to Avoid

• Treating repeated measures as independent: inflates Type I error and corrupts p-values.
• Confusing wide and long formats: long-format functions like lme() reject wide data.
• Ignoring the source of correlation: assuming pure heterogeneity when serial correlation is also present misspecifies the model.
• Conflating correlation with causation: temporal ordering helps, but it is not proof.
• Thinking ignored correlation biases \(\hat{\beta}\): \(\hat{\beta}\) stays consistent; it is the standard errors that break.

Canonical Resources

• Author slides: content.sph.harvard.edu/fitzmaur/ala2e/
• Datasets: content.sph.harvard.edu/fitzmaur/ala2e/datasets.html
• Sample R code: content.sph.harvard.edu/fitzmaur/ala2e/SampleR.html

For Next Time

Read Chapter 3 of Fitzmaurice, Laird & Ware (2011)

• Lecture 3 (Ch. 3): the general linear model for longitudinal data