Repeated measures in R

Advanced

Each row of a repeated-measures table is a visit, and the independent unit is the person. A random intercept and a marginal model are two commands for that dependence, and they answer different questions when the link is not the identity.

What you need first. Read Random intercepts and Fitted models as objects.

Part of Using R · 1 question on this topic · What links here

Repeated measures in R

Part Using R
Level Advanced
Chapter 10 of 12 in this part
Previous Generalised additive models in R
Next Time series in R
Questions 1
Definitions Exchangeable correlation

Orthodont records a craniofacial distance, in millimetres, for 27 children at ages 8, 10, 12, and 14. Each row is a visit. Each child is a cluster of four visits. A regression that treats all 108 visits as independent observations uses the wrong unit: it pretends there are 108 children, and its standard errors are then too small for any contrast that varies between children. Random intercepts, and the difference between a fixed block and a random intercept, are in Random intercepts. A marginal mean with a working correlation and a sandwich variance is in Marginal models. This page is the two commands, lme4::lmer and geepack::geeglm, and the constraint on an exchangeable correlation.

data("Orthodont", package = "nlme")
teeth <- as.data.frame(Orthodont) |>
  dplyr::mutate(Subject = factor(Subject), age_from_8 = age - 8) |>
  dplyr::arrange(Subject, age)
teeth |>
  dplyr::summarise(visits = dplyr::n(), children = dplyr::n_distinct(Subject))

The summary should report 108 visits and 27 children. age_from_8 puts the intercept at age 8, the first occasion, so the intercept is a distance at the start of the study rather than a distance at age 0, which is an extrapolation. arrange(Subject, age) puts each child’s visits in a contiguous block. The estimating-equation command below requires that order.

ggplot2::ggplot(teeth, ggplot2::aes(age, distance, group = Subject, colour = Sex)) +
  ggplot2::geom_line(alpha = 0.42) +
  ggplot2::geom_point(alpha = 0.55, size = 1.2) +
  ggplot2::labs(x = "Age (years)", y = "Distance (mm)", colour = "Sex") +
  ggplot2::theme_minimal(base_size = 12)

One line is one child. The picture is the reason a random intercept is plausible: children differ in level by more than the visit-to-visit wiggle. It is not yet a fitted model.

A random intercept and a random slope

The random-intercept model gives each child a level \alpha_j, drawn from a common distribution, and a mean that is linear in age and sex. Adding a random slope lets the age trend differ by child as well. The fixed effects describe the mean structure. The variance components describe the estimated spread of those child-specific departures. They are not standard errors of the fixed effects.

random_intercept <- lme4::lmer(
  distance ~ age_from_8 * Sex + (1 | Subject),
  data = teeth, REML = TRUE
)
random_slope <- lme4::lmer(
  distance ~ age_from_8 * Sex + (age_from_8 | Subject),
  data = teeth, REML = TRUE
)
lme4::fixef(random_slope)
lme4::VarCorr(random_slope)
lme4::isSingular(random_slope)

(1 | Subject) is an intercept that varies by child. (age_from_8 | Subject) is an intercept and a slope, both varying by child, with a covariance between them. fixef is the estimated mean structure. VarCorr is the estimated child-level covariance. isSingular is true when that covariance was estimated on the boundary, typically a variance or a correlation pinned at zero. A singular fit means the richer random-effects model is not supported by these 27 children. The fit has then collapsed toward a simpler one, and the printed variance of zero is a boundary, not a precise estimate of a tiny variance. Prefer the simpler model, or treat the singular fit as a warning, rather than quoting the zero as a finding.

Restricted maximum likelihood is the right default for the variance components when the mean structure is held fixed. A likelihood-ratio comparison of two models that differ in their fixed effects should not use restricted likelihoods. The two restricted likelihoods are not likelihoods of the same contrasts. Refit both with REML = FALSE and then compare them. A test of a variance component against zero has its own problem: the null is on the boundary of the parameter space, so a naive chi-squared reference distribution from anova is not the null distribution. A significant-looking p-value from that comparison can be too small.

