STAT 516 Lec 10 (accessible)

Randomized complete block split-plot design

Author

Karl Gregory

Published

April 17, 2026

Alfalfa data from Dr. Longnecker’s notes

Six fields; three plots in each field; four sub-plots in each plot.
Each plot randomly assigned a type of alfalfa.
Each sub-plot randomly assigned a cutting date.
Response for each subplot is yield in tons/acre in the following year.

alfalfa <- data.frame(yield = c(2.17,1.88,1.62,2.34,1.58,1.66,
                                1.56,1.26,1.22,1.59,1.25,0.94,
                                2.29,1.60,1.67,1.91,1.39,1.12,
                                2.23,2.01,1.82,2.10,1.66,1.10,
                                2.33,2.01,1.70,1.78,1.42,1.35,
                                1.38,1.30,1.85,1.09,1.13,1.06,
                                1.86,1.70,1.81,1.54,1.67,0.88,
                                2.27,1.81,2.01,1.40,1.31,1.06,
                                1.75,1.95,2.13,1.78,1.31,1.30,
                                1.52,1.47,1.80,1.37,1.01,1.31,
                                1.55,1.61,1.82,1.56,1.23,1.13,
                                1.56,1.72,1.99,1.55,1.51,1.33),
                      variety = as.factor(c(rep("ladak",24),
                                            rep("cossack",24),
                                            rep("ranger",24))),
                      date = as.factor(rep(c(rep("none",6),
                                             rep("sep01",6),
                                             rep("sep20",6),
                                             rep("oct07",6)),
                                           3)),
                      field = as.factor(rep(1:6,12)))

head(alfalfa,n = 26)

   yield variety  date field
1   2.17   ladak  none     1
2   1.88   ladak  none     2
3   1.62   ladak  none     3
4   2.34   ladak  none     4
5   1.58   ladak  none     5
6   1.66   ladak  none     6
7   1.56   ladak sep01     1
8   1.26   ladak sep01     2
9   1.22   ladak sep01     3
10  1.59   ladak sep01     4
11  1.25   ladak sep01     5
12  0.94   ladak sep01     6
13  2.29   ladak sep20     1
14  1.60   ladak sep20     2
15  1.67   ladak sep20     3
16  1.91   ladak sep20     4
17  1.39   ladak sep20     5
18  1.12   ladak sep20     6
19  2.23   ladak oct07     1
20  2.01   ladak oct07     2
21  1.82   ladak oct07     3
22  2.10   ladak oct07     4
23  1.66   ladak oct07     5
24  1.10   ladak oct07     6
25  2.33 cossack  none     1
26  2.01 cossack  none     2

Randomized complete block split plot design

Each EU randomly assigned to one level of whole-plot factor.
Each EU receives all levels of split-plot factor in random order.
Groups of EUs over which this is replicated are treated as blocks.

Treatment effects model for RCB split-plot design

Assume $Y_{i j k} = μ + τ_{i} + γ_{j} + (τ γ)_{i j} + C_{k} + (τ C)_{i k} + ε_{i j k},$ for $i = 1, \dots, a$ , $j = 1, \dots, b$ and $k = 1, \dots, c$ , where

$Y_{i j k}$ is the response of sub-plot $j$ from whole plot $i$ in block $k$ .
the $τ_{i}$ are treatment effects for the whole-plot factor.
the $γ_{j}$ are treatment effects for the split-plot factor.
the $(τ γ)_{i j}$ are interaction effects between the factors.
the $C_{k}$ are independent $Normal (0, σ_{C}^{2})$ block effects.
the $(τ C)_{i k}$ are indep $Normal (0, σ_{A C}^{2})$ whole-plot $\times$ block effects.
the $ε_{i j k}$ are independent $Normal (0, σ_{ε}^{2})$ error terms.

Define the cell means as $μ_{i j} = μ + τ_{i} + γ_{j} + (τ γ)_{i j}, i = 1, \dots, a, j = 1, \dots, b .$

Goals in the RCB split-plot design

In the RCB split-plot model we will focus on how to:

Decompose the variability in the $Y_{i j k}$ into its sources.
Test for significance of main effects and interaction.
Estimate the variance components $σ_{C}^{2}$ , $σ_{A C}^{2}$ , and $σ_{ε}^{2}$ .
Make comparisons between treatment means.
Check whether the model assumptions are satisfied.

Sums of squares for the RCB split-plot design

SS	Symbol	Formula
Total	${SS}_{Tot}$	$\sum_{i = 1}^{a} \sum_{j = 1}^{b} \sum_{k = 1}^{c} (Y_{i j k} - {\bar{Y}}_{. . .})^{2}$
A	${SS}_{A}$	$b c \sum_{i = 1}^{a} ({\bar{Y}}_{i . .} - {\bar{Y}}_{. . .})^{2}$
C	${SS}_{C}$	$a b \sum_{k = 1}^{c} ({\bar{Y}}_{. . k} - {\bar{Y}}_{. . .})^{2}$
AC	${SS}_{AC}$	$b \sum_{i = 1}^{a} \sum_{k = 1}^{c} ({\bar{Y}}_{i . k} - ({\bar{Y}}_{i . .} + {\bar{Y}}_{. . k} - {\bar{Y}}_{. . .}))^{2}$
B	${SS}_{B}$	$a c \sum_{j = 1}^{b} ({\bar{Y}}_{. j .} - {\bar{Y}}_{. . .})^{2}$
AB	${SS}_{AB}$	$c \sum_{i = 1}^{a} \sum_{j = 1}^{b} (Y_{i j .} - ({\bar{Y}}_{i . .} + {\bar{Y}}_{. j .} - {\bar{Y}}_{. . .}))^{2}$
Error	${SS}_{Error}$	${SS}_{Tot} - ({SS}_{A} + {SS}_{B} + {SS}_{AB} + {SS}_{C} + {SS}_{AC})$

We have ${SS}_{Tot} = {SS}_{A} + {SS}_{B} + {SS}_{AB} + {SS}_{C} + {SS}_{AC}$ .

Expected mean squares in RCB split plot design

Source	Df	Expected mean square
A	$a - 1$	$b c θ_{A}^{2} + b σ_{A C}^{2} + σ_{ε}^{2}$
C	$c - 1$	$a b σ_{C}^{2} + b σ_{A C}^{2} + σ_{ε}^{2}$
AC	$(a - 1) (c - 1)$	$b σ_{A C}^{2} + σ_{ε}^{2}$
B	$b - 1$	$a c θ_{B}^{2} + σ_{ε}^{2}$
AB	$(a - 1) (b - 1)$	$c θ_{A B}^{2} + σ_{ε}^{2}$
Error	$a (b - 1) (c - 1)$	$σ_{ε}^{2}$

In the above

$θ_{A}^{2} = (a - 1)^{- 1} \sum_{i = 1}^{a} ({\bar{μ}}_{i .} - {\bar{μ}}_{. .})^{2}$
$θ_{B}^{2} = (b - 1)^{- 1} \sum_{j = 1}^{b} ({\bar{μ}}_{. j} - {\bar{μ}}_{. .})^{2}$
$θ_{A B}^{2} = [(a - 1) (b - 1)]^{- 1} \sum_{i = 1}^{a} \sum_{j = 1}^{b} (μ_{i j} - ({\bar{μ}}_{i .} + {\bar{μ}}_{. j} - {\bar{μ}}_{. .}))^{2}$

ANOVA table for RCB split-plot design

Source	Df	SS	MS	F value
A	$a - 1$	${SS}_{A}$	${MS}_{A}$	$F_{A} = {MS}_{A} / {MS}_{AC}$
C	$c - 1$	${SS}_{C}$	${MS}_{C}$	$F_{C} = {MS}_{C} / {MS}_{AC}$
AC	$(a - 1) (c - 1)$	${SS}_{AC}$	${MS}_{AC}$	$F_{AC} = {MS}_{AC} / {MS}_{Error}$
B	$b - 1$	${SS}_{B}$	${MS}_{B}$	$F_{B} = {MS}_{B} / {MS}_{Error}$
AB	$(a - 1) (b - 1)$	${SS}_{AB}$	${MS}_{AB}$	$F_{AB} = {MS}_{AB} / {MS}_{Error}$
Error	$a (b - 1) (c - 1)$	${SS}_{Error}$	${MS}_{Error}$
Total	$a b c - 1$	${SS}_{Tot}$

Reject $H_{0}$ : ${\bar{μ}}_{1.} = \dots = {\bar{μ}}_{a .}$ if $F_{A} > F_{a - 1, (a - 1) (c - 1), α}$ .
Reject $H_{0}$ : $σ_{C}^{2} = 0$ if $F_{C} > F_{c - 1, (a - 1) (c - 1), α}$ .
Reject $H_{0}$ : $σ_{A C}^{2} = 0$ if $F_{AC} > F_{(a - 1) (c - 1), a (b - 1) (c - 1), α}$ .
Reject $H_{0}$ : ${\bar{μ}}_{.1} = \dots = {\bar{μ}}_{. b}$ if $F_{B} > F_{b - 1, a (b - 1) (c - 1), α}$ .
R. $H_{0}$ : $μ_{i j} = {\bar{μ}}_{i .} + {\bar{μ}}_{. j} - {\bar{μ}}_{. .} \forall i j$ if $F_{AB} > F_{(a - 1) (b - 1), a (b - 1) (c - 1), α}$ .

Alfalfa data (cont)

anova() on lm() output gives wrong p-vals for variety and field.

lm_out <- lm(yield ~ variety + date + variety:date + field + field:variety, 
             data = alfalfa)
anova(lm_out)

Analysis of Variance Table

Response: yield
              Df Sum Sq Mean Sq F value    Pr(>F)    
variety        2 0.1753 0.08763  3.1198 0.0538447 .  
date           3 1.9727 0.65758 23.4123 2.789e-09 ***
field          5 4.1388 0.82775 29.4710 3.798e-13 ***
variety:date   6 0.2147 0.03579  1.2742 0.2883071    
variety:field 10 1.3574 0.13574  4.8330 0.0001022 ***
Residuals     45 1.2639 0.02809                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

y <- alfalfa$yield
y... <- predict(lm(yield ~ 1,data = alfalfa))
yi.. <- predict(lm(yield ~ variety,data = alfalfa))
y.j. <- predict(lm(yield ~ date,data = alfalfa))
y..k <- predict(lm(yield ~ field,data = alfalfa))
yij. <- predict(lm(yield ~ variety + date + variety:date,data = alfalfa))
yi.k <- predict(lm(yield ~ variety + field + variety:field,data = alfalfa))

SST <- sum((y - y...)^2)
SSA <- sum((yi.. - y...)^2)
SSC <- sum((y..k - y...)^2)
SSAC <- sum((yi.k - (yi.. + y..k - y...))^2)
SSB <- sum((y.j. - y...)^2)
SSAB <- sum((yij. - (yi.. + y.j. - y...))^2)
SSE <- SST - (SSA + SSC + SSAC + SSB + SSAB)

a <- 3
b <- 4
c <- 6
  
MSA <- SSA / (a-1)
MSC <- SSC / (c-1)
MSAC <- SSAC / ((c-1)*(a-1))
MSB <- SSB / (b-1)
MSAB <- SSAB / ((a-1)*(b-1))
MSE <- SSE / (a*(b-1)*(c-1))

FA_incorrect <- MSA / MSE
FC_incorrect <- MSC / MSE
FA <- MSA / MSAC
FC <- MSC / MSAC
FAC <- MSAC / MSE
FB <- MSB / MSE
FAB <- MSAB / MSE

pA_incorrect <- 1 - pf(FA_incorrect,a-1,a*(b-1)*(c-1))
pC_incorrect <- 1 - pf(FC_incorrect,c-1,a*(b-1)*(c-1))
pA <- 1 - pf(FA,a-1,(c-1)*(a-1))
pC <- 1 - pf(FC,c-1,(c-1)*(a-1))
pAC  <- 1 - pf(FAC,(c-1)*(a-1),a*(b-1)*(c-1))
pB <- 1 - pf(FB,b-1,a*(b-1)*(c-1))
pAB <- 1 - pf(FAB,(a-1)*(b-1),a*(b-1)*(c-1))

Correct ANOVA table:

Source	Df	SS	MS	F value	p value
A	2	0.175	0.088	0.646	0.5449
C	5	4.139	0.828	6.098	0.0076
AC	10	1.357	0.136	4.833	0.0001
B	3	1.973	0.658	23.412	0.0000
AB	6	0.215	0.036	1.274	0.2883
Error	45	1.264	0.028
Total	71	9.123

Correct p values for fixed effects with lmerTest package:

library(lmerTest) # first time run install.packages("lmerTest")

Loading required package: lme4

Warning: package 'lme4' was built under R version 4.4.1

Loading required package: Matrix

Warning: package 'Matrix' was built under R version 4.4.1


Attaching package: 'lmerTest'

The following object is masked from 'package:lme4':

    lmer

The following object is masked from 'package:stats':

    step

lmer_out <- lmer(yield ~ variety + date + variety:date + (1|field) + (1|field:variety), 
                 data = alfalfa)
anova(lmer_out)

Type III Analysis of Variance Table with Satterthwaite's method
              Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
variety      0.03626 0.01813     2    10  0.6455    0.5449    
date         1.97274 0.65758     3    45 23.4123 2.789e-09 ***
variety:date 0.21473 0.03579     6    45  1.2742    0.2883    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Rescales ${SS}_{A}$ and ${MS}_{A}$ by the factor ${MS}_{Error} / {MS}_{AC}$ .

Note: If ${\hat{σ}}_{C}^{2} = 0$ , it will compute the whole-plot effect p value differently!

interaction.plot(alfalfa$date,alfalfa$variety,alfalfa$yield)

MoMs for variance components in RCB split-plot design

The method of moments estimators for $σ_{C}^{2}$ , $σ_{A C}^{2}$ , and $σ_{ε}^{2}$ are

${\dot{σ}}_{ε}^{2} = {MS}_{Error}$
${\dot{σ}}_{A C}^{2} = \frac{{MS}_{AC} - {MS}_{Error}}{b}$
${\dot{σ}}_{C}^{2} = \frac{{MS}_{C} - {MS}_{AC}}{a b}$

It is possible to obtain negative values for ${\dot{σ}}_{A C}^{2}$ and ${\dot{σ}}_{C}^{2}$ . Use REML.

Alfalfa data (cont)

Obtain REML estimates of $σ_{C}^{2}$ , $σ_{A C}^{2}$ , and $σ_{ε}^{2}$ from lmer().

summary(lmer_out)$varcor

 Groups        Name        Std.Dev.
 field:variety (Intercept) 0.16406 
 field         (Intercept) 0.24014 
 Residual                  0.16759

Variances of some means and difference in means

Contrast	Variance	MoM variance estimator
${\bar{Y}}_{i . .} - {\bar{Y}}_{i^{'} . .}$	$\frac{2}{b c} (b σ_{A C}^{2} + σ_{ε}^{2})$	$\frac{2}{b c} {MS}_{AC}$

${\bar{Y}}_{. j .} - {\bar{Y}}_{. j^{'} .}$	$\frac{2}{a c} σ_{ε}^{2}$	$\frac{2}{a c} {MS}_{Error}$

${\bar{Y}}_{i j .} - {\bar{Y}}_{i^{'} j .}$	$\frac{2}{c} (σ_{A C}^{2} + σ_{ε}^{2})$	$\frac{2}{c} [{MS}_{AC} + (b - 1) {MS}_{Error}]$

${\bar{Y}}_{i j .} - {\bar{Y}}_{i j^{'} .}$	$\frac{2}{c} σ_{ε}^{2}$	$\frac{2}{c} {MS}_{Error}$

${\bar{Y}}_{i j .} - {\bar{Y}}_{i^{'} j^{'} .}$	$\frac{2}{c} (σ_{A C}^{2} + σ_{ε}^{2})$	$\frac{2}{c} [{MS}_{AC} + (b - 1) {MS}_{Error}]$

${\bar{Y}}_{i . .}$	$\frac{1}{b c} [b σ_{C}^{2} + b σ_{A C}^{2} + σ_{ε}^{2}]$	$\frac{1}{a b c} [{MS}_{C} + (a - 1) {MS}_{AC}]$

${\bar{Y}}_{. j .}$	$\frac{1}{b c} [a σ_{C}^{2} + σ_{A C}^{2} + σ_{ε}^{2}]$	$\frac{1}{a b c} [{MS}_{C} + (b - 1) {MS}_{AC}]$

${\bar{Y}}_{i j .}$	$\frac{1}{c} [σ_{C}^{2} + σ_{A C}^{2} + σ_{ε}^{2}]$	$\frac{1}{a b c} [{MS}_{C} + (a - 1) {MS}_{AC} + a (b - 1) {MS}_{Error}]$

Some (unadjusted) CIs in RCB split plot design

Target	$(1 - α) 100 %$ confidence interval
${\bar{μ}}_{i .} - {\bar{μ}}_{i^{'} .}$	${\bar{Y}}_{i . .} - {\bar{Y}}_{i^{'} . .} \pm t_{(a - 1) (c - 1), α / 2} \sqrt{{MS}_{AC}} \sqrt{\frac{2}{b c}}$
${\bar{μ}}_{. j} - {\bar{μ}}_{. j^{'}}$	${\bar{Y}}_{. j .} - {\bar{Y}}_{. j^{'} .} \pm t_{a (b - 1) (c - 1), α / 2} \sqrt{{MS}_{Error}} \sqrt{\frac{2}{a c}}$
$μ_{i j} - μ_{i^{'} j}$	${\bar{Y}}_{i j .} - {\bar{Y}}_{i^{'} j .} \pm t_{ν^{*}, α / 2} \sqrt{{MS}_{AC} + (b - 1) {MS}_{Error}} \sqrt{\frac{2}{c}}$
$μ_{i j} - μ_{i j^{'}}$	${\bar{Y}}_{i j .} - {\bar{Y}}_{i j^{'} .} \pm t_{a (b - 1) (c - 1), α / 2} \sqrt{{MS}_{Error}} \sqrt{\frac{2}{c}}$
$μ_{i j} - μ_{i^{'} j^{'}}$	${\bar{Y}}_{i j .} - {\bar{Y}}_{i^{'} j^{'} .} \pm t_{ν^{*}, α / 2} \sqrt{{MS}_{AC} + (b - 1) {MS}_{Error}} \sqrt{\frac{2}{c}}$

In the above $ν^{*} = \frac{{MS}_{AC} + (b - 1) {MS}_{Error}}{\frac{{MS}_{AC}^{2}}{(a - 1) (c - 1)} + \frac{(b - 1)^{2} {MS}_{Error}^{2}}{a (b - 1) (c - 1)}}$ à la Satterthwaite.

Alfalfa data (cont)

Dunnett’s comparison of marginal date means with no cutting as baseline: ${\bar{Y}}_{. j .} - {\bar{Y}}_{.1 .} \pm d_{b - 1, a (b - 1) (c - 1), α} \sqrt{{MS}_{Error}} \sqrt{\frac{2}{a c}}$ Replace, in previous slide, $t_{a (b - 1) (c - 1), α / 2}$ with $d_{b - 1, a (b - 1) (c - 1), α}$ .

In Dunnett’s table we cannot find $d_{4, 45, 0.05}$ , so take $d_{4, 40, 0.05} = 2.44$ .

Table A.5 from Mohr, Wilson, and Freund (2021)

y.1.bar <- mean(alfalfa$yield[alfalfa$date=="none"])
y.2.bar <- mean(alfalfa$yield[alfalfa$date=="sep01"])
y.3.bar <- mean(alfalfa$yield[alfalfa$date=="sep20"])
y.4.bar <- mean(alfalfa$yield[alfalfa$date=="oct07"])


se <- sqrt(MSE) * sqrt(2/(a*c))
me <- 2.44 * se

dtab <- rbind(c(y.2.bar - y.1.bar - me,y.2.bar - y.1.bar + me),
              c(y.3.bar - y.1.bar - me,y.3.bar - y.1.bar + me),
              c(y.4.bar - y.1.bar - me,y.4.bar - y.1.bar + me))

colnames(dtab) <- c("lower","upper")
rownames(dtab) <- c("sep01 - none",
                    "sep20 - none",
                    "oct07 - none")
round(dtab,3)

              lower  upper
sep01 - none -0.578 -0.305
sep20 - none -0.343 -0.070
oct07 - none -0.226  0.046

Unadjusted CIs with ls_means() from R package lmerTest:

ls_means(lmer_out,which="date")

Least Squares Means table:

          Estimate Std. Error  df t value   lower   upper  Pr(>|t|)    
datenone   1.78111    0.11255 6.1  15.825 1.50641 2.05581 3.684e-06 ***
dateoct07  1.69111    0.11255 6.1  15.026 1.41641 1.96581 5.010e-06 ***
datesep01  1.33944    0.11255 6.1  11.901 1.06474 1.61415 1.977e-05 ***
datesep20  1.57444    0.11255 6.1  13.989 1.29974 1.84915 7.646e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  Confidence level: 95%
  Degrees of freedom method: Satterthwaite

ls_means(lmer_out, pairwise = TRUE, which = "date")

Least Squares Means table:

                        Estimate Std. Error df t value      lower      upper
datenone - dateoct07   0.0900000  0.0558639 45  1.6111 -0.0225156  0.2025156
datenone - datesep01   0.4416667  0.0558639 45  7.9061  0.3291511  0.5541823
datenone - datesep20   0.2066667  0.0558639 45  3.6995  0.0941511  0.3191823
dateoct07 - datesep01  0.3516667  0.0558639 45  6.2951  0.2391511  0.4641823
dateoct07 - datesep20  0.1166667  0.0558639 45  2.0884  0.0041511  0.2291823
datesep01 - datesep20 -0.2350000  0.0558639 45 -4.2067 -0.3475156 -0.1224844
                       Pr(>|t|)    
datenone - dateoct07  0.1141598    
datenone - datesep01  4.723e-10 ***
datenone - datesep20  0.0005859 ***
dateoct07 - datesep01 1.138e-07 ***
dateoct07 - datesep20 0.0424494 *  
datesep01 - datesep20 0.0001219 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  Confidence level: 95%
  Degrees of freedom method: Satterthwaite

References

Mohr, Donna L, William J Wilson, and Rudolf J Freund. 2021. Statistical Methods. Academic Press.