Mixed Models in R: lme4 and nlme

Introduction

Clinical trials frequently collect repeated measurements from the same patient. Examples include blood pressure measured at multiple visits, tumor burden assessed repeatedly, laboratory biomarkers collected over time, pulmonary function measurements, pain scores, and patient-reported outcomes.

These observations create an important statistical problem: measurements from the same patient are generally correlated.

A conventional linear regression model assumes independent observations. Applying that assumption to repeated measurements can produce incorrect standard errors and therefore misleading statistical inference.

Mixed-effects models provide one of the most flexible frameworks for addressing this problem.

Why Mixed Models Matter in Clinical Trials

Consider a randomized clinical trial comparing an experimental treatment with placebo. Suppose the primary continuous endpoint is measured at Weeks 4, 8, 12, 24, and 36.

The data might look conceptually like this:

The three measurements from Patient 001 are not three independent people. They belong to the same individual.

A mixed model explicitly accounts for this structure.

The Basic Linear Mixed Model

A useful starting point is:

where:

• \(Y_{ij}\) is the outcome for subject \(i\) at measurement \(j\).
• \(\boldsymbol{\beta}\) contains the fixed-effect parameters.
• \(\mathbf{b}_i\) contains subject-specific random effects.
• \(\mathbf{x}_{ij}\) describes the fixed-effect design.
• \(\mathbf{z}_{ij}\) describes the random-effect design.
• \(\epsilon_{ij}\) is the residual error.

The random effects are commonly assumed to follow:

and the residual errors may follow:

The distinction between \(\mathbf{D}\) and \(\mathbf{R}_i\) becomes extremely important when deciding whether to use lme4 or nlme.

Fixed Effects Versus Random Effects

A Clinical Example

Suppose a trial measures a continuous biomarker at baseline and Weeks 4, 8, 12, and 24.

The scientific question is:

This naturally suggests a model containing:

• Treatment
• Visit
• Treatment × Visit
• A subject-level random effect

A typical model might be written:

Here, \(\beta_3\) represents the treatment-by-time interaction.

Why the Treatment-by-Time Interaction Matters

Suppose the treatment effect is zero at baseline but becomes increasingly negative after treatment begins.

A model containing only treatment and time would force the treatment effect to be constant unless additional structure were introduced.

The interaction allows treatment effects to vary by visit.

For a categorical visit variable, the model can estimate a treatment contrast at each post-baseline visit.

Visualizing Patient-Level Variation

The individual lines represent subjects. The gold line represents a population-level mean trajectory.

A mixed model allows both levels of information to coexist mathematically: population-level effects are represented by fixed effects, while subject-specific deviations are represented by random effects.

Random Intercept Models

The simplest useful longitudinal mixed model gives each subject their own intercept.

The subject-specific intercept is:

Therefore, patients can start at different levels while sharing the same average slope.

Clinical Interpretation

Suppose the endpoint is systolic blood pressure.

One patient may have a baseline level of 150 mmHg and another 125 mmHg. A random-intercept model acknowledges that these patients begin at different levels.

The model does not require every patient to have the same baseline outcome.

Random Intercept in lme4

library(lme4)

fit_ri <- lmer(
  outcome ~ treatment + week + (1 | USUBJID),
  data = dat,
  REML = TRUE
)

summary(fit_ri)

The term:

(1 | USUBJID)

means that each subject receives a subject-specific random intercept.

Random Slopes

Patients may differ not only in their baseline levels but also in their rates of change.

A random-slope model can be written:

The individual slope becomes:

Thus, different patients can follow different trajectories.

Random Intercept and Random Slope in lme4

fit_rs <- lmer(
  outcome ~ treatment * week +
    (1 + week | USUBJID),
  data = dat,
  REML = TRUE
)

The expression:

(1 + week | USUBJID)

allows each subject to have:

• A random intercept
• A random slope for week
• A covariance between the random intercept and slope

What Does the Random-Effect Correlation Mean?

Suppose the model estimates a negative correlation between random intercepts and random slopes.

This could indicate that patients with higher baseline outcomes tend to have more negative longitudinal slopes.

The correlation is not automatically a treatment effect.

It describes the relationship between subject-specific deviations.

Removing the Random-Effect Correlation

In lme4, the notation:

(1 + week | USUBJID)

estimates a correlated random intercept and slope.

The notation:

(1 + week || USUBJID)

removes the intercept-slope covariance.

This can sometimes improve model stability.

Centering Time

Time is often easier to interpret when centered.

For example:

dat$week_c <- dat$week - 12

Now the model intercept corresponds to the expected outcome at Week 12 rather than at Week 0.

Centering can be particularly helpful when interpreting random intercepts.

Why Time Coding Matters

Consider:

outcome ~ treatment * week

versus:

