Random, nested and split-plot factors

Fixed or random?

The difference lies not in the factor itself but in the question. A fixed factor – the levels in the experiment are exactly those that interest us, and "the inference space is the specific set of levels tested" (Montgomery, 2020, p. 493). A random factor – the levels "were chosen at random from a larger population of levels, and the researcher wants to draw conclusions about the whole population" (ibid.). Maxwell, Delaney and Kelley (2018) put it in terms of replication: a fixed factor is one whose same levels would also be included in a replication of the study; for a random one the goal is "to generalize not to the effects of a few particular levels, but to a population of levels of which the selected levels are merely representative" (p. 548). Such factors are "sampled" rather than "manipulated" (p. 549): schools, classes, therapists, operators, batches, days – and the most common random factor of all, the participants (p. 547).

Even when the schools are not a true random sample – for example those that agreed to take part – the question is whether the conclusion is meant to apply only to them. If so, the factor is fixed; if one wants to speak about "schools in general" – random (Maxwell et al., 2018, p. 549). And for groups chosen for convenience, Murray and Blitstein (2003, as cited in Maxwell et al., 2018, p. 607) write that it is not only permissible but "necessary" to define the group as a random effect. Ignoring a random factor, or defining it as fixed, "can lead to some of the tests being positively biased" (p. 551).

One random factor: variance components and intraclass correlation

In the one-way random model, y_ij = μ + τ_i + ε_ij, the effects τ_i are a random variable with variance σ²_τ, and the variance of an observation is σ²_τ + σ² – the "variance components" (Montgomery, 2020, p. 101). The computation is identical to an ordinary analysis of variance; "the conclusions, however, are quite different" (p. 103): the hypothesis is H₀: σ²_τ = 0, about the whole population of levels. From the expectations E(MS_Tr) = σ² + nσ²_τ and E(MS_E) = σ² follow the estimates σ̂² = MS_E and σ̂²_τ = (MS_Tr − MS_E)/n (Equations 3.48-3.51, p. 103).

σ²_τ itself has no exact confidence interval (p. 106), but the ratio σ²_τ/(σ²_τ + σ²) – the intraclass correlation (ICC) – does (Equations 3.56-3.57, p. 106), and the page computes it. The ICC "indicates the proportion of variance attributable to the random factor" (Maxwell et al., 2018, p. 547), and it is analogous to ω² for fixed effects (p. 565). In Montgomery's looms example (Example 3.10) – an ICC between .38 and .98 (p. 106). In school achievement data a typical ICC is .15-.25 (Hedges & Hedberg, 2007, as cited in Maxwell et al., 2018, p. 608).

Expected mean squares: how the denominator is chosen

In an ordinary analysis of variance every test is against the residual. When there is a random factor this is no longer correct: "when there is a random factor in the design, any interaction of it with another factor affects the mean square of the main effect of the other factor" (Maxwell et al., 2018, p. 555). For example, in curriculum (fixed) × schools (random): if the curriculum works better in some of the schools, the difference between the curricula in the sample depends on which schools were sampled – and a test against the residual "gives F values that are too large, unless the interaction in the population is exactly zero" (p. 558).

The tool is the expected mean square E(MS) of each term. Montgomery (2020, section 13.4, pp. 505-507) gives rules for any balanced design – crossed, nested or combined: the expectation of a term is σ², plus its own component, plus the components of all the random terms that contain it; "the coefficient of each component is the number of observations at each distinct value of it" (p. 506). F is built so that the numerator and denominator "differ only by the component or the fixed factor of interest" (p. 505). The page applies the rules exactly this way, and shows the expectations in Table 2 – and these can be compared with Tables 13.10, 14.1, 14.11 and 14.17 in the book. Maxwell et al. (2018) summarize with a rule of thumb: "besides the effect itself, the only contributors to the expectation are interactions with random factors or nested random effects" (p. 593).

Mixed model: restricted or unrestricted?

When a fixed factor A is crossed with a random factor B, there are two versions of the interaction. In the restricted model the sum of the interaction effects across the levels of A is zero at each level of B (Montgomery, 2020, Equation 13.8, p. 500); then E(MS_B) = σ² + anσ²_β, and B is tested against the residual. In the unrestricted model there is no such restriction (p. 502); then E(MS_B) also includes the interaction, and B is tested against MS_AB (p. 503). The test of the fixed factor is identical in both.

