Common Mistakes in Factorial Mixed Models

Why Small Modeling Choices Matter

A factorial mixed model combines systematic effects of experimental factors with variation associated with groups such as hospitals, clinics, schools, or sites. Its conclusions depend on more than the formula’s list of terms. A model can be computationally well behaved and still answer the wrong question if sampled groups are treated as the entire target, plausible treatment variation is suppressed, or a coefficient is interpreted without regard to its coding.

This tutorial focuses on four diagnostic mistakes. The goal is not to include every possible random term or to prefer one fitting method in every setting. It is to make the model’s assumptions and the intended interpretation agree.

A Model to Use for Diagnosis

Suppose a response $Y_{gijr}$ is recorded in group $g$, at levels $i$ and $j$ of two factorial factors, with replicate $r$. For illustration, code factor A as $x_{gi}\in\{0,1\}$ and factor B as $z_{gj}\in\{0,1\}$. A model with fixed factorial effects and a group-specific intercept and A slope is

Y_{gijr}=\beta_0+\beta_Ax_{gi}+\beta_Bz_{gj}+\beta_{AB}x_{gi}z_{gj} +u_{0g}+u_{Ag}x_{gi}+\varepsilon_{gijr}.

Here the $\beta$ terms describe the fixed-effects mean structure. The group effects $u_{0g}$ and $u_{Ag}$ allow groups to differ in baseline response and in their A effect. A standard Gaussian version assumes independent groups, group effects with mean zero and a specified covariance matrix, and residual errors with mean zero that are independent of the group effects and of one another, with a common residual variance. Normality of the group effects and residuals is an additional assumption for the usual Gaussian likelihood and small-sample procedures. The intended sampling scheme and scientific question must support treating the groups as draws from a wider population; a random-effects formula cannot establish that on its own.

These assumptions are a starting model, not a guarantee that every term belongs in every analysis. In particular, a random slope is meaningful only when the corresponding factor varies across groups in a way that can inform group-to-group differences and when that variation is part of the scientific target.

Mistake 1: Treating Sampled Groups as Fixed by Default

A group term is not automatically fixed just because each group has a name. The key question is the target of inference. If the study concerns only a prespecified set of particular hospitals, hospital-specific fixed effects may be appropriate. If hospitals were sampled to learn about variation across a broader population of hospitals, a random-effects description may better reflect that target. These choices give different interpretations: a fixed effect describes the groups included in the model, while a random-effects distribution describes group variation in a population.

Worked Example: What Does the Group Term Represent?

Suppose a hypothetical study samples four clinics and estimates their baseline mean response, in points, as 12, 14, 11, and 15. The intended question is about clinics like these in a wider clinic population, not only the four named clinics.

The average of these four estimates is

\bar{x}=\frac{12+14+11+15}{4}=\frac{52}{4}=13\text{ points}.

Their sample standard deviation is

s=\sqrt{\frac{(12-13)^2+(14-13)^2+(11-13)^2+(15-13)^2}{4-1}} =\sqrt{\frac{1+1+4+4}{3}} =\sqrt{\frac{10}{3}} \approx 1.83\text{ points}.

These calculations describe the four displayed estimates; by themselves, they are not a reliable estimate of a population distribution from only four clinics. But they clarify the modeling target. A random clinic effect represents clinic-to-clinic variation around a population mean, whereas four fixed clinic indicators describe differences among these particular clinics. If the scientific question were restricted to the four named clinics, fixed indicators could be sensible. Calling clinics random simply because they are clinics would not be enough; the sampling target and study design must justify generalization.

Mistake 2: Omitting a Scientifically Justified Random Slope

A random-intercept-only model allows groups to start at different levels but assumes that the modeled treatment effect is the same in every group. That assumption can be too restrictive when the treatment varies within groups and the scientific question concerns how its effect varies across the group population. A random slope represents that variation. It is not a rule to add a slope for every factor: the design must provide information about within-group treatment differences, and the slope must have a scientific interpretation.