outcome ~ treatment * factor(week)

These are different models.

Categorical Visit Models

A common longitudinal clinical-trial model uses visit as a factor:

dat$AVISIT <- factor(
  dat$AVISIT,
  levels = c(
    "Baseline",
    "Week 4",
    "Week 8",
    "Week 12",
    "Week 24"
  )
)

fit_visit <- lmer(
  outcome ~ treatment * AVISIT +
    (1 | USUBJID),
  data = dat,
  REML = TRUE
)

This estimates treatment differences at the individual visits rather than assuming that the treatment effect changes linearly over time.

The Baseline Visit Requires Care

Clinical-trial longitudinal models frequently include a baseline observation, but baseline handling requires careful thought.

For example, an analysis might model:

• Change from baseline as the dependent variable
• Observed post-baseline outcome with baseline adjustment
• Repeated post-baseline outcomes with baseline as part of the covariance structure

These are not interchangeable formulations.

Change From Baseline

Suppose:

A mixed model can then use change from baseline as the outcome.

fit_chg <- lmer(
  CHG ~ treatment * AVISIT +
    (1 | USUBJID),
  data = postbaseline,
  REML = TRUE
)

This can be appropriate when the estimand and analysis plan are defined in terms of change from baseline.

Baseline Adjustment

Another approach is to model the post-baseline outcome while adjusting for baseline:

fit_adj <- lmer(
  outcome ~ baseline +
    treatment * AVISIT +
    (1 | USUBJID),
  data = postbaseline,
  REML = TRUE
)

This is conceptually closer to an ANCOVA-style analysis extended across multiple post-baseline visits.

Entering the World of nlme

The nlme package is especially important because it provides explicit facilities for modeling the residual covariance structure.

This is one of the most important practical differences between lme4 and nlme.

Basic nlme Syntax

library(nlme)

fit_nlme <- lme(
  fixed = outcome ~ treatment * AVISIT,
  random = ~ 1 | USUBJID,
  data = dat,
  method = "REML"
)

summary(fit_nlme)

The syntax is somewhat different from lme4, but the underlying mixed-model concepts are the same.

Random Slopes in nlme

fit_nlme_rs <- lme(
  fixed = outcome ~ treatment * week,
  random = ~ week | USUBJID,
  data = dat,
  method = "REML"
)

The expression:

random = ~ week | USUBJID

corresponds to a subject-specific random intercept and random slope.

Random Effects Versus Residual Correlation

This distinction is fundamental.

There are two broad ways repeated observations can become correlated:

1. Through shared random effects.
2. Through a structured residual covariance matrix.

Mixed models can use either or both mechanisms.

Compound Symmetry

Under compound symmetry, observations from the same patient have equal pairwise correlation regardless of how far apart they are in time.

Conceptually:

This may be reasonable when correlation is approximately constant across visits.

AR(1) Correlation

An autoregressive AR(1) structure assumes that observations closer together in time are more strongly correlated.

The correlation between two observations can be represented as:

for equally spaced time points under the simplest formulation.

Thus:

• Adjacent visits have correlation \(\rho\).
• Visits two time units apart have correlation \(\rho^2\).
• Visits three time units apart have correlation \(\rho^3\).

AR(1) in nlme

fit_ar1 <- lme(
  fixed = outcome ~ treatment * AVISIT,
  random = ~ 1 | USUBJID,
  correlation = corAR1(
    form = ~ week | USUBJID
  ),
  data = dat,
  method = "REML"
)

summary(fit_ar1)

This is one of the major reasons nlme remains particularly useful for clinical longitudinal analyses.

Continuous-Time Correlation

When assessments occur at irregular times, the assumption that visits are equally spaced may be inappropriate.

In such settings, a continuous-time correlation structure may be preferable.

corCAR1(
  form = ~ day | USUBJID
)

The key idea is that correlation is related to elapsed time rather than merely the visit number.

Heterogeneous Variances

Clinical-trial outcomes may have different residual variances at different visits.

For example, variability may increase as patients are followed longer.

The nlme package provides variance functions for this situation.

fit_het <- lme(
  fixed = outcome ~ treatment * AVISIT,
  random = ~ 1 | USUBJID,
  weights = varIdent(
    form = ~ 1 | AVISIT
  ),
  data = dat,
  method = "REML"
)

This allows residual variance to differ by visit.

Combining Heterogeneous Variance and AR(1)

A particularly useful model for longitudinal clinical data can include both:

• A correlation structure
• Visit-specific residual variances

fit_ar1_het <- lme(
  fixed = outcome ~ treatment * AVISIT,
  random = ~ 1 | USUBJID,
  correlation = corAR1(
    form = ~ week | USUBJID
  ),
  weights = varIdent(
    form = ~ 1 | AVISIT
  ),
  data = dat,
  method = "REML"
)