Montgomery "prefers the restricted model" and notes that it is "actually slightly more general" (p. 503); Maxwell et al. (2018) "reject the unrestricted approach", because in it the null hypothesis about the random factor no longer means that the marginal means are equal (p. 606). On the other hand, JMP and Minitab (by default) work with the unrestricted model (Montgomery, 2020, p. 504), as does lmer in R; and according to Maxwell et al. (2018) – so do SPSS and SAS (pp. 605-606). Therefore the page's default is the restricted model, with a switch to the unrestricted model; "Interpreting the result" notes which tests change when switching. The relationship between the components in the two models: σ²_β = σ²_γ + σ²_αγ/a – in the book (p. 503) the formula is printed the other way round (finding E3, see below); in Example 13.2: 10.2798 + (−0.1399)/3 = 10.2332.

When there is no exact test: approximate F

With two or more random factors there is sometimes no single mean square whose expectation fits. Then a combination is built: MS' = MS_r + … + MS_s in the numerator and MS'' = MS_u + … + MS_v in the denominator, so that the difference in expectations is exactly the component being tested, and F = MS'/MS'' is approximately F-distributed with Satterthwaite's degrees of freedom (1946, as cited in Montgomery, 2020): p = (ΣMS)²/Σ(MS²/f) for each side (Equations 13.14-13.18, p. 509). Mean squares can also be subtracted, but Montgomery prefers a combination that "has no mean squares with a negative coefficient" (p. 512) – and the page always builds such a combination, and chooses the shortest one. In Example 13.5 (pressure drop) the test of temperature is (MS_A + MS_ABC)/(MS_AB + MS_AC) = 2.22 with 2.11 and 7.88 degrees of freedom (pp. 511-512), and Minitab, which subtracts, gets 2.30. This test is "only approximately F-distributed" (Maxwell et al., 2018, p. 560).

Pooling. An alternative is to remove a negligible interaction from the model and pool it with the residual. But deciding about an interaction from the data "can be dangerous", because it is exposed to errors of both types (Montgomery, 2020, p. 508). The practical rule: "if the original error mean square has six or more degrees of freedom – do not pool"; and otherwise only if the F of the term is not significant at a high α, such as 0.25 (p. 509). Maxwell et al. (2018) add that the conservative approach is to keep using the interaction as the denominator regardless of the result of the preliminary test (p. 562), and that the most practical way to increase power is simply to add levels to the random factor (p. 563). On the page, terms to pool are checked in the settings.

Nested factors

"A factor is nested in another factor if each of its levels appears together with only one level of the other factor" (Maxwell et al., 2018, p. 573): classes in schools, batches from suppliers, therapists in treatment methods. Montgomery's practical test: "if the levels of the factor can be renumbered arbitrarily" – batch 1 from supplier A and batch 1 from supplier B are not the same batch – "the factor is nested" (2020, p. 519). Therefore there is no interaction between a nested factor and the factor it is nested in (p. 519), and its sum of squares is the sum of the differences between its levels within each parent (p. 521). On the page the nested levels can be numbered afresh in each parent (1-4 for each supplier) or consecutively (1-12); when the numbering is consecutive, the page detects the nesting by itself.

The nested factor is usually random (Maxwell et al., 2018, p. 574), and then it is the denominator of the factor above it: the comparison between treatment methods is tested against the variance between the therapists, not between the participants – "the question of whether male trainees rate differently from female trainees is properly answered relative to the variance among the trainees" (p. 583). The same applies to studies in which whole classes are assigned to conditions (p. 575). Wampold and Serlin (2000, as cited in Maxwell et al., 2018, p. 590) show that ignoring a random nested factor with an ICC of .3 can inflate the Type I error rate threefold or more. Analyzing a nested design as if it were crossed produces a meaningless "interaction"; the SS_B + SS_A×B of the crossed analysis is exactly SS_B(A) (Montgomery, 2020, pp. 522-524). Combined designs – some factors crossed and some nested, as in Example 14.2 (fixtures × layouts, operators nested in layout) – are analyzed with the same rules (pp. 530-534).

Split-plot

Sometimes "it is not possible to run the experiment in a completely random order" (Montgomery, 2020, p. 534). In the example of the tensile strength of paper (section 14.4) three pulp preparation methods are tested at four temperatures: a batch is prepared with one method, divided into four parts, and each is baked at a different temperature – "nine batches in all" instead of 36 (p. 535). The batch is a whole plot, the method – the whole-plot factor, and each part – a subplot (p. 535). The principle: "the hard-to-change factors form the whole plots, and the easy-to-change ones – the subplots" (p. 539). The design has two errors: the whole-plot factor is tested against replicate × whole-plot factor (the whole-plot error), and the subplot factor and the interaction – against subplot errors, which are usually smaller, "because the subplots are more homogeneous" (p. 537). Therefore "it is advisable to put the factor of greater interest in the subplots" (p. 535).

