V The V Lab

Longitudinal Data Analysis in Depth

From repeated measurements and covariance structures to marginal models, mixed-effects trajectories, missing data, prediction, and responsible reporting

The V Lab Home Biostatistics Foundations Medical Statistical Tests Linear Regression Logistic Regression Design of Experiments Clinical Trials Poisson and Zero-Inflated Models Multilevel Modelling Longitudinal Data Analysis Causal Inference Propensity Score Matching Survival Analysis AFT and Cox PH Meta-Analysis 中文版

AudienceLearners in public health, epidemiology, medicine, psychology, and health data science

Study timeApproximately 210–270 minutes

PrerequisitesLinear regression, confidence intervals, model matrices, and basic R

About the simulated cohort This module follows 360 simulated adults assigned to usual care or a blood-pressure intervention and scheduled at months 0, 3, 6, 9, and 12. Outcomes contain participant-specific intercepts and slopes, serially correlated measurement error, and monotone dropout. No real people or clinical records are represented.

How to use this tutorial

Define the longitudinal estimand and time origin → audit people, visits, and missingness → visualize individual and mean trajectories → choose a mean and covariance model → estimate interpretable contrasts → diagnose fit and missing-data assumptions → report uncertainty and limits.

Learning objectives

After completing this tutorial, you should be able to:

• distinguish longitudinal, repeated cross-sectional, clustered, and simple pre–post data;
• move safely between long and wide layouts while preserving IDs and visit times;
• describe within-person dependence using random effects and residual covariance structures;
• explain why ordinary least-squares standard errors are generally invalid for repeated rows;
• distinguish population-averaged from subject-specific questions;
• fit independence-working marginal models, gls() covariance models, MMRM, and lme() models;
• interpret treatment main effects, treatment-by-time interactions, and nonlinear time terms;
• standardize model-matrix predictions into visit means and change contrasts with confidence intervals;
• distinguish fixed-part predictions from empirical Bayes/BLUP subject predictions;
• evaluate residuals, covariance assumptions, dropout, and sensitivity-analysis boundaries;
• decompose time-varying covariates into within- and between-person components; and
• plan prospective power and report a reproducible longitudinal analysis.

1 What makes data longitudinal?

1.1 The person and occasion are both analysis units

Longitudinal data repeatedly observe the same units through time. Rows from one participant share biology, history, environment, and measurement processes, so they are not independent replicates. The structure creates two kinds of information:

• between-person information: why people have different typical outcome levels;
• within-person information: how an outcome changes for the same person as time or exposure changes.

A repeated cross-sectional study samples different people at each wave. It can estimate population means over calendar time, but cannot identify individual change without linking the same people. A two-time-point pre–post study is longitudinal, although it offers little evidence about trajectory shape or serial covariance.

Design	Repeated people?	Typical target	Main dependence issue
Longitudinal cohort/trial	Yes	Mean change, trajectory, within-person association	Correlated outcomes and dropout
Repeated cross-section	No	Population mean at each wave	Sampling design, not person-level change
Clustered cross-section	Usually no	Group-adjusted association	Shared cluster context
Intensive longitudinal study	Yes, many occasions	Dynamics and short-term covariation	Dense serial dependence and time-varying confounding

Let $Y_{ij}$ be systolic blood pressure (SBP) for participant $i$ at time $t_{ij}$. A mean model alone might be

$$E(Y_{ij}\mid X_i,t_{ij})= \beta_0+\beta_1\text{Program}_i+\beta_2t_{ij} +\beta_3(\text{Program}_i t_{ij})+\beta_4t_{ij}^2+\beta_X^TX_i.$$

It does not yet say how repeated residuals correlate. Mean structure, covariance structure, target population, and missing-data assumptions are separate parts of the analysis.

Check your understanding: Do five visits from 360 people produce 1,800 independent observations?

No. There are up to 1,800 rows, but only 360 independently sampled or randomized participants. Repeated rows add information about trajectories and within-person variability; they do not create new independent treatment assignments.

2 Build and audit the longitudinal data

2.1 The simulated data-generating mechanism

Program is randomized once per person. The true average trajectory is mildly curved; the intervention has no baseline effect by construction and lowers the monthly slope by 0.55 mmHg. Each participant has a correlated random intercept and random slope. Measurement errors follow an AR(1) process with lag-one correlation 0.55.

Dropout is monotone: baseline is always present, but once a follow-up is missed, subsequent values are missing. While a participant remains active, retention depends on the previous observed SBP and program. The full outcomes are retained only to validate this simulation; real analysts never observe counterfactual missing values.

audit <- data.frame(
  Item = c(
    "Participants", "Scheduled rows", "Observed rows", "Scheduled visits",
    "Duplicate ID-time pairs", "Missing baseline outcomes",
    "Minimum observed visits per person", "Maximum observed visits per person"
  ),
  Value = c(
    nlevels(long[["id"]]), nrow(long), nrow(observed_long),
    length(unique(long[["time"]])),
    sum(duplicated(long[c("id", "time")])),
    sum(is.na(long[["sbp"]][long[["time"]] == 0])),
    min(as.vector(table(observed_long[["id"]]))),
    max(as.vector(table(observed_long[["id"]])))
  )
)
knitr::kable(audit, caption = "Structural and missingness audit of the simulated cohort")

Structural and missingness audit of the simulated cohort

Item	Value
Participants	360
Scheduled rows	1800
Observed rows	1642
Scheduled visits	5
Duplicate ID-time pairs	0
Missing baseline outcomes	0
Minimum observed visits per person	1
Maximum observed visits per person	5

Always check that one row means one person-occasion, time units are explicit, IDs do not recycle across sites, and visits are ordered within person. A numeric visit number is not automatically elapsed time: irregular schedules need the actual elapsed scale.

2.2 Visit retention and monotone missingness

retention <- aggregate(
  !is.na(sbp) ~ time + program, data = long, FUN = mean
)
names(retention)[3] <- "Retention"

pattern_matrix <- xtabs(~ id + visit, data = observed_long) > 0
pattern_text <- apply(pattern_matrix, 1, paste0, collapse = "")
pattern_table <- sort(table(pattern_text), decreasing = TRUE)

knitr::kable(retention, digits = 3,
  caption = "Observed proportion at each scheduled visit and in each randomized group")

Observed proportion at each scheduled visit and in each randomized group

time	program	Retention
0	Usual care	1.000
3	Usual care	0.939
6	Usual care	0.894
9	Usual care	0.878
12	Usual care	0.811
0	Intervention	1.000
3	Intervention	0.950
6	Intervention	0.922
9	Intervention	0.872
12	Intervention	0.856

knitr::kable(
  data.frame(Pattern = names(pattern_table), Participants = as.vector(pattern_table)),
  caption = "Visit patterns, where 1 is observed and 0 is missing"
)

Visit patterns, where 1 is observed and 0 is missing

Pattern	Participants
TRUETRUETRUETRUETRUE	300
TRUEFALSEFALSEFALSEFALSE	20
TRUETRUETRUETRUEFALSE	15
TRUETRUEFALSEFALSEFALSE	13
TRUETRUETRUEFALSEFALSE	12

The patterns should be interpreted in the scheduled order 0, 3, 6, 9, 12 months. Here they are monotone by design. Intermittent missingness would produce patterns such as 10111 and may reflect different processes.

Do not calculate retention only among people who return later At each visit, the denominator should be the people expected under the study protocol, with deaths, withdrawals, administrative censoring, and ineligibility described separately. Conditioning on later attendance can hide dropout.

2.3 Long and wide formats

Long format is usually best for regression: one row per person-visit, one outcome column, and an explicit time column. Wide format is useful for data entry, visit-specific summaries, and changes between named visits.

wide_sbp <- reshape(
  long[c("id", "time", "sbp")], idvar = "id", timevar = "time",
  direction = "wide"
)
wide_sbp <- wide_sbp[order(as.integer(wide_sbp[["id"]])), ]
names(wide_sbp) <- sub("sbp\\.", "sbp_month_", names(wide_sbp))

long_preview <- long[long[["id"]] %in% levels(long[["id"]])[1:3],
                     c("id", "program", "time", "sbp")]
knitr::kable(long_preview, digits = 2,
  caption = "Long-format records for the first three participants")

Long-format records for the first three participants

id	program	time	sbp
1	Intervention	0	146
1	Intervention	3	143
1	Intervention	6	141
1	Intervention	9	138
1	Intervention	12	135
2	Usual care	0	151
2	Usual care	3	149
2	Usual care	6	147
2	Usual care	9	143
2	Usual care	12	150
3	Usual care	0	151
3	Usual care	3	155
3	Usual care	6	149
3	Usual care	9	148
3	Usual care	12	156