This is much closer to the type of covariance modeling often encountered in serious longitudinal clinical-trial work than a simple random-intercept model.

Random Effects and Residual Covariance Are Not the Same Thing

Consider a model with a random intercept:

Because all observations from a patient share \(b_i\), they become correlated.

However, this induced correlation has a particular structure.

A random-intercept model therefore should not automatically be interpreted as an arbitrary repeated-measures covariance model.

lme4 Versus nlme: A Practical Decision

A Complete Simulated Clinical Trial

We will now build a small simulated longitudinal trial from the ground up.

Suppose:

• 200 patients are randomized.
• 100 receive Active treatment.
• 100 receive Placebo.
• The endpoint is a continuous biomarker.
• Measurements occur at Weeks 0, 4, 8, 12, and 24.

Generate the Data

set.seed(2026)

n_subjects <- 200

subjects <- data.frame(
  USUBJID = sprintf("SUBJ-%03d", 1:n_subjects),
  treatment = rep(
    c("Placebo", "Active"),
    each = n_subjects / 2
  )
)

visits <- c(0, 4, 8, 12, 24)

dat <- merge(
  subjects,
  data.frame(week = visits)
)

dat <- dat[
  order(dat$USUBJID, dat$week),
]

dat$USUBJID <- factor(dat$USUBJID)

dat$treatment <- factor(
  dat$treatment,
  levels = c("Placebo", "Active")
)

Simulate Patient-Level Random Effects

subject_re <- data.frame(
  USUBJID = levels(dat$USUBJID),
  random_intercept = rnorm(
    n_subjects,
    mean = 0,
    sd = 8
  )
)

dat <- merge(
  dat,
  subject_re,
  by = "USUBJID"
)

dat <- dat[
  order(dat$USUBJID, dat$week),
]

Generate the Outcome

dat$outcome <-
  100 +
  0.05 * dat$week +
  ifelse(
    dat$treatment == "Active",
    -0.30 * dat$week,
    0
  ) +
  dat$random_intercept +
  rnorm(nrow(dat), 0, 5)

The active treatment therefore has a more favorable longitudinal trajectory.

Fit the Basic lme4 Model

library(lme4)

fit_lme4 <- lmer(
  outcome ~ treatment * week +
    (1 | USUBJID),
  data = dat,
  REML = TRUE
)

summary(fit_lme4)

Interpreting the Fixed Effects

Because treatment uses Placebo as the reference level:

The treatment-by-week interaction is therefore the key parameter if the scientific question concerns differential longitudinal change.

Extracting Fixed Effects

fixef(fit_lme4)

To obtain confidence intervals:

confint(
  fit_lme4,
  parm = "beta_",
  method = "Wald"
)

Different confidence-interval methods can have different computational and small-sample properties. The method should be chosen deliberately rather than automatically.

Extracting Random Effects

ranef(fit_lme4)

The result contains estimated subject-specific deviations from the population fixed effects.

These are often called empirical Bayes or conditional-mode estimates, depending on the modeling framework and terminology.

Random Effects Are Not Ordinary Fixed-Effect Estimates

A common conceptual error is to interpret every estimated random effect as though it were a precisely estimated individual parameter.

Random effects are shrunk toward the population mean.

Patients with limited information tend to have stronger shrinkage toward zero than patients with extensive information.

Adding a Random Slope

fit_lme4_rs <- lmer(
  outcome ~ treatment * week +
    (1 + week | USUBJID),
  data = dat,
  REML = TRUE
)

summary(fit_lme4_rs)

This model allows patients to differ in both baseline level and rate of change.

Comparing Random-Effects Structures

When comparing models that differ only in random-effects structure, likelihood comparisons can be informative.

anova(
  fit_lme4,
  fit_lme4_rs
)

However, mixed-model comparison requires careful attention to estimation method, nesting, boundary parameters, and the scientific purpose of the comparison.

REML Versus Maximum Likelihood

Mixed models are commonly estimated using either:

• Restricted maximum likelihood (REML)
• Maximum likelihood (ML)

REML estimates variance components while accounting for the loss of degrees of freedom associated with estimating fixed effects.

Maximum likelihood estimates all parameters under the likelihood formulation.

When Should REML Be Used?

REML is often preferred for final estimation of variance components when the fixed-effects structure is already established.

lmer(
  outcome ~ treatment * week +
    (1 + week | USUBJID),
  data = dat,
  REML = TRUE
)

When Should ML Be Used?

When comparing models with different fixed-effects structures using likelihood criteria, ML is generally used rather than REML.

fit_ml_1 <- lmer(
  outcome ~ treatment + week +
    (1 | USUBJID),
  data = dat,
  REML = FALSE
)