In the full model (Equation 14.15) the replicates are random, and the replicates and replicate × subplot have no test (p. 536); if the interactions with the replicates in the subplot are assumed to be negligible, they are pooled and a simpler model is obtained (Equation 14.16, p. 537) – this is the model JMP uses (ibid.), and also aov(y ~ A*B + Error(rep/A)) in R. On the page: in the "Split-plot design" tab – the run order and the expected table; in the analysis – replicate random, both factors fixed, and for pooling check replicate × subplot and replicate × the interaction. In the behavioral sciences the same structure is also called a "mixed design": a between-subjects factor and a within-subjects factor, with the participants nested in the between-subjects factor (Maxwell et al., 2018, pp. 689-690) – on the page: the participant is random and nested in the group, and crossed with the within-subjects factor.

Split-plot "by accident". When a factorial experiment is run in an order that minimizes the number of changes of a hard-to-change factor, and then analyzed as if the order had been random, "one tends to make too many Type I errors on the whole-plot factors, and too many Type II errors on the subplot factors" (Montgomery, 2020, p. 538). And a rule of thumb: "if all the effects and interactions are significant, it is worth checking carefully how the experiment was actually run" (p. 540).

Variance components: estimates, confidence intervals and negative values

The page estimates the variance components by the ANOVA method: solving the expectation equations of the random terms "after eliminating the mean squares containing fixed factors" (Montgomery, 2020, p. 501). This is a method of moments, which "does not require the normality assumption" (p. 103). The standard error – √(Σ2c_i²MS_i²/f_i), as in the JMP output in the book (Tables 3.20 and 14.7, pp. 109, 527). The error variance has an exact chi-square confidence interval (Equation 13.19, p. 512); the rest – the Satterthwaite approximation (Equations 13.20-13.23, pp. 514-515), which "can be used only when the estimate is positive" (p. 515), and works well when the degrees of freedom are large and the coefficients are positive (p. 516). In Example 13.6: σ̂²_τβ = 31.36, r = 5.19, and the interval 12.39 to 179.83 (in the book 180.84, because of the rounding of the quantile – finding E25).

A negative estimate. The difference between two mean squares can come out negative. Montgomery lists three ways: accept it as a sign that the component is zero – but "zero in place of the negative estimate may disturb the properties of other estimates"; re-estimate with a method that always gives non-negative values; or see it as a sign that the model is wrong (p. 104). Another practical way: when the p of the interaction is large, remove it and reanalyze – "a relatively easy approach that often works nearly as well as more sophisticated methods" (p. 497). Maxwell et al. (2018) note that the convention is to set the estimate to zero (p. 566), and that it is usually a sign of a small sample or a component close to zero (p. 964). The page shows the negative estimate as it is (in orange), and computes the share of the total variance as if it were zero.

ANOVA or REML? REML is "the accepted version of maximum likelihood for estimating variance components" (Montgomery, 2020, p. 108), and its advantage is that "it produces standard errors, so confidence intervals are easy to find" (p. 513). In a balanced design "the two methods agree" (pp. 109, 499) – but only when REML is not restricted to non-negative values (as in JMP). lmer in R restricts, so when one of the ANOVA estimates is negative, it sets it to zero and changes the others too: in Example 13.1 – parts 10.2513 instead of 10.2798 (finding E9). In the page's example – method × school: −0.089 with ANOVA, 0 with REML. In an unbalanced design a mixed model is used – "the method of choice when sample sizes are unequal" (Maxwell et al., 2018, p. 954); the page analyzes only balanced designs, and for an unbalanced design offers an R script with lmer.

How to read the output?

Table 1 – the structure of the design as the page understood it: the type of each factor, the number of levels, and what it is nested in. It is worth checking it first – most errors in such an analysis are errors in the structure.

Table 2 – the analysis of variance: SS, df and MS as usual; the E(MS) column – what each mean square estimates (σ²_X – the variance component of a random term, θ²_X – the effect size of a fixed term, as in Maxwell et al., 2018, p. 557); "test denominator" – what each term is tested against. F with ^a – approximate F, with non-integer degrees of freedom. "No test" – there is no suitable combination; "serves as residual" – the highest-order term when there are no replicates in the cells.

