Reading Mixed-Model Software Output for Crossover Data

Why Software Output Needs Interpretation

A mixed-effects model can produce a long table of coefficients, variance estimates, tests, and warnings. The numbers are useful only when connected to the fitted model and the question being asked. In a two-period crossover analysis, the central question is usually the estimated difference between formulations on the log exposure scale. A coefficient labelled “Treatment” is not automatically that difference: its meaning depends on the model’s coding and reference levels.

This tutorial considers a balanced, complete two-sequence crossover, with subjects assigned to sequence AB or BA and each subject observed in both periods. The response $Y$ is the log-transformed exposure measure. The aim is to read typical mixed-model output, not to choose among model structures or make an equivalence decision.

The Model Behind the Output

Model: A common two-period model includes fixed effects for sequence, period, and treatment, a random subject effect nested within sequence, and a residual error: $$ Y_{ijk}=\mu+\text{Sequence}_i+\text{Subject}_{j(i)}+\text{Period}_k+\text{Treatment}_{\ell(ijk)}+\varepsilon_{ijk}. $$ Here $i$ indexes sequence, $j$ indexes subjects within sequence, $k$ indexes period, and $\ell(ijk)$ identifies the formulation received. The random subject effect has mean zero and variance $\sigma_S^2$; the residual errors have mean zero and variance $\sigma_e^2$. In the usual model, subject effects and residual errors are independent, and residual errors are independent across observations. The model also assumes that the fixed-effects structure adequately describes mean log exposure and that the residual variance is common across formulations.

The output’s coefficient table describes the fixed-effects part of this model. Its variance-components table describes estimated variability associated with the random subject effect and residual error. A test or interval for treatment requires the estimated treatment contrast, its standard error, and an appropriate denominator degrees of freedom. These are related outputs, but they answer different questions.

Software often uses reference-level coding. If reference formulation R is compared with formulation T, the treatment coefficient may represent T minus R, R minus T, or a contrast coded another way. Check the stated reference levels and contrast definitions before interpreting the sign. On the log scale, a treatment contrast $\widehat{\theta}$ is a difference in mean log exposures; exponentiating it gives a ratio-scale estimate.

Worked Example: Reading a Treatment Coefficient

Suppose hypothetical software output uses R as the reference formulation and reports the coefficient for T as $-0.0350$, with standard error $0.0280$, denominator degrees of freedom $28$, and a two-sided $p$-value of $0.221$. The model’s contrast is T minus R.

The estimated log-scale contrast is $\widehat{\theta}=-0.0350$. Its test statistic for a zero contrast is $$ t=\frac{\widehat{\theta}}{SE}=\frac{-0.0350}{0.0280}=-1.25. $$

Exponentiating gives the estimated geometric mean ratio: $$ \exp(\widehat{\theta})=\exp(-0.0350)=0.9656. $$ Thus, the fitted geometric mean exposure for T is estimated to be 0.9656 times that for R, or about 3.44% lower, for this hypothetical analysis. The $p$-value of 0.221 is not strong evidence against a zero log-scale contrast under the specified test. It does not establish that the formulations are equivalent.

For illustration, a 95% $t$-interval uses $t_{0.975,28}\approx2.048$: $$ -0.0350\pm2.048(0.0280)=[-0.0923,\ 0.0223]. $$ Exponentiating the endpoints gives an interval of approximately $[0.9118,\ 1.0226]$ for the geometric mean ratio. In context, the interval procedure has 95% coverage under repeated sampling when its model and assumptions hold; it is not a 95% probability statement about this fixed ratio.

Variance Components and What They Describe

A variance-components table may report a variance estimate for subjects and another for residual error. These are on the squared log-outcome scale. Taking square roots expresses the corresponding standard deviations in log units. A larger subject variance indicates more variation among subjects after accounting for the fixed effects; a larger residual variance indicates more within-subject unexplained variation under the fitted model.

The subject variance is not the standard error of the treatment estimate. Nor should the residual variance be read as the variability of the original, untransformed exposure in its original units. Both refer to the fitted log-scale model. For a random-intercept model, the ratio $\sigma_S^2/(\sigma_S^2+\sigma_e^2)$ describes the model-implied correlation between two observations from the same subject, conditional on the fixed effects and under the stated covariance assumptions. That correlation is not itself a measure of treatment effect.

Worked Example: Translating Variance Estimates

Consider hypothetical output with an estimated subject variance of $0.0121$ and residual variance of $0.0180$. Both are variances in squared log units.

The estimated subject standard deviation is $$ \sqrt{0.0121}=0.1100 $$ log units, and the estimated residual standard deviation is $$ \sqrt{0.0180}=0.1342 $$ log units. These values describe different sources of variation in the fitted model; they are not standard errors for the treatment contrast.

The implied within-subject correlation for this random-intercept model is $$ \frac{0.0121}{0.0121+0.0180} =\frac{0.0121}{0.0301} =0.4020. $$ Under this model, two observations from the same subject have an estimated correlation of about 0.402 after accounting for the fixed effects. This is a model-based summary, not a direct correlation calculated without adjustment.

Why Denominator Degrees of Freedom Can Be Fractional

A mixed-model test statistic for a fixed effect is commonly reported in the form $t=\widehat{\theta}/SE$, but its reference distribution depends on how uncertainty in the variance components is handled. Software may approximate the denominator degrees of freedom using methods such as Satterthwaite or Kenward–Roger. These methods account, in different ways, for the fact that variance components are estimated rather than known. The resulting degrees of freedom need not be an integer.

