From Linear Regression to Mixed-Effects Models: Understanding Fixed and Random Effects from First Principles

อ่านฉบับภาษาไทย (Thai version)

Abstract

Mixed-effects models become easier to understand when each component is introduced before the next. This article begins with ordinary linear regression and one blood-pressure example, then distinguishes shared intercepts and slopes from personal departures. It explains why repeated observations are correlated, how random-effect variation differs from residual variation, and why a treatment interaction is different from a random slope. Practical sections introduce covariance choices, estimation, diagnostics, and alternatives without treating any method as universally best. Anatomical sections show how categorical fixed contrasts can coexist with one numeric random position slope. A short WSS model specification separates the mean question from covariance choices. The closing chapter extends the same reasoning to wards, patients, anatomical segments and crossed raters. All numerical clinical examples are illustrative rather than empirical study findings. The aim is a self-contained path from the first regression equation to an informed interpretation of repeated measurements and clusters.

Jump to: 1 · 2 · 3 · 4 · 5 · 6 · 7 · 8 · 9 · 10 · 11 · 12 · 13 · 14 · 15 · 16

🔍 Click to enlarge Click anywhere to close

Part I. Understand Regression Before Mixed Models

Chapter 1. Linear Regression: What Are Intercept and Slope?

Start with systolic blood pressure, abbreviated BP, measured in millimetres of mercury (mmHg). We want to describe how average BP changes with months since baseline. A straight line is a useful first description, provided its limitations are checked.

Let \(Y\) denote an observed BP and \(X\) denote months since baseline. Let \(\beta_0\) be the mean BP at month zero, \(\beta_1\) the mean change per month, and \(\varepsilon\) the observation's departure from its mean. Then the regression equation is:

\[Y=\beta_0+\beta_1X+\varepsilon.\]

For illustration, set the intercept to 130 mmHg and the slope to -3.5 mmHg per month. The predicted mean is 130 at baseline, 126.5 after one month, and 119.5 after three months. Multiplying the slope by time converts a rate into an amount of change.

Suppose a month-one reading is 129 mmHg. Its residual relative to this illustrative line is 2.5 mmHg. A residual is an observed discrepancy, whereas the error term represents the unobserved departure in the statistical model.

What the line does and does not say

The intercept is meaningful only at the chosen zero of the predictor. If time zero represents baseline, it describes baseline mean BP. If zero lies far outside the observed range, the intercept may be mostly an extrapolation.

The slope describes an association in this model. A decreasing BP line alone does not establish that follow-up time causes the decrease. Treatment changes, selection, and other factors can also contribute.

A mean line is not a promise that every observed BP follows that line exactly. Measurement noise and unmeasured influences create departures. Before adding complexity, distinguish the location of the mean from the scatter around it.

The same arithmetic also prevents a common language error: a negative slope means BP falls as time increases, not that BP itself must become negative. We use the line within a clinically sensible observation window. We do not extend it indefinitely.

• Key idea: An intercept locates the mean line; a slope describes its change; a residual is the observed departure from that line.

Chapter 2. Fixed Intercept and Fixed Slope: One Regression Equation for Everyone

A fixed intercept and fixed slope provide shared coefficients for the population mean. In our simple example, everyone contributes to estimating the same baseline level and the same monthly change. Different measured values remain possible because observations depart from the mean.

The word fixed does not mean the coefficient is known, exact, or free of sampling uncertainty. A fitted slope has a standard error and a confidence interval. Those describe uncertainty in estimating the shared relationship from a finite sample.

Continuous and categorical predictors

Months are continuous in meaning even if visits occur at whole months. A one-unit change means one month, so the coefficient has units of mmHg per month. Rescaling time to years changes the numerical coefficient and its unit, without changing the underlying fitted line.

A categorical predictor needs a different interpretation. Let \(T\) equal zero for the reference BP treatment group and one for the comparison group, and let \(\alpha_0\) and \(\alpha_T\) denote its intercept and group coefficient. A model without time would be:

The notation \(E(Y\mid T)\) means the expected BP for the specified group.

\[E(Y\mid T)=\alpha_0+\alpha_TT.\]

The notation \(E(Y\mid T)\) means the expected BP for the specified group. Here \(\alpha_T\) is a difference between group means, not a rate of physical movement between two groups. Coding a category numerically does not make it a continuous quantity.

A terminology panel

In regression, fixed effects are shared mean coefficients. In panel-data discussions, a fixed-effects estimator often introduces person-specific intercepts as unrestricted nuisance parameters or removes them by within-person transformation. That approach can allow association between those intercepts and the predictors.

A conventional random-effects regression instead needs assumptions about how unobserved person effects relate to included predictors. If those assumptions are unsuitable, person-level means of changing predictors or a different analysis may be needed. The phrase fixed effects therefore requires its setting.

Meta-analysis uses fixed-effect or common-effect language for another question: whether study estimates target a common underlying effect. That is a different modelling context. None of these meanings says that a fixed coefficient must be identical to every individual's biological response.

• Key idea: Fixed coefficients describe shared population relationships and still have estimation uncertainty.

Chapter 3. Repeated Measurements: Why Four Observations Are Not Four Independent Patients

Now measure each person's BP at months zero, one, two, and three. The data contain four rows per person. Four rows describe four measurements, but they do not create four independent people.

The same person carries persistent influences across visits. Their physiology, baseline health, and measurement setting may make their BP values resemble one another. A high reading at one visit can therefore be informative about another reading from the same person.

Describe the rows before choosing a model

Use a person identifier, visit time, BP outcome, and any prespecified covariates. Each row should represent one measurement at a clearly defined time. Check whether duplicate person-time rows are intended replicates or recording errors.

Continuous describes the outcome's scale. Repeated describes how observations are organized. These are separate properties: repeated measurements can also record binary events, counts, or ordered categories.

A Gaussian model uses a normal distribution for errors remaining after the modelled mean and any personal departures are specified. It can suit continuous BP measurements when that approximation is reasonable. Other outcome types need other distributions and often other links between predictors and the mean. Repetition alone does not determine the outcome model.

Why independence matters

An analysis that treats every BP row as an unrelated person may misstate precision. The direction and size of the error depend on the design and the comparison being estimated. Such a comparison, for example a difference between treatment means, is called a contrast. There is no universal rule that every repeated-data standard error must increase by the same factor.

A within-person change uses the pairing between visits. A between-treatment comparison also depends on how many independent people were assigned to each treatment. The same data can therefore contain different amounts of information for different questions.

Keeping one average per person can sometimes answer a narrow average-BP question. It discards the shape of follow-up and may fail to answer a time-trend question. Our goal is to retain the repeated measurements while representing what they share.

• Key idea: Repeated observations share a person, so their information cannot generally be counted as independent patients.

Part II. Building a Mixed-Effects Model

Chapter 4. What Makes a Model Mixed? Fixed Effects + Random Effects

A mixed-effects model combines shared mean coefficients with person-specific random departures. Mixed refers to those two types of effect. It does not mean a mixture of continuous outcomes and repeated measurements.

Let \(i\) identify a person and \(j\) identify a measurement. Let \(Y_{ij}\) be BP and \(t_{ij}\) be months since baseline. Keep \(\beta_0\) and \(\beta_1\) for the shared intercept and time slope, introduce \(b_{0i}\) for person \(i\)'s intercept departure, and use \(\varepsilon_{ij}\) for the remaining error.

\[Y_{ij}=\beta_0+\beta_1t_{ij}+b_{0i}+\varepsilon_{ij}.\]

The fixed part describes the population mean line. The random intercept describes a person who stays above or below that mean line across visits. The remaining error describes departures not captured by either component.

What random means here

We model the person departures as draws from a population distribution, commonly a normal distribution with mean zero. Zero mean defines the fixed intercept as the population centre. It does not require every sampled person's departure to cancel exactly.

The random-effects distribution allows information to be shared across people. A person with few noisy observations has a less certain individual departure than someone with many informative observations. Predictions therefore tend to move toward the shared relationship, an effect often called partial pooling or shrinkage.

This is different from fitting an unrelated line for every person. Separate lines do not use a shared distribution to stabilize the individual estimates. The mixed model estimates common coefficients and variation across people together.

For a standard Gaussian mixed model, random effects and remaining errors are usually assumed independent of one another, with appropriate assumptions concerning the predictors. These assumptions support interpretation and estimation. They deserve scrutiny rather than being treated as automatic truths.

The model can describe clustered data beyond longitudinal follow-up. Our BP example stays with repeated visits so that only one new idea is introduced at a time. The mathematical meaning of mixed remains the same.

• Key idea: Mixed means fixed effects and random effects occur in the same model.

Chapter 5. Fixed Intercept vs Random Intercept

The fixed intercept locates the population mean at baseline. A random intercept moves one person's whole mean line up or down. The person's baseline is the sum of the shared intercept and their departure.

With the shared baseline of 130 mmHg, a departure of +10 gives a baseline of 140. A departure of -10 gives a baseline of 120. If both people share the same slope, their mean trajectories remain parallel.

One offset, repeated across rows

The important feature is that each person has one intercept departure reused at every visit. It is not a new independent draw for each BP reading. Reusing the same offset creates a shared source of variation across that person's observations.

Suppose one person has a positive departure. Even before considering measurement-specific errors, all their visits tend to sit above the population line. Another person with a negative departure tends to sit below it.

Variance measures spread as the mean squared departure from a mean. Its square root is the standard deviation (SD). Correlation ranges from -1 to +1: positive values indicate a tendency to vary together. Zero indicates no linear association, while negative values indicate opposite movement.

Let \(\tau_0^2\) denote the population variance of the intercept departures and \(\sigma^2\) the variance of independent remaining errors. Assume these components are independent. Two distinct visits from the same person then share covariance \(\tau_0^2\).

Covariance describes how two quantities vary together. It has squared outcome units; correlation rescales it to a unitless measure. A positive shared intercept variance induces positive covariance between visits under these assumptions.