Table 3 and Figure 1 – the variance components: estimate, standard error, share of the total variance (the ICC of each component), and confidence interval. The "Estimate" column shows the formula – it can be reproduced by hand from Table 2. ICC (in a one-way design) – with an exact confidence interval. Means – for each fixed factor, with a standard error from the denominator of its test (Montgomery, 2020, p. 501).

In the example (three teaching methods × four schools, two classes per method in each school, three pupils per class): the method is tested against method × school, F(2, 6) = 5.18, p = .049 – barely significant, with only six degrees of freedom; the school, in the restricted model, is tested against the classes, F(3, 12) = 2.14, p = .148; the differences between the classes are significant, F(12, 48) = 2.48, p = .013, and the class accounts for about 30% of the variance. The method × school interaction – a small negative estimate.

Assumptions and warnings

Balance. The expectation rules apply to balanced designs only – the same number of levels in each parent and the same number of observations in each cell (Montgomery, 2020, p. 505). For a random nested factor with unequal sizes "there is no exact test" (Maxwell et al., 2018, p. 586).

Distribution. The random effects and the errors are normally distributed, independent, with homogeneous variance – also within each level of the fixed factor (Maxwell et al., 2018, p. 585). Approximate F and Satterthwaite intervals – approximations only.

Power. "Power depends both on the number of participants and on the number of levels of the random factor – and the number of levels is usually the more important consideration" (Maxwell et al., 2018, p. 572). With four schools, the test of the method has only six degrees of freedom, even if each school has hundreds of pupils.

Many crossed random factors. When there is more than one random factor besides those in the term, there is usually no exact test (Maxwell et al., 2018, p. 559), and the approximate test requires many levels.

What this page does not compute. Unbalanced designs (for them – an lmer script), staggered nested designs (staggered; Montgomery, 2020, Section 14.1.4), unreplicated fractional split-plot analyzed with half-normal plots (section 14.5.1), MLS intervals (Equations 13.24-13.27), random slopes and group-level predictors (Maxwell et al., 2018, Chapter 16), and ICC for rater reliability (chapter 11). Split-split-plot and strip-plot designs (sections 14.5.2-14.5.3) are analyzed on the page as a design with a random replicate and three or two fixed factors, with the appropriate terms pooled. On the "Design of experiments" page – factorial experiments with fixed factors, 2^k and fractions; and on the "Response surface methodology" page – quantitative factors and the optimum.

And how do you do it in SPSS?

In Analyze › General Linear Model › Univariate put the random factors into Random Factor(s); a nested term is written in syntax, /DESIGN = A B(A). SPSS computes expected mean squares and builds tests, including approximate ones – in the unrestricted model (Maxwell et al., 2018, pp. 605-606), so for comparison choose "unrestricted" on the page. The tests of the restricted model are added with /TEST = method VS method*school (p. 568); the page's syntax writes them by itself. Variance components – with VARCOMP and /METHOD=SSTYPE(3) (ANOVA estimates, unrestricted). The menu paths are general information about the software, and the syntax was not run in SPSS. PSPP, the free software that reads SPSS syntax, has no UNIANOVA and VARCOMP, and its GLM command does not accept nested terms; we checked the sums of squares of a crossed design in it, and they are identical to the page.

And how do you do it in R?

Two ways. The classical one – lm for the sums of squares, then tests and components according to the expectations, by hand; this is what the page's script does, with the denominators the page chose. aov with Error() gives the tests of the unrestricted model (for example aov(y ~ A*B + Error(rep/A)) for split-plot). The modern one – lmer from the lme4 package: (1 | school) for a random factor, (1 | method:school) for its interaction, and (1 | school:class) for a nested factor; lmerTest adds F tests. lmer is an unrestricted model with non-negative components – see "ANOVA or REML". The script below was run as is on the English sample file, and the values in the comments are identical to the page. The site rule for comments in R: a UTF-8 file, every comment line starts with an English word, and variable names in English.

# Random and nested factors – Professor Stat. This file is saved as UTF-8; every comment line starts with an English word
library(readxl)   # Reading Excel
library(lme4)     # Mixed models (install.packages("lme4"))

# Loading: long format – row = pupil; the sheet "Data"
d <- as.data.frame(read_excel("C:/Rdata/DATA_mix_example_eng.xlsx", sheet = "Data"))
for (v in c("method", "school", "class")) d[[v]] <- factor(d[[v]])