fit_ml_2 <- lmer(
  outcome ~ treatment * week +
    (1 | USUBJID),
  data = dat,
  REML = FALSE
)

anova(
  fit_ml_1,
  fit_ml_2
)

AIC and BIC

Information criteria can help compare competing models:

AIC(
  fit_ml_1,
  fit_ml_2
)

BIC(
  fit_ml_1,
  fit_ml_2
)

Lower AIC or BIC generally indicates a preferred balance between fit and complexity under the corresponding criterion.

These criteria should not replace scientific judgment.

Modeling Time as a Factor

For many confirmatory longitudinal analyses, visit is categorical.

dat$visit <- factor(
  dat$week,
  levels = c(0, 4, 8, 12, 24)
)

fit_factor <- lmer(
  outcome ~ treatment * visit +
    (1 | USUBJID),
  data = dat,
  REML = TRUE
)

This allows the mean trajectory to take a different shape at each scheduled visit.

Why Factor Visit Is Often Attractive in Clinical Trials

Suppose the observed mean trajectory is:

A linear time effect would impose a straight-line relationship.

A factor-time model does not.

It can estimate the Week 4, Week 8, Week 12, and Week 24 effects independently, subject to the model structure.

Estimated Marginal Means

In clinical reporting, the raw model coefficients are often not the most useful quantities to present.

Estimated marginal means are frequently more interpretable.

library(emmeans)

emm <- emmeans(
  fit_factor,
  ~ treatment | visit
)

emm

This provides estimated treatment means at each visit.

Treatment Comparisons at Each Visit

contrast(
  emm,
  method = "revpairwise"
)

A common clinical question is:

The emmeans framework makes such contrasts much easier to obtain than manually reconstructing them from coefficient tables.

Confidence Intervals for Treatment Differences

pairs(
  emm,
  adjust = "none"
)

Multiplicity adjustment should follow the statistical analysis plan rather than being selected solely because it produces a preferred presentation.

Example of a Clinical Summary Table

The numbers above are illustrative.

Why This Is More Useful Than the Raw Coefficient Table

A coefficient such as:

treatmentActive:visit24

depends on the model parameterization and reference categories.

The estimated marginal mean contrast directly answers the clinical question:

Covariance Structures in Clinical Longitudinal Data

One of the most important advanced decisions is how to represent within-patient correlation.

Common choices include:

Unstructured Covariance

With four visits, an unstructured covariance matrix can be represented as:

Every variance and covariance is estimated.

This is flexible but can require many parameters.

Why Unstructured Covariance Can Become Difficult

With \(m\) repeated measurements, an unstructured covariance matrix requires:

variance and covariance parameters.

For 5 visits:

parameters are required.

For 10 visits:

parameters are required.

This illustrates why covariance selection becomes increasingly important as the number of repeated assessments grows.

Using nlme for More Explicit Covariance Modeling

fit_cs <- lme(
  fixed = outcome ~ treatment * AVISIT,
  random = ~ 1 | USUBJID,
  correlation = corCompSymm(
    form = ~ week | USUBJID
  ),
  data = dat,
  method = "REML"
)

fit_ar1 <- lme(
  fixed = outcome ~ treatment * AVISIT,
  random = ~ 1 | USUBJID,
  correlation = corAR1(
    form = ~ week | USUBJID
  ),
  data = dat,
  method = "REML"
)

Comparing Covariance Structures

A common development workflow is:

Missing Data

Longitudinal clinical-trial datasets commonly contain missing post-baseline measurements.

Mixed models are attractive partly because they can use available observations without requiring every patient to have complete follow-up.

Under standard likelihood-based estimation, this can provide valid inference under a missing-at-random framework when the model is correctly specified and the assumptions underlying the missingness mechanism are appropriate.

MAR Does Not Mean "Missing Data Are Ignored"

Suppose patients with worsening disease are more likely to discontinue.

The probability of missingness may depend on observed prior outcomes.

A likelihood-based mixed model can potentially accommodate such a mechanism under MAR.

But the model does not magically make missingness irrelevant.

Missing Not at Random

If the probability of missingness depends on an unobserved outcome after conditioning on observed information, the missingness mechanism may be missing-not-at-random (MNAR).

Sensitivity analyses may then be needed.

Potential approaches include:

• Pattern-mixture models
• Selection models
• Reference-based imputation
• Tipping-point analyses

The appropriate strategy depends on the estimand and trial context.

Dropout Is Not the Same as a Missing Value

In clinical-trial reporting, it is useful to distinguish:

• Intermittent missing assessments
• Permanent discontinuation
• Death
• Administrative censoring
• Missing data caused by treatment failure

These may have different implications for the estimand and sensitivity analysis.

Diagnostics

A mixed model should not be accepted simply because the software returned an estimate.

