Estimating Crossover Model Effects by REML

Why Estimate Variance Components by REML?

In a crossover analysis, observations from the same subject are related: a subject-level tendency to have higher or lower exposure affects both periods. A mixed-effects model represents that tendency with a random subject effect, while fixed effects describe features such as treatment, period, and sequence. Estimating the variance of the random effect matters because it describes between-subject variability and contributes to the covariance of the observations used to estimate fixed effects.

Restricted maximum likelihood, or REML, estimates variance components after accounting for the fixed effects. This distinction matters in finite samples: ordinary maximum likelihood estimates variance parameters using a likelihood that also estimates the fixed-effect means, which can produce downward-biased variance estimates. REML instead bases estimation on information that does not depend on the unknown fixed-effect values.

The Random-Intercept Crossover Model

Consider a complete two-period, two-treatment crossover with log-transformed pharmacokinetic outcomes. Let \(y_{ij}\) be the log outcome for subject \(j\) in observation \(i\), where \(i\) indexes that subject’s two periods. A convenient model is

y_{ij} = \mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta} + b_j + \varepsilon_{ij}.

The vector \(\boldsymbol{\beta}\) contains fixed effects, including the intercept and the estimable sequence, period, and treatment effects. The design vector \(\mathbf{x}_{ij}\) specifies which fixed effects apply to observation \(ij\). The random intercept \(b_j\) represents the subject’s persistent departure from the population level. The residual \(\varepsilon_{ij}\) represents remaining within-subject variation.

Model assumptions: The subject effects are independent with \(b_j\sim N(0,\sigma_s^2)\). The residuals are independent with \(\varepsilon_{ij}\sim N(0,\sigma_e^2)\), and are independent of the subject effects. The fixed-effects design matrix \(X\) has full column rank. The variance components satisfy \(\sigma_s^2\geq 0\) and \(\sigma_e^2>0\). In this model, the two observations from one subject have covariance \(\sigma_s^2\), while observations from different subjects have covariance zero.

In matrix form, \(y=X\boldsymbol{\beta}+Zb+\varepsilon\), so the marginal distribution is \(y\sim N(X\boldsymbol{\beta},V)\), with \(V=ZGZ^{\mathsf T}+R\). For the random-intercept model, \(G=\sigma_s^2 I\), and \(R=\sigma_e^2 I\). The model therefore imposes a specific covariance pattern; it is not a distribution-free way to handle repeated measurements.

What REML Maximizes

If the variance parameters are collected in \(\theta\), the generalized least-squares estimate of the fixed effects for a given \(V(\theta)\) is

\hat{\boldsymbol{\beta}}(\theta) = (X^{\mathsf T}V^{-1}X)^{-1}X^{\mathsf T}V^{-1}y.

REML does not simply substitute this estimate into the ordinary likelihood and treat the result as a variance likelihood. Instead, it constructs likelihood information from linear combinations of the data that eliminate the fixed effects. The resulting restricted log likelihood, up to a constant that does not depend on \(\theta\), is

-2\ell_R(\theta) = \log|V| + \log|X^{\mathsf T}V^{-1}X| + y^{\mathsf T}Py,

where

P = V^{-1} - V^{-1}X(X^{\mathsf T}V^{-1}X)^{-1}X^{\mathsf T}V^{-1}.

The REML estimates are the nonnegative variance-component values that minimize this expression. The second determinant term is an adjustment for estimating the fixed effects. Omitting it does not give the restricted likelihood.

Justification: Choose a matrix \(K\) whose columns span the directions orthogonal to the columns of \(X\), so \(K^{\mathsf T}X=0\). Then \(K^{\mathsf T}y=K^{\mathsf T}X\boldsymbol{\beta}+K^{\mathsf T}e=K^{\mathsf T}e\), where \(e\) has covariance \(V\). Thus these contrasts have distribution \(N(0,K^{\mathsf T}VK)\), independent of the unknown fixed-effect values. Their likelihood provides the restricted likelihood. The matrix expression above is an equivalent form of that likelihood, with terms that account for the covariance scale and the fixed-effects design.

A Balanced-Design Calculation

In a balanced, complete two-period crossover, the REML calculation has a useful transparent form. For each subject define the average and half-difference of their two log observations:

m_j=\frac{y_{j1}+y_{j2}}{2}, \qquad d_j=\frac{y_{j2}-y_{j1}}{2}.

Under the stated model, the subject effect contributes to \(m_j\) but cancels from \(d_j\). The residual contributions to the two quantities are independent, each with variance \(\sigma_e^2/2\). Therefore,