Moving the time origin

Centred time means subtracting a chosen reference month from each time value. If we centre at month two, the intercept now describes mean BP at month two. The model's predicted BP values remain unchanged when the coefficients are transformed consistently.

In a random-intercept-only model, centring changes the meaning of the fixed intercept while the shared person offset stays the same. Once person-specific slopes are introduced, changing the origin also changes which person-specific level is called the random intercept. This matters when interpreting intercept-slope association.

A random intercept is not a standard error for a person's baseline. Its population spread describes heterogeneity. The uncertainty in estimating that spread, or in predicting one person's departure, is another quantity.

• Key idea: A random intercept is one person-specific offset shared across that person's measurements.

Chapter 6. Fixed Slope vs Random Slope

Parallel BP trajectories may be too restrictive. Some people improve rapidly, while others change slowly. We need a way to retain an average monthly change while allowing individual changes to differ.

Let \(b_{1i}\) denote person \(i\)'s departure from the fixed time slope \(\beta_1\). Their total slope is \(\beta_1+b_{1i}\). Adding the random slope gives:

\[Y_{ij}=(\beta_0+b_{0i})+(\beta_1+b_{1i})t_{ij}+\varepsilon_{ij}.\]

The fixed slope remains in the equation. It describes the population's mean change, while the random slope describes how one person differs from that mean. Replacing the fixed slope with a mean-zero random departure would impose a zero population slope unless another term supplied it.

Read the units

An intercept departure has units of mmHg. A time-slope departure has units of mmHg per month. Multiplying the latter by months produces an outcome-scale departure.

At baseline, time equals zero, so a slope departure contributes nothing to BP at that instant. Its effect grows in magnitude as we move from the reference time. This explains why people can become more or less separated during follow-up.

People who start high might decline faster, decline slower, or show no systematic relation between starting level and change. An intercept-slope covariance can represent these possibilities. Its interpretation depends on the time origin.

What the data must support

Learning about slope variation requires within-person variation in time and enough independent people. Two measurements can define a noisy personal change, but disentangling several variance components can still be difficult. More visits help describe trajectories; more people help characterize the population distribution.

A random time slope also assumes a particular shape for individual departures: linear in the chosen time scale. It cannot automatically represent every person's curved BP course. Curvature may require a different mean function, additional random terms, or a more suitable scientific model.

Individual trajectory plots are therefore useful before fitting. They show whether straight lines are plausible and whether people appear to differ mainly in level, change, or both. They guide a model question rather than prove a variance structure.

• Key idea: A random slope adds a person-specific deviation to the fixed slope.

Chapter 7. Four Combinations of Intercepts and Slopes

From components to R and Stata formulas

The R formulas show model structure: A uses lm; B–D use lme4::lmer. The term (1 + X | id) adds personal deviations; it does not remove the fixed intercept and slope already in Y ~ X. In C, noconstant removes only the random intercept from the id equation; the fixed intercept remains. The lesson adds reml to Stata mixed commands to match R estimation.

In the generic WSS specification, WSS ~ section * group + covariates + (1 | id) is Model B with a fixed section-by-group interaction. The alternative WSS ~ section * group + covariates + (1 + position | id) is Model D. The random position slope adds personal deviations while the fixed effects still estimate the mean group contrasts at each section. Here covariates is a placeholder for adjustment variables justified by the research question, not executable code by itself.

Use \(Y_{ij}\) for person \(i\)'s BP at visit \(j\), \(t_{ij}\) for months, \(\beta_0,\beta_1\) for shared coefficients, \(b_{0i},b_{1i}\) for personal departures, and \(\varepsilon_{ij}\) for remaining error. The fully annotated equation is:

\[Y_{ij}=\underbrace{\beta_0}_{\text{average intercept}}+\underbrace{b_{0i}}_{\text{person's intercept departure (random effect)}}+\left(\underbrace{\beta_1}_{\text{average slope}}+\underbrace{b_{1i}}_{\text{person's slope departure (random effect)}}\right)t_{ij}+\underbrace{\varepsilon_{ij}}_{\text{residual}}.\]

# These synthetic distributions are invented for teaching, never study data.
# Run from this termcast directory. CSVs are shared unchanged by R and Stata.
options(digits = 17)
save_once <- function(d, example, filename = "simulated.csv") {
  dir.create(example, showWarnings = FALSE, recursive = TRUE)
  p <- file.path(example, filename)
  if (!file.exists(p)) write.csv(d, p, row.names = FALSE)
  invisible(p)
}
set.seed(1067001)
n <- 80; k <- 6
id <- rep(seq_len(n), each = k)
X <- rep(seq(-1, 1, length.out = k), n)
u0 <- rnorm(n, 0, 1)
u1 <- 0.2 * u0 + rnorm(n, 0, 0.75)
Y <- 4 + 1.5 * X + u0[id] + u1[id] * X + rnorm(n*k, 0, 0.6)
save_once(data.frame(id, X, Y), "ch07-four-models")

set.seed(1061101)
n <- 64; k <- 5
id <- rep(seq_len(n), each = k)
month <- rep(0:4, n)
treatment <- rep(rep(0:1, length.out = n), each = k)
u0 <- rnorm(n, 0, 4)
u1 <- 0.08 * u0 + rnorm(n, 0, 1)
BP <- 120 - month - 3*treatment - 0.8*month*treatment +
  u0[id] + u1[id]*month + rnorm(n*k, 0, 2)
save_once(data.frame(id, month, treatment, BP), "ch11-month")

set.seed(1061301)
n <- 72; k <- 4
id <- rep(seq_len(n), each = k)
position <- rep(0:3, n)
section <- position + 1
group <- rep(rep(0:1, length.out = n), each = k)
u0 <- rnorm(n, 0, 1)
u1 <- 0.15*u0 + rnorm(n, 0, 0.45)
y <- 5 + c(0, 0.4, 1, 1.2)[section] +
  group*c(0.4, 0.7, 0.2, 1)[section] +
  u0[id] + u1[id]*position + rnorm(n*k, 0, 0.5)
save_once(data.frame(id, position, section, group, y),
          "ch13-position")

set.seed(1061401)
n <- 96; k <- 4
id <- rep(seq_len(n), each = k)
section <- rep(1:4, n)
group <- rep(rep(0:1, length.out = n), each = k)
u0 <- rnorm(n, 0, 0.8)
y <- 5 + c(0, 0.4, 1, 1.2)[section] + 0.5*group +
  u0[id] + rnorm(n*k, 0, c(0.35, 0.6, 0.9, 1.2)[section])
save_once(data.frame(id, section, group, y), "ch14-residuals")

set.seed(1061601)
nw <- 32; np <- 8; nr <- 3
ward <- rep(seq_len(nw), each = np*nr)
patient <- rep(rep(seq_len(np), each = nr), nw)
treatment <- (patient - 1) %% 2
x <- rep(c(-1, 0, 1), nw*np)
wi <- rnorm(nw, 0, 1)
ws <- 0.2*wi + rnorm(nw, 0, 0.7)
pi <- rnorm(nw*np, 0, 0.6)
pid <- (ward - 1)*np + patient
y <- 5 + 0.7*x + 0.5*treatment + wi[ward] +
  ws[ward]*treatment + pi[pid] + rnorm(length(x), 0, 0.4)
save_once(data.frame(ward, patient, treatment, x, y),
          "ch16-grouping", "simulated-nested.csv")
set.seed(1061602)
patient <- rep(seq_len(50), each = 12)
rater <- rep(seq_len(12), 50)
pi <- rnorm(50, 0, 1)
ri <- rnorm(12, 0, 0.6)
y <- 5 + pi[patient] + ri[rater] + rnorm(600, 0, 0.5)
save_once(data.frame(patient, rater, y), "ch16-grouping",
          "simulated-crossed.csv")
cat("Six seeded synthetic CSV files ready; existing files preserved.\n")

setwd("ch07-four-models")
suppressPackageStartupMessages(library(lme4))
d <- read.csv("simulated.csv")

# Fixed = population average; random = deviation.
# These are synthetic observations, not study data.
cat(nrow(d), "simulated rows\n")

print(head(d, 6))

A <- lm(Y ~ X, data = d)
f <- data.frame(Estimate = coef(A))
cat("A: population coefficients fitted\n")

print(f)

v <- data.frame(residual_variance = sigma(A)^2)
cat("Variance scale; covariance may be negative.\n")

print(v)

B <- lmer(Y ~ X + (1 | id), data = d)
f <- data.frame(Estimate = fixef(B))
cat("B: population coefficients fitted\n")

print(f)

v <- as.data.frame(VarCorr(B))
v <- v[, c("grp", "var1", "var2", "vcov")]
v[is.na(v)] <- ""
cat("Variance scale; covariance may be negative.\n")

print(v)

C <- lmer(Y ~ X + (0 + X | id), data = d)
f <- data.frame(Estimate = fixef(C))
cat("C: population coefficients fitted\n")

print(f)

v <- as.data.frame(VarCorr(C))
v <- v[, c("grp", "var1", "var2", "vcov")]
v[is.na(v)] <- ""
cat("Variance scale; covariance may be negative.\n")

print(v)

D <- lmer(Y ~ X + (1 + X | id), data = d)
f <- data.frame(Estimate = fixef(D))
cat("D: population coefficients fitted\n")

print(f)

v <- as.data.frame(VarCorr(D))
v <- v[, c("grp", "var1", "var2", "vcov")]
v[is.na(v)] <- ""
cat("Variance scale; covariance may be negative.\n")

print(v)

setwd("..")

cd "ch07-four-models"
clear all
set more off
quietly import delimited ///
  "simulated.csv", clear case(preserve) asdouble
* Fixed = population average; random = deviation.
* These are synthetic observations, not study data.
quietly mkmat id X Y in 1/6, matrix(preview)
display _N " simulated rows"

matrix list preview

quietly regress Y X
estimates store A
matrix f = (_b[_cons] \ ///
  _b[X])