knitr::kable(head(wide_sbp, 6), digits = 2,
  caption = "The same outcomes reshaped to one row per participant")

The same outcomes reshaped to one row per participant

	id	sbp_month_0	sbp_month_3	sbp_month_6	sbp_month_9	sbp_month_12
1	1	146	143	141	138	135
6	2	151	149	147	143	150
11	3	151	155	149	148	156
16	4	129	128	135	128	128
21	5	140	141	136	131	131
26	6	139	140	129	135	135

After every reshape, verify participant count, distinct ID-time count, expected column names, and the number of nonmissing outcomes. Never rely on row position alone to merge time-varying records.

3 Visualize before modeling

3.1 Individual trajectories

set.seed(20261007)
shown_ids <- sample(levels(long[["id"]]), 30)
plot_data <- observed_long[observed_long[["id"]] %in% shown_ids, ]
plot(
  NA, xlim = range(times),
  ylim = c(min(plot_data[["sbp"]]),
           max(plot_data[["sbp"]]) + 0.10 * diff(range(plot_data[["sbp"]]))),
  xlab = "Months since baseline", ylab = "Systolic blood pressure (mmHg)",
  main = "Observed individual trajectories", xaxt = "n"
)
axis(1, at = times)
for (person in shown_ids) {
  d <- plot_data[plot_data[["id"]] == person, ]
  col_i <- if (d[["program_num"]][1] == 1) palette_long["orange"] else palette_long["blue"]
  lines(d[["time"]], d[["sbp"]], col = grDevices::adjustcolor(col_i, 0.45), lwd = 1)
  points(d[["time"]], d[["sbp"]], col = col_i, pch = 16, cex = 0.45)
}
legend("top", c("Usual care", "Intervention"), horiz = TRUE,
       col = c(palette_long["blue"], palette_long["orange"]), lwd = 2, bty = "n")

Observed systolic blood-pressure trajectories for 30 randomly selected participants. Thin lines connect repeated measurements from the same participant; colors indicate randomized program. Individual intercepts, slopes, noise, and unequal follow-up are all visible.

A spaghetti plot reveals heterogeneity and missing tails, but hundreds of opaque lines can become ink without information. Show a reproducible subset, state how it was selected, and pair it with summaries of the full cohort.

3.2 Group mean trajectories

mean_rows <- do.call(rbind, lapply(split(observed_long, list(observed_long[["program"]], observed_long[["time"]])), function(d) {
  n_d <- nrow(d)
  m_d <- mean(d[["sbp"]])
  se_d <- sd(d[["sbp"]]) / sqrt(n_d)
  data.frame(program = d[["program"]][1], time = d[["time"]][1],
             n = n_d, mean = m_d, lower = m_d - qt(0.975, n_d - 1) * se_d,
             upper = m_d + qt(0.975, n_d - 1) * se_d)
}))
mean_rows <- mean_rows[order(mean_rows[["program"]], mean_rows[["time"]]), ]

mean_ylim <- range(mean_rows[c("lower", "upper")])
mean_ylim[2] <- mean_ylim[2] + 0.14 * diff(mean_ylim)
plot(
  NA, xlim = range(times), ylim = mean_ylim,
  xlab = "Months since baseline", ylab = "Observed mean SBP (mmHg)",
  main = "Available-case group means", xaxt = "n"
)
axis(1, at = times)
for (g in levels(long[["program"]])) {
  d <- mean_rows[mean_rows[["program"]] == g, ]
  col_g <- if (g == "Intervention") palette_long["orange"] else palette_long["blue"]
  arrows(d[["time"]], d[["lower"]], d[["time"]], d[["upper"]],
         angle = 90, code = 3, length = 0.04, col = col_g)
  lines(d[["time"]], d[["mean"]], col = col_g, lwd = 2)
  points(d[["time"]], d[["mean"]], col = col_g, pch = 16)
}
legend("top", levels(long[["program"]]), horiz = TRUE,
       col = c(palette_long["blue"], palette_long["orange"]),
       lwd = 2, pch = 16, bty = "n")

Observed mean systolic blood pressure by randomized group and visit with pointwise 95% t intervals. These descriptive means use available observations and therefore may change composition as participants drop out.

These intervals describe visit-specific means as if each visit were summarized separately. They do not adjust for repeated-measure covariance, baseline covariates, or informative changes in who remains observed.

3.3 Change scores are useful but incomplete

complete_change <- merge(
  long[long[["time"]] == 0, c("id", "program", "age_c", "female", "sbp")],
  long[long[["time"]] == 12 & !is.na(long[["sbp"]]), c("id", "sbp")],
  by = "id", suffixes = c("_0", "_12")
)
complete_change[["change_0_12"]] <- complete_change[["sbp_12"]] - complete_change[["sbp_0"]]
boxplot(
  change_0_12 ~ program, data = complete_change,
  col = c(grDevices::adjustcolor(palette_long["blue"], 0.25),
          grDevices::adjustcolor(palette_long["orange"], 0.25)),
  border = c(palette_long["blue"], palette_long["orange"]),
  xlab = "Randomized group", ylab = "Month 12 minus baseline SBP (mmHg)",
  main = "Complete-pair change distributions"
)
stripchart(change_0_12 ~ program, data = complete_change, vertical = TRUE,
           method = "jitter", add = TRUE, pch = 16, cex = 0.45,
           col = grDevices::adjustcolor(palette_long["gray"], 0.45))
abline(h = 0, lty = 2, col = palette_long["gray"])

Distribution of baseline-to-month-12 SBP change among participants observed at both visits. Negative values indicate reduced blood pressure. This complete-pair display is descriptive and excludes participants who dropped out before month 12.

Change scores answer a meaningful question, but one baseline and one endpoint discard intermediate trajectory information. Complete-pair analysis additionally changes the population to people with both values. A longitudinal model can use all observed visits under explicit assumptions.

4 Represent time deliberately

4.1 Continuous, categorical, polynomial, and spline time

Time coding determines the estimand and shape. Month 0 makes the program main effect a baseline difference and program:time a difference in monthly slope. Centering time at month 12 instead makes the program main effect the month-12 difference.

time_forms <- data.frame(
  Representation = c("Linear", "Quadratic", "Categorical", "Natural spline", "Piecewise linear"),
  R_syntax = c(
    "program * time", "program * time + I(time^2)",
    "program * visit_factor", "program * splines::ns(time, df = 3)",
    "program * (time + pmax(time - 6, 0))"
  ),
  Interpretation = c(
    "Constant slope and one slope difference",
    "Smooth curvature; hierarchy retains the linear term",
    "Separate mean at every scheduled visit",
    "Flexible smooth curve with basis-dependent coefficients",
    "Different slopes before and after a prespecified knot"
  )
)
knitr::kable(time_forms, caption = "Common representations of longitudinal time")

Common representations of longitudinal time

Representation	R_syntax	Interpretation
Linear	program * time	Constant slope and one slope difference
Quadratic	program * time + I(time^2)	Smooth curvature; hierarchy retains the linear term
Categorical	program * visit_factor	Separate mean at every scheduled visit
Natural spline	program * splines::ns(time, df = 3)	Flexible smooth curve with basis-dependent coefficients
Piecewise linear	program * (time + pmax(time - 6, 0))	Different slopes before and after a prespecified knot

Categorical time is robust for a small number of scheduled visits and supports visit-specific contrasts, but cannot directly interpolate. Continuous time is parsimonious and useful with irregular timing, but a wrong shape biases contrasts. Plot fitted trajectories on the outcome scale; do not interpret spline basis coefficients one by one.

4.2 Baseline belongs to the estimand strategy

For a randomized trial, common choices include:

1. ANCOVA/MMRM: analyze postbaseline outcomes and adjust for baseline outcome. This typically improves precision and permits treatment differences after baseline.
2. Constrained longitudinal data analysis: include baseline as an outcome and constrain randomized-group baseline means to equality.
3. Unconstrained repeated-outcome model: include baseline as an outcome with a group main effect, estimating the realized baseline difference.
4. Change-score analysis: model endpoint minus baseline, usually less flexible with multiple follow-ups.

Do not both condition on baseline as a covariate and also include the same baseline record as an ordinary response without deriving the model. In observational studies, baseline adjustment can change the scientific estimand and should follow a causal diagram rather than a ritual.

5 Dependence and covariance

5.1 Four useful covariance structures

For five visits, a residual covariance matrix can be modeled as:

• independence: all off-diagonal covariances are zero;
• compound symmetry (CS): one variance and one common correlation;
• AR(1): correlation declines as $\rho^{|j-k|}$ with visit lag;
• unstructured: every variance and covariance is estimated separately.