Diagnostics should be part of the analysis workflow.

Residuals Versus Fitted Values

plot(
  fitted(fit_lme4),
  resid(fit_lme4)
)

abline(
  h = 0,
  lty = 2
)

This can reveal:

• Nonlinearity
• Heteroscedasticity
• Outliers
• Model misspecification

Normal Q-Q Plot

qqnorm(
  resid(fit_lme4)
)

qqline(
  resid(fit_lme4)
)

The normality assumption is primarily relevant to the inferential framework and residual distribution; modest deviations do not automatically invalidate a mixed model.

Checking Random Effects

qqnorm(
  ranef(fit_lme4)$USUBJID[[1]]
)

qqline(
  ranef(fit_lme4)$USUBJID[[1]]
)

Strong departures from the assumed random-effects distribution can indicate model inadequacy, although diagnostics should be interpreted in context.

Singular Fits in lme4

One common warning is:

boundary (singular) fit

This can occur when the estimated random-effects covariance matrix is on the boundary of the parameter space.

For example, a random-slope variance may be estimated essentially as zero.

isSingular(
  fit_lme4_rs,
  tol = 1e-4
)

A singular fit is not merely an annoying warning to suppress.

Convergence Problems

Mixed models are nonlinear optimization problems and can sometimes fail to converge cleanly.

Potential causes include:

• Overly complex random-effects structures
• Highly correlated predictors
• Poor scaling of continuous variables
• Too few observations per subject
• Too little information to estimate covariance parameters
• Boundary solutions

Scaling Time

If time is measured in days:

outcome ~ treatment * day +
  (1 + day | USUBJID)

the random slope can become numerically difficult if the time scale is large.

Using weeks:

dat$week <- dat$day / 7

may make the parameters easier to interpret and sometimes improve numerical behavior.

Alternative Optimizers

For difficult lme4 models, alternative optimizers can sometimes help diagnose optimization issues.

fit_bobyqa <- lmer(
  outcome ~ treatment * week +
    (1 + week | USUBJID),
  data = dat,
  REML = TRUE,
  control = lmerControl(
    optimizer = "bobyqa"
  )
)

Optimizer changes should not be used to conceal a fundamentally unidentified model.

Model Validation Workflow

Nested Random Effects

Clinical data sometimes contain hierarchical structures such as:

• Patients within sites
• Patients within investigators
• Measurements within patients
• Patients within treatment centers

For example:

outcome ~ treatment * visit +
  (1 | site/USUBJID)

The notation represents nested random effects.

This can be expanded conceptually as:

(1 | site) +
(1 | site:USUBJID)

Crossed Random Effects

Not all hierarchical structures are nested.

If observations are associated with multiple non-nested grouping factors, crossed random effects may be appropriate.

outcome ~ treatment +
  (1 | patient) +
  (1 | assessor)

This could be relevant in some measurement studies where multiple assessors evaluate multiple patients.

Multiple Random Effects

A model can contain multiple sources of random variation:

lmer(
  outcome ~ treatment * visit +
    (1 | site) +
    (1 | USUBJID),
  data = dat
)

Whether such a structure is justified depends on the design and amount of information available for each variance component.

Random Slopes by Treatment?

A common question is whether treatment should appear in the random-effects structure.

Usually, treatment is a fixed effect when the scientific objective is to estimate the treatment effect.

The random slope typically represents subject-specific variation in a continuous covariate such as time.

(1 + week | USUBJID)

is therefore fundamentally different from:

(1 + treatment | USUBJID)

The latter requires a meaningful subject-level replication structure for the random treatment coefficient and is generally not the standard representation for a randomized treatment effect.

Repeated Measures Versus Longitudinal Random Effects

These concepts overlap but are not identical.

A repeated-measures model can represent within-subject covariance directly. A random-effects model represents subject-specific deviations that induce correlation.

In practice, the two approaches can produce similar covariance patterns in some cases.

But the implied covariance matrices are not generally identical.

Clinical-Trial Example: Blood Pressure

Suppose a trial measures systolic blood pressure at:

• Baseline
• Week 2
• Week 4
• Week 8
• Week 12

A plausible model might be:

fit_bp <- lmer(
  SBP_change ~ treatment * visit +
    baseline_SBP +
    (1 | USUBJID),
  data = bp_dat,
  REML = TRUE
)

The treatment difference at Week 12 can then be obtained with:

emmeans(
  fit_bp,
  pairwise ~ treatment | visit
)

Clinical-Trial Example: Oncology Biomarker

Suppose a biomarker is measured repeatedly and the scientific objective is to compare longitudinal treatment trajectories.

fit_bio <- lmer(
  log_biomarker ~ treatment * visit +
    baseline_log_biomarker +
    (1 + week | USUBJID),
  data = biomarker_dat,
  REML = TRUE
)