matrix colnames f = Estimate
matrix rownames f = Intercept ///
  X
display "A: population coefficients fitted"

matrix list f

matrix v = (e(rmse)^2)
matrix colnames v = residual_variance
display "Variance scale; covariance may be negative."

matrix list v

quietly mixed Y X || id:, reml
estimates store B
matrix f = (_b[_cons] \ ///
  _b[X])
matrix colnames f = Estimate
matrix rownames f = Intercept ///
  X
display "B: population coefficients fitted"

matrix list f

quietly estat recovariance
matrix v = r(Cov2)
display "Variance scale; covariance may be negative."

matrix list v

scalar residual_var = exp(2*_b[lnsig_e:_cons])
display "Residual variance = " %9.6f residual_var

display residual_var

quietly mixed Y X || id: X, noconstant reml
estimates store C
matrix f = (_b[_cons] \ ///
  _b[X])
matrix colnames f = Estimate
matrix rownames f = Intercept ///
  X
display "C: population coefficients fitted"

matrix list f

quietly estat recovariance
matrix v = r(Cov2)
display "Variance scale; covariance may be negative."

matrix list v

scalar residual_var = exp(2*_b[lnsig_e:_cons])
display "Residual variance = " %9.6f residual_var

display residual_var

quietly mixed Y X || id: X, ///
  covariance(unstructured) reml
estimates store D
matrix f = (_b[_cons] \ ///
  _b[X])
matrix colnames f = Estimate
matrix rownames f = Intercept ///
  X
display "D: population coefficients fitted"

matrix list f

quietly estat recovariance
matrix v = r(Cov2)
display "Variance scale; covariance may be negative."

matrix list v

scalar residual_var = exp(2*_b[lnsig_e:_cons])
display "Residual variance = " %9.6f residual_var

display residual_var

cd ".."

Reconstruct two synthetic trajectories

Choose a shared intercept of 130 and a shared slope of -3.5. Give A an intercept departure of +10 and a slope departure of -1.5. Give B an intercept departure of -10 and a slope departure of +1.5.

A starts at 140 and changes by -5 per month. B starts at 120 and changes by -2 per month. At month three, A reaches 125 and B reaches 114, exactly reconstructing the table.

This is a synthetic decomposition of two straight lines. It is not an estimate of a population random-effect variance, a fitted software result, or evidence that two people are enough to fit a reliable mixed model. We have deliberately set the residuals to zero to isolate the arithmetic.

What each restriction allows

Model A supplies one mean line, so it cannot represent persistent personal levels or slopes in its mean structure. Model B allows different levels but keeps trajectories parallel. It can express A starting above B, but cannot reproduce their different rates exactly without residual departures.

Model C permits different rates but requires a common personal mean level at time zero because it has no random intercept. This is a substantive restriction tied to the chosen origin. It cannot reproduce the two distinct baseline levels in this table as exact personal mean trajectories.

Model D allows both personal levels and personal rates, so the illustrative decomposition fits its structure. That flexibility is useful only when supported by the study design and enough information. More flexible does not automatically mean more credible.

The A–D category table compares model components, and the two-person BP table is a hand-built illustration; distinguish these from the recorded fits to simulated data in the lesson. Actual fitted fixed coefficients and residuals would depend on the data and estimation procedure. We should never present this hand-built illustration as a model comparison with estimated likelihoods or significance tests.

• Key idea: Random classifies a varying coefficient whose components are a fixed population average and an individual deviation. Model D still contains both \(\beta\) and \(b\).

Chapter 8. Random-Effect Variance vs Residual Variance

Random-effect variance describes how personal mean trajectories differ. Residual variance describes remaining variation around those trajectories under the chosen model. A large amount of one does not imply a large amount of the other.

Keep \(\tau_0^2\) for random-intercept variance and \(\sigma^2\) for independent residual variance. In a random-intercept model, the marginal variance of a BP observation is their sum. Marginal means that we average over the population distribution of personal offsets.

\[\operatorname{Var}(Y_{ij})=\tau_0^2+\sigma^2.\]

For two distinct visits with the same residual variance and no residual correlation, define the intraclass correlation coefficient, abbreviated ICC, as their common within-person correlation. It is:

\[\operatorname{ICC}=\frac{\tau_0^2}{\tau_0^2+\sigma^2}.\]

For a synthetic variance of 64 between people and 36 remaining within-person variance, the total is 100 and the ICC is 0.64. The corresponding standard deviations are 8 and 6 mmHg. These numbers illustrate the definition and are not estimates from the two trajectories.

Conditional and marginal views

Conditional means holding a person's random effects fixed. Under independent residual errors, their remaining measurement departures are then independent. Marginally, those same measurements are correlated because they share an unknown personal offset across the population.

There is no contradiction. Knowing a person's persistent level removes one source of dependence. Averaging over uncertainty in that shared level restores it.

With a random slope

Let \(\tau_1^2\) be slope-departure variance and \(\tau_{01}\) the intercept-slope covariance. At time \(t\), the marginal BP variance becomes:

\[\operatorname{Var}\{Y_i(t)\}=\tau_0^2+2t\tau_{01}+t^2\tau_1^2+\sigma^2.\]

Between-person spread can now change with time. Correlations between visits can depend on both visit times, so one constant ICC no longer summarizes every pair. Reporting a single number without defining the pair or reference time can conceal this structure.

A random-effect standard deviation measures heterogeneity in outcome or slope units. A standard error measures uncertainty in an estimated coefficient or statistic. Confusing them makes a model appear to answer a different question from the one it actually answers.

• Key idea: Between-person variation and remaining measurement variation answer different questions.

Documentation supporting the model and software details: lme4 model formulation.

Part III. From Model Equations to Real Repeated-Measures Analysis

Chapter 9. Fixed Interaction Is Not a Random Slope

Suppose the BP study contains two treatment groups. One question asks whether their average BP changes differ. Another asks whether people within the same treatment group still have different changes.

Let \(T_i\) equal one for the comparison treatment and zero for reference treatment. Let \(\beta_T\) describe their baseline mean difference and \(\beta_{tT}\) describe their difference in monthly change. A model with both an interaction and individual trajectory variation is:

Here \(\beta_t\) names the fixed time slope in the reference treatment group, previously written as \(\beta_1\).

\[Y_{ij}=\beta_0+\beta_tt_{ij}+\beta_TT_i+\beta_{tT}t_{ij}T_i+b_{0i}+b_{1i}t_{ij}+\varepsilon_{ij}.\]

The reference group's mean slope is \(\beta_t\). The comparison group's mean slope is \(\beta_t+\beta_{tT}\). Within either group, person \(i\)'s slope also includes the departure \(b_{1i}\).

Systematic difference versus remaining heterogeneity

A fixed treatment-by-time interaction represents a systematic difference between group mean trends. A random time slope represents person-to-person variation left after the fixed part has been specified. Both can belong in the same model because they answer different questions.

For illustration, the reference group might average -3.5 mmHg per month while the comparison group averages -4.5. The interaction is then -1 mmHg per month. Individuals may still depart above or below their own group's average slope.

Adding random slopes cannot substitute for a scientifically required treatment interaction. Without the interaction, the model imposes the same mean slope across groups, even if it allows individuals to vary. Conversely, an interaction cannot guarantee that all people within a group have parallel trajectories.

Separate between-person and within-person information

Time varies within a person, while assigned treatment generally does not. A person's treatment indicator therefore cannot directly reveal how that same person's BP would have evolved under both assignments. Causal interpretation also depends on allocation, follow-up, and the target question.

Other changing predictors can carry both within-person and between-person information. For example, a higher BP-related covariate within one person need not have the same association as a higher average level across people. Including person means alongside within-person deviations can make that distinction explicit when relevant.

Interpret interactions at stated reference values. A baseline treatment coefficient applies at the chosen time origin; after centring time at month two, it becomes the mean treatment difference at month two. Centring changes the question attached to coefficients without changing fitted values when done consistently.

• Key idea: A fixed interaction describes systematic group differences; a random slope describes remaining person-to-person differences.

Chapter 10. How Mixed Models Handle Correlation and Variance

A random effect is one source of dependence. Remaining errors may also be correlated, or have different variances at different visits. These choices describe different parts of the same covariance model.

Let \(\mathbf Y_i\) collect a person's BP measurements. Let \(X_i\) contain their fixed-effect predictor rows and \(\beta\) the shared coefficients. Let \(Z_i\) contain the predictor rows attached to the person's random-effect vector \(\mathbf b_i\), and let \(\varepsilon_i\) collect the remaining errors.

Write \(G\) for the covariance matrix of the random effects and \(R_i\) for the residual covariance matrix. A covariance matrix records variances on its diagonal and covariances between distinct quantities off the diagonal. Assuming independence between random effects and residual errors gives:

The symbol \(\top\) means transpose: turn a row into a column, or interchange the rows and columns of a matrix. For an intercept and time slope, each row of \(Z_i\) is \((1,t_{ij})\). Matrix multiplication adds the shared random contributions to each pair of measurements.

\[\mathbf Y_i=X_i\beta+Z_i\mathbf b_i+\varepsilon_i.\]

\[\operatorname{Var}(\mathbf Y_i)=Z_iGZ_i^\top+R_i.\]

Compare the covariance choices

The equal-variance random-intercept model has compound symmetry: the same marginal variance at each visit and the same covariance between every distinct pair. With unequal residual variances, its common covariance survives but its correlations generally differ. Correlation divides covariance by the two marginal standard deviations.

For AR(1), adjacent equally spaced visits have one residual correlation parameter; correlations at longer lags are powers of that parameter. Unequal actual intervals need a suitable time representation or continuous-time correlation model. A visit number is not automatically a substitute for elapsed time.

Alternatives and redundancy

A marginal model can estimate the population mean and model its covariance directly without personal random-effect predictions. Generalized estimating equations, or GEE, target mean relationships using a working correlation and appropriate uncertainty estimates. Cluster-robust standard errors can protect mean-coefficient inference against some covariance misspecification when there are enough independent clusters and suitable corrections.