Random intercepts induce a CS-like persistent correlation. Random slopes produce correlations that vary with time. Residual AR(1) captures extra serial similarity after random effects.

cs_cor <- matrix(0.50, 5, 5); diag(cs_cor) <- 1
ar1_cor <- 0.55^abs(outer(0:4, 0:4, "-"))
rownames(cs_cor) <- colnames(cs_cor) <- paste0("M", times)
rownames(ar1_cor) <- colnames(ar1_cor) <- paste0("M", times)
knitr::kable(cs_cor, digits = 2, caption = "Compound-symmetry correlation example")

Compound-symmetry correlation example

	M0	M3	M6	M9	M12
M0	1.0	0.5	0.5	0.5	0.5
M3	0.5	1.0	0.5	0.5	0.5
M6	0.5	0.5	1.0	0.5	0.5
M9	0.5	0.5	0.5	1.0	0.5
M12	0.5	0.5	0.5	0.5	1.0

knitr::kable(ar1_cor, digits = 2, caption = "AR(1) correlation example")

AR(1) correlation example

	M0	M3	M6	M9	M12
M0	1.00	0.55	0.30	0.17	0.09
M3	0.55	1.00	0.55	0.30	0.17
M6	0.30	0.55	1.00	0.55	0.30
M9	0.17	0.30	0.55	1.00	0.55
M12	0.09	0.17	0.30	0.55	1.00

“Unstructured” is not assumption-free: it still assumes a shared covariance matrix, needs enough people and visit overlap, and can become unstable with many times. Information criteria, convergence, residual patterns, design knowledge, and sensitivity analyses should all inform covariance choice.

Check your understanding: Is AR(1) automatically appropriate whenever time is ordered?

No. Standard discrete AR(1) treats correlation as a function of visit lag and assumes equally spaced occasions. For irregular elapsed times, a continuous-time correlation such as corCAR1(form = ~ actual_time | id) may be more plausible. Random slopes can also explain correlation that AR(1) alone cannot.

6 Why naive ordinary regression fails

6.1 An independence-working mean model

Ordinary least squares can estimate the intended mean coefficients when its mean model and exogeneity assumptions are correct, but the usual covariance formula treats every residual as independent. That generally understates or otherwise distorts uncertainty.

working_lm <- lm(
  sbp ~ program * time + I(time^2) + age_c + female,
  data = observed_long
)
working_vcov_cluster <- cluster_sandwich(working_lm, observed_long[["id"]])

working_terms <- c("programIntervention", "time", "programIntervention:time")
working_cluster_se <- sqrt(diag(working_vcov_cluster))[working_terms]
working_cluster_df <- length(unique(observed_long[["id"]])) - 1
working_t_critical <- qt(0.975, df = working_cluster_df)
working_table <- data.frame(
  Term = working_terms,
  Estimate = coef(working_lm)[working_terms],
  `Naive OLS SE` = sqrt(diag(vcov(working_lm)))[working_terms],
  `Simple CR1 participant-cluster SE` = working_cluster_se,
  `Cluster df (G - 1)` = working_cluster_df,
  `CR1 t-based 95% CI lower` =
    coef(working_lm)[working_terms] - working_t_critical * working_cluster_se,
  `CR1 t-based 95% CI upper` =
    coef(working_lm)[working_terms] + working_t_critical * working_cluster_se,
  check.names = FALSE
)
knitr::kable(working_table, digits = 3,
  caption = "Naive OLS uncertainty versus a simple participant-cluster CR1 teaching calculation")

Naive OLS uncertainty versus a simple participant-cluster CR1 teaching calculation

	Term	Estimate	Naive OLS SE	Simple CR1 participant-cluster SE	Cluster df (G - 1)	CR1 t-based 95% CI lower	CR1 t-based 95% CI upper
programIntervention	programIntervention	0.052	0.701	0.937	359	-1.791	1.895
time	time	-0.253	0.177	0.093	359	-0.436	-0.070
programIntervention:time	programIntervention:time	-0.496	0.098	0.071	359	-0.635	-0.356

The sandwich calculation above is a transparent teaching implementation: it sums participant-level score contributions, applies a simple CR1 multiplier, and uses $G-1=359$ cluster degrees of freedom for the displayed $t$ intervals. It is not CR2 and has no leverage-specific adjustment. Production analyses should use a validated GEE or cluster-robust package with prespecified small-sample, leverage, and degrees-of-freedom methods. A sandwich estimator does not repair a misspecified mean model, informative dropout, confounding, or too few independent clusters.

6.2 Population-averaged and subject-specific questions

A marginal model describes the average outcome in a population, for example the average program difference at month 12. A mixed model describes outcomes conditional on latent participant random effects and can predict an observed person’s trajectory.

For Gaussian identity-link models, fixed-effect coefficients often have the same numerical conditional and marginal mean interpretation because random effects have mean zero. For logistic and other nonlinear links, conditional and marginal effects differ. Decide the target before selecting software.

7 Marginal covariance models with gls()

7.1 Compound symmetry versus AR(1)

nlme::gls() directly models a marginal covariance while keeping the same fixed mean. Here CS assumes one within-person correlation, whereas AR(1) allows correlation to decay across scheduled lags.

gls_cs <- nlme::gls(
  sbp ~ program * time + I(time^2) + age_c + female,
  correlation = nlme::corCompSymm(form = ~ 1 | id),
  data = observed_long, method = "REML", na.action = na.omit
)
gls_ar1 <- nlme::gls(
  sbp ~ program * time + I(time^2) + age_c + female,
  correlation = nlme::corAR1(form = ~ visit_index | id),
  data = observed_long, method = "REML", na.action = na.omit
)

gls_compare <- data.frame(
  Model = c("GLS with compound symmetry", "GLS with AR(1)"),
  Parameters = c(attr(logLik(gls_cs), "df"), attr(logLik(gls_ar1), "df")),
  AIC = c(AIC(gls_cs), AIC(gls_ar1)),
  BIC = c(BIC(gls_cs), BIC(gls_ar1)),
  Correlation = c(
    coef(gls_cs[["modelStruct"]][["corStruct"]], unconstrained = FALSE),
    coef(gls_ar1[["modelStruct"]][["corStruct"]], unconstrained = FALSE)
  )
)
knitr::kable(gls_compare, digits = 3,
  caption = "Marginal covariance models with the same fixed-effect structure")

Marginal covariance models with the same fixed-effect structure

	Model	Parameters	AIC	BIC	Correlation
Rho	GLS with compound symmetry	9	10095	10143	0.806
Phi	GLS with AR(1)	9	9869	9917	0.877

Because the fixed effects and data are identical, REML information criteria can compare these covariance structures. Compare fixed-effect structures with maximum likelihood, then refit the selected model by REML for final variance estimation. A lower AIC is evidence about relative predictive fit under the candidate set, not proof that the covariance mechanism is true.

7.2 Fixed-effect estimates under AR(1)

gls_tab <- summary(gls_ar1)[["tTable"]]
gls_results <- data.frame(
  Term = rownames(gls_tab), Estimate = gls_tab[, "Value"],
  `Standard error` = gls_tab[, "Std.Error"],
  `Large-sample 95% CI lower` =
    gls_tab[, "Value"] - large_sample_z * gls_tab[, "Std.Error"],
  `Large-sample 95% CI upper` =
    gls_tab[, "Value"] + large_sample_z * gls_tab[, "Std.Error"],
  `p-value` = gls_tab[, "p-value"], check.names = FALSE
)
knitr::kable(gls_results, digits = 3,
  caption = "AR(1) marginal GLS estimates for the simulated SBP trajectory")

AR(1) marginal GLS estimates for the simulated SBP trajectory

	Term	Estimate	Standard error	Large-sample 95% CI lower	Large-sample 95% CI upper	p-value
(Intercept)	(Intercept)	143.369	0.788	141.824	144.915	0.000
programIntervention	programIntervention	-0.105	0.904	-1.876	1.666	0.907
time	time	-0.253	0.087	-0.424	-0.082	0.004
I(time^2)	I(time^2)	0.016	0.006	0.004	0.028	0.008
age_c	age_c	1.212	0.403	0.423	2.001	0.003
femaleYes	femaleYes	-2.559	0.825	-4.177	-0.941	0.002
programIntervention:time	programIntervention:time	-0.466	0.072	-0.607	-0.326	0.000

With time coded in months from baseline:

• programIntervention is the adjusted randomized-group difference at month 0;
• time is the usual-care instantaneous linear component at month 0;
• I(time^2) captures shared curvature;
• programIntervention:time is the intervention-minus-usual-care difference in monthly linear slope.

The total derivative for usual care is $\beta_{time}+2\beta_{time^2}t$; the intervention derivative adds $\beta_{program\times time}$. Report contrasts at meaningful visits rather than relying only on component coefficients.