Derivation: For a treatment coded $x\in\{0,1\}$, the group-specific random contribution in the model above is $u_{0g}+xu_{Ag}$. The difference in that contribution between treatment levels 1 and 0 is $$ (u_{0g}+u_{Ag})-(u_{0g}+0u_{Ag})=u_{Ag}. $$ Thus the random intercept cancels from the group-specific treatment contrast, while the random slope remains. If $\operatorname{Var}(u_{Ag})=\tau_{AA}$, then the between-group variance in this contrast is $\tau_{AA}$. An intercept-only model sets this source of treatment-effect variation to zero.

Worked Example: Checking Whether a Shared Treatment Effect Is Plausible

Imagine four hypothetical sites where the observed treatment-minus-control differences, in outcome units, are 2, 0, 4, and 2. Their average difference is

\bar{x}=\frac{2+0+4+2}{4}=2\text{ units}.

The sample variance of these four differences is

s^2=\frac{(2-2)^2+(0-2)^2+(4-2)^2+(2-2)^2}{4-1} =\frac{0+4+4+0}{3} =\frac{8}{3}\approx 2.67\text{ units}^2.

The sites do not show identical observed differences: one is 0 units and another is 4 units. This pattern makes it important to consider whether the scientific model should allow site-to-site treatment variation. It does not prove that the population random-slope variance is nonzero: estimates from four sites are imprecise, and observed differences also contain sampling error. Nor does it mean a random slope is always estimable or warranted. The diagnostic question is whether a shared-effect assumption is defensible for the study’s target, given how treatment varies across sites and the available data.

A useful covariance check makes the assumption concrete. Under the stated model, two distinct observations in the same group at factor settings $x$ and $x'$ share random effects. If residual errors are independent, their covariance from those group effects is