The log transformation may be appropriate when the biomarker distribution is strongly right-skewed or multiplicative effects are scientifically meaningful.

Clinical-Trial Example: Patient-Reported Outcome

Patient-reported outcomes often have repeated assessments and substantial between-patient heterogeneity.

fit_pro <- lmer(
  PRO_change ~ treatment * visit +
    baseline_PRO +
    (1 | USUBJID),
  data = pro_dat,
  REML = TRUE
)

If the instrument has unusual distributional properties, the appropriateness of a Gaussian mixed model should be evaluated rather than assumed.

Generalized Mixed Models

Although this tutorial focuses on linear mixed models, the same conceptual framework extends to non-Gaussian outcomes.

For binary outcomes, lme4 provides glmer().

fit_binary <- glmer(
  response ~ treatment * visit +
    (1 | USUBJID),
  data = dat_binary,
  family = binomial
)

This is a generalized linear mixed model rather than a linear mixed model.

What lme4 Does Especially Well

• Clear random-effects formula syntax
• Efficient estimation for many mixed-model applications
• Random intercepts and slopes
• Nested and crossed random effects
• Generalized mixed models through glmer()
• Useful extraction methods for fitted models

What nlme Does Especially Well

• Explicit residual correlation structures
• AR(1) modeling
• Continuous-time correlation structures
• Heterogeneous residual variances
• Flexible variance functions
• Nonlinear mixed-effects models

A More Advanced nlme Model

Consider:

fit_advanced <- lme(
  fixed =
    outcome ~ treatment * AVISIT + baseline,
  random =
    ~ 1 | USUBJID,
  correlation =
    corAR1(
      form = ~ week | USUBJID
    ),
  weights =
    varIdent(
      form = ~ 1 | AVISIT
    ),
  data = dat,
  method = "REML"
)

This model includes:

• Fixed treatment effects
• Fixed visit effects
• Treatment-by-visit interaction
• Baseline adjustment
• Subject-specific random intercept
• AR(1) residual correlation
• Visit-specific residual variance

That is a substantially richer representation of longitudinal dependence than a basic random-intercept model.

Model Complexity Must Be Earned

A more complicated covariance model is not automatically a better model.

Every additional variance or covariance parameter requires information.

A useful principle is:

Common Mistake: Random Intercept by Default

A very common workflow is:

lmer(
  outcome ~ treatment * visit +
    (1 | USUBJID),
  data = dat
)

followed by assuming that the repeated-measures problem has been completely resolved.

That may be adequate.

But it should be evaluated rather than assumed.

The subject-level random intercept imposes a particular covariance structure.

Common Mistake: Random Slope by Default

The opposite mistake is fitting:

(1 + week | USUBJID)

to every dataset simply because longitudinal observations exist.

A random slope requires enough within-subject information to estimate subject-specific changes.

With only two observations per patient, a complex random-effects structure may be poorly identified.

Common Mistake: Treating Every Visit as Independent

A model such as:

lm(
  outcome ~ treatment * visit,
  data = dat
)

does not account for repeated observations from the same patient.

The standard errors can therefore be inappropriate.

Common Mistake: Ignoring Treatment Assignment in the Longitudinal Structure

A model such as:

outcome ~ treatment + visit

assumes the treatment effect is constant across visits.

If the scientific question concerns whether trajectories differ, the treatment-by-visit interaction should be considered.

Common Mistake: Automatically Using Numeric Visit

This:

treatment * week

is a linear-time model.

This:

treatment * factor(week)

is a categorical-visit model.

They answer different questions.

Common Mistake: Comparing REML Fits With Different Fixed Effects

Suppose:

fit1 <- lmer(
  outcome ~ treatment + visit +
    (1 | USUBJID),
  data = dat,
  REML = TRUE
)

fit2 <- lmer(
  outcome ~ treatment * visit +
    (1 | USUBJID),
  data = dat,
  REML = TRUE
)

Comparing their REML likelihoods for fixed-effect selection is generally not the appropriate approach.

Use ML for fixed-effect model comparison:

fit1_ml <- update(
  fit1,
  REML = FALSE
)

fit2_ml <- update(
  fit2,
  REML = FALSE
)

anova(
  fit1_ml,
  fit2_ml
)

Common Mistake: Confusing Random Effects With Repeated-Measures Covariance

A random intercept captures subject-specific heterogeneity.

An AR(1) residual structure describes how residual correlation changes with time separation.

These are related concepts but not interchangeable.

Common Mistake: Reporting Only the p-Value

A clinical-trial mixed-model result should generally include:

• Estimated treatment difference
• Confidence interval
• P-value where appropriate
• Analysis visit/time point
• Analysis population
• Model specification
• Covariance/random-effects structure