\operatorname{Var}(m_j)=\tau^2=\sigma_s^2+\frac{\sigma_e^2}{2}, \qquad \operatorname{Var}(d_j)=q=\frac{\sigma_e^2}{2}.

In the standard two-sequence design, the fixed effects allow a separate mean for the \(m_j\) values in each sequence and a separate mean for the \(d_j\) values in each sequence. With \(n\) subjects in total, removing those two means leaves \(n-2\) residual degrees of freedom for each set of transformed observations. Let \(SS_m\) and \(SS_d\) be the within-sequence sums of squares for the averages and half-differences. If the unconstrained estimates satisfy \(\hat{\tau}^2\geq\hat{q}\), then

\hat{\tau}^2=\frac{SS_m}{n-2}, \qquad \hat{q}=\frac{SS_d}{n-2}, \qquad \hat{\sigma}_e^2=2\hat{q}, \qquad \hat{\sigma}_s^2=\hat{\tau}^2-\hat{q}.

These formulas follow from the restricted likelihoods of independent normal residual contrasts: for a normal sample with an estimated mean, the restricted-likelihood variance estimate is its residual sum of squares divided by its residual degrees of freedom. The last two equations then convert the estimates of \(\tau^2\) and \(q\) back to the model’s residual and subject variances. If the difference estimate is negative, the nonnegative constraint on \(\sigma_s^2\) must be enforced by maximizing the likelihood at the boundary; simply reporting a negative variance is not permissible.

Worked Example: Estimating Both Variance Components

Suppose a hypothetical balanced crossover has four subjects in each sequence, so \(n=8\). After calculating \(m_j\) and \(d_j\), subtract the mean within each sequence. Assume the resulting centered \(m_j\) values are \(-0.35,-0.15,0.15,0.35\) in each sequence, and the centered \(d_j\) values are \(-0.10,-0.05,0.05,0.10\) in each sequence. These are illustrative log-scale data, not results from an actual study.

For the averages, the sum of squares in one sequence is \(0.35^2+0.15^2+0.15^2+0.35^2=0.29\). Across the two sequences, \(SS_m=0.58\). Similarly, each sequence contributes \(0.10^2+0.05^2+0.05^2+0.10^2=0.025\) to the half-difference sum of squares, so \(SS_d=0.05\). The residual degrees of freedom for each are \(n-2=6\).

The estimates are

\hat{\tau}^2=\frac{0.58}{6}=0.09667, \qquad \hat{q}=\frac{0.05}{6}=0.008333.

Since \(\hat{\tau}^2>\hat{q}\), the subject variance estimate is positive. Converting to the model components gives

\hat{\sigma}_e^2=2(0.008333)=0.01667, \qquad \hat{\sigma}_s^2=0.09667-0.008333=0.08833.

Both estimates are in squared log-outcome units. The fitted subject-intercept standard deviation is \(\sqrt{0.08833}=0.2972\) log units, and the residual standard deviation is \(\sqrt{0.01667}=0.1291\) log units. These summarize different sources of variation: persistent differences between subjects and remaining variation between a subject’s observations.

Boundary Estimates and Interpretation

Worked Example: A Subject Variance at the Boundary

Consider another hypothetical set of transformed data with \(n=8\), \(SS_m=0.12\), and \(SS_d=0.30\). The unconstrained estimates are \(\hat{\tau}^2=0.12/6=0.0200\) and \(\hat{q}=0.30/6=0.0500\). Their difference would imply \(\hat{\sigma}_s^2=-0.0300\), which is not a valid variance.

The constrained model requires \(\tau^2\geq q\). Its maximum in this case occurs at the boundary \(\tau^2=q=v\), which corresponds to \(\sigma_s^2=0\). The restricted log-likelihood contribution that depends on \(v\) is equivalent to minimizing

12\log(v)+\frac{SS_m+SS_d}{v}.

Differentiating with respect to \(v\) and setting the derivative to zero gives \(12/v-0.42/v^2=0\), hence \(v=0.42/12=0.0350\). Therefore \(\hat{q}=0.0350\), \(\hat{\sigma}_e^2=2(0.0350)=0.0700\), and \(\hat{\sigma}_s^2=0\). The data and model do not support a positive subject-intercept variance estimate under this fit; zero is a boundary estimate, not proof that subjects are identical.

Worked Example: Reading Fitted Variances on the Log Scale