8 MMRM with categorical visits

8.1 Postbaseline outcomes adjusted for baseline SBP

A common trial MMRM treats postbaseline visit as categorical, includes program-by-visit interactions, adjusts for baseline, and uses an unstructured within-person covariance with visit-specific variances. It uses no random effects; the covariance is marginal.

baseline_lookup <- long[
  long[["time"]] == 0, c("id", "sbp", "age_c", "female")
]
names(baseline_lookup)[2] <- "baseline_sbp"
baseline_lookup[["baseline_centered"]] <-
  baseline_lookup[["baseline_sbp"]] - mean(baseline_lookup[["baseline_sbp"]])
post <- merge(
  long[long[["time"]] > 0 & !is.na(long[["sbp"]]), ],
  baseline_lookup[c("id", "baseline_sbp", "baseline_centered")],
  by = "id", sort = FALSE
)
post <- post[order(as.integer(post[["id"]]), post[["time"]]), ]
post[["post_visit"]] <- factor(post[["time"]], levels = c(3, 6, 9, 12))
post[["post_index"]] <- match(post[["time"]], c(3, 6, 9, 12))

mmrm_fit <- nlme::gls(
  sbp ~ program * post_visit + baseline_centered + age_c + female,
  correlation = nlme::corSymm(form = ~ post_index | id),
  weights = nlme::varIdent(form = ~ 1 | post_visit),
  data = post, method = "REML", na.action = na.omit,
  control = nlme::glsControl(maxIter = 200, msMaxIter = 300)
)

mmrm_info <- data.frame(
  Quantity = c("Observed postbaseline rows", "Participants represented", "AIC", "Residual SD at reference visit"),
  Value = c(nrow(post), nlevels(droplevels(post[["id"]])), AIC(mmrm_fit), sigma(mmrm_fit))
)
knitr::kable(mmrm_info, digits = 3,
  caption = "Fitted categorical-visit MMRM with unstructured correlation and visit-specific variance")

Fitted categorical-visit MMRM with unstructured correlation and visit-specific variance

Quantity	Value
Observed postbaseline rows	1282.00
Participants represented	340.00
AIC	7290.87
Residual SD at reference visit	4.07

This likelihood analysis is valid under its mean/covariance model and a missing at random (MAR) assumption conditional on variables included in the likelihood. MAR is not guaranteed by using MMRM. Include strong predictors of missingness and outcome when scientifically justified, and plan sensitivity analyses for plausible departures.

8.2 Standardized MMRM program differences

Rather than report reference-person coefficients, average design rows over the study’s baseline covariate distribution. This is model-matrix standardization.

mmrm_beta <- coef(mmrm_fit)
mmrm_V <- vcov(mmrm_fit)
mmrm_formula <- ~ program * post_visit + baseline_centered + age_c + female

mmrm_design <- function(program_level, visit_value) {
  base_people <- baseline_lookup
  nd <- data.frame(
    program = factor(program_level, levels = levels(long[["program"]])),
    post_visit = factor(visit_value, levels = c(3, 6, 9, 12)),
    baseline_centered = base_people[["baseline_centered"]],
    age_c = base_people[["age_c"]], female = base_people[["female"]]
  )
  X <- model.matrix(mmrm_formula, nd)
  colMeans(X[, names(mmrm_beta), drop = FALSE])
}

mmrm_contrasts <- do.call(rbind, lapply(c(3, 6, 9, 12), function(month) {
  L <- mmrm_design("Intervention", month) - mmrm_design("Usual care", month)
  linear_estimate(L, mmrm_beta, mmrm_V, paste0("Intervention minus usual care at month ", month))
}))
knitr::kable(mmrm_contrasts, digits = 3,
  caption = "Covariate-standardized MMRM program differences by postbaseline visit")

Covariate-standardized MMRM program differences by postbaseline visit

Contrast	Estimate	Standard error	Large-sample 95% CI lower	Large-sample 95% CI upper
Intervention minus usual care at month 3	-1.41	0.442	-2.27	-0.543
Intervention minus usual care at month 6	-2.50	0.579	-3.64	-1.370
Intervention minus usual care at month 9	-4.39	0.633	-5.63	-3.151
Intervention minus usual care at month 12	-5.52	0.705	-6.90	-4.134

Pointwise intervals answer four separate questions. If all visits are confirmatory, prespecify multiplicity control or a joint test. A visit-specific $p$-value pattern is not evidence that effects suddenly appear or disappear between visits.

9 Linear mixed-effects trajectories

9.1 Random intercept, random slope, and residual AR(1)

A subject-specific model adds random intercept $b_{0i}$ and slope $b_{1i}$:

$$Y_{ij}=X_{ij}^T\beta+b_{0i}+b_{1i}t_{ij}+\varepsilon_{ij}, \qquad \begin{pmatrix}b_{0i}\\b_{1i}\end{pmatrix}\sim N(0,D), \quad \operatorname{Cor}(\varepsilon_{ij},\varepsilon_{ik})=\rho^{|j-k|}.$$

Random effects explain stable trajectory heterogeneity; AR(1) residuals explain remaining short-range serial dependence. Avoid adding both automatically: the data must support their distinct roles.

lme_ri_ml <- nlme::lme(
  sbp ~ program * time + I(time^2) + age_c + female,
  random = ~ 1 | id, data = observed_long, method = "ML",
  na.action = na.omit, control = nlme::lmeControl(opt = "optim")
)
lme_rs_ml <- nlme::lme(
  sbp ~ program * time + I(time^2) + age_c + female,
  random = ~ time | id,
  correlation = nlme::corAR1(form = ~ visit_index | id),
  data = observed_long, method = "ML", na.action = na.omit,
  control = nlme::lmeControl(opt = "optim")
)
lme_fit <- nlme::lme(
  sbp ~ program * time + I(time^2) + age_c + female,
  random = ~ time | id,
  correlation = nlme::corAR1(form = ~ visit_index | id),
  data = observed_long, method = "REML", na.action = na.omit,
  control = nlme::lmeControl(opt = "optim")
)

lme_fit_audit <- function(fit) {
  ap_var <- fit[["apVar"]]
  data.frame(
    `Returned without convergence error` = inherits(fit, "lme"),
    `Finite log likelihood` = is.finite(as.numeric(logLik(fit))),
    `Finite fixed effects` = all(is.finite(nlme::fixef(fit))),
    `Finite approximate variance (apVar)` =
      is.matrix(ap_var) && length(ap_var) > 0 && all(is.finite(ap_var)),
    Iterations = if (is.null(fit[["numIter"]])) NA_integer_ else fit[["numIter"]],
    check.names = FALSE
  )
}
lme_audit <- rbind(lme_fit_audit(lme_ri_ml), lme_fit_audit(lme_rs_ml))
lme_compare <- data.frame(
  Model = c("Random intercept", "Random intercept + slope + residual AR(1)"),
  AIC_ML = c(AIC(lme_ri_ml), AIC(lme_rs_ml)),
  BIC_ML = c(BIC(lme_ri_ml), BIC(lme_rs_ml)),
  Parameters = c(attr(logLik(lme_ri_ml), "df"), attr(logLik(lme_rs_ml), "df")),
  lme_audit,
  check.names = FALSE
)
knitr::kable(lme_compare, digits = 3,
  caption = "ML comparison and numerical audit of two candidate dependence structures")

ML comparison and numerical audit of two candidate dependence structures

Model	AIC_ML	BIC_ML	Parameters	Returned without convergence error	Finite log likelihood	Finite fixed effects	Finite approximate variance (apVar)	Iterations
Random intercept	10079	10128	9	TRUE	TRUE	TRUE	TRUE	NA
Random intercept + slope + residual AR(1)	9838	9903	12	TRUE	TRUE	TRUE	TRUE	NA

The second candidate changes two components at once: it adds a random slope and residual AR(1) correlation. Its AIC difference therefore cannot be attributed to either component alone. Fit random-intercept-plus-AR(1) and random-slope-without-AR(1) candidates when that separation matters, and judge all candidates using design knowledge, parameter estimates, diagnostics, stability, and sensitivity—not AIC alone. A successful lme() return plus finite likelihood, fixed effects, and apVar is a practical numerical audit, not proof of a unique optimum or correct covariance model.

Use ML, not REML, when comparing fixed-effect structures. REML likelihoods depend on the fixed-effect design and are not comparable when fixed effects differ. Boundary likelihood-ratio tests for random effects need special calibration; AIC and scientific plausibility are useful but not sufficient.

9.2 Fixed effects and interpretation

