Stats

Is a GLMM appropriate for a split-split-plot experiment with very limited replication and repeated measures?

Step-by-step statistics solution: Is a GLMM appropriate for a split-split-plot experiment with very limited replication and repeated measures?

As an Amazon Associate, I earn from qualifying purchases. For more practice problems like this, see Schaum’s Outline of Statistics, 6th Edition.

1. Restating the problem in plain language

You have a split‑split‑plot experiment that was run at three sites.

Level Factor (treatment) Number of experimental units
Whole‑plot canopy (3 levels) 8 whole plots (2 per site)
Sub‑plot fencing (2 levels) 16 sub‑plots (2 per whole plot)
Sub‑sub‑plot planting (2 levels) 32 sub‑sub‑plots (2 per sub‑plot)

Each sub‑sub‑plot was counted once per year (3 years at two sites, 4 years at the third).
The response is a count of plants, so a negative‑binomial GLMM was tried:

glmer.nb( plant_count ~ canopy + fencing + planting +
          (1|trial_year) + (1|site) +
          (1|wholeplot_id/subplot_id/subsubplot_id),
          data = count_data )

A colleague warned that the design is “too small” for such a model and that R should have complained.

You want to know:

  1. Is the design hopeless for estimating the three treatment effects?
  2. How can you tell whether the model is inappropriate (wrong structure) versus merely imprecise because of limited data?
  3. What diagnostics should be examined?

Below is a step‑by‑step worked solution that addresses these questions and shows how to proceed responsibly.

2. Step‑by‑step solution

Step 1️⃣ Clarify the statistical goals and the hierarchy of the experiment

Hierarchy Random component (potential) Fixed component (interest)
Site Site (3 levels) – represents a sample of the population of sites – (often treated as fixed because only 3 levels)
Year Trial year (3‑4 levels) – repeated measures within a site – (could be fixed if you want to test year effects)
Whole‑plot Whole‑plot ID (nested in site) – the experimental unit on which canopy was applied Canopy treatment
Sub‑plot Sub‑plot ID (nested in whole‑plot) – the unit on which fencing was applied Fencing treatment
Sub‑sub‑plot Sub‑sub‑plot ID (nested in sub‑plot) – the unit on which planting was applied Planting treatment
Observation Observation within a sub‑sub‑plot in a given year – (counts)

Because the three treatments are applied at different levels, the random‑effects structure must respect that nesting. The model should allow for correlation within the same whole‑plot, sub‑plot, and sub‑sub‑plot across years.

Step 2️⃣ Realise the limitations imposed by the few levels of the highest random factors

  • Site has only 3 levels. Variance components estimated from ≤ 5 levels are notoriously unstable; standard software will still give a numeric estimate, but the standard error of that variance is huge and the likelihood surface is flat.
  • Whole‑plot has only 8 levels (2 per site). Again, the between‑whole‑plot variance is very poorly identified.

When a random factor has so few levels, the usual remedy is to treat it as a fixed effect (or to drop it if it is not of primary interest). This does not invalidate the experiment; it just changes the way we account for that source of variation.

Step 3️⃣ Choose a parsimonious mixed‑model that matches the data

A sensible starting point is:

library(lme4)
mod1 <- glmer.nb(
  plant_count ~ canopy * fencing * planting +      # all main & interaction effects
                year +                             # optional fixed year effect
                (1|site) +                         # treat site as random – keep an eye on convergence
                (1|site:wholeplot) +               # whole‑plot variance (nested in site)
                (1|site:wholeplot:subplot) +       # subplot variance (nested)
                (1|site:wholeplot:subplot:subsub), # sub‑sub‑plot variance (nested)
  data = count_data,
  control = glmerControl(optimizer = "bobyqa")
)

Why this formulation?

  • canopy, fencing, and planting are fixed because they are the treatments of scientific interest.
  • The three random terms represent the three hierarchical “plot” levels. The syntax site:wholeplot is equivalent to “whole‑plot nested in site”.
  • Adding year as a fixed effect (or as a random slope if you think the treatment response changes over years) handles the repeated‑measure aspect without inflating the random‑effects structure.

If the model fails to converge or gives singular fits (variance estimates at zero), simplify:

Situation Action
Site variance ≈ 0 Drop (1|site) and instead include site as a fixed factor (+ site)
Whole‑plot variance ≈ 0 Drop (1|site:wholeplot) – the whole‑plot level is then absorbed into the residual
Singular fit Remove the highest‑level random term first (site), then whole‑plot, etc.

Step 4️⃣ Check identifiability and model adequacy

4.1 Convergence & singularity

summary(mod1)          # look for “Variance … = 0 (singular)”
isSingular(mod1)       # TRUE → simplify

If isSingular() returns TRUE, the random structure is over‑parameterised for the data.