\operatorname{Cov}(Y_g(x),Y_g(x'))= \tau_{00}+(x+x')\tau_{0A}+xx'\tau_{AA},

where $\tau_{00}=\operatorname{Var}(u_{0g})$, $\tau_{AA}=\operatorname{Var}(u_{Ag})$, and $\tau_{0A}=\operatorname{Cov}(u_{0g},u_{Ag})$. Omitting a scientifically justified slope removes the $\tau_{AA}$ contribution and the intercept–slope covariance term. This can affect uncertainty and the dependence structure, not just the wording of the fitted coefficients.

Mistake 3: Treating a Coded Coefficient as a Universal Effect

In factorial models, a coefficient describes a particular contrast implied by the coding and the other factors’ reference values. With treatment coding for both factors, $\beta_A$ in the model above is the A difference when B is at its reference level $z=0$. If the interaction is present, the A difference can change with B. The fitted cell means and explicitly defined contrasts are generally more informative than an isolated coefficient.

Worked Example: The A Effect Depends on B and on Coding

Suppose four hypothetical cell means, in outcome units, are $10$ for A0/B0, $14$ for A1/B0, $12$ for A0/B1, and $20$ for A1/B1. With treatment coding, the fitted values give

\beta_0=10,\qquad \beta_A=14-10=4,\qquad \beta_B=12-10=2,

and the interaction is the difference between the two A effects:

\beta_{AB}=(20-12)-(14-10)=8-4=4.

Therefore the A effect is 4 units when B is at level 0, but $4+4=8$ units when B is at level 1. Calling $\beta_A=4$ “the A effect in every setting” would be incorrect. If the two B levels are equally weighted, the average of the two A differences is $(4+8)/2=6$ units.

Now code each factor with values $-1$ and $+1$, using a model with terms for the intercept, A, B, and their interaction. The coefficient of A is one quarter of the sum of the cell values multiplied by the A code:

\beta_A^{*}=\frac{-10+14-12+20}{4}=\frac{12}{4}=3.

Under this coding, changing A from $-1$ to $+1$, averaged equally over B, corresponds to a difference of $2\beta_A^{*}=6$ units. The coefficient is 3, not 6, because the coded change in A is two units. Both codings describe the same cell means; their coefficients have different meanings. State the coding, reference values, and contrast before interpreting a coefficient.

Mistake 4: Comparing REML Fits With Different Fixed Effects

Restricted maximum likelihood (REML) uses linear combinations of the response that remove the fixed-effects mean structure when estimating covariance parameters. The combinations available depend on the fixed-effects design matrix. Consequently, restricted likelihood values from models with different fixed-effects structures are not directly comparable as if they were likelihoods calculated from the same response representation. A smaller or larger reported REML value does not, by itself, select between those fixed-effects models.

For a comparison of nested fixed-effects structures, use maximum likelihood (ML), provided the models use the same observations and response and the nesting and regularity conditions for the likelihood-ratio comparison are satisfied. For a one-parameter addition under the usual large-sample reference, the likelihood-ratio statistic is compared with a chi-square distribution with one degree of freedom. This is not a claim that every model comparison is a simple chi-square test; the comparison and its assumptions must be checked.

Worked Example: Why Two REML Numbers Do Not Decide a Fixed Effect

Suppose a hypothetical analysis uses the same 6 observations to compare a model with an intercept only against a model with an intercept and one treatment coefficient. Their fixed-effects design matrices have ranks 1 and 2, respectively. The corresponding residual contrast spaces have dimensions $6-1=5$ and $6-2=4$. If the reported restricted log likelihoods are $-8.1$ and $-6.9$, their numerical difference is

$$ -6.9-(-8.1)=1.2. $$

That subtraction is not a valid likelihood-ratio statistic for the fixed-effect comparison: the two restricted likelihoods are based on different residual contrast spaces. The apparent improvement of 1.2 cannot be interpreted as evidence for the treatment term.

Suppose instead the ML log likelihoods for these same nested models are $-10.2$ and $-8.0$. The ML likelihood-ratio statistic is

2\{\ell_{\text{larger}}-\ell_{\text{smaller}}\} =2\{-8.0-(-10.2)\} =2(2.2)=4.4.

For a one-parameter addition, the usual large-sample reference is $\chi^2_1$; a statistic of 4.4 gives an approximate $p$-value of 0.036. Under the stated comparison conditions, this is evidence against the smaller fixed-effects model at the 0.05 level. It is not a 3.6% probability that the null model is true, and the result does not remove the need to interpret the treatment effect in context.

A Practical Diagnostic Sequence

Before interpreting a factorial mixed-model result, check the modeling choices in a consistent order. This sequence helps expose mistakes without treating any one formula as universally correct.

Name the target.
Decide whether the conclusions concern the particular groups observed or a broader population represented by sampled groups.

Check group-level variation.
Ask whether the treatment or factor varies within groups and whether the scientific question includes variation in its effect across groups.

Translate coefficients into contrasts.
Write down the factor coding and reference levels, then use fitted means or explicit contrasts to state the effect being described.

Match the fitting method to the comparison.
Do not compare REML likelihoods across different fixed-effects structures; use an appropriate ML comparison when its conditions hold.

Common Mistakes

“The groups have names, so they must be fixed.” Names do not determine the model. The sampling target determines whether the analysis is about those particular groups or a broader group population.
“A random intercept handles all group differences.” It handles baseline differences, but it assumes a shared treatment slope unless the model includes slope variation.
“The A coefficient is the A effect.” With interactions, it may be the A effect only at a reference level of another factor. Coding can also change the coefficient’s scale and interpretation.
“The REML value improved, so the added fixed effect is supported.” REML values from different fixed-effects spaces are not directly comparable for that decision.

These are not merely reporting details. They determine which population is described, what treatment variation is allowed, which contrast a number represents, and whether a comparison is statistically meaningful. A clear interpretation should identify the target, the model’s group-level assumptions, the relevant coding, and the basis for any model comparison.

Check Your Understanding

For each question, explain the modeling issue rather than giving only a yes-or-no answer.

A study includes six randomly sampled clinics and aims to describe treatment response across clinics like those sampled. What is the key concern in deciding whether clinic should be fixed or random?
A model contains a random intercept for site. What does this allow to vary, and what treatment-effect assumption remains if no random slope is included?
With treatment coding and an A-by-B interaction, what must be stated before interpreting the coefficient for A?
Why can’t REML log likelihoods from models with different fixed-effects structures be compared directly as a test of the added fixed effect?
In the worked cell-mean example, what is the A difference at B1, and why is it not the same as the treatment-coded coefficient $\beta_A$?