lme_tab <- summary(lme_fit)[["tTable"]]
lme_results <- data.frame(
  Term = rownames(lme_tab), Estimate = lme_tab[, "Value"],
  `Standard error` = lme_tab[, "Std.Error"],
  `Large-sample 95% CI lower` =
    lme_tab[, "Value"] - large_sample_z * lme_tab[, "Std.Error"],
  `Large-sample 95% CI upper` =
    lme_tab[, "Value"] + large_sample_z * lme_tab[, "Std.Error"],
  `p-value` = lme_tab[, "p-value"], check.names = FALSE
)
knitr::kable(lme_results, digits = 3,
  caption = "REML linear mixed-model fixed effects")

REML linear mixed-model fixed effects

	Term	Estimate	Standard error	Large-sample 95% CI lower	Large-sample 95% CI upper	p-value
(Intercept)	(Intercept)	143.324	0.809	141.740	144.909	0.000
programIntervention	programIntervention	-0.090	0.928	-1.909	1.730	0.923
time	time	-0.257	0.082	-0.417	-0.096	0.002
I(time^2)	I(time^2)	0.016	0.006	0.004	0.027	0.007
age_c	age_c	1.185	0.408	0.385	1.985	0.004
femaleYes	femaleYes	-2.463	0.837	-4.103	-0.823	0.003
programIntervention:time	programIntervention:time	-0.465	0.064	-0.590	-0.340	0.000

The program main effect is the month-0 difference; randomization implies its population value is zero, although the realized sample need not be exactly balanced. The program-by-time coefficient is the adjusted difference in monthly slope. The quadratic coefficient is common to both groups because no program-by-quadratic interaction was specified.

Fixed effects describe the average trajectory conditional on measured covariates and with random effects set to their mean of zero. Random slopes describe unexplained person-to-person deviations; they are not estimates of a causal treatment-effect distribution.

10 Turn coefficients into useful estimands

10.1 Standardized group means at every visit

We average model-matrix rows over all 360 participants’ age and sex distribution, assigning everyone in turn to each program. Because program was randomized, this targets the study population without relying on an arbitrary reference age or sex.

lme_beta <- nlme::fixef(lme_fit)
lme_V <- vcov(lme_fit)
lme_formula <- ~ program * time + I(time^2) + age_c + female
people <- unique(long[c("id", "age_c", "female")])

lme_design <- function(program_level, month) {
  nd <- data.frame(
    program = factor(program_level, levels = levels(long[["program"]])),
    time = month, age_c = people[["age_c"]], female = people[["female"]]
  )
  X <- model.matrix(lme_formula, nd)
  colMeans(X[, names(lme_beta), drop = FALSE])
}

standardized_means <- do.call(rbind, lapply(levels(long[["program"]]), function(g) {
  do.call(rbind, lapply(times, function(month) {
    L <- lme_design(g, month)
    out <- linear_estimate(L, lme_beta, lme_V, paste(g, "at month", month))
    data.frame(Program = g, Month = month, out[-1], check.names = FALSE)
  }))
}))
knitr::kable(standardized_means, digits = 3,
  caption = "Covariate-standardized mixed-model mean SBP at each scheduled visit")

Covariate-standardized mixed-model mean SBP at each scheduled visit

Program	Month	Estimate	Standard error	Large-sample 95% CI lower	Large-sample 95% CI upper
Usual care	0	142	0.659	141	143
Usual care	3	141	0.613	140	142
Usual care	6	141	0.601	140	142
Usual care	9	141	0.600	140	142
Usual care	12	141	0.639	140	142
Intervention	0	142	0.659	140	143
Intervention	3	140	0.613	139	141
Intervention	6	138	0.600	137	139
Intervention	9	137	0.598	135	138
Intervention	12	135	0.635	134	137

These are fixed-part population mean predictions. Adding random effects would instead predict particular participants represented in the fitted data.

10.2 Month-12 difference and changes from baseline

L_u0 <- lme_design("Usual care", 0)
L_u12 <- lme_design("Usual care", 12)
L_i0 <- lme_design("Intervention", 0)
L_i12 <- lme_design("Intervention", 12)

target_contrasts <- rbind(
  linear_estimate(L_i12 - L_u12, lme_beta, lme_V,
                  "Month 12: intervention minus usual care"),
  linear_estimate(L_u12 - L_u0, lme_beta, lme_V,
                  "Usual care: month 12 minus baseline"),
  linear_estimate(L_i12 - L_i0, lme_beta, lme_V,
                  "Intervention: month 12 minus baseline"),
  linear_estimate((L_i12 - L_i0) - (L_u12 - L_u0), lme_beta, lme_V,
                  "Difference in baseline-to-month-12 changes")
)
knitr::kable(target_contrasts, digits = 3,
  caption = "Prespecified visit and change contrasts from the full coefficient covariance matrix")

Prespecified visit and change contrasts from the full coefficient covariance matrix

Contrast	Estimate	Standard error	Large-sample 95% CI lower	Large-sample 95% CI upper
Month 12: intervention minus usual care	-5.672	0.894	-7.42	-3.920
Usual care: month 12 minus baseline	-0.801	0.544	-1.87	0.264
Intervention: month 12 minus baseline	-6.383	0.538	-7.44	-5.328
Difference in baseline-to-month-12 changes	-5.582	0.765	-7.08	-4.083

Every interval uses $L\widehat V L^T$, preserving covariance among coefficients. Subtracting two independently reported standard errors would be wrong. The reusable linear_estimate() function and manually constructed intervals on this page use a large-sample Wald normal reference (qnorm(0.975)); they do not inherit nlme denominator degrees of freedom. For studies with fewer independent people or complex covariance estimation, prespecify a validated finite-sample approach—for example, an appropriate nlme approximate $t$ procedure or Kenward–Roger/Satterthwaite inference in software that supports the fitted model—rather than choosing a method after seeing results. In this model the difference-in-changes equals 12 times the program-by-time coefficient because curvature is shared, but explicit contrast code remains auditable if the model changes.

10.3 Display the standardized trajectories

grid_month <- seq(0, 12, length.out = 121)
curve_rows <- do.call(rbind, lapply(levels(long[["program"]]), function(g) {
  do.call(rbind, lapply(grid_month, function(month) {
    L <- lme_design(g, month)
    est <- drop(L %*% lme_beta)
    se <- sqrt(drop(L %*% lme_V %*% L))
    data.frame(program = g, time = month, fit = est,
               lower = est - large_sample_z * se,
               upper = est + large_sample_z * se)
  }))
}))
curve_ylim <- range(curve_rows[c("lower", "upper")])
curve_ylim[2] <- curve_ylim[2] + 0.14 * diff(curve_ylim)
plot(
  NA, xlim = c(0, 12), ylim = curve_ylim,
  xlab = "Months since baseline", ylab = "Standardized mean SBP (mmHg)",
  main = "Model-based mean trajectories", xaxt = "n"
)
axis(1, at = times)
for (g in levels(long[["program"]])) {
  d <- curve_rows[curve_rows[["program"]] == g, ]
  col_g <- if (g == "Intervention") palette_long["orange"] else palette_long["blue"]
  polygon(c(d[["time"]], rev(d[["time"]])), c(d[["lower"]], rev(d[["upper"]])),
          col = grDevices::adjustcolor(col_g, 0.16), border = NA)
  lines(d[["time"]], d[["fit"]], col = col_g, lwd = 2.5)
}
legend("topright", levels(long[["program"]]), inset = 0.02,
       col = c(palette_long["blue"], palette_long["orange"]),
       lwd = 2.5, bty = "n", cex = 0.88)

Covariate-standardized fixed-part SBP trajectories from the random-intercept/random-slope model. Shaded pointwise 95% confidence bands reflect fixed-effect estimation uncertainty, not the much wider distribution of individual future outcomes.

11 Subject-specific prediction and covariance

11.1 Fixed-part versus BLUP predictions

predict(..., level = 0) uses only fixed effects. level = 1 adds empirical Bayes estimates of each observed participant’s random intercept and slope. These BLUPs partially pool noisy individual trajectories toward the population distribution.

complete_ids <- names(which(table(observed_long[["id"]]) == 5))[1:9]
pred_data <- observed_long[observed_long[["id"]] %in% complete_ids, ]
pred_data[["fixed_prediction"]] <- predict(lme_fit, newdata = pred_data, level = 0)
pred_data[["subject_prediction"]] <- predict(lme_fit, newdata = pred_data, level = 1)

old_par <- par(mfrow = c(3, 3), mar = c(2.4, 2.2, 2.2, 0.6),
               oma = c(4.2, 4.4, 1, 0))