None is universally the best method. Choose according to the target quantity, such as a treatment difference at month three. Also consider cluster count, missingness assumptions, design, and whether individual predictions are needed. Robust uncertainty does not repair a wrong mean model or informative missingness.

Do not combine unrestricted covariance pieces indiscriminately. A random intercept plus an unrestricted residual covariance can assign the same covariance to either component, preventing separate identification. A model must distinguish the information available in the total covariance from the decomposition we wish to interpret.

• Key idea: Model covariance to answer the scientific question; do not rank methods by a universal best-to-worst ladder.

Documentation supporting the model and software details: nlme: covariance and unequal residual variances; nlme: AR(1) time covariate.

Chapter 11. Fit, Interpret, and Check a Mixed Model

Start with the intended contrast, the person identifier, actual visit times, and a defensible mean function. Plot the BP observations and missing visits before fitting. Decide which coefficients describe the scientific question and which covariance components the design can support.

Read real R and Stata output on simulated data

The lesson uses one seeded simulated long-format dataset. R creates the CSV and Stata reads the same file, with the same treatment reference group. Compare the random-intercept model with a random month slope, then inspect covariance, singularity and residuals. The displayed results come from real computation on simulated data.

# These synthetic distributions are invented for teaching, never study data.
# Run from this termcast directory. CSVs are shared unchanged by R and Stata.
options(digits = 17)
save_once <- function(d, example, filename = "simulated.csv") {
  dir.create(example, showWarnings = FALSE, recursive = TRUE)
  p <- file.path(example, filename)
  if (!file.exists(p)) write.csv(d, p, row.names = FALSE)
  invisible(p)
}
set.seed(1067001)
n <- 80; k <- 6
id <- rep(seq_len(n), each = k)
X <- rep(seq(-1, 1, length.out = k), n)
u0 <- rnorm(n, 0, 1)
u1 <- 0.2 * u0 + rnorm(n, 0, 0.75)
Y <- 4 + 1.5 * X + u0[id] + u1[id] * X + rnorm(n*k, 0, 0.6)
save_once(data.frame(id, X, Y), "ch07-four-models")

set.seed(1061101)
n <- 64; k <- 5
id <- rep(seq_len(n), each = k)
month <- rep(0:4, n)
treatment <- rep(rep(0:1, length.out = n), each = k)
u0 <- rnorm(n, 0, 4)
u1 <- 0.08 * u0 + rnorm(n, 0, 1)
BP <- 120 - month - 3*treatment - 0.8*month*treatment +
  u0[id] + u1[id]*month + rnorm(n*k, 0, 2)
save_once(data.frame(id, month, treatment, BP), "ch11-month")

set.seed(1061301)
n <- 72; k <- 4
id <- rep(seq_len(n), each = k)
position <- rep(0:3, n)
section <- position + 1
group <- rep(rep(0:1, length.out = n), each = k)
u0 <- rnorm(n, 0, 1)
u1 <- 0.15*u0 + rnorm(n, 0, 0.45)
y <- 5 + c(0, 0.4, 1, 1.2)[section] +
  group*c(0.4, 0.7, 0.2, 1)[section] +
  u0[id] + u1[id]*position + rnorm(n*k, 0, 0.5)
save_once(data.frame(id, position, section, group, y),
          "ch13-position")

set.seed(1061401)
n <- 96; k <- 4
id <- rep(seq_len(n), each = k)
section <- rep(1:4, n)
group <- rep(rep(0:1, length.out = n), each = k)
u0 <- rnorm(n, 0, 0.8)
y <- 5 + c(0, 0.4, 1, 1.2)[section] + 0.5*group +
  u0[id] + rnorm(n*k, 0, c(0.35, 0.6, 0.9, 1.2)[section])
save_once(data.frame(id, section, group, y), "ch14-residuals")

set.seed(1061601)
nw <- 32; np <- 8; nr <- 3
ward <- rep(seq_len(nw), each = np*nr)
patient <- rep(rep(seq_len(np), each = nr), nw)
treatment <- (patient - 1) %% 2
x <- rep(c(-1, 0, 1), nw*np)
wi <- rnorm(nw, 0, 1)
ws <- 0.2*wi + rnorm(nw, 0, 0.7)
pi <- rnorm(nw*np, 0, 0.6)
pid <- (ward - 1)*np + patient
y <- 5 + 0.7*x + 0.5*treatment + wi[ward] +
  ws[ward]*treatment + pi[pid] + rnorm(length(x), 0, 0.4)
save_once(data.frame(ward, patient, treatment, x, y),
          "ch16-grouping", "simulated-nested.csv")
set.seed(1061602)
patient <- rep(seq_len(50), each = 12)
rater <- rep(seq_len(12), 50)
pi <- rnorm(50, 0, 1)
ri <- rnorm(12, 0, 0.6)
y <- 5 + pi[patient] + ri[rater] + rnorm(600, 0, 0.5)
save_once(data.frame(patient, rater, y), "ch16-grouping",
          "simulated-crossed.csv")
cat("Six seeded synthetic CSV files ready; existing files preserved.\n")

setwd("ch11-month")
suppressPackageStartupMessages(library(lme4))
d <- read.csv("simulated.csv")
d$treatment <- factor(d$treatment)
# Fixed = population average; random = deviation.
# These are synthetic observations, not study data.
cat(nrow(d), "simulated rows\n")

print(head(d, 6))

M0 <- lmer(BP ~ month*treatment + (1 | id),
  data = d)
f <- data.frame(Estimate = fixef(M0))
cat("M0: population coefficients fitted\n")

print(f)

v <- as.data.frame(VarCorr(M0))
v <- v[, c("grp", "var1", "var2", "vcov")]
v[is.na(v)] <- ""
cat("Variance scale; covariance may be negative.\n")

print(v)

M1 <- lmer(BP ~ month*treatment +
  (1 + month | id), data = d)
f <- data.frame(Estimate = fixef(M1))
cat("M1: population coefficients fitted\n")

print(f)

v <- as.data.frame(VarCorr(M1))
v <- v[, c("grp", "var1", "var2", "vcov")]
v[is.na(v)] <- ""
cat("Variance scale; covariance may be negative.\n")

print(v)

singular <- isSingular(M1, tol = 1e-4)
print(singular)

print(singular)

plot(fitted(M1), resid(M1), pch = 16, cex = 0.5,
  xlab = "Fitted / ค่าทำนาย",
  ylab = "Residual / residual",
  main = "Simulated / ข้อมูลจำลอง")
abline(h = 0, lty = 2)

setwd("..")

cd "ch11-month"
clear all
set more off
quietly import delimited ///
  "simulated.csv", clear case(preserve) asdouble
* Fixed = population average; random = deviation.
* These are synthetic observations, not study data.
quietly mkmat id month treatment BP in 1/6, matrix(preview)
display _N " simulated rows"

matrix list preview