A p-value without the estimated treatment effect provides limited clinical information.

Reporting a Mixed-Model Result

A concise report might state:

The exact wording and estimand should follow the statistical analysis plan.

Clinical Study Report Considerations

Model Specification Example

A full statistical-programming specification might read conceptually:

Traceability in Clinical Programming

For regulated clinical-trial analyses, the model result should be traceable from the source analysis dataset through derived variables to the final reported table or figure.

Reproducibility

A clinical-trial mixed-model program should record:

• R version
• Package versions
• Model specification
• Analysis dataset version
• Derivation programs
• Random seed when simulation or stochastic procedures are involved
• Model convergence status
• Warnings
• Output-generation code

Package Versions

sessionInfo()

packageVersion("lme4")
packageVersion("nlme")
packageVersion("emmeans")

For regulated environments, package versions should be controlled and documented as part of the validated or appropriately qualified computational environment.

Extracting a Clean Results Dataset

For reporting, it is often useful to transform model-derived results into a structured dataset.

emm_results <- as.data.frame(
  emmeans(
    fit_factor,
    ~ treatment | visit
  )
)

contrast_results <- as.data.frame(
  pairs(
    emmeans(
      fit_factor,
      ~ treatment | visit
    )
  )
)

The resulting datasets can feed standardized table-generation programs.

Predicted Trajectories

Mixed models can also generate predicted trajectories.

newdat <- expand.grid(
  treatment = levels(dat$treatment),
  week = seq(0, 24, by = 1)
)

newdat$pred <- predict(
  fit_lme4,
  newdata = newdat,
  re.form = NA
)

The argument:

re.form = NA

requests predictions based on the population-level fixed effects rather than subject-specific random effects.

Population-Level Versus Subject-Level Prediction

The distinction matters greatly when creating clinical figures.

Visualization of Estimated Treatment Trajectories

library(ggplot2)

ggplot(
  newdat,
  aes(
    x = week,
    y = pred,
    color = treatment
  )
) +
  geom_line(
    linewidth = 1
  ) +
  labs(
    x = "Week",
    y = "Estimated Mean Outcome",
    color = "Treatment"
  ) +
  theme_classic()

For a publication or CSR, confidence intervals around the estimated marginal means should generally be included where they add useful inferential context.

Model-Based Versus Raw Means

A frequent reporting mistake is to call model-adjusted estimates "observed means."

They are not the same.

Unbalanced Data

One advantage of mixed models is their ability to work with unbalanced longitudinal data.

For example:

A mixed model can use the available observations under its modeling assumptions rather than automatically discarding every patient with an incomplete trajectory.

But Unbalanced Does Not Mean Unproblematic

Suppose patients who deteriorate rapidly are more likely to discontinue.

The observed data could then become increasingly enriched for patients doing well.

The mixed model can account for observed longitudinal information under an appropriate missingness assumption, but the analyst must still investigate the dropout process.

Sensitivity Analyses

For important confirmatory analyses, sensitivity analyses may examine whether conclusions change under alternative assumptions.

Potential sensitivity analyses include:

• Alternative covariance structures
• Alternative missing-data assumptions
• Alternative baseline definitions where scientifically justified
• Alternative time representations
• Alternative handling of post-treatment events
• Pattern-mixture or reference-based approaches

Post-Baseline Intercurrent Events

Clinical-trial mixed models must be considered in the context of the estimand framework.

Examples of intercurrent events include:

• Treatment discontinuation
• Rescue medication
• Switching treatment
• Death
• Protocol-defined treatment failure

The model should not be selected independently of the estimand.

A Practical End-to-End Analysis Script

library(lme4)
library(emmeans)
library(ggplot2)

#--------------------------------------------------
# 1. Prepare analysis variables
#--------------------------------------------------

dat$USUBJID <- factor(dat$USUBJID)

dat$treatment <- factor(
  dat$treatment,
  levels = c("Placebo", "Active")
)

dat$visit <- factor(
  dat$visit,
  levels = c(
    "Baseline",
    "Week 4",
    "Week 8",
    "Week 12",
    "Week 24"
  )
)

#--------------------------------------------------
# 2. Fit the longitudinal mixed model
#--------------------------------------------------

fit <- lmer(
  outcome ~ baseline +
    treatment * visit +
    (1 | USUBJID),
  data = dat,
  REML = TRUE
)

#--------------------------------------------------
# 3. Review model
#--------------------------------------------------

summary(fit)

#--------------------------------------------------
# 4. Check singularity
#--------------------------------------------------

isSingular(fit)

#--------------------------------------------------
# 5. Extract fixed effects
#--------------------------------------------------

fixef(fit)

#--------------------------------------------------
# 6. Estimated marginal means
#--------------------------------------------------