4.2 Over‑dispersion & distributional fit

  • Over‑dispersion is already handled by the negative‑binomial variance function, but you should still verify that the estimated shape parameter (theta) is reasonable.

  • Use the DHARMa package for residual diagnostics (simulation‑based, works for GLMMs):

library(DHARMa)
sim_res <- simulateResiduals(mod1, n = 250)
plot(sim_res)                # uniformity, outliers, heteroscedasticity
testDispersion(sim_res)      # residual over‑dispersion
testZeroInflation(sim_res)   # zero‑inflation check

If the residuals are uniform and no serious patterns appear, the model’s distributional assumption is acceptable.

4.3 Confidence intervals & precision

Obtain profile or bootstrapped confidence intervals for the fixed‑effect coefficients:

confint(mod1, parm = "beta_", method = "profile")   # profile CI
confint(mod1, parm = "beta_", method = "boot", nsim = 500) # bootstrap CI

Wide intervals (often spanning zero) signal that the data simply do not contain enough information to estimate the effect precisely. That is imprecision, not model misspecification.

Step 5️⃣ Distinguish model‑inadequacy from limited‑information

Indicator Interpretation
Convergence warnings / singular fit Model is too complex for the data → drop/re‑parameterise random terms.
Residual patterns (non‑uniform, heteroscedastic) Mis‑specified distribution (e.g., zero‑inflation, omitted covariate) → try zero‑inflated NB, add covariates, or transform response.
Very large standard errors for fixed effects while the model converges and residuals look good Limited replication – the design is capable of estimating the effect, but the data are insufficient for precision.
Variance estimates on higher‑level random terms are zero The corresponding source of variation is not estimable; either treat it as fixed or remove it.

Thus, diagnostics (convergence, residuals) tell you if the model is unsuitable. Precision measures (standard errors, confidence‑interval width) tell you if the data are insufficient.

Step 6️⃣ Practical recommendations for this particular dataset

  1. Treat “site” as a fixed factor (only 3 levels).
  2. Keep the whole‑plot and subplot random effects if the model does not become singular; otherwise drop the highest‑level random term.

  3. Model the repeated‑measure structure by adding a random intercept for subsubplot_id and, optionally, a random slope for year within each sub‑sub‑plot:

    mod2 <- glmer.nb(
  plant_count ~ canopy * fencing * planting + year +
    site +                                   # fixed
    (1|site:wholeplot) +                     # whole‑plot random
    (1|site:wholeplot:subplot) +             # subplot random
    (1|subsub_id) +                          # sub‑sub‑plot random (repeated measures)
    (0 + year | subsub_id),                  # random slope for year (optional)
  data = count_data
)
  4. Report:
    • estimated fixed‑effect coefficients with profile/bootstrapped CI,
    • the variance components that could be estimated (with their SE),
    • the fact that site was fitted as a fixed effect because only three sites were sampled,
    • the diagnostics (DHARMa plots, over‑dispersion test).
  5. Interpret cautiously: if a CI for a treatment effect includes zero, the data do not provide evidence for an effect; it is not a flaw of the model but a limitation of the experimental replication.

Step 7️⃣ Optional: Bayesian alternative

Bayesian hierarchical models (e.g., with brms or rstanarm) can incorporate weakly‑informative priors on variance components, which stabilises estimation when the number of levels is small. The workflow is the same, but you obtain posterior credible intervals that naturally reflect the uncertainty due to few sites.

library(brms)
bayes_mod <- brm(
  plant_count ~ canopy * fencing * planting + year + site +
    (1|site:wholeplot) + (1|site:wholeplot:subplot) + (1|subsub_id),
  data = count_data,
  family = negbinomial(),
  prior = c(set_prior("normal(0,5)", class = "b"),
            set_prior("cauchy(0,2)", class = "sd")),
  iter = 4000, warmup = 1000, chains = 4
)

The posterior draws give a full picture of the uncertainty, which is particularly valuable with only three sites.

3. Final answer (summary)

  • The design does not rule out a meaningful analysis, but the very small number of sites (3) and whole‑plots (8) means that the highest‑level variance components are poorly estimated.
  • A reasonable GLMM can be fitted if the random‑effects structure is simplified to match the data: treat site as a fixed effect, keep whole‑plot and subplot random effects only if they are not singular, and model the repeated measurements with a random intercept (and optionally a random slope) for each sub‑sub‑plot.
  • Model adequacy is checked by:
    1. Convergence and singular‑fit diagnostics (summary(), isSingular()),
    2. Residual diagnostics (DHARMa plots, over‑dispersion, zero‑inflation tests),
    3. Confidence‑interval width for fixed effects (imprecision vs. misspecification).
  • If the model converges, residuals look good, and the only problem is wide confidence intervals, the model is appropriate but the data simply do not provide precise estimates.
  • If convergence

Original question: Is a GLMM appropriate for a split-split-plot experiment with very limited replication and repeated measures? on Cross Validated (Stats Stack Exchange), licensed CC BY-SA.

More Stats Problems