quietly mixed BP c.month##i.treatment || id:, reml
estimates store M0
matrix f = (_b[_cons] \ ///
  _b[month] \ ///
  _b[1.treatment] \ ///
  _b[1.treatment#c.month])
matrix colnames f = Estimate
matrix rownames f = Intercept ///
  month ///
  treatment1 ///
  month_treatment1
display "M0: population coefficients fitted"

matrix list f

quietly estat recovariance
matrix v = r(Cov2)
display "Variance scale; covariance may be negative."

matrix list v

scalar residual_var = exp(2*_b[lnsig_e:_cons])
display "Residual variance = " %9.6f residual_var

display residual_var

quietly mixed BP c.month##i.treatment || id: month, ///
  covariance(unstructured) reml
estimates store M1
matrix f = (_b[_cons] \ ///
  _b[month] \ ///
  _b[1.treatment] \ ///
  _b[1.treatment#c.month])
matrix colnames f = Estimate
matrix rownames f = Intercept ///
  month ///
  treatment1 ///
  month_treatment1
display "M1: population coefficients fitted"

matrix list f

quietly estat recovariance
matrix v = r(Cov2)
display "Variance scale; covariance may be negative."

matrix list v

scalar residual_var = exp(2*_b[lnsig_e:_cons])
display "Residual variance = " %9.6f residual_var

display residual_var

estimates restore M1
quietly predict double residual, residuals
quietly predict double fitted, fitted
quietly scatter residual fitted, yline(0) ///
  title("Simulated / ข้อมูลจำลอง") ///
  xtitle("Fitted / ค่าทำนาย") ///
  ytitle("Residual / residual")

cd ".."

The Stata covariance option concerns the random intercept and random month slope in this command. It does not request an unrestricted residual covariance across visits. Matching software requires matching the mean model, random terms, residual assumptions, estimation method, data, and reference coding.

Read output in layers

Fixed coefficients describe baseline, change, group differences, and interactions at stated reference values. Their standard errors and confidence intervals describe estimation uncertainty. A random-effect SD describes population heterogeneity, while residual SD describes the remaining outcome-scale spread.

A random-effect correlation describes association between personal intercept and slope departures. Its meaning changes with time centring. A correlation close to either boundary may signal insufficient information for the chosen covariance structure.

Estimation and checks

Maximum likelihood (ML) selects parameters that make the observed data most compatible with the assumed probability model. Likelihood measures that compatibility as the parameters change. Restricted maximum likelihood (REML) estimates covariance after accounting for the fixed effects that must also be estimated. Likelihood comparisons of different fixed-effect structures generally use ML on identical observations; REML comparisons require the same fixed-effect structure.

A converged optimizer has met its numerical stopping criteria. It does not certify that the mean, covariance, or missingness assumptions are adequate. A singular fit means the estimated random-effect covariance has reduced rank, for example a zero variance or a perfectly dependent combination of random effects.

Check residuals against fitted values, time, treatment, and person. Examine departures from linearity, changing spread, unusual individuals, and remaining serial patterns. Normality plots are diagnostic aids, not a replacement for assessing whether the scientific contrast is stable under plausible alternatives.

For variance components on a boundary, an ordinary likelihood-ratio test reference distribution may be inappropriate. Use a justified procedure rather than reading every model comparison as a routine chi-squared test. Report sensitivity analyses that materially change uncertainty or interpretation.

• Key idea: Recorded simulated results make the code concrete; convergence, singularity, assumptions, and uncertainty still need separate assessment.

Documentation supporting the model and software details: lme4: singular fits; lme4: estimation and model specification; Stata: mixed; Stata: mixed postestimation.

Part IV. Advanced Mixed Models: From Time to Aortic Sections

Chapter 12. Repeated Measures Without Time: Anatomical Sections

We now move explicitly from BP follow-up to repaired tetralogy of Fallot, abbreviated TOF. TOF is a congenital heart condition; repaired TOF describes people who have undergone surgical repair. Four-dimensional flow magnetic resonance imaging, or 4D flow MRI, can characterize blood flow across three spatial dimensions over the cardiac cycle.

Wall shear stress, abbreviated WSS, describes tangential stress associated with flowing blood at the vessel wall. An imaging-derived WSS measure needs a stated definition, unit, processing method, and location. We use illustrative WSS in pascals (Pa) below, without claiming any observed clinical finding.

Four locations in one person

The four teaching labels are AAO, STJ, Before branching, and DAO. AAO denotes the ascending aorta, STJ the sinotubular junction, and DAO the descending aorta. Before branching is the study's chosen location label and requires an explicit anatomical definition in a real imaging protocol.

These measurements share a person even when obtained in the same examination. Their dependence does not require follow-up time. Repeated measurements can describe locations, conditions, or visits.

Section is categorical here. Numbering the four locations one through four and fitting one numerical slope would impose equally spaced linear mean changes across those codes. Anatomical order alone does not justify that restriction.

Reference coding

Take AAO as reference. Let \(\beta_0\) be its mean, and define indicator variables \(S_{\mathrm{STJ}}\), \(S_{\mathrm{Before}}\), \(S_{\mathrm{DAO}}\), each equal to one at its named section and zero otherwise. Let the corresponding coefficients describe each section's mean difference from AAO.

\[E(Y\mid\text{section})=\beta_0+\beta_{\mathrm{STJ}}S_{\mathrm{STJ}}+\beta_{\mathrm{Before}}S_{\mathrm{Before}}+\beta_{\mathrm{DAO}}S_{\mathrm{DAO}}.\]

There are four means but only three differences from a chosen reference. A reference change alters the coefficient labels and numerical contrasts, while the four fitted means remain the same. The intercept is now a reference-section mean, not a baseline-time mean.

A random intercept adds one personal offset to all four section means. It represents a person whose WSS tends to be higher or lower throughout the profile after fixed predictors are considered. It does not by itself permit every person's section contrasts to differ.

• Key idea: Repeated measurements need a shared observational unit, not necessarily a time axis.

Chapter 13. Random Slopes for Categorical Predictors: What Can We Estimate?

1. Start with a binary section contrast

Temporarily retain only AAO and DAO. Let \(D\) be zero for AAO and one for DAO. Let \(\beta_D\) be the population mean DAO-minus-AAO difference; it is a fixed contrast.

With a random intercept \(b_{0i}\), a person's AAO mean is \(\beta_0+b_{0i}\), and their DAO mean is \(\beta_0+b_{0i}+\beta_D\). The personal offset cancels when taking DAO minus AAO. Everyone therefore has the same modelled mean contrast.

2. Allow that contrast to vary

Introduce \(b_{Di}\), a person-specific departure from the DAO contrast. The model becomes:

\[Y_{ij}=\beta_0+b_{0i}+(\beta_D+b_{Di})D_{ij}+\varepsilon_{ij}.\]

For a synthetic illustration, set the AAO population mean to 1.0 Pa and the DAO contrast to -0.2 Pa. Give one person a +0.1 personal level offset and a +0.1 contrast departure: their two means are 1.1 and 1.0. Give another person a -0.1 level offset and a -0.1 contrast departure: their means are 0.9 and 0.6.

Their DAO-minus-AAO contrasts are -0.1 and -0.3 Pa. A categorical random slope means a varying contrast. It does not require physical time or a continuous predictor.

3. Expand to four levels

Four section levels require three independent contrasts. A full (1 + section | id) structure contains four random coefficients per person: an intercept and three contrast departures. With an unrestricted covariance among these coefficients, it allows a fully flexible random profile.

Allowing only one prespecified contrast departure is a smaller model. Such a restricted structure may be identifiable from the available covariance patterns. The failure of a saturated structure is not a prohibition against all categorical random slopes.

4. Explain what cannot be separated

For the two-section case, let \(H\) be the covariance of the person's latent AAO and DAO mean departures, with diagonal entries \(h_{11},h_{22}\) and covariance \(h_{12}\). With one reading per section and independent equal residual variance \(\sigma^2\), the observed covariance \(V\) is:

\[V=\begin{pmatrix}h_{11}+\sigma^2&h_{12}\\h_{12}&h_{22}+\sigma^2\end{pmatrix}.\]

Two variances and one covariance provide three observable covariance quantities. The latent covariance contributes three unknowns, and free residual variance adds a fourth. Counting suggests a problem; an explicit construction proves it.

Let \(\delta\) be a small positive variance amount below the residual variance, and let \(I\) be the identity matrix.

The identity matrix has ones on its diagonal and zeros elsewhere, so adding \(\delta I\) changes only the diagonal variances.

Replace \(H\) with \(H+\delta I\) and replace \(\sigma^2\) with \(\sigma^2-\delta\). Both covariance components remain valid for a suitable \(\delta\), while their sum \(V\) is unchanged.

For example, observed diagonal variances of 0.10 and covariance 0.04 can arise from latent diagonal variances of 0.06 plus residual variance 0.04. They can also arise from latent diagonal variances of 0.08 plus residual variance 0.02, retaining the same latent covariance 0.04. Both latent covariance matrices are mathematically valid and yield the same observed covariance.

A four-section covariance matrix has four diagonal variances and six distinct off-diagonal covariances: ten quantities. A full latent random profile already uses ten covariance parameters; adding free residual variance gives eleven. The same diagonal-shifting construction leaves the observed four-by-four covariance unchanged.

The intercept-plus-three-contrasts representation is an invertible recoding of the four latent section levels. It therefore has the same identifiability limitation. More people estimate the total covariance more precisely, but do not resolve this structural ambiguity by themselves.

Within-section replicate readings, justified covariance restrictions, or externally known measurement error can add identifying information or assumptions. An unrestricted marginal covariance estimates the total profile covariance without claiming its decomposition. A zero residual variance constraint changes the interpretation and needs explicit scientific justification.

Use the same \(H\), \(V\), and \(h_{ab}\) notation for all four section levels, where \(a,b\) index sections. The diagonal entries of \(H\) are latent variances; its off-diagonal entries are latent covariances. The observed covariance is:

\[V=\begin{pmatrix}h_{11}+\sigma^2&h_{12}&h_{13}&h_{14}\\h_{12}&h_{22}+\sigma^2&h_{23}&h_{24}\\h_{13}&h_{23}&h_{33}+\sigma^2&h_{34}\\h_{14}&h_{24}&h_{34}&h_{44}+\sigma^2\end{pmatrix}.\]

5. One random slope along a numeric position

The limit in subsection 4 concerns the full unrestricted four-section random profile plus free residual variance. It does not prohibit every random slope. Keep the fixed section-by-group interaction categorical, so that the group contrast at each section remains the question, and add one numeric position slope:

WSS ~ section * group + covariates + (1 + position | id)

# These synthetic distributions are invented for teaching, never study data.
# Run from this termcast directory. CSVs are shared unchanged by R and Stata.
options(digits = 17)
save_once <- function(d, example, filename = "simulated.csv") {
  dir.create(example, showWarnings = FALSE, recursive = TRUE)
  p <- file.path(example, filename)
  if (!file.exists(p)) write.csv(d, p, row.names = FALSE)
  invisible(p)
}
set.seed(1067001)
n <- 80; k <- 6
id <- rep(seq_len(n), each = k)
X <- rep(seq(-1, 1, length.out = k), n)
u0 <- rnorm(n, 0, 1)
u1 <- 0.2 * u0 + rnorm(n, 0, 0.75)
Y <- 4 + 1.5 * X + u0[id] + u1[id] * X + rnorm(n*k, 0, 0.6)
save_once(data.frame(id, X, Y), "ch07-four-models")

set.seed(1061101)
n <- 64; k <- 5
id <- rep(seq_len(n), each = k)
month <- rep(0:4, n)
treatment <- rep(rep(0:1, length.out = n), each = k)
u0 <- rnorm(n, 0, 4)
u1 <- 0.08 * u0 + rnorm(n, 0, 1)
BP <- 120 - month - 3*treatment - 0.8*month*treatment +
  u0[id] + u1[id]*month + rnorm(n*k, 0, 2)
save_once(data.frame(id, month, treatment, BP), "ch11-month")

set.seed(1061301)
n <- 72; k <- 4
id <- rep(seq_len(n), each = k)
position <- rep(0:3, n)
section <- position + 1
group <- rep(rep(0:1, length.out = n), each = k)
u0 <- rnorm(n, 0, 1)
u1 <- 0.15*u0 + rnorm(n, 0, 0.45)
y <- 5 + c(0, 0.4, 1, 1.2)[section] +
  group*c(0.4, 0.7, 0.2, 1)[section] +
  u0[id] + u1[id]*position + rnorm(n*k, 0, 0.5)
save_once(data.frame(id, position, section, group, y),
          "ch13-position")

set.seed(1061401)
n <- 96; k <- 4
id <- rep(seq_len(n), each = k)
section <- rep(1:4, n)
group <- rep(rep(0:1, length.out = n), each = k)
u0 <- rnorm(n, 0, 0.8)
y <- 5 + c(0, 0.4, 1, 1.2)[section] + 0.5*group +
  u0[id] + rnorm(n*k, 0, c(0.35, 0.6, 0.9, 1.2)[section])
save_once(data.frame(id, section, group, y), "ch14-residuals")

set.seed(1061601)
nw <- 32; np <- 8; nr <- 3
ward <- rep(seq_len(nw), each = np*nr)
patient <- rep(rep(seq_len(np), each = nr), nw)
treatment <- (patient - 1) %% 2
x <- rep(c(-1, 0, 1), nw*np)
wi <- rnorm(nw, 0, 1)
ws <- 0.2*wi + rnorm(nw, 0, 0.7)
pi <- rnorm(nw*np, 0, 0.6)
pid <- (ward - 1)*np + patient
y <- 5 + 0.7*x + 0.5*treatment + wi[ward] +
  ws[ward]*treatment + pi[pid] + rnorm(length(x), 0, 0.4)
save_once(data.frame(ward, patient, treatment, x, y),
          "ch16-grouping", "simulated-nested.csv")
set.seed(1061602)
patient <- rep(seq_len(50), each = 12)
rater <- rep(seq_len(12), 50)
pi <- rnorm(50, 0, 1)
ri <- rnorm(12, 0, 0.6)
y <- 5 + pi[patient] + ri[rater] + rnorm(600, 0, 0.5)
save_once(data.frame(patient, rater, y), "ch16-grouping",
          "simulated-crossed.csv")
cat("Six seeded synthetic CSV files ready; existing files preserved.\n")

setwd("ch13-position")
suppressPackageStartupMessages(library(lme4))
d <- read.csv("simulated.csv")
d$section <- factor(d$section)
d$group <- factor(d$group)
# Fixed = population average; random = deviation.
# These are synthetic observations, not study data.
cat(nrow(d), "simulated rows\n")

print(head(d, 6))

P0 <- lmer(y ~ section*group + (1 | id), data = d)
f <- data.frame(Estimate = fixef(P0))
cat("P0: population coefficients fitted\n")

print(f)

v <- as.data.frame(VarCorr(P0))
v <- v[, c("grp", "var1", "var2", "vcov")]
v[is.na(v)] <- ""
cat("Variance scale; covariance may be negative.\n")

print(v)

P1 <- lmer(y ~ section*group +
  (1 + position | id), data = d)
f <- data.frame(Estimate = fixef(P1))
cat("P1: population coefficients fitted\n")

print(f)

v <- as.data.frame(VarCorr(P1))
v <- v[, c("grp", "var1", "var2", "vcov")]
v[is.na(v)] <- ""
cat("Variance scale; covariance may be negative.\n")

print(v)

setwd("..")

cd "ch13-position"
clear all
set more off
quietly import delimited ///
  "simulated.csv", clear case(preserve) asdouble
* Fixed = population average; random = deviation.
* These are synthetic observations, not study data.
quietly mkmat id position section group y in 1/6, matrix(preview)
display _N " simulated rows"

matrix list preview

quietly mixed y i.section##i.group || id:, reml
estimates store P0
matrix f = (_b[_cons] \ ///
  _b[2.section] \ ///
  _b[3.section] \ ///
  _b[4.section] \ ///
  _b[1.group] \ ///
  _b[2.section#1.group] \ ///
  _b[3.section#1.group] \ ///
  _b[4.section#1.group])
matrix colnames f = Estimate
matrix rownames f = Intercept ///
  section2 ///
  section3 ///
  section4 ///
  group1 ///
  section2_group1 ///
  section3_group1 ///
  section4_group1
display "P0: population coefficients fitted"

matrix list f

quietly estat recovariance
matrix v = r(Cov2)
display "Variance scale; covariance may be negative."

matrix list v

scalar residual_var = exp(2*_b[lnsig_e:_cons])
display "Residual variance = " %9.6f residual_var

display residual_var

quietly mixed y i.section##i.group || id: position, ///
  covariance(unstructured) reml
estimates store P1
matrix f = (_b[_cons] \ ///
  _b[2.section] \ ///
  _b[3.section] \ ///
  _b[4.section] \ ///
  _b[1.group] \ ///
  _b[2.section#1.group] \ ///
  _b[3.section#1.group] \ ///
  _b[4.section#1.group])
matrix colnames f = Estimate
matrix rownames f = Intercept ///
  section2 ///
  section3 ///
  section4 ///
  group1 ///
  section2_group1 ///
  section3_group1 ///
  section4_group1
display "P1: population coefficients fitted"

matrix list f

quietly estat recovariance
matrix v = r(Cov2)
display "Variance scale; covariance may be negative."

matrix list v

scalar residual_var = exp(2*_b[lnsig_e:_cons])
display "Residual variance = " %9.6f residual_var

display residual_var

cd ".."

For example, assign anatomical order along the flow path as STJ = 0, AAO = 1, arch (the article's “Before branching” section) = 2, and DAO = 3. Alternatively, use measured centreline distance. The random intercept gives each person an offset; the random position slope gives each person a departure that changes along this position axis. The fixed section means need not lie on a straight line.

With one reading at each of four distinct positions, the outcome vector has ten observable variance/covariance entries. This restricted model requests four covariance parameters: intercept variance, slope variance, their covariance, and one common independent residual variance. Four unknowns against ten entries passes a necessary counting check; counting alone is not proof of identifiability.

To see why the parameters can be separated here, write the position scores as a column vector \(p\), the column of ones as \(\mathbf 1\), and the identity matrix as \(I\). Then:

\[V=\tau_{00}\mathbf 1\mathbf 1^\top+\tau_{01}(\mathbf 1p^\top+p\mathbf 1^\top)+\tau_{11}pp^\top+\sigma^2 I.\]

Here \(\tau_{00}\) is intercept variance, \(\tau_{11}\) slope variance, and \(\tau_{01}\) their covariance. For four distinct position scores, these four covariance basis matrices are linearly independent. In particular, the off-diagonal entries identify the constant, sum-of-positions, and product-of-positions terms, and the remaining diagonal contribution identifies residual variance. This establishes structural identifiability under the stated residual assumptions. These are population covariance parameters estimated across independent people, not four parameters fitted separately from one person's four readings. Identifiability does not guarantee precise estimates with a small number of people, a stable fit, or a nonzero slope variance.

Unequal spacing is allowed, just as visits at 0, 3, 6, and 12 months are allowed: the slope uses the actual numeric values. Assigned scores 0, 1, 2, 3 define change per position-score step; centreline distance in cm defines change per cm. Equal score increments need not represent equal anatomical distances. The substantive assumption is that each person's departure from the mean section profile is linear in the chosen position, not that visits or anatomical locations must be equally spaced.

Check that assumption against the profiles, inspect convergence and singular fits, and compare with the random-intercept model. For covariance comparisons, keep the same observations and fixed-effects design; REML can then be used. If the fixed mean model changes, use ML rather than comparing REML likelihoods. Testing zero slope variance is a boundary problem, with the intercept–slope covariance constrained as well; an ordinary chi-squared likelihood-ratio reference is not automatically valid. Use a justified boundary-aware procedure, such as a suitable parametric bootstrap.

A random slope for one prespecified contrast, such as DAO versus the other sections, is another restricted option. It still needs an identifiable covariance specification and enough independent people. These are possible model specifications. The lesson reports simulated fits of the random-intercept and numeric-position-slope models only, not unpublished study results or a fitted full random profile.

• Key idea: The full random-profile limit does not prohibit one random position slope; identifiability depends on the covariance structure and design.

Chapter 14. Common Variance vs Section-Specific Variance

A random intercept can coexist with either one residual variance or section-specific residual variances. Those choices change remaining measurement spread. They do not introduce personal section-contrast slopes.

Let \(\tau^2\) be the random-intercept variance and let \(\sigma_s^2\) be the residual variance for section \(s\). Assume independent residual errors across sections. The marginal section variance is \(\tau^2+\sigma_s^2\), and every distinct pair shares covariance \(\tau^2\).

If the residual variances differ, the pairwise correlations generally differ because their denominators differ. A shared offset alone does not guarantee equal correlations. Equal correlations additionally require the relevant symmetry in marginal variances.

Why standard errors can be equal

With the same complete subjects at every section, equal group sizes across sections, a common covariance structure, and symmetric mean contrasts, standard errors can be exactly equal. The sections have different means but equivalent information under the model. Equal standard errors are therefore compatible with legitimate model calculations.

They do not establish equal biological variability. Observed SD, model residual SD, random-intercept SD, and SE of a mean or contrast are distinct quantities. Each answers a different question.

Let \(n\) be the number of complete independent subjects in one group. In the simple common-variance random-intercept model, the variance of a section sample mean is \((\tau^2+\sigma^2)/n\). There is no section label in that expression, explaining equal section mean SEs in this setting.

The variance of an estimated AAO-minus-DAO mean difference is instead \(2\sigma^2/n\) under the same assumptions. The shared person intercept cancels in the within-person difference. Equal SEs for the individual means do not imply that those means are independent.

Investigate assumptions, not appearance

Inspect residual spread by section and group, missingness, influential people, and fitted covariance adequacy. Consider a section-specific variance model when the observed pattern and scientific design justify it. Do not change a model merely to make a column of SEs look less repetitive.

In R's nlme framework, weights = varIdent(form = ~1 | section) permits different residual scales across section levels. It does not add random slopes. The model can keep the same random intercept and fixed group-by-section mean structure.

# These synthetic distributions are invented for teaching, never study data.
# Run from this termcast directory. CSVs are shared unchanged by R and Stata.
options(digits = 17)
save_once <- function(d, example, filename = "simulated.csv") {
  dir.create(example, showWarnings = FALSE, recursive = TRUE)
  p <- file.path(example, filename)
  if (!file.exists(p)) write.csv(d, p, row.names = FALSE)
  invisible(p)
}
set.seed(1067001)
n <- 80; k <- 6
id <- rep(seq_len(n), each = k)
X <- rep(seq(-1, 1, length.out = k), n)
u0 <- rnorm(n, 0, 1)
u1 <- 0.2 * u0 + rnorm(n, 0, 0.75)
Y <- 4 + 1.5 * X + u0[id] + u1[id] * X + rnorm(n*k, 0, 0.6)
save_once(data.frame(id, X, Y), "ch07-four-models")

set.seed(1061101)
n <- 64; k <- 5
id <- rep(seq_len(n), each = k)
month <- rep(0:4, n)
treatment <- rep(rep(0:1, length.out = n), each = k)
u0 <- rnorm(n, 0, 4)
u1 <- 0.08 * u0 + rnorm(n, 0, 1)
BP <- 120 - month - 3*treatment - 0.8*month*treatment +
  u0[id] + u1[id]*month + rnorm(n*k, 0, 2)
save_once(data.frame(id, month, treatment, BP), "ch11-month")

set.seed(1061301)
n <- 72; k <- 4
id <- rep(seq_len(n), each = k)
position <- rep(0:3, n)
section <- position + 1
group <- rep(rep(0:1, length.out = n), each = k)
u0 <- rnorm(n, 0, 1)
u1 <- 0.15*u0 + rnorm(n, 0, 0.45)
y <- 5 + c(0, 0.4, 1, 1.2)[section] +
  group*c(0.4, 0.7, 0.2, 1)[section] +
  u0[id] + u1[id]*position + rnorm(n*k, 0, 0.5)
save_once(data.frame(id, position, section, group, y),
          "ch13-position")

set.seed(1061401)
n <- 96; k <- 4
id <- rep(seq_len(n), each = k)
section <- rep(1:4, n)
group <- rep(rep(0:1, length.out = n), each = k)
u0 <- rnorm(n, 0, 0.8)
y <- 5 + c(0, 0.4, 1, 1.2)[section] + 0.5*group +
  u0[id] + rnorm(n*k, 0, c(0.35, 0.6, 0.9, 1.2)[section])
save_once(data.frame(id, section, group, y), "ch14-residuals")

set.seed(1061601)
nw <- 32; np <- 8; nr <- 3
ward <- rep(seq_len(nw), each = np*nr)
patient <- rep(rep(seq_len(np), each = nr), nw)
treatment <- (patient - 1) %% 2
x <- rep(c(-1, 0, 1), nw*np)
wi <- rnorm(nw, 0, 1)
ws <- 0.2*wi + rnorm(nw, 0, 0.7)
pi <- rnorm(nw*np, 0, 0.6)
pid <- (ward - 1)*np + patient
y <- 5 + 0.7*x + 0.5*treatment + wi[ward] +
  ws[ward]*treatment + pi[pid] + rnorm(length(x), 0, 0.4)
save_once(data.frame(ward, patient, treatment, x, y),
          "ch16-grouping", "simulated-nested.csv")
set.seed(1061602)
patient <- rep(seq_len(50), each = 12)
rater <- rep(seq_len(12), 50)
pi <- rnorm(50, 0, 1)
ri <- rnorm(12, 0, 0.6)
y <- 5 + pi[patient] + ri[rater] + rnorm(600, 0, 0.5)
save_once(data.frame(patient, rater, y), "ch16-grouping",
          "simulated-crossed.csv")
cat("Six seeded synthetic CSV files ready; existing files preserved.\n")

setwd("ch14-residuals")
suppressPackageStartupMessages(library(lme4))
d <- read.csv("simulated.csv")
d$section <- factor(d$section)
d$group <- factor(d$group)
# Fixed = population average; random = deviation.
# These are synthetic observations, not study data.
cat(nrow(d), "simulated rows\n")

print(head(d, 6))

H <- nlme::lme(y ~ section*group,
  random = ~1 | id, data = d, method = "REML",
  weights = nlme::varIdent(form = ~1 | section))
f <- data.frame(Estimate = nlme::fixef(H))
cat("H: population coefficients fitted\n")

print(f)

ratios <- coef(H$modelStruct$varStruct,
  unconstrained = FALSE, allCoef = TRUE)
v <- data.frame(section = names(ratios),
  residual_variance = sigma(H)^2 * ratios^2)
cat("Variance scale; covariance may be negative.\n")

print(v)

random_var <- as.numeric(nlme::VarCorr(H)[1, "Variance"])
print(random_var)

print(random_var)

setwd("..")

cd "ch14-residuals"
clear all
set more off
quietly import delimited ///
  "simulated.csv", clear case(preserve) asdouble
* Fixed = population average; random = deviation.
* These are synthetic observations, not study data.
quietly mkmat id section group y in 1/6, matrix(preview)
display _N " simulated rows"

matrix list preview

quietly mixed y i.section##i.group || id:, ///
  reml residuals(independent, by(section))
estimates store H
matrix f = (_b[_cons] \ ///
  _b[2.section] \ ///
  _b[3.section] \ ///
  _b[4.section] \ ///
  _b[1.group] \ ///
  _b[2.section#1.group] \ ///
  _b[3.section#1.group] \ ///
  _b[4.section#1.group])
matrix colnames f = Estimate
matrix rownames f = Intercept ///
  section2 ///
  section3 ///
  section4 ///
  group1 ///
  section2_group1 ///
  section3_group1 ///
  section4_group1
display "H: population coefficients fitted"

matrix list f

matrix v = (exp(2*_b[lnsig_e:_cons]) \ ///
  exp(2*(_b[lnsig_e:_cons]+_b[r_lns2ose:_cons])) \ ///
  exp(2*(_b[lnsig_e:_cons]+_b[r_lns3ose:_cons])) \ ///
  exp(2*(_b[lnsig_e:_cons]+_b[r_lns4ose:_cons])))
matrix colnames v = residual_variance
matrix rownames v = section1 section2 section3 section4
display "Variance scale; covariance may be negative."

matrix list v

scalar random_var = exp(2*_b[lns1_1_1:_cons])
display "Random-intercept variance = " %9.6f random_var

display random_var

cd ".."

Changing covariance assumptions can change SEs, confidence intervals, and P values, even when point estimates remain similar. Missingness and adjustment can also change estimates. Any sensitivity that changes a clinically relevant inference belongs in the report, including results near an arbitrary significance threshold.

• Key idea: Equal precision can follow from a symmetric model and design; covariance adequacy still needs independent assessment.

Documentation supporting the model and software details: nlme: section-specific variance.

Chapter 15. A Short Worked Model Specification

The question

For a continuous WSS outcome measured across anatomical sections, ask whether the mean difference between two groups varies by section. Keep section categorical to estimate each section's group contrast.

The random-intercept model

WSS ~ section * group + covariates + (1 | id)

The fixed interaction describes group contrasts by section. The shared person offset describes within-person dependence, with a separately specified residual structure.

The random-position-slope alternative

WSS ~ section * group + covariates + (1 + position | id)

This keeps the same mean question while allowing personal departures to change linearly along the numeric position axis defined in Chapter 13, subsection 5.

Use the R/Stata lesson and full code in Chapter 13, subsection 5 to run through these two structures on the same simulated four-section dataset, with position 0–3. Chapters 13 and 15 share that lesson so the identifiability explanation and model specification refer to the same example. The notation covariates is a generic placeholder; the executable teaching models contain no adjustment-variable list from an unpublished study.

Choosing between them

Choose a covariance structure that the design can identify and the profiles support. Examine fit stability and singularity; compare the models using the estimation and variance-boundary cautions in Chapter 13. A simpler model is a defensible choice only with its assumptions understood. The recorded comparison is a simulated demonstration, not a comparison from a patient study.

• Key idea: Keep the scientific mean question fixed while considering identifiable random-intercept and random-position-slope structures.

Chapter 16. Beyond Time: Clusters Such as Wards, Patients and Anatomical Segments

A random effect belongs to a grouping factor, not to time. Follow one continuous teaching example: a continuous patient outcome in a hospital service. Patients in the same ward may share care conditions even when each patient contributes only one observation. A model such as outcome ~ treatment + age + (1 | ward) gives each ward a shared departure from the overall mean.

Follow the grouping structure

If those patients remain in one ward and have repeated measurements, (1 | ward/patient) expands to (1 | ward) + (1 | ward:patient): ward offsets plus patient offsets within wards. Patient identifiers and the data structure must represent that nesting correctly; transfers between wards require a different specification. Repeated observations within one person may be indexed by time or location. Eyes within patients and anatomical segments within a vessel use the same idea; patients within hospitals are another clustering level.

Now suppose raters in this service each assess several patients and each patient can be assessed by several raters. Patient and rater are crossed rather than nested. (1 | patient) + (1 | rater) represents both shared sources of variation. It does not by itself add a patient-by-rater interaction or account for every possible remaining correlation.

A slope can also vary between wards

If treatment varies among patients within each ward and there are enough wards with treatment contrasts, outcome ~ treatment + age + (1 + treatment | ward) allows the treatment association to differ between wards, around its fixed mean association. This requires within-ward treatment variation and adequate replication. If each ward has only one treatment, a ward intercept and a ward-specific treatment slope cannot be separated by this model; use a model appropriate to the cluster-level treatment design.

Ward staffing is a cluster-level predictor: it varies between wards. Patient age is an individual-level predictor and can vary within a ward. More patients do not create more independent ward-level units. For an individual predictor, distinguish its within-ward association from the association of ward means when these represent different questions. Neither a random effect nor a fixed adjustment alone makes an observational treatment association causal.

# These synthetic distributions are invented for teaching, never study data.
# Run from this termcast directory. CSVs are shared unchanged by R and Stata.
options(digits = 17)
save_once <- function(d, example, filename = "simulated.csv") {
  dir.create(example, showWarnings = FALSE, recursive = TRUE)
  p <- file.path(example, filename)
  if (!file.exists(p)) write.csv(d, p, row.names = FALSE)
  invisible(p)
}
set.seed(1067001)
n <- 80; k <- 6
id <- rep(seq_len(n), each = k)
X <- rep(seq(-1, 1, length.out = k), n)
u0 <- rnorm(n, 0, 1)
u1 <- 0.2 * u0 + rnorm(n, 0, 0.75)
Y <- 4 + 1.5 * X + u0[id] + u1[id] * X + rnorm(n*k, 0, 0.6)
save_once(data.frame(id, X, Y), "ch07-four-models")

set.seed(1061101)
n <- 64; k <- 5
id <- rep(seq_len(n), each = k)
month <- rep(0:4, n)
treatment <- rep(rep(0:1, length.out = n), each = k)
u0 <- rnorm(n, 0, 4)
u1 <- 0.08 * u0 + rnorm(n, 0, 1)
BP <- 120 - month - 3*treatment - 0.8*month*treatment +
  u0[id] + u1[id]*month + rnorm(n*k, 0, 2)
save_once(data.frame(id, month, treatment, BP), "ch11-month")

set.seed(1061301)
n <- 72; k <- 4
id <- rep(seq_len(n), each = k)
position <- rep(0:3, n)
section <- position + 1
group <- rep(rep(0:1, length.out = n), each = k)
u0 <- rnorm(n, 0, 1)
u1 <- 0.15*u0 + rnorm(n, 0, 0.45)
y <- 5 + c(0, 0.4, 1, 1.2)[section] +
  group*c(0.4, 0.7, 0.2, 1)[section] +
  u0[id] + u1[id]*position + rnorm(n*k, 0, 0.5)
save_once(data.frame(id, position, section, group, y),
          "ch13-position")

set.seed(1061401)
n <- 96; k <- 4
id <- rep(seq_len(n), each = k)
section <- rep(1:4, n)
group <- rep(rep(0:1, length.out = n), each = k)
u0 <- rnorm(n, 0, 0.8)
y <- 5 + c(0, 0.4, 1, 1.2)[section] + 0.5*group +
  u0[id] + rnorm(n*k, 0, c(0.35, 0.6, 0.9, 1.2)[section])
save_once(data.frame(id, section, group, y), "ch14-residuals")

set.seed(1061601)
nw <- 32; np <- 8; nr <- 3
ward <- rep(seq_len(nw), each = np*nr)
patient <- rep(rep(seq_len(np), each = nr), nw)
treatment <- (patient - 1) %% 2
x <- rep(c(-1, 0, 1), nw*np)
wi <- rnorm(nw, 0, 1)
ws <- 0.2*wi + rnorm(nw, 0, 0.7)
pi <- rnorm(nw*np, 0, 0.6)
pid <- (ward - 1)*np + patient
y <- 5 + 0.7*x + 0.5*treatment + wi[ward] +
  ws[ward]*treatment + pi[pid] + rnorm(length(x), 0, 0.4)
save_once(data.frame(ward, patient, treatment, x, y),
          "ch16-grouping", "simulated-nested.csv")
set.seed(1061602)
patient <- rep(seq_len(50), each = 12)
rater <- rep(seq_len(12), 50)
pi <- rnorm(50, 0, 1)
ri <- rnorm(12, 0, 0.6)
y <- 5 + pi[patient] + ri[rater] + rnorm(600, 0, 0.5)
save_once(data.frame(patient, rater, y), "ch16-grouping",
          "simulated-crossed.csv")
cat("Six seeded synthetic CSV files ready; existing files preserved.\n")

setwd("ch16-grouping")
suppressPackageStartupMessages(library(lme4))
d <- read.csv("simulated-nested.csv")

# Fixed = population average; random = deviation.
# These are synthetic observations, not study data.
cat(nrow(d), "simulated rows\n")

print(head(d, 6))

d <- subset(d, x == 0)
cat(nrow(d), "rows: one row per patient\n")

print(nrow(d))

W0 <- lmer(y ~ treatment + (1 | ward), data = d)
f <- data.frame(Estimate = fixef(W0))
cat("W0: population coefficients fitted\n")

print(f)

v <- as.data.frame(VarCorr(W0))
v <- v[, c("grp", "var1", "var2", "vcov")]
v[is.na(v)] <- ""
cat("Variance scale; covariance may be negative.\n")

print(v)

W1 <- lmer(y ~ treatment +
  (1 + treatment | ward), data = d)
f <- data.frame(Estimate = fixef(W1))
cat("W1: population coefficients fitted\n")

print(f)

v <- as.data.frame(VarCorr(W1))
v <- v[, c("grp", "var1", "var2", "vcov")]
v[is.na(v)] <- ""
cat("Variance scale; covariance may be negative.\n")

print(v)

d <- read.csv("simulated-nested.csv")
cat(nrow(d), "repeated observations\n")

print(head(d, 6))

N <- lmer(y ~ x + (1 | ward/patient), data = d)
f <- data.frame(Estimate = fixef(N))
cat("N: population coefficients fitted\n")

print(f)

v <- as.data.frame(VarCorr(N))
v <- v[, c("grp", "var1", "var2", "vcov")]
v[is.na(v)] <- ""
cat("Variance scale; covariance may be negative.\n")

print(v)

d <- read.csv("simulated-crossed.csv")
cat(nrow(d), "simulated ratings\n")

print(head(d, 6))

K <- lmer(y ~ 1 + (1 | patient) + (1 | rater),
  data = d)
f <- data.frame(Estimate = fixef(K))
cat("K: population coefficients fitted\n")

print(f)

v <- as.data.frame(VarCorr(K))
v <- v[, c("grp", "var1", "var2", "vcov")]
v[is.na(v)] <- ""
cat("Variance scale; covariance may be negative.\n")

print(v)

setwd("..")

cd "ch16-grouping"
clear all
set more off
quietly import delimited ///
  "simulated-nested.csv", clear case(preserve) asdouble
* Fixed = population average; random = deviation.
* These are synthetic observations, not study data.
quietly mkmat ward patient treatment x y in 1/6, matrix(preview)
display _N " simulated rows"

matrix list preview

quietly keep if x == 0
display _N " rows: one row per patient"

display _N

quietly mixed y treatment || ward:, reml
estimates store W0
matrix f = (_b[_cons] \ ///
  _b[treatment])
matrix colnames f = Estimate
matrix rownames f = Intercept ///
  treatment
display "W0: population coefficients fitted"

matrix list f

quietly estat recovariance
matrix v = r(Cov2)
display "Variance scale; covariance may be negative."

matrix list v

scalar residual_var = exp(2*_b[lnsig_e:_cons])
display "Residual variance = " %9.6f residual_var

display residual_var

quietly mixed y treatment || ward: treatment, ///
  covariance(unstructured) reml
estimates store W1
matrix f = (_b[_cons] \ ///
  _b[treatment])
matrix colnames f = Estimate
matrix rownames f = Intercept ///
  treatment
display "W1: population coefficients fitted"

matrix list f

quietly estat recovariance
matrix v = r(Cov2)
display "Variance scale; covariance may be negative."

matrix list v

scalar residual_var = exp(2*_b[lnsig_e:_cons])
display "Residual variance = " %9.6f residual_var

display residual_var

quietly import delimited ///
  "simulated-nested.csv", clear asdouble
quietly mkmat ward patient x y in 1/6, matrix(preview)
display _N " repeated observations"

matrix list preview

quietly mixed y x || ward: || patient:, reml
estimates store N
matrix f = (_b[_cons] \ ///
  _b[x])
matrix colnames f = Estimate
matrix rownames f = Intercept ///
  x
display "N: population coefficients fitted"

matrix list f

quietly estat recovariance
matrix v = (r(Cov2)[1,1] \ ///
  r(Cov3)[1,1] \ exp(2*_b[lnsig_e:_cons]))
matrix colnames v = Variance
matrix rownames v = ward patient residual
display "Variance scale; covariance may be negative."

matrix list v

quietly import delimited ///
  "simulated-crossed.csv", clear asdouble
quietly mkmat patient rater y in 1/6, matrix(preview)
display _N " simulated ratings"

matrix list preview

quietly mixed y || _all: R.rater || patient:, reml
estimates store K
matrix f = (_b[_cons])
matrix colnames f = Estimate
matrix rownames f = Intercept
display "K: population coefficients fitted"

matrix list f

quietly estat recovariance
matrix v = (r(Cov2)[1,1] \ ///
  r(Cov3)[1,1] \ exp(2*_b[lnsig_e:_cons]))
matrix colnames v = Variance
matrix rownames v = rater patient residual
display "Variance scale; covariance may be negative."

matrix list v

cd ".."

What the ICC describes, and what few clusters change

In the simple Gaussian ward random-intercept model with common independent residual variance and no additional random effects, the conditional-on-covariates intraclass correlation is:

\[\mathrm{ICC}=\frac{\tau_{\mathrm{ward}}^2}{\tau_{\mathrm{ward}}^2+\sigma^2}.\]

It is the share of remaining outcome variance due to the ward and the correlation of two distinct patients in the same ward under those assumptions. Additional patient or rater effects, unequal residual variances, or random slopes change the variance decomposition and may make correlations depend on the pair and predictor values. One constant ICC is not a universal summary.

With few clusters, for example five wards, the ward variance and any slope heterogeneity are estimated poorly, and inference can be fragile. Fixed ward effects or a marginal approach with appropriate small-sample inference may better match the question and design. Cluster-robust standard errors are also unreliable with very few independent clusters; they are not an automatic remedy. Fixed ward effects cannot separately identify predictors that are constant within ward, such as ward-level treatment.

Mixed models handle repeated measurements and clusters. The grouping structure, together with the scientific question and available variation, determines the random effects; the word “time” does not.

• Key idea: Identify which observations share a ward, patient, vessel or rater, then specify random effects and uncertainty that those independent groups can support.

Related in the wiki: [[describe-data-structure-models]] [[mixed-model-random-slopes]] [[correlated-data-model-guide]] [[mixed-model-random-intercept-shared-offset]]