Why the Randomization Structure Matters
In a split-plot experiment, treatments are assigned at more than one level. One treatment is randomized to relatively large experimental units, called whole plots; another is randomized to smaller units within them, called subplots. Responses from subplots in the same whole plot are generally dependent because those observations share the whole-plot conditions. A generalized linear mixed model (GLMM) represents that dependence with a whole-plot random effect.
The random-effects structure follows the units randomized, not merely the layout of the data table. A treatment assigned to a whole plot has whole-plot-level replication. Recording many subplot responses within one whole plot improves information about that plot, but does not create additional independent whole plots for estimating the whole-plot treatment effect.
For a GLMM, the response distribution and link are specified conditionally on the random effects. For binary outcomes, a common choice is a Bernoulli distribution with a logit link. The same design logic applies to other conditional response distributions and links.
From Randomization to a GLMM
Suppose blocks are indexed by \(i\), whole plots within blocks by \(j\), and subplots within whole plots by \(k\). Treatment \(A\) is assigned to whole plots, and treatment \(B\) is assigned to subplots. If there are multiple blocks, a random block intercept can represent variation among the sampled blocks; a separate random intercept for each whole plot represents variation among whole plots within blocks. One binary-outcome model is
Here \(a_j\) and \(b_k\) are contrast-coded treatment values, and the fixed effects describe the treatment means on the conditional log-odds scale. The block effects \(u_i\) and whole-plot effects \(v_{ij}\) are commonly modeled as independent normal random variables with mean zero and variances \(\sigma_u^2\) and \(\sigma_v^2\), respectively. Conditional on these effects, responses are assumed independent. If there is no block structure, omit \(u_i\); if whole plots are not grouped into blocks, index and model them at the appropriate level.
In this model, \(v_{ij}\) is nested within block \(i\), and all subplot observations in whole plot \(ij\) share it. The fixed treatment factors \(A\) and \(B\), however, are crossed when the design includes their combinations. A treatment factor is not itself a random effect merely because its units are nested. Whether a factor is fixed or random depends on the scientific target and sampling assumptions; randomization level determines which experimental units provide independent treatment assignments.
Worked Example: Translating a Split-Plot Layout
Hypothetical design. Six greenhouses each contain two beds. Within every greenhouse, one bed is randomly assigned standard irrigation and the other reduced irrigation. Each bed contains three subplots, and the three fertilizer treatments are randomly assigned once each within that bed. A binary outcome records whether a plant meets a prespecified growth criterion. Identify the treatment and random-effects structure.
There are \(6\times2\times3=36\) subplot observations. Irrigation is a whole-plot treatment because it is assigned to beds; fertilizer is a subplot treatment because it is assigned to subplots within beds. Each fertilizer level occurs with each irrigation level, so the treatment factors are crossed. But each subplot belongs to only one bed, so the experimental units are nested.
Let \(Y_{ijk}\) denote the outcome in greenhouse \(i\), bed \(j\), and subplot \(k\). Use \(a_j=-1\) for standard irrigation and \(a_j=1\) for reduced irrigation. Use two contrasts for the three fertilizer levels, denoted together by the vector \(b_k\). A suitable model has a greenhouse intercept \(u_i\), a bed-within-greenhouse intercept \(v_{ij}\), fixed irrigation and fertilizer contrasts, and their interaction:
The greenhouse effect accounts for responses sharing a greenhouse; the bed effect accounts for additional dependence among subplots sharing a bed. The model assumes independent normal random intercepts at those levels and conditional independence of the binary responses. If the six greenhouses are the only blocks and are treated as fixed design controls rather than a sample from a population of blocks, the block term may instead be handled as a fixed blocking effect. The bed random intercept remains tied to the whole-plot assignment structure.
There are only 12 bed-level assignments, not 36 independent irrigation assignments. Thus, the number of subplots does not substitute for replication of irrigation across whole plots. The model's coefficient for irrigation must be informed by variation across beds assigned to different irrigation levels, accounting for the blocking and clustering.
Why the Shared Random Effect Creates Dependence
The random intercept gives a direct mathematical representation of the dependence induced by shared whole-plot conditions. Consider two distinct subplots \(k\) and \(k'\) in the same bed and block. Their linear predictors contain the same \(u_i\) and \(v_{ij}\). If all other random effects are independent, their covariance on the latent linear-predictor scale is
For subplots in different beds of the same greenhouse, the shared block effect remains but the bed effects differ, giving covariance \(\sigma_u^2\) on this scale. This follows by expanding the covariance: independent random effects contribute zero cross-covariance, while each shared mean-zero random effect contributes its variance. These are covariances of the latent predictors; in a nonlinear GLMM they are not the observed binary-response covariances or correlations.
The distinction matters because treating the 36 responses in the example as independent would discard the dependence represented by the shared bed and greenhouse conditions. It could also give misleading uncertainty for treatment contrasts, especially contrasts assigned at the whole-plot level.
Worked Example: Interpreting Conditional Treatment Effects
Consider a hypothetical binary-outcome split-plot GLMM with two treatments at each level. Let \(A=0\) or \(1\) be assigned to whole plots and \(B=0\) or \(1\) to subplots. Suppose the model is
For a whole plot with \(u=0.3\) and \(v=-0.2\), calculate the conditional probability when \(A=1\) and \(B=1\). Substitution gives
Therefore,
This is the model-based probability for a subplot with those treatment assignments, conditional on the stated random effects. It is not the probability averaged over the population distribution of whole plots.
At \(B=1\), changing \(A\) from 0 to 1 changes the conditional log odds by \(\beta_A+\beta_{AB}=0.6+0.2=0.8\). The corresponding conditional odds ratio is \(\exp(0.8)\approx2.226\). Thus, at this subplot treatment level and for the same random effects, the modeled odds are about 2.226 times as high under \(A=1\) as under \(A=0\). This is an association implied by the model; a causal interpretation additionally depends on the design and its assumptions.
Nested Treatment Factors Are Not the Same as Nested Units
A separate question is whether the treatment factor levels are crossed or nested. If the same levels of \(B\) are used within every level of \(A\), the treatment factors are crossed, even if \(B\) is assigned to units nested within the \(A\)-level units. A factor is nested when its levels are specific to levels of another factor, so that the same named lower-level treatment levels do not recur meaningfully across the higher-level factor.
For example, suppose a hypothetical experiment compares two delivery routes. Each route has two route-specific formulations, and the formulations are meaningful only within their route. Formulation is then nested within route, not crossed with route. The model should represent the route-specific treatment combinations rather than imply that a formulation used for one route is also used for the other.
Worked Example: Recognizing Route-Specific Treatments
Suppose patients are assigned to one of two routes, oral or injection. The oral route has formulations O1 and O2; the injection route has formulations I1 and I2. Each patient is assigned to a treatment combination, and patients are grouped in clinics. The binary outcome is response. Assume, hypothetically, that clinics are sampled from a population of clinics and that treatment assignment is randomized within the study.
The four meaningful treatment combinations are oral-O1, oral-O2, injection-I1, and injection-I2. It would be misleading to define one common formulation factor with levels O1, O2, I1, and I2 and then treat those levels as interchangeable across routes. A fixed-effect cell-means specification can represent the four combinations, with a random clinic intercept to account for patients sharing a clinic. Route comparisons can be constructed from the appropriate combination means; comparisons of O1 directly with I1 are comparisons of treatment combinations, not of a formulation level common to both routes.
If each route has \(n_r\) patients per formulation, the total number of patients is \(2(n_{\text{oral}}+n_{\text{injection}})\). That patient count does not change the fact that formulation levels are route-specific. The random clinic effect is justified by the shared clinic membership, not by the formulation nesting. If treatments were instead randomized at the clinic level, treatment replication would have to be counted across clinics, not across patients.
Assumptions, Diagnostics, and Alternatives
A split-plot GLMM relies on assumptions about both the experiment and the conditional model. Check that the recorded assignments match the actual randomization: whole-plot treatments must be represented at the whole-plot level, and subplot treatments at the subplot level. The conditional response distribution and link must suit the outcome, and the fixed-effects structure must include the treatment contrasts relevant to the scientific question. The usual random-intercept formulation assumes mean-zero random effects with specified distributions, commonly independent normal effects, and conditional independence of responses after those effects and predictors are included.
Diagnostics should include checks for coding errors in unit identifiers and treatment assignments, convergence problems, and whether the fitted random-effect variances are plausible. A variance estimate near zero can indicate that a component is weakly supported, but it can also arise from limited information or a design that cannot separate components. Few whole plots can make whole-plot variation and whole-plot treatment effects difficult to estimate reliably, even when there are many subplot observations. Examine conditional fit and residual patterns appropriate to the response distribution; unexplained clustering or overdispersion can signal that the assumed dependence structure or response model is inadequate.
The main alternative is to fit a GLM as though every observation were independent. That omits the dependence from shared units and is not a faithful analysis of the split-plot randomization. Aggregating to whole-plot summaries may lose subplot-level treatment information and is not generally equivalent to fitting the GLMM. A design-based analysis that follows the actual randomization can be valuable, especially when the number of randomized whole plots is small, but it answers questions through a different inferential framework and does not automatically provide the same model-based estimates.
A random-intercept GLMM is also limited: it assumes a particular distributional form for unit heterogeneity and may be poorly supported when there are very few units at a level. It does not repair confounding, missing treatment combinations, or treatment assignment that was not randomized. Most importantly, many subplots cannot compensate for too few independently randomized whole plots.
Common Pitfall
Check Your Understanding
Answer each question, then compare with the brief solutions.
- A treatment is randomized to pens, and a second treatment is randomized to animals within each pen. Which level must be replicated to estimate the first treatment effect?
- In a split-plot design, every whole-plot treatment is paired with every subplot treatment. Are the treatment factors crossed, nested, or both?
- Two subplots share the same whole plot and block. If the model includes independent random intercepts for block and whole plot, which random effects do their linear predictors share?
- In a nested treatment design, two formulations are unique to each delivery route. Is it appropriate to interpret a formulation contrast as the same contrast across routes?
Solutions. (1) The pens: the first treatment is assigned to pens, so independent replication for that effect is at the pen level. (2) The treatment factors are crossed because all combinations occur; the experimental units are still nested because subplots belong to whole plots. (3) They share both the block and whole-plot random intercepts. (4) No. The formulations are route-specific, so a contrast within one route does not automatically define the same contrast in another.