for (person in complete_ids) {
  d <- pred_data[pred_data[["id"]] == person, ]
  plot(d[["time"]], d[["sbp"]], pch = 16, col = palette_long["navy"],
       xlim = c(0, 12), ylim = range(pred_data[["sbp"]]), xaxt = "n",
       xlab = "", ylab = "", main = paste("Participant", person), cex = 0.65)
  axis(1, at = times, cex.axis = 0.75)
  lines(d[["time"]], d[["fixed_prediction"]], lty = 2, lwd = 1.7,
        col = palette_long["gray"])
  lines(d[["time"]], d[["subject_prediction"]], lwd = 2,
        col = palette_long["teal"])
}
mtext("Month", side = 1, outer = TRUE, line = 2.2)
mtext("Systolic blood pressure (mmHg)", side = 2, outer = TRUE, line = 2.4)

Observed, fixed-part, and participant-specific BLUP trajectories for nine complete participants. The shared fixed-part curves depend on program and covariates; BLUP curves adapt to each participant while shrinking toward the population model.

par(old_par)

BLUPs are estimates, not directly observed personal traits. Ranking participants by BLUP can be unstable, especially with few visits. Prediction for a new participant integrates over unknown random effects and needs a prediction interval containing random-effect and residual variation, not merely a fixed-effect confidence band.

11.2 Time-varying ICC and within-person correlation

With random slopes, the random-effect variance at time $t$ is $z_t^TDz_t$ for $z_t=(1,t)^T$. The proportion of marginal variance attributable to random effects therefore changes with time. Residual AR(1) adds lag-specific covariance.

D_hat <- as.matrix(nlme::getVarCov(lme_fit, type = "random.effects"))
sigma2_hat <- sigma(lme_fit)^2
rho_hat <- as.numeric(coef(
  lme_fit[["modelStruct"]][["corStruct"]], unconstrained = FALSE
))

random_variance <- function(t) drop(c(1, t) %*% D_hat %*% c(1, t))
total_variance <- function(t) random_variance(t) + sigma2_hat
random_share <- vapply(times, function(t) random_variance(t) / total_variance(t), numeric(1))
cor_with_baseline <- vapply(seq_along(times), function(k) {
  t <- times[k]
  re_cov <- drop(c(1, 0) %*% D_hat %*% c(1, t))
  residual_cov <- sigma2_hat * rho_hat^(k - 1)
  (re_cov + residual_cov) / sqrt(total_variance(0) * total_variance(t))
}, numeric(1))

old_par <- par(mfrow = c(1, 2), mar = c(4.2, 4.2, 3, 1))
plot(times, random_share, type = "b", pch = 16, lwd = 2,
     col = palette_long["teal"], ylim = c(0, 1), xaxt = "n",
     xlab = "Month", ylab = "Random-effect variance / total variance",
     main = "Time-varying variance share")
axis(1, at = times)
plot(times, cor_with_baseline, type = "b", pch = 16, lwd = 2,
     col = palette_long["purple"], ylim = c(0, 1), xaxt = "n",
     xlab = "Month", ylab = "Model-implied correlation with baseline",
     main = "Correlation with month 0")
axis(1, at = times)

Model-implied random-effect variance proportion and correlation with baseline across follow-up time. Random slopes make between-person heterogeneity change with time, while residual AR(1) correlation decays with visit separation.

par(old_par)

Calling the first quantity an “ICC” is common, but in longitudinal random-slope models there is no single ICC. The correlation between two measurements also depends on both times, not only one variance share.

12 Diagnostics and sensitivity

12.1 Examine the mean, tails, variance, and people

diag_data <- observed_long
diag_data[["fitted_conditional"]] <- fitted(lme_fit, level = 1)
diag_data[["resid_standard"]] <- residuals(lme_fit, type = "normalized")
person_resid <- aggregate(abs(resid_standard) ~ id, diag_data, mean)
person_resid <- person_resid[order(person_resid[["abs(resid_standard)"]]), ]

old_par <- par(mfrow = c(2, 2), mar = c(4.1, 4.3, 3, 1))
plot(diag_data[["fitted_conditional"]], diag_data[["resid_standard"]],
     pch = 16, cex = 0.4,
     col = grDevices::adjustcolor(palette_long["navy"], 0.35),
     xlab = "Conditional fitted SBP", ylab = "Normalized residual",
     main = "Residuals versus fitted")
abline(h = 0, lty = 2, col = palette_long["vermillion"])
qqnorm(diag_data[["resid_standard"]], pch = 16, cex = 0.4,
       col = grDevices::adjustcolor(palette_long["navy"], 0.4), main = "Normal Q-Q")
qqline(diag_data[["resid_standard"]], col = palette_long["vermillion"], lwd = 2)
boxplot(resid_standard ~ visit, data = diag_data, las = 2,
        col = grDevices::adjustcolor(palette_long["sky"], 0.3),
        xlab = "", ylab = "Normalized residual", main = "Residuals by visit")
plot(seq_len(nrow(person_resid)), person_resid[["abs(resid_standard)"]],
     pch = 16, cex = 0.45, col = palette_long["teal"],
     xlab = "Participants ordered by residual magnitude",
     ylab = "Mean absolute normalized residual", main = "Person-level residual screen")

Four diagnostics for the mixed-effects trajectory model: conditional residuals versus fitted values, a normal Q-Q plot, residual distributions by visit, and participant-level mean absolute standardized residuals. These panels assess mean shape, tails, visit-specific variance, and influential people.

par(old_par)

Also inspect random-effect distributions, fitted covariance parameters near boundaries, convergence messages, sensitivity to covariance structure, and influence by deleting whole participants—not individual rows. Residual normality matters differently for mean estimates, small-sample inference, and prediction; it should not become an automatic hypothesis-test ritual.

12.2 A diagnostic checklist

• Plot observed and fitted trajectories, not only coefficient tables.
• Check whether curvature or treatment-specific curvature is omitted.
• Compare residual spread across visits and groups.
• Assess normalized residual tails and outlying participant trajectories.
• Inspect random-effect variances, correlation, and near-singular covariance.
• Verify time ordering and the scale assumed by AR(1).
• Refit scientifically plausible covariance structures as sensitivity analyses.
• Report convergence controls and any model-selection decisions.

13 Missing outcomes and dropout

13.1 MCAR, MAR, and MNAR are assumptions about distributions

• MCAR: missingness is independent of observed and unobserved outcomes and covariates. This is strong and rarely established by a test.
• MAR: conditional on included observed history, missingness does not additionally depend on the missing outcome.
• MNAR: missingness still depends on unobserved outcomes after conditioning; sensitivity assumptions are required.

Likelihood-based GLS/MMRM/mixed models can use all observed outcomes under MAR when their outcome and covariance models are correct. They do not impute a literal value into every empty cell, and MAR is not a property conferred by software.

dropout_summary <- do.call(rbind, lapply(levels(long[["id"]]), function(person) {
  d <- long[long[["id"]] == person, ]
  first_missing <- which(is.na(d[["sbp"]]))[1]
  data.frame(
    id = person, program = d[["program"]][1], baseline_sbp = d[["sbp"]][1],
    observed_visits = sum(!is.na(d[["sbp"]])),
    last_observed_month = max(d[["time"]][!is.na(d[["sbp"]])]),
    dropout_month = if (is.na(first_missing)) "Completed" else as.character(d[["time"]][first_missing])
  )
}))
dropout_by_group <- aggregate(
  cbind(observed_visits, baseline_sbp) ~ program, dropout_summary, mean
)
knitr::kable(dropout_by_group, digits = 2,
  caption = "Baseline SBP and observed-visit count by randomized group")

Baseline SBP and observed-visit count by randomized group

program	observed_visits	baseline_sbp
Usual care	4.52	142
Intervention	4.60	142

knitr::kable(with(dropout_summary, addmargins(table(program, dropout_month))),
  caption = "First missed visit, or completion, by randomized group")

First missed visit, or completion, by randomized group

	12	3	6	9	Completed	Sum
Usual care	12	11	8	3	146	180
Intervention	3	9	5	9	154	180
Sum	15	20	13	12	300	360

Record reasons and timing of missingness. Death, treatment discontinuation, missed measurement, and administrative loss can correspond to different estimands. A treatment-policy estimand may continue outcomes after treatment discontinuation; a while-on-treatment estimand answers another question.

13.2 Why complete-case analysis changes the target sample

complete_fit <- lm(
  change_0_12 ~ program + age_c + female,
  data = complete_change
)
complete_term <- summary(complete_fit)[["coefficients"]]["programIntervention", ]
mixed_change <- target_contrasts[target_contrasts[["Contrast"]] ==
                                   "Difference in baseline-to-month-12 changes", ]