The fixed effect of age_from_8 is a change in mean distance per year of age, for the baseline sex, inside this model. The interaction age_from_8:Sex is the difference of that yearly change between the two sexes. Both are conditional on the random-effects model. They are not the same numbers, in general, as the coefficients of a marginal model that does not condition on the child.

A marginal model

A generalised estimating equation specifies the marginal mean and a working correlation inside each child. The mean here uses an identity link and a Gaussian family, so a coefficient is a change in millimetres, and it can be compared in unit with the linear mixed model’s fixed effects. For a nonlinear link the comparison is not that direct. A mixed-model coefficient is conditional on the random effect. A marginal coefficient averages over it. On a logit scale those are different estimands, and a larger conditional coefficient is not evidence that the marginal fit is biased.

gee_fit <- geepack::geeglm(
  distance ~ age_from_8 * Sex,
  id = Subject,
  data = teeth,
  family = gaussian(),
  corstr = "exchangeable"
)
summary(gee_fit)$coefficients

id = Subject says which rows share a working correlation. geeglm requires the rows of one id to be contiguous. The arrange above is that requirement, not a cosmetic sort. corstr = "exchangeable" uses one common pairwise correlation inside a child. An independence working correlation fits the same mean and uses a different set of estimating equations for the coefficients. The sandwich standard errors are meant to stay valid when the working correlation is wrong, under conditions that include a large number of independent clusters and a correct mean. Twenty-seven children is a modest number of clusters. Treat the asymptotic standard errors as approximate. The sandwich does not repair a wrong mean, dropout that depends on the unobserved distance, or dependence between children who were treated as independent clusters.

gee_fit is fit on the complete rows of the variables in the formula. A child with a missing visit is not handled by an imputation inside this command. If missing visits were dropped before the call, the sample is the visits that remain, and the sandwich does not put the dropped visits back.

The bound on an equicorrelation

Definition 1 (Exchangeable correlation) An exchangeable correlation matrix for a cluster of K observations, with common pairwise correlation \rho, is R(\rho)=(1-\rho)I_K+\rho\,\mathbf{1}\mathbf{1}^{\top}. Its eigenvalues are 1-\rho, with multiplicity K-1, and 1+(K-1)\rho, with multiplicity one. R(\rho) is positive semidefinite if and only if -1/(K-1)\le\rho\le 1. At either endpoint it is singular.

Article: Exchangeable correlation →

The eigenvalues are worth deriving once, because the bound is easy to miss in the output. R(\rho)\mathbf{1}=(1-\rho)\mathbf{1}+\rho K\mathbf{1}=(1+(K-1)\rho)\mathbf{1}, so 1+(K-1)\rho is an eigenvalue with eigenvector \mathbf{1}. Any vector orthogonal to \mathbf{1} is sent to (1-\rho) times itself, which gives the other eigenvalue and the multiplicity K-1. Both eigenvalues are nonnegative exactly on -1/(K-1)\le\rho\le 1.

For four visits, the lower bound is -1/3. A working correlation of -0.5 cannot be the correlation of four visits inside one child: no such correlation matrix exists. For a cluster of size 6 the bound is -0.2. geeglm can report a working correlation at or near the boundary when the data ask for a value the matrix cannot hold. That report is a constraint of linear algebra, not a precise estimate. Whether an exchangeable correlation is plausible at all, for distances that drift smoothly with age, is a separate scientific question. A child whose distance rises steadily will have adjacent ages more alike than ages six years apart, which is a patterned correlation, not one common \rho.

Key idea

The independent unit is the child, not the visit. A singular random-effects fit is a boundary. An exchangeable correlation for K visits cannot lie below -1/(K-1). Conditional and marginal coefficients agree in unit here because the link is the identity; they do not agree in general.

Questions on this topic

Open a question to attempt it. Its worked answer is folded away at the bottom of the page.

1 question.

Question Level Answer
A correlation four visits cannot have Advanced ✓

See also