emm <- emmeans(
  fit,
  ~ treatment | visit
)

#--------------------------------------------------
# 7. Treatment comparisons
#--------------------------------------------------

contrasts <- pairs(
  emm,
  adjust = "none"
)

#--------------------------------------------------
# 8. Convert to reporting datasets
#--------------------------------------------------

emm_results <- as.data.frame(emm)

contrast_results <- as.data.frame(
  contrasts
)

#--------------------------------------------------
# 9. Diagnostics
#--------------------------------------------------

plot(
  fitted(fit),
  resid(fit)
)

qqnorm(
  resid(fit)
)

qqline(
  resid(fit)
)

#--------------------------------------------------
# 10. Record environment
#--------------------------------------------------

sessionInfo()

How to Choose Between lme4 and nlme

Important Clinical-Trial Distinction

There is no universal rule that says:

or:

They provide different modeling capabilities.

The correct package depends on the model required by the scientific question and analysis plan.

Checklist for Clinical-Trial Mixed Models

Quick Reference: lme4

library(lme4)

# Random intercept
lmer(
  y ~ x + (1 | subject),
  data = dat
)

# Random intercept + slope
lmer(
  y ~ x + (1 + x | subject),
  data = dat
)

# Uncorrelated intercept and slope
lmer(
  y ~ x + (1 + x || subject),
  data = dat
)

# Generalized mixed model
glmer(
  y ~ treatment + (1 | subject),
  data = dat,
  family = binomial
)

Quick Reference: nlme

library(nlme)

# Random intercept
lme(
  fixed = y ~ x,
  random = ~ 1 | subject,
  data = dat
)

# Random intercept + slope
lme(
  fixed = y ~ x,
  random = ~ x | subject,
  data = dat
)

# AR(1)
lme(
  fixed = y ~ x,
  random = ~ 1 | subject,
  correlation =
    corAR1(
      form = ~ time | subject
    ),
  data = dat
)

# Heterogeneous residual variance
lme(
  fixed = y ~ x,
  random = ~ 1 | subject,
  weights =
    varIdent(
      form = ~ 1 | visit
    ),
  data = dat
)

Final Conceptual Framework

The easiest way to remember mixed models is to separate three questions.

This distinction is the key to understanding why lme4 and nlme can produce different-looking model specifications even when they are addressing the same general longitudinal problem.

The Most Important Clinical-Trial Example

Suppose a randomized trial has a continuous endpoint measured repeatedly after baseline.

A strong starting conceptual model is:

The programming implementation then depends on how subject variation and residual dependence are intended to be represented.

With lme4, this might be:

lmer(
  outcome ~ baseline +
    treatment * visit +
    (1 | USUBJID),
  data = dat,
  REML = TRUE
)

With nlme, a richer covariance specification might be:

lme(
  fixed =
    outcome ~ baseline +
      treatment * visit,

  random =
    ~ 1 | USUBJID,

  correlation =
    corAR1(
      form = ~ week | USUBJID
    ),

  data = dat,

  method = "REML"
)

The second model explicitly describes residual correlation in addition to subject-level heterogeneity.

What to Remember

Summary

Mixed-effects models are among the most important statistical tools for analyzing longitudinal clinical-trial data.

They provide a principled way to combine population-level treatment effects with subject-level heterogeneity while accounting for the correlation created by repeated measurements.

In R, lme4 provides an intuitive and powerful framework for random-effects modeling, including random intercepts, random slopes, nested effects, crossed effects, and generalized mixed models.

nlme adds particularly important capabilities for explicit residual covariance and variance modeling, including AR(1), continuous-time correlation, compound symmetry, and heterogeneous residual variance.

For clinical-trial statisticians, the most important conceptual distinction is not simply "which package should I use?"

The more important questions are:

• What is the estimand?
• What are the fixed effects?
• How much subject-level heterogeneity should be modeled?
• What residual covariance structure is plausible?
• How should missing observations be handled?
• What treatment contrast is clinically relevant?
• How will the result be validated and reported?

Once those questions are answered, lme4 and nlme become implementation tools for expressing the statistical model in R.

References

Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1–48.

Pinheiro, J.C., & Bates, D.M. (2000). Mixed-Effects Models in S and S-PLUS. Springer.

Pinheiro, J., Bates, D., DebRoy, S., Sarkar, D., & R Core Team. nlme: Linear and Nonlinear Mixed Effects Models.

Verbeke, G., & Molenberghs, G. (2000). Linear Mixed Models for Longitudinal Data. Springer.

Fitzmaurice, G.M., Laird, N.M., & Ware, J.H. (2011). Applied Longitudinal Analysis. Wiley.

McCulloch, C.E., Searle, S.R., & Neuhaus, J.M. (2008). Generalized, Linear, and Mixed Models. Wiley.