comparison_missing <- data.frame(
  Analysis = c("Complete-pair change score", "Mixed model using all observed visits"),
  Participants_or_records = c(nrow(complete_change), nrow(observed_long)),
  Estimate = c(complete_term["Estimate"], mixed_change[["Estimate"]]),
  `Standard error` = c(complete_term["Std. Error"], mixed_change[["Standard error"]]),
  check.names = FALSE
)
knitr::kable(comparison_missing, digits = 3,
  caption = "Complete-pair and likelihood analyses target data through different assumptions")

Complete-pair and likelihood analyses target data through different assumptions

	Analysis	Participants_or_records	Estimate	Standard error
Estimate	Complete-pair change score	300	-5.17	0.788
	Mixed model using all observed visits	1642	-5.58	0.765

The numerical estimates need not be dramatically different in every simulation. The problem is structural: complete-case analysis conditions on remaining observed through month 12 and is unbiased only under restrictive missingness conditions.

13.3 Tools and their boundaries

Strategy	What it can do	Important boundary
Likelihood under MAR	Use incomplete outcome vectors without single imputation	Relies on outcome/covariance model and MAR
Multiple imputation	Include auxiliary variables and propagate imputation uncertainty	Imputation must honor multilevel/time structure and analysis compatibility
Inverse-probability weighting	Reweight observed histories using modeled attendance probabilities	Sensitive to positivity, extreme weights, and missingness-model error
Pattern-mixture sensitivity	Shift imputed or modeled missing outcomes by plausible deltas	Delta values require clinical justification
Selection or joint model	Couple dropout/survival and outcome processes	Strong unverifiable structure; not automatically superior

MNAR cannot be identified from observed data alone without assumptions or external information. Prespecify a plausible range, show how conclusions change, and separate primary MAR analysis from sensitivity analyses.

14 Time-varying covariates

14.1 Separate within-person and between-person associations

Suppose weekly activity is measured at every visit. Entering raw activity alone mixes whether more-active people have lower SBP with whether a person’s more-active-than-usual visit coincides with lower SBP.

# Construct the deterministic teaching covariate only on observed records,
# using the observed concurrent outcome; it does not change the shared DGM.
activity_observed <- observed_long
activity_observed[["activity_minutes"]] <- with(
  activity_observed, 35 + 3 * program_num - 0.25 * time - 0.18 * (sbp - 140) +
    4 * sin(as.integer(id) * 0.7 + time)
)
activity_mean <- ave(
  activity_observed[["activity_minutes"]], activity_observed[["id"]], FUN = mean
)
activity_observed[["activity_between"]] <- activity_mean - mean(activity_mean)
activity_observed[["activity_within"]] <-
  activity_observed[["activity_minutes"]] - activity_mean

activity_fit <- nlme::lme(
  sbp ~ program * time + I(time^2) + age_c + female +
    activity_between + activity_within,
  random = ~ time | id, data = activity_observed, method = "REML",
  na.action = na.omit, control = nlme::lmeControl(opt = "optim")
)
activity_terms <- summary(activity_fit)[["tTable"]][c("activity_between", "activity_within"), ]
activity_results <- data.frame(
  Component = c("Between-person activity mean", "Within-person deviation from own mean"),
  Estimate = activity_terms[, "Value"], `Standard error` = activity_terms[, "Std.Error"],
  check.names = FALSE
)
knitr::kable(activity_results, digits = 3,
  caption = "Within-between decomposition of a simulated time-varying covariate")

Within-between decomposition of a simulated time-varying covariate

	Component	Estimate	Standard error
activity_between	Between-person activity mean	-4.610	0.110
activity_within	Within-person deviation from own mean	-0.122	0.027

The within coefficient is an adjusted association, not automatically the effect of increasing activity. Here the teaching variable is deliberately constructed from concurrently observed SBP, making its endogeneity explicit; no unobserved postdropout SBP is used. In real studies, time-varying activity may respond to prior SBP, symptoms, or treatment. If treatment changes activity, adjusting for it in the primary treatment model changes the total-effect estimand and may introduce postrandomization bias. Mediation or time-varying confounding questions need causal methods and explicit assumptions.

15 Outcomes beyond Gaussian SBP

15.1 Binary and count repeated outcomes

The dependence problem persists when outcomes change type:

Outcome	Population-averaged approach	Subject-specific approach	Typical effect scale
Binary symptom indicator	GEE with binomial family	Logistic GLMM	Marginal or conditional odds ratio; preferably probabilities/risks too
Count of events per interval	Poisson/negative-binomial GEE with offset	Count GLMM with offset	Marginal or conditional rate ratio
Ordinal severity	Ordinal marginal model	Ordinal mixed model	Cumulative odds or standardized probabilities
Time until recurrent event	Recurrent-event marginal model	Frailty model	Rate, hazard, or mean cumulative function contrast

For a binary GLMM one might write lme4::glmer(event ~ program * time + (time | id), family = binomial). This page does not execute that model, keeping dependencies to nlme. With nonlinear links, conditional GLMM coefficients are not population-average effects. Standardize predicted probabilities or integrate over random effects when the estimand is marginal.

Overdispersion, zero inflation, exposure time, and event dependence require separate checks for counts. An offset is a known log exposure coefficient of one, not an ordinary covariate.

16 Prospective power and precision

16.1 Power belongs before outcomes are observed

Longitudinal power depends on the target contrast, visit schedule, covariance, heterogeneity, dropout, allocation, and analysis. More visits help only to the extent that they add information rather than highly redundant measurements.

# Rough endpoint-change approximation for planning discussion only.
var_change_0_12 <- 12^2 * sd_slope^2 +
  2 * residual_sd^2 * (1 - rho^4)
sd_change_0_12 <- sqrt(var_change_0_12)
planned_effect <- 0.55 * 12

power_rows <- do.call(rbind, lapply(c(100, 140, 180), function(per_group) {
  do.call(rbind, lapply(c(0, 0.10, 0.20), function(attrition) {
    retained <- floor(per_group * (1 - attrition))
    pwr <- power.t.test(
      n = retained, delta = planned_effect, sd = sd_change_0_12,
      sig.level = 0.05, type = "two.sample", alternative = "two.sided"
    )[["power"]]
    data.frame(Planned_per_group = per_group, Attrition = attrition,
               Approximate_analyzed_per_group = retained, Approximate_power = pwr)
  }))
}))
knitr::kable(power_rows, digits = 3,
  caption = "Rough prospective endpoint-change power across enrollment and attrition scenarios")

Rough prospective endpoint-change power across enrollment and attrition scenarios

Planned_per_group	Attrition	Approximate_analyzed_per_group	Approximate_power
100	0.0	100	1
100	0.1	90	1
100	0.2	80	1
140	0.0	140	1
140	0.1	126	1
140	0.2	112	1
180	0.0	180	1
180	0.1	162	1
180	0.2	144	1

This approximation ignores intermediate visits, nonlinear time, covariance estimation, MAR likelihood, and unequal dropout. A defensible plan simulates the complete prespecified design and primary analysis across plausible covariance and missingness scenarios. Use clinically meaningful effects, not the effect observed in the completed study. “Observed power” is a transformation of the estimate and standard error and adds no evidence beyond them.

17 A reproducible analysis and reporting workflow

17.1 Selected results from the simulated cohort

case_summary <- data.frame(
  Result = c(
    "Randomized participants", "Scheduled visits", "Observed outcome records",
    "Month-12 retention", "Mixed-model program-by-time estimate",
    "Mixed-model program-by-time SE", "AR(1) residual correlation",
    "Difference in baseline-to-month-12 changes"
  ),
  Value = c(
    n, length(times), nrow(observed_long),
    mean(!is.na(long[["sbp"]][long[["time"]] == 12])),
    lme_beta["programIntervention:time"],
    sqrt(diag(lme_V))["programIntervention:time"], rho_hat,
    mixed_change[["Estimate"]]
  )
)
knitr::kable(case_summary, digits = 3,
  caption = "Selected design and model results from the simulated longitudinal study")

Selected design and model results from the simulated longitudinal study

Result	Value
Randomized participants	360.000
Scheduled visits	5.000
Observed outcome records	1642.000
Month-12 retention	0.833
Mixed-model program-by-time estimate	-0.465
Mixed-model program-by-time SE	0.064
AR(1) residual correlation	0.529
Difference in baseline-to-month-12 changes	-5.582

17.1.1 Auditable reporting template

We analyzed [number] participants scheduled at [times and units], with [number/proportion] outcomes observed at each visit. The primary estimand was [population, outcome, treatment strategy, contrast, summary measure]. Time was represented as [categorical/continuous/spline] and baseline was handled by [strategy]. We fitted [marginal/mixed model] with fixed effects [list], participant dependence represented by [covariance/random effects], and estimation by [ML/REML/GEE]. We report [standardized visit means/change contrast/effect scale] with 95% confidence intervals. Likelihood inference assumed missing at random conditional on [variables/history]; sensitivity analyses used [methods and ranges]. We assessed [mean form, residuals, variance, covariance, influence, convergence]. Limitations include [missingness, timing, measurement, model dependence, causal scope, generalizability].