# Sums of squares: method, school, the interaction, and class nested in method × school
options(contrasts = c("contr.sum", "contr.poly"))
a <- anova(lm(score ~ method + school + method:school + method:school:class, data = d))
MS <- setNames(a[["Mean Sq"]], c("method", "school", "method:school", "class", "Residuals")); round(MS, 3)
# Result: 511.097 212.759 98.690 99.222 40.056

# Tests in the restricted model (the page's default): school is tested against class, not against method × school
c(method = MS[["method"]] / MS[["method:school"]], school = MS[["school"]] / MS[["class"]],
  method.school = MS[["method:school"]] / MS[["class"]], class = MS[["class"]] / MS[["Residuals"]])
# Result: 5.179 (2, 6), 2.144 (3, 12), 0.995 (6, 12), 2.477 (12, 48)

# Variance components by the ANOVA method, in the restricted model
c(school = (MS[["school"]] - MS[["class"]]) / 18, method.school = (MS[["method:school"]] - MS[["class"]]) / 6,
  class = (MS[["class"]] - MS[["Residuals"]]) / 3, residual = MS[["Residuals"]])
# Result: 6.308, -0.089, 19.722, 40.056

# For comparison: REML – unrestricted model, and non-negative components (method × school is set to zero)
m <- lmer(score ~ method + (1 | school) + (1 | method:school) + (1 | method:school:class), data = d)
print(VarCorr(m), comp = "Variance")
# Result: school 6.318, method × school 0.000, class 19.663, residual 40.056

Reporting the results

The report should make it possible to reproduce the design and the analysis: which factors are fixed and which random and why, what is nested in what, how many levels and how many observations in each cell; which model (restricted or unrestricted) and which terms were pooled; for each test – F with the degrees of freedom of the actual denominator (not of the residual, when the denominator is different) and p; and the variance components or the ICC, with confidence intervals. APA (2020, section 6.43) asks for "sufficient information for the reader to corroborate the analyses". The formatting rules (section 6.36): p and ICC without a zero before the decimal point; F and variance components – with it. The results paragraph in the output is built this way, in English.

Common mistakes

Testing everything against the residual. When there is a random factor, this inflates the F of the fixed factors crossed with it or above it (Maxwell et al., 2018, pp. 558, 590).

Analyzing a nested design as if it were crossed. The result is an "interaction" that has no interpretation (Montgomery, 2020, pp. 522-524). Ask: is batch 1 from supplier A the same batch as batch 1 from supplier B?

Ignoring the whole plots. A split-plot analyzed as a completely randomized design – too many Type I errors on the whole-plot factor and Type II errors on the subplots (p. 538).

Few levels of the random factor. Three or four schools – and the test of the method has almost no power, no matter how many pupils (Maxwell et al., 2018, p. 572).

Comparing variance components across software without checking the model. Restricted or unrestricted, ANOVA or REML, with or without setting to zero – each gives a different number for the same component. In the paper strength example, JMP reports 2.4757 for the replicates (Table 14.19, p. 538 – unrestricted, pooled), and the restricted model gives 2.90.

Pooling an interaction just because it is not significant. With six or more degrees of freedom for error – do not pool; otherwise – only if the p is above 0.25 (Montgomery, 2020, p. 509).

References

American Psychological Association. (2020). Publication manual of the American Psychological Association (7th ed.). https://doi.org/10.1037/0000165-000

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2018). Designing experiments and analyzing data: A model comparison perspective (3rd ed.). Routledge.

Montgomery, D. C. (2020). Design and analysis of experiments (10th ed.). Wiley.

Satterthwaite (1946), Murray and Blitstein (2003), Hedges and Hedberg (2007) and Wampold and Serlin (2000) – as cited in Montgomery (2020) and in Maxwell et al. (2018). Under the site rule, the reference list includes only sources that were read in the original; citations are by printed page. Read: section 3.9 and chapters 13-14 of Montgomery, and chapter 10, sections of chapters 11-12 and chapter 16 of Maxwell et al. The page's computation is verified against R 4.3.3 (lm, aov with Error, lme4, lmerTest), and against all the numbers printed in Examples 3.10, 13.1-13.6, 14.1-14.2 and the paper strength example (section 14.4) of Montgomery and in Tables 10.5, 10.9 and 12.19 of Maxwell et al.