Suppose REML software reports hypothetical estimates \(\hat{\sigma}_s^2=0.0400\) and \(\hat{\sigma}_e^2=0.0900\) for a natural-log AUC model. The corresponding standard deviations are \(\sqrt{0.0400}=0.200\) log units for subject effects and \(\sqrt{0.0900}=0.300\) log units for residuals.

Under the normal random-effect assumption, a subject effect one standard deviation above zero is \(0.200\) log units above the population level. On the original multiplicative scale, \(\exp(0.200)=1.221\), so this is a factor of about 1.22 in the subject-level scale relative to the population level. It is not a 22.1% estimate of a treatment effect. Likewise, the residual standard deviation describes unexplained observation-level variation, not the uncertainty of the treatment estimate. The fixed treatment estimate and its uncertainty must be interpreted separately.

Diagnostics, Alternatives, and Limitations

Before interpreting REML estimates, check whether the fitted model is a reasonable description of the data. Confirm that the analysis uses the intended log endpoint and that sequence, period, and treatment are represented consistently in the fixed-effects design. Inspect residuals for markedly unusual patterns or observations, and check whether the fitted model’s assumptions of normality and a common residual variance are plausible. A variance estimate at zero, an estimate very close to zero, or substantial sensitivity to a small number of subjects deserves attention. These findings do not by themselves establish a model failure, but they can signal weak information or a mismatch between the assumed and observed variability.

For a balanced complete two-period design under this random-intercept model, the transformed-sums-of-squares calculation gives a direct way to understand the REML estimates. In an unbalanced design, the simple sums-of-squares decomposition generally does not apply; REML instead optimizes the matrix likelihood using the actual design. The ordinary ANOVA approach can be convenient when its balance and error-term conditions hold, while REML handles unequal information through the specified covariance model. Maximum likelihood is another likelihood-based option, but unlike REML it does not adjust its variance likelihood for the number of fixed effects estimated. These approaches may therefore produce different variance estimates, especially with limited data.

REML does not repair an incorrect model. Its estimates depend on the fixed-effects design, the random-effects structure, and the residual assumptions. With few subjects, variance estimates can be imprecise or lie on a boundary. A random-intercept model also makes a particular covariance assumption; REML alone does not determine whether that assumption is scientifically adequate. In addition, restricted likelihoods from models with different fixed-effects designs are not directly comparable as if they were calculated from the same restricted data space.

Common Mistakes

Confusing the two variance components. The subject variance describes persistent between-subject variation; the residual variance describes remaining observation-level variation under the model.
Interpreting a variance as a standard deviation. A variance is in squared log units. Take its square root to obtain a standard deviation in log units.
Allowing a negative variance estimate. Variance components are constrained to be nonnegative. A negative unconstrained difference indicates that the constrained likelihood must be optimized at or near a boundary.
Calling a fitted variance a treatment effect. Variance components describe variability, not the estimated treatment contrast or its confidence interval.
Assuming REML makes the model assumption-free. REML adjusts for fixed-effect estimation, but it still relies on the chosen covariance structure and distributional assumptions.

Key takeaway: REML estimates variance components from information remaining after fixed effects are accounted for. In a balanced complete two-period crossover, transformed averages and half-differences reveal how the subject and residual variances contribute; in every design, interpret the estimates in the units and assumptions of the fitted model.

Check Your Understanding

Answer each question, then compare your reasoning with the brief solutions.

In the random-intercept model, what is the covariance between two observations from the same subject?
With \(n=10\), \(SS_m=0.64\), and \(SS_d=0.16\), calculate \(\hat{\sigma}_s^2\) and \(\hat{\sigma}_e^2\), assuming the unconstrained estimates are feasible.
Why does REML include an adjustment for the fixed-effects design?
A fitted residual variance is \(0.0625\) squared log units. What is the residual standard deviation, and in what units?
What should be done if the unconstrained estimate of \(\sigma_s^2\) is negative?

Solutions. (1) The covariance is \(\sigma_s^2\), because the two observations share the same random intercept and have independent residuals. (2) Here \(n-2=8\), so \(\hat{\tau}^2=0.64/8=0.0800\) and \(\hat{q}=0.16/8=0.0200\). Thus \(\hat{\sigma}_e^2=2(0.0200)=0.0400\) and \(\hat{\sigma}_s^2=0.0800-0.0200=0.0600\), in squared log units. (3) REML uses contrasts that eliminate the fixed-effect means, so its likelihood must account for the information consumed in estimating those effects. (4) The standard deviation is \(\sqrt{0.0625}=0.250\) log units. (5) Enforce the nonnegative constraint and maximize the restricted likelihood at the boundary as appropriate; do not report the negative value as a variance.