17.2 Analysis checklist

1. Define the population, outcome, time origin, treatment/exposure strategy, and target contrast.
2. Draw the visit schedule and distinguish scheduled, actual, and analysis time.
3. Verify IDs, duplicate person-times, units, ordering, impossible values, and protocol windows.
4. Tabulate observed records, retention, reasons, and visit patterns by important groups.
5. Plot individual trajectories, visit distributions, group means, and missing tails.
6. Prespecify baseline handling and the fixed time shape.
7. Choose marginal or subject-specific inference to match the question.
8. Specify random effects and/or residual covariance from design and scientific knowledge.
9. Fit using the correct ML/REML comparison strategy and record convergence.
10. Translate coefficients into standardized visit means and prespecified contrasts.
11. Diagnose mean shape, variance, tails, covariance, and participant influence.
12. State missingness assumptions and conduct clinically interpretable sensitivity analyses.
13. Distinguish association, randomized total effect, mediation, and prediction.
14. Report estimates, intervals, analysis population, software, and limitations.

18 Common errors at a glance

Error	Why it is wrong	Better approach
Treat repeated rows as independent	Standard OLS uncertainty assumes zero within-person covariance	Use an appropriate covariance, GEE, or mixed model
Call a visit number elapsed time	Visits may be delayed or irregular	Store scheduled and actual time explicitly
Use only endpoint completers	Selection and information loss can change the target population	Use all observed outcomes under stated assumptions; assess sensitivity
Select time shape after inspecting significance	Data-dependent selection distorts inference	Prespecify a plausible shape and test sensitivity
Interpret the program main effect as month 12	With baseline coded zero, it is the month-0 difference	Estimate an explicit month-12 contrast
Interpret a linear term despite quadratic time	The slope changes with time	Report derivatives or visit contrasts
Choose unstructured covariance automatically	It can be unstable and data-hungry	Match complexity to visits, overlap, and sample size
Add random slope and AR(1) mechanically	They can compete to explain the same dependence	Justify each layer and inspect estimates/sensitivity
Treat BLUPs as observed personal effects	They are shrunken estimates with uncertainty	Use them cautiously for prediction and show uncertainty
Say MMRM “handles all missing data”	Its likelihood still relies on MAR and model correctness	Include predictors, document assumptions, assess MNAR sensitivity
Adjust for a postrandomization covariate in the primary model	It can change the estimand or introduce bias	Define total, direct, or mediation targets before adjustment
Report observed power	It duplicates information in the result	Plan prospective power under plausible scenarios

19 Exercises and answers

1. Why is a repeated cross-section unable to estimate individual change?
2. What does programIntervention mean when time 0 is baseline and an interaction is included?
3. What does programIntervention:time mean in the quadratic model used here?
4. Why can an AR(1) covariance be inappropriate for irregular visit times?
5. Contrast the targets of a marginal model and a mixed model.
6. Why is an unstructured covariance not “assumption-free”?
7. What assumptions allow likelihood-based mixed models to use incomplete outcome vectors?
8. Why should activity be split into within- and between-person components?
9. Which uncertainty components are absent from a fixed-part confidence band for a future individual?
10. Why is complete-pair change analysis not automatically protected by randomization?

Show exercise answers

1. Different people are observed at different waves, so a change in the sample mean cannot be decomposed into within-person change versus changing composition.
2. It is the adjusted intervention-minus-usual-care difference at month 0, the reference time—not the overall or month-12 difference.
3. It is the intervention-minus-usual-care difference in the linear monthly component. Both groups share the quadratic curvature in this specification, so their slope difference is constant.
4. Discrete AR(1) indexes visit lags. Two adjacent visits are assigned the same correlation even if their elapsed gaps differ; continuous-time correlation may better match irregular timing.
5. A marginal model targets the population-average response. A mixed model conditions on latent participant effects and also supports subject-specific prediction. Under nonlinear links their coefficients differ numerically and conceptually.
6. It assumes one common positive-definite covariance matrix and estimates many parameters. Sparse visit overlap or too few participants can make it unstable.
7. Correctly specified mean and covariance models plus MAR conditional on variables represented in the likelihood, along with regularity and distinct-parameter assumptions.
8. The raw coefficient mixes differences between typically active people with visit-to-visit deviations from a person’s own activity. The decomposition makes these associations explicit.
9. Future random-effect heterogeneity, future residual variation, and sometimes covariance-parameter uncertainty. A confidence band describes uncertainty in the mean, not individual prediction.
10. Randomization protects assigned-group comparisons at baseline, but conditioning on being observed at month 12 is postrandomization selection. Dropout can depend on prognosis or treatment.

20 Quick reference

20.1 Core models and R entry points

quick_reference <- data.frame(
  Goal = c(
    "Independence-working mean", "Participant-cluster teaching sandwich",
    "Marginal CS/AR(1)", "Categorical-visit MMRM", "Random intercept and slope",
    "Fixed-part prediction", "Observed-person BLUP prediction",
    "Continuous-time residual correlation", "Prospective two-group approximation"
  ),
  R_entry = c(
    "lm(y ~ group * time + covariates)",
    "Sum participant score outer products",
    "gls(..., correlation = corCompSymm()/corAR1())",
    "gls(y ~ group * visit + baseline, corSymm(), weights = varIdent())",
    "lme(..., random = ~ time | id, correlation = corAR1())",
    "predict(fit, level = 0)", "predict(fit, level = 1)",
    "corCAR1(form = ~ actual_time | id)",
    "power.t.test(...) or design-aligned simulation"
  )
)
knitr::kable(quick_reference, caption = "Longitudinal analysis quick reference")

Longitudinal analysis quick reference

Goal	R_entry
Independence-working mean	lm(y ~ group * time + covariates)
Participant-cluster teaching sandwich	Sum participant score outer products
Marginal CS/AR(1)	gls(…, correlation = corCompSymm()/corAR1())
Categorical-visit MMRM	gls(y ~ group * visit + baseline, corSymm(), weights = varIdent())
Random intercept and slope	lme(…, random = ~ time | id, correlation = corAR1())
Fixed-part prediction	predict(fit, level = 0)
Observed-person BLUP prediction	predict(fit, level = 1)
Continuous-time residual correlation	corCAR1(form = ~ actual_time | id)
Prospective two-group approximation	power.t.test(…) or design-aligned simulation

20.2 Final interpretation checklist

• Is the effect population-averaged or conditional on participant random effects?
• What time point does the intercept and group main effect reference?
• Does the reported slope account for interactions and nonlinear terms?
• Are visit means standardized to an explicit covariate distribution?
• Do confidence intervals use the full covariance of coefficient contrasts?
• Does prediction concern the population mean, an observed person, or a new person?
• Are random effects distinguished from residual serial correlation?
• Is baseline handled once in a coherent estimand strategy?
• Are dropout reasons, MAR assumptions, and MNAR sensitivity ranges transparent?
• Are time-varying covariates treated as possible postexposure variables?
• Does prospective power match the scheduled visits, covariance, attrition, and primary model?
• Are causal claims supported by design and assumptions beyond repeated measurement?

Next directions for study

Next topics include generalized estimating equations, small-sample robust inference, nonlinear mixed models, generalized linear mixed models, intensive longitudinal and state-space models, functional data, joint longitudinal–survival models, multiple imputation for multilevel data, pattern-mixture and tipping-point analyses, causal inference with time-varying confounding, estimands under intercurrent events, dynamic treatment regimes, and transportability.

sessionInfo()

## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.5.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
## 
## time zone: America/Edmonton
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.39   R6_2.6.1        fastmap_1.2.0   xfun_0.60       lattice_0.22-9 
##  [6] cachem_1.1.0    knitr_1.51      htmltools_0.5.9 rmarkdown_2.31  stats4_4.6.1   
## [11] lifecycle_1.0.5 cli_3.6.6       grid_4.6.1      sass_0.4.10     jquerylib_0.1.4
## [16] compiler_4.6.1  tools_4.6.1     nlme_3.1-169    evaluate_1.0.5  bslib_0.12.0   
## [21] yaml_2.3.12     rlang_1.3.0     jsonlite_2.0.0

For teaching only: This module introduces longitudinal data structure, marginal covariance models, mixed-effects trajectories, model-matrix contrasts, missingness, diagnostics, and reporting. It does not replace a protocol-specific statistical analysis plan, clinical judgment, data-monitoring procedures, or consultation with a biostatistician. All data and numerical results are simulated.