The denominator degrees of freedom determine the reference $t$-distribution used for a test or confidence interval. They are not the number of subjects, the number of observations, or a count of residuals. In a simple balanced crossover, familiar hand calculations can provide a useful check on the treatment estimate and its uncertainty. Mixed-model procedures can nevertheless report different degrees of freedom depending on the method and model. Read the method named in the output rather than assuming all software uses the same calculation.

Worked Example: Interpreting Approximate Degrees of Freedom

Suppose hypothetical output reports a treatment contrast of $-0.0420$ log units, $SE=0.0310$, denominator $df=18.6$, and a two-sided $p$-value of about 0.192. The software identifies the degrees-of-freedom method as Satterthwaite.

The test statistic is $$ t=\frac{-0.0420}{0.0310}=-1.355. $$ The fractional $df=18.6$ is an approximation used to choose the reference $t$-distribution; it is not an error in the subject count. Using $t_{0.975,18.6}\approx2.096$, the illustrative 95% interval is $$ -0.0420\pm2.096(0.0310) =-0.0420\pm0.0650 =[-0.1070,\ 0.0230]. $$

Because this interval includes zero and the two-sided $p$-value is 0.192, the output does not provide strong evidence of a nonzero log-scale treatment contrast at the 0.05 level. The result concerns evidence against a zero contrast; it does not settle a separate equivalence question.

A Short Justification for the Treatment Contrast

The treatment coefficient can be understood by forming each subject’s within-subject difference, T minus R. Under the stated model, the subject random intercept cancels from that difference. Let $\pi$ be the period 2 effect relative to period 1 and $\theta$ the T-minus-R treatment effect. In sequence AB, R is given in period 1 and T in period 2, so the expected within-subject difference is $\theta+\pi$. In sequence BA, T is given in period 1 and R in period 2, so it is $\theta-\pi$. Therefore, with equal numbers in the two sequences, the average of the two sequence-specific mean differences targets the treatment effect:

\frac{(\theta+\pi)+(\theta-\pi)}{2}=\theta.

This cancellation explains why the treatment estimate is not simply the mean of the AB subjects’ differences: that mean also includes the period effect. Equal sequence weighting cancels the period contribution in this balanced setting. The argument relies on the stated fixed-effects model and complete observations in both periods; it should not be carried over automatically to other designs or missing-data patterns.

Convergence Messages and Diagnostics

Convergence is a numerical statement: the fitting algorithm has met its specified stopping criteria for the optimization problem. It is necessary for trusting the reported estimates, but it does not prove that the model is scientifically appropriate, that its assumptions hold, or that the treatment result is reliable. Conversely, a warning should not be ignored merely because the treatment coefficient looks plausible.

Worked Example: Responding to a Boundary Warning

Suppose hypothetical output reports a subject variance estimate of $0.00004$, a residual variance estimate of $0.0144$, and a message that the fitted random-effects covariance is singular or near a boundary. The estimated subject standard deviation is $$ \sqrt{0.00004}=0.00632 $$ log units, while the residual standard deviation is $$ \sqrt{0.0144}=0.1200 $$ log units. The fitted subject variance is very small relative to the residual variance, which is consistent with a boundary warning.

The warning does not, by itself, show that the treatment estimate is biased or that the model should be changed. Check whether the optimizer reports convergence, whether the gradient and Hessian diagnostics indicate a numerical problem, and whether estimates are stable under reasonable starting values or fitting settings. Then consider whether the near-zero component is plausible for the design and model. Do not treat a boundary estimate as proof of no subject variability, and do not remove a model component solely to obtain a warning-free table.

Reading Output Without Overclaiming

A useful reading order is to confirm the fitted model and coding, identify the treatment contrast, inspect its standard error and denominator degrees of freedom, then examine variance components and convergence messages. This keeps a common confusion in check: the fixed-effect table concerns estimated mean differences, while the variance table describes modeled variability. Neither table substitutes for the other.

The main alternatives are not interchangeable output formats but different inferential calculations and model specifications. In a complete balanced crossover, the hand-calculated ANOVA provides a valuable benchmark for the treatment estimate and standard error. A mixed-model procedure can express the same fixed treatment contrast while using an explicit random-effects covariance structure and an approximate degrees-of-freedom method. If those analyses differ, compare the model terms, variance assumptions, coding, contrast, and degrees-of-freedom method before concluding that one output is wrong.

Common failure signs include an unexpected sign caused by reversed reference levels, a reported test for an overall factor when the desired result is a specific contrast, variance estimates at or near boundaries, unstable estimates, and optimizer or Hessian warnings. Residual plots and other model diagnostics assess assumptions that a convergence flag cannot check. A converged fit can still be poorly supported by the data or incompatible with the intended analysis.

Common pitfall: “The software says the model converged, so the treatment result is valid.” Convergence only indicates that the numerical algorithm met its stopping rule. The correction is to verify the model, contrast coding, uncertainty calculation, assumptions, and diagnostic messages separately.

Check Your Understanding

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

If the output coefficient is T minus R and equals $0.020$ log units, what is the estimated geometric mean ratio T/R?
What does a subject variance estimate describe, and why is it not the treatment standard error?
Why might a mixed-model output report denominator degrees of freedom of 16.7?
What does a convergence message establish, and what does it not establish?

Solutions. (1) The ratio estimate is $\exp(0.020)=1.020$, so the fitted geometric mean for T is about 1.020 times that for R. (2) It estimates between-subject variance after accounting for fixed effects, on the squared log scale; treatment uncertainty also depends on the design, residual variation, and the contrast’s standard error calculation. (3) An approximation such as Satterthwaite or Kenward–Roger can produce fractional degrees of freedom when variance components are estimated. (4) Convergence indicates that the optimizer met its numerical stopping criteria. It does not verify model assumptions, appropriateness, or the scientific validity of the treatment conclusion.