STAT 516 Lec 09
Randomized complete block designs
2026-10-26
Nitrogen fertilization data from Kuehl (2000)
Six nitrogen timing schedules in each of four plots in a field of wheat. Response is the nitrate content in stem tissue samples.
Nitrogen fertilization
Randomized complete block design (RCBD)
- EUs belong to blocks—groups of EUs homogeneous in some way.
- Each EU in a block is randomly assigned to a fixed treatment.
- All treatments appear exactly once in each block.
- Purpose is to capture the between-block variability among the EUs.
- This helps us detect treatment effects with greater power.
Treatment effects model for the RCBD
Assume
\[ Y_{ij} = \mu + \tau_i + B_j + \varepsilon_{ij}, \quad \text{ for } i=1,\dots,a, \quad j = 1,\dots,b, \] where
- \(Y_{ij}\) is the response of the EU in block \(j\) receiving treatment \(i\).
- the \(\tau_i\) are the fixed effects of the treatment.
- the \(B_j\) are independent \(\mathcal{N}(0,\sigma_B^2)\) random block effects.
- the \(\varepsilon_{ij}\) are independent \(\mathcal{N}(0,\sigma_\varepsilon^2)\) error/interaction terms.
- \(\mu\) is an overall or baseline mean.
Define the cell means as \[ \mu_i = \mu + \tau_i, \quad i = 1,\dots,a. \]
Goals in the RCBD
In the randomized complete block design we wish to
- Visualize the data.
- Decompose the variability in the \(Y_{ij}\) into its sources.
- Estimate the variance components \(\sigma_B^2\) and \(\sigma_\varepsilon^2\).
- Test whether the treatment has any effect.
- Make comparisons between treatment means.
- Check whether the model assumptions are satisfied.
Nitrogen fertilization
Nitrogen fertilization
Nitrogen fertilization
Nitrogen fertilization
Sums of squares for the RCBD
| Sum of squares | Symbol | Formula |
|---|---|---|
| Total | \(\operatorname{SS}_{\operatorname{Tot}}\) | \(\sum_{i=1}^a\sum_{j=1}^b(Y_{ij} - \bar Y_{..})^2\) |
| Treatment | \(\operatorname{SS}_{\operatorname{A}}\) | \(b\sum_{i=1}^a (\bar Y_{i.} - \bar Y_{..})^2\) |
| Block | \(\operatorname{SS}_{\operatorname{B}}\) | \(a\sum_{j=1}^b (\bar Y_{.j} - \bar Y_{..})^2\) |
| Error | \(\operatorname{SS}_{\operatorname{Error}}\) | \(\sum_{i=1}^a \sum_{j=1}^b (Y_{ij} - (\bar Y_{i.} + \bar Y_{.j} - \bar Y_{..}))^2\) |
- We can make the decomposition \(\operatorname{SS}_{\operatorname{Tot}}= \operatorname{SS}_{\operatorname{A}}+ \operatorname{SS}_{\operatorname{B}}+ \operatorname{SS}_{\operatorname{Error}}\).
- The \(\operatorname{SS}_{\operatorname{Error}}\) is really the interaction sum of squares \(\operatorname{SS}_{\operatorname{AB}}\).
- But without replication, we cannot estimate an interaction.
- So the interaction serves as the error term.
ANOVA table for RCBD
| Source | Df | SS | MS | F value |
|---|---|---|---|---|
| A | \(a-1\) | \(\operatorname{SS}_{\operatorname{A}}\) | \(\operatorname{MS}_{\operatorname{A}}\) | \(F_{\operatorname{A}} = \operatorname{MS}_{\operatorname{A}}/ \operatorname{MS}_{\operatorname{Error}}\) |
| B | \(b-1\) | \(\operatorname{SS}_{\operatorname{B}}\) | \(\operatorname{MS}_{\operatorname{B}}\) | \(F_{\operatorname{B}} = \operatorname{MS}_{\operatorname{B}}/ \operatorname{MS}_{\operatorname{Error}}\) |
| Error | \((a-1)(b-1)\) | \(\operatorname{SS}_{\operatorname{Error}}\) | \(\operatorname{MS}_{\operatorname{Error}}\) | |
| Total | \(ab-1\) | \(\operatorname{SS}_{\operatorname{Tot}}\) |
- Reject \(H_0\): \(\mu_1 = \dots = \mu_a\) if \(F_{\operatorname{A}} > F_{a-1,(a-1)(b-1),\alpha}\).
- Reject \(H_0\): \(\sigma_B^2 = 0\) if \(F_{\operatorname{B}} > F_{b-1,(a-1)(b-1),\alpha}\).
Nitrogen fertilization
Obtain ANOVA table “manually”:
y <- ntr
y.. <- predict(lm(ntr ~ 1))
yi. <- predict(lm(ntr ~ trt))
y.j <- predict(lm(ntr ~ block))
SST <- sum((y - y..)^2)
SSA <- sum((yi. - y..)^2)
SSB <- sum((y.j - y..)^2)
SSE <- SST - (SSA + SSB)
a <- 6
b <- 4
MSA <- SSA/(a-1)
MSB <- SSB/(b-1)
MSE <- SSE/((a-1)*(b-1))
FA <- MSA / MSE
FB <- MSB / MSE
pA <- 1 - pf(FA,a-1,(a-1)*(b-1))
pB <- 1 - pf(FB,b-1,(a-1)*(b-1))Nitrogen fertilization
| Source | Df | SS | MS | F value | p value |
|---|---|---|---|---|---|
| A | 5 | 201.32 | 40.26 | 5.5917 | 0.0042 |
| B | 3 | 197.00 | 65.67 | 9.1198 | 0.0011 |
| Error | 15 | 108.01 | 7.20 | ||
| Total | 23 | 506.33 |
Nitrogen fertilization
Obtain ANOVA table with anova() function on lm() output:
Analysis of Variance Table
Response: ntr
Df Sum Sq Mean Sq F value Pr(>F)
trt 5 201 40.3 5.59 0.0042 **
block 3 197 65.7 9.12 0.0011 **
Residuals 15 108 7.2
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Expected mean squares in the RCBD
| Source | Df | Expected mean square |
|---|---|---|
| A | \(a-1\) | \(b\theta_A^2 + \sigma^2_\varepsilon\) |
| B | \(b-1\) | \(a\sigma_B^2 + \sigma^2_\varepsilon\) |
| Error | \((a-1)(b-1)\) | \(\sigma_\varepsilon^2\) |
In the above \(\theta_A^2 = (a-1)^{-1}\sum_{i=1}^a(\mu_i - \bar\mu_.)^2\).
Method of moments for variance components in RCBD
Equating \(\operatorname{MS}_{\operatorname{B}}\) and \(\operatorname{MS}_{\operatorname{Error}}\) with their expectations gives
- \(\dot{\sigma}_{\varepsilon}^2 = \operatorname{MS}_{\operatorname{Error}}\).
- \(\displaystyle \dot \sigma_B^2 = \frac{\operatorname{MS}_{\operatorname{B}}- \operatorname{MS}_{\operatorname{Error}}}{a}\)
May obtain \(\dot \sigma_B^2 < 0\), so one should use REML estimation.
Nitrogen fertilization
Obtain REML estimators of \(\sigma_B^2\) and \(\sigma_{\varepsilon}^2\) on the fertilization data.
Linear mixed model fit by REML ['lmerModLmerTest']
Formula: ntr ~ trt + (1 | block)
REML criterion at convergence: 102
Random effects:
Groups Name Std.Dev.
block (Intercept) 3.12
Residual 2.68
Number of obs: 24, groups: block, 4
Fixed Effects:
(Intercept) trt2 trt3 trt4 trt5 trt6
38.28 5.75 8.49 2.34 1.23 4.95 Nitrogen fertilization
Obtain MoMs estimators for \(\sigma_B^2\) and \(\sigma_{\varepsilon}^2\) on the fertilization data.
We have \(\dot \sigma_B = 3.122\) and \(\dot \sigma_{\varepsilon}= 2.683\).
Variances of means and difference in means
| Contrast | Variance | MoM variance estimator |
|---|---|---|
| \(\bar Y_{i.}\) | \(\frac{1}{b}(\sigma_B^2 + \sigma_{\varepsilon}^2)\) | \(\frac{2}{ab}[\operatorname{MS}_{\operatorname{B}}+ (b-1)\operatorname{MS}_{\operatorname{Error}}]\) |
| \(\bar Y_{i.} - \bar Y_{i'.}\) | \(\frac{2}{b}\sigma_{\varepsilon}^2\) | \(\frac{2}{b}\operatorname{MS}_{\operatorname{Error}}\) |
Some (unadjusted) CIs in RCB split plot design
| Target | \((1-\alpha)100\%\) confidence interval |
|---|---|
| \(\mu_i\) | \(\bar Y_{i.} \pm t_{\nu^*,\alpha/2}\sqrt{\operatorname{MS}_{\operatorname{B}}+ (b-1)\operatorname{MS}_{\operatorname{Error}}}\sqrt{\frac{2}{ab}}\) |
| \(\mu_i - \mu_{i'}\) | \(\bar Y_{i.} - \bar Y_{i'.} \pm t_{(a-1)(b-1),\alpha/2}\sqrt{\operatorname{MS}_{\operatorname{Error}}}\sqrt{\frac{2}{b}}\) |
In the above \(\displaystyle \nu^* = \frac{\operatorname{MS}_{\operatorname{B}}+ (b-1)\operatorname{MS}_{\operatorname{Error}}}{\frac{\operatorname{MS}_{\operatorname{B}}^2}{(b-1)} + \frac{(b-1)^2\operatorname{MS}_{\operatorname{Error}}^2}{(a-1)(b-1)}}\) à la Satterthwaite.
Nitrogen fertilization
Unadjusted CIs with ls_means() from R package lmerTest
Least Squares Means table:
Estimate Std. Error df t value lower upper Pr(>|t|)
trt1 38.28 2.06 6.8 18.6 33.38 43.18 4.5e-07 ***
trt2 44.03 2.06 6.8 21.4 39.13 48.93 1.8e-07 ***
trt3 46.77 2.06 6.8 22.7 41.87 51.67 1.2e-07 ***
trt4 40.62 2.06 6.8 19.7 35.72 45.51 3.0e-07 ***
trt5 39.51 2.06 6.8 19.2 34.61 44.41 3.6e-07 ***
trt6 43.23 2.06 6.8 21.0 38.33 48.12 2.0e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Confidence level: 95%
Degrees of freedom method: Satterthwaite Nitrogen fertilization
Least Squares Means table:
Estimate Std. Error df t value lower upper Pr(>|t|)
trt1 - trt2 -5.755 1.897 15 -3.03 -9.799 -1.711 0.00839 **
trt1 - trt3 -8.492 1.897 15 -4.48 -12.537 -4.448 0.00044 ***
trt1 - trt4 -2.337 1.897 15 -1.23 -6.382 1.707 0.23694
trt1 - trt5 -1.232 1.897 15 -0.65 -5.277 2.812 0.52580
trt1 - trt6 -4.947 1.897 15 -2.61 -8.992 -0.903 0.01980 *
trt2 - trt3 -2.738 1.897 15 -1.44 -6.782 1.307 0.16965
trt2 - trt4 3.418 1.897 15 1.80 -0.627 7.462 0.09182 .
trt2 - trt5 4.522 1.897 15 2.38 0.478 8.567 0.03080 *
trt2 - trt6 0.808 1.897 15 0.43 -3.237 4.852 0.67646
trt3 - trt4 6.155 1.897 15 3.24 2.111 10.199 0.00545 **
trt3 - trt5 7.260 1.897 15 3.83 3.216 11.304 0.00165 **
trt3 - trt6 3.545 1.897 15 1.87 -0.499 7.589 0.08138 .
trt4 - trt5 1.105 1.897 15 0.58 -2.939 5.149 0.56897
trt4 - trt6 -2.610 1.897 15 -1.38 -6.654 1.434 0.18916
trt5 - trt6 -3.715 1.897 15 -1.96 -7.759 0.329 0.06911 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Confidence level: 95%
Degrees of freedom method: Satterthwaite Multiple comparisons in the RCBD
Tukey’s for comparing all pairs of means among \(\mu_1,\dots,\mu_a\): \[ \bar Y_{i.} - \bar Y_{i'.} \pm q_{a,(a-1)(b-1),\alpha}\sqrt{\operatorname{MS}_{\operatorname{Error}}} \sqrt{\frac{1}{b}}, \quad 1 \leq i < i' \leq a. \]
Dunnett’s for comparing \(\mu_2,\dots,\mu_a\) to a baseline \(\mu_1\): \[ \bar Y_{i.} - \bar Y_{1.} \pm d_{a,(a-1)(b-1),\alpha}\sqrt{\operatorname{MS}_{\operatorname{Error}}}\sqrt{\frac{2}{b}}, \quad i=2,\dots,a. \]
Nitrogen fertilization
Compare all pairs of fertilizers with Tukey’s CIs for mean differences.
alpha <- 0.05
a <- 6
b <- 4
se <- sqrt(MSE) * sqrt(2/b) # have compute MSE previously
me <- qtukey(1-alpha,a,(a-1)*(b-1)) / sqrt(2) * se
ntr_means <- aggregate(ntr, by = list(trt), mean)$x
CIs <- matrix(NA,choose(a,2),2)
comp <- numeric(choose(a,2))
k <- 1
for(i in 1:(a-1))
for(j in (i+1):a){ # double loop takes us through all pairs
dij <- ntr_means[i] - ntr_means[j]
CIs[k,] <- c(dij - me, dij + me)
comp[k] <- paste(i,"-",j)
k <- k + 1
}
colnames(CIs) <- c("lower","upper")
rownames(CIs) <- compNitrogen fertilization
Nitrogen fertilization
Compare all fertilizers to fertilizer 1 using Dunnett’s method.
alpha <- 0.05
a <- 6
b <- 4
me <- 2.82 * sqrt(MSE) * sqrt(2/b) # value 2.82 from Dunnett's table
ntr_means <- aggregate(ntr, by = list(trt), mean)$x
CIs <- matrix(NA,a-1,2)
k <- 1
for(i in 2:a){
di <- ntr_means[i] - ntr_means[1]
CIs[k,] <- c(di - me,di + me)
k <- k + 1
}
colnames(CIs) <- c("lower","upper")
rownames(CIs) <- paste(2:a,"- 1")Nitrogen fertilization
Nitrogen fertilization
Nitrogen fertilization
Ignoring the blocks in the nitrogen data
If we ignore the blocks, the design looks like a one-way ANOVA:
| Treatment | 1 | 2 | 3 | 4 | 5 | 6 |
|---|---|---|---|---|---|---|
| 34.98 | 40.89 | 42.07 | 37.18 | 37.99 | 34.89 | |
| 41.22 | 46.69 | 49.42 | 45.85 | 41.99 | 50.15 | |
| 36.94 | 46.65 | 52.68 | 40.23 | 37.61 | 44.57 | |
| 39.97 | 41.90 | 43.29 | 39.20 | 40.45 | 43.29 |
Suppose we fit \(Y_{ij} = \mu + \tau_i + \varepsilon_{ij}\), \(i = 1,2,3,4,5,6\), \(j = 1,2,3,4\).
Nitrogen fertilization
We lose power to detect a treatment effect!
Cabbage counts versus nitrogen, Kuehl (2000)
Heads of cabbage in subplots of two field plots.
Subplots in each plot randomly assigned to nitrogen levels.
Cabbage counts
Cabbage counts
Treatment effects model for RCBD with replication
Assume \[ Y_{ijk} + \mu + \tau_i + B_j + (\tau B)_{ij} + \varepsilon_{ijk} \] for \(i=1,\dots,a\), \(j=1,\dots,b\), and \(k = 1,\dots,n_{ij}\), where
- \(Y_{ijk}\) is response of EU \(k\) in block \(j\) receiving treatment \(i\).
- the \(\tau_i\) are the fixed effects of the treatment.
- the \(B_j\) are independent \(\mathcal{N}(0,\sigma_B^2)\) random block effects.
- the \((\tau B)_{ij}\) are indep. \(\mathcal{N}(0,\sigma_{AB}^2)\) random interaction effects.
- the \(\varepsilon_{ijk}\) are independent \(\mathcal{N}(0,\sigma_\varepsilon^2)\) error terms.
- \(\mu\) is an overall or baseline mean.
Define the cell means as \(\mu_i = \mu + \tau_i, \quad i = 1,\dots,a.\)
Assume for now a balanced design: \(n_{ij} = n\) for all \(i,j\).
Sums of squares for RCBD with replication
| Sum of squares | Symbol | Formula |
|---|---|---|
| Total | \(\operatorname{SS}_{\operatorname{Tot}}\) | \(\sum_{i=1}^a\sum_{j=1}^b \sum_{k=1}^n (Y_{ijk} - \bar Y_{...})^2\) |
| Treatment | \(\operatorname{SS}_{\operatorname{A}}\) | \(nb\sum_{i=1}^a (\bar Y_{i..} - \bar Y_{...})^2\) |
| Block | \(\operatorname{SS}_{\operatorname{B}}\) | \(na\sum_{j=1}^b (\bar Y_{.j.} - \bar Y_{...})^2\) |
| Block \(\times\) Treatment | \(\operatorname{SS}_{\operatorname{AB}}\) | \(n\sum_{i=1}^a \sum_{j=1}^b (\bar Y_{ij.} - (\bar Y_{i..} + \bar Y_{.j.} - \bar Y_{...}))^2\) |
| Error | \(\operatorname{SS}_{\operatorname{Error}}\) | \(\sum_{i=1}^a \sum_{j=1}^b \sum_{k=1}^n(Y_{ijk} - \bar Y_{ij.})^2\) |
- Can make the decomposition \(\operatorname{SS}_{\operatorname{Tot}}= \operatorname{SS}_{\operatorname{A}}+ \operatorname{SS}_{\operatorname{B}}+ \operatorname{SS}_{\operatorname{AB}}+ \operatorname{SS}_{\operatorname{Error}}\).
ANOVA table for RCBD with replication
| Source | Df | SS | MS | F value |
|---|---|---|---|---|
| A | \(a-1\) | \(\operatorname{SS}_{\operatorname{A}}\) | \(\operatorname{MS}_{\operatorname{A}}\) | \(F_{\operatorname{A}} = \operatorname{MS}_{\operatorname{A}}/ \textcolor{red}{\operatorname{MS}_{\operatorname{AB}}}\) |
| B | \(b-1\) | \(\operatorname{SS}_{\operatorname{B}}\) | \(\operatorname{MS}_{\operatorname{B}}\) | \(F_{\operatorname{B}} = \operatorname{MS}_{\operatorname{B}}/ \textcolor{red}{\operatorname{MS}_{\operatorname{AB}}}\) |
| AB | \((a-1)(b-1)\) | \(\operatorname{SS}_{\operatorname{AB}}\) | \(\operatorname{MS}_{\operatorname{AB}}\) | \(F_{\operatorname{AB}} = \operatorname{MS}_{\operatorname{AB}}/ \operatorname{MS}_{\operatorname{Error}}\) |
| Error | \(ab(n-1)\) | \(\operatorname{SS}_{\operatorname{Error}}\) | \(\operatorname{MS}_{\operatorname{Error}}\) | |
| Total | \(abn-1\) | \(\operatorname{SS}_{\operatorname{Tot}}\) |
- Reject \(H_0\): \(\mu_1 = \dots = \mu_a\) if \(F_{\operatorname{A}} > F_{a-1,\textcolor{red}{(a-1)(b-1)},\alpha}\).
- Reject \(H_0\): \(\sigma_B^2 = 0\) if \(F_{\operatorname{B}} > F_{b-1,\textcolor{red}{(a-1)(b-1)},\alpha}\).
- Reject \(H_0\): \(\sigma_{AB}^2 = 0\) if \(F_{\operatorname{AB}} > F_{(a-1)(b-1),ab(n-1),\alpha}\).
Expected MS in RCBD with replication
| Source | Df | Expected mean square |
|---|---|---|
| A | \(a-1\) | \(nb\theta_A^2 + n \sigma_{AB}^2 + \sigma^2_\varepsilon\) |
| B | \(b-1\) | \(na\sigma_B^2 + n \sigma_{AB}^2 + \sigma^2_\varepsilon\) |
| AB | \((a-1)(b-1)\) | \(n \sigma_{AB}^2 + \sigma^2_\varepsilon\) |
| Error | \(ab(n-1)\) | \(\sigma_\varepsilon^2\) |
In the above \(\theta_A^2 = (a-1)^{-1}\sum_{i=1}^a(\mu_i - \bar\mu_.)^2\).
Check expected MS values under each \(H_0\) on previous slide.
Cabbage counts
a <- nlevels(nitr)
b <- nlevels(blk)
n <- 2
y... <- predict(lm(y ~ 1))
yi.. <- predict(lm(y ~ nitr))
y.j. <- predict(lm(y ~ blk))
yij. <- predict(lm(y ~ nitr + blk + nitr:blk))
SST <- sum((y - y...)**2)
SSA <- sum((yi.. - y...)**2)
SSB <- sum((y.j. - y...)**2)
SSAB <- sum((yij. - (yi.. + y.j. - y...))**2)
SSE <- sum((y - yij.)**2)
MSA <- SSA / (a-1)
MSB <- SSB / (b-1)
MSAB <- SSAB / ((a-1)*(b-1))
MSE <- SSE / (a*b*(n-1))
FA <- MSA / MSAB
FB <- MSB / MSAB
FAB <- MSAB / MSE
pA <- 1 - pf(FA,a-1,(a-1)*(b-1))
pB <- 1 - pf(FB,b-1,(a-1)*(b-1))
pAB <- 1 - pf(FAB,(a-1)*(b-1),a*b*(n-1))Cabbage counts
| Source | Df | SS | MS | F value | p value |
|---|---|---|---|---|---|
| A | 4 | 4813.00 | 1203.25 | 16.7234 | 0.0092 |
| B | 1 | 1022.45 | 1022.45 | 14.2106 | 0.0196 |
| AB | 4 | 287.80 | 71.95 | 1.7030 | 0.2253 |
| Error | 10 | 422.50 | 42.25 | ||
| Total | 19 | 6545.75 |
Cabbage counts
Note that anova() on lm() output does not give desired output:
Method of moments for estimating \(\sigma_B\), \(\sigma_{AB}^2\), and \(\sigma_\varepsilon^2\)
Set the mean squares equal to their expectations and solve: \[ \begin{aligned} \operatorname{MS}_{\operatorname{B}}&\overset{\text{set}}{=} na\sigma_B^2 + n \sigma_{AB}^2 + \sigma_\varepsilon^2\\ \operatorname{MS}_{\operatorname{AB}}&\overset{\text{set}}{=} n \sigma_{AB}^2 + \sigma_\varepsilon^2\\ \operatorname{MS}_{\operatorname{Error}}&\overset{\text{set}}{=} \sigma_\varepsilon^2 \end{aligned} \] Solve for \(\sigma_B^2\), \(\sigma_{AB}^2\), and \(\sigma_\varepsilon^2\): \[ \begin{aligned} \dot \sigma_B^2 & =\frac{1}{na}(\operatorname{MS}_{\operatorname{B}}- \operatorname{MS}_{\operatorname{AB}})\\ \dot \sigma_{AB}^2 &=\frac{1}{n}(\operatorname{MS}_{\operatorname{AB}}- \operatorname{MS}_{\operatorname{Error}})\\ \dot \sigma_\varepsilon^2 &=\operatorname{MS}_{\operatorname{Error}} \end{aligned} \]
Possible to obtain negative values for \(\dot \sigma_B^2\) and \(\dot \sigma_{AB}^2\). Then use REML.
Cabbage counts
Obtain the method of moments estimates for the variance components.
We obtain \(\dot \sigma_B = 9.749\), \(\dot \sigma_{AB} = 3.854\), and \(\dot \sigma_\varepsilon = 6.500\).
Cabbage counts
Obtain the REML estimates for the variance components.
Linear mixed model fit by REML ['lmerModLmerTest']
Formula: y ~ nitr + (1 | blk) + (1 | nitr:blk)
REML criterion at convergence: 111
Random effects:
Groups Name Std.Dev.
nitr:blk (Intercept) 3.85
blk (Intercept) 9.75
Residual 6.50
Number of obs: 20, groups: nitr:blk, 10; blk, 2
Fixed Effects:
(Intercept) nitr50 nitr100 nitr150 nitr200
112.8 32.7 36.2 44.7 36.2 CIs (unadjusted) for diffs in means in RCBD with rep.
A \((1-\alpha)100\%\) confidence interval for \(\mu_i - \mu_{i'}\) is \[ \bar Y_{i..} - \bar Y_{i'..} \pm t_{(a-1)(b-1),\alpha/2}\sqrt{\operatorname{MS}_{\operatorname{AB}}}\sqrt{\frac{2}{nb}} \] for \(1 \leq i < i' \leq a\).
If \((a-1)(b-1) = 1\), this interval will be frightfully wide.
In this case one should drop the interaction term from the model!
CI for difference in means for cabbage count data
Compare the means of the nitrogen level \(50\) and \(0\) groups:
CIs (unadjusted) for parameters in RCBD with replication
Use confint() on output of lmer():
2.5 % 97.5 %
.sig01 0.00 6.25
.sig02 2.72 29.86
.sigma 4.42 9.06
(Intercept) 95.10 130.40
nitr50 23.25 42.25
nitr100 26.75 45.75
nitr150 35.25 54.25
nitr200 26.75 45.75These CIs are not adjusted to achieve a familywise coverage probability!
RCBD with replication without an interaction term
If we omit the interaction term in the RCBD with replication, we have
\[ Y_{ijk} + \mu + \tau_i + B_j + \varepsilon_{ijk} \] for \(i=1,\dots,a\), \(j=1,\dots,b\), and \(k = 1,\dots,n_{ij}\), where
- \(Y_{ijk}\) is response of EU \(k\) in block \(j\) receiving treatment \(i\).
- the \(\tau_i\) are the fixed effects of the treatment.
- the \(B_j\) are independent \(\mathcal{N}(0,\sigma_B^2)\) random block effects.
- the \(\varepsilon_{ijk}\) are independent \(\mathcal{N}(0,\sigma_\varepsilon^2)\) error terms.
- \(\mu\) is an overall or baseline mean.
Assume for now a balanced design: \(n_{ij} = n\) for all \(i,j\).
Use this model if \(a = b = 2\), under which \((a-1)(b-1)=1\).
Sums of squares for RCBD with rep., no interaction
| SS | Symbol | Formula |
|---|---|---|
| Total | \(\operatorname{SS}_{\operatorname{Tot}}\) | \(\sum_{i=1}^a\sum_{j=1}^b \sum_{k=1}^n (Y_{ijk} - \bar Y_{...})^2\) |
| Treatment | \(\operatorname{SS}_{\operatorname{A}}\) | \(nb\sum_{i=1}^a (\bar Y_{i..} - \bar Y_{...})^2\) |
| Block | \(\operatorname{SS}_{\operatorname{B}}\) | \(na\sum_{j=1}^b (\bar Y_{.j.} - \bar Y_{...})^2\) |
| Error | \(\operatorname{SS}_{\operatorname{Error}}\) | \(\sum_{i=1}^a \sum_{j=1}^b \sum_{k=1}^n(Y_{ijk} - (\bar Y_{i..} + \bar Y_{.j.} - \bar Y_{...}))^2\) |
- Can make the decomposition \(\operatorname{SS}_{\operatorname{Tot}}= \operatorname{SS}_{\operatorname{A}}+ \operatorname{SS}_{\operatorname{B}}+ \operatorname{SS}_{\operatorname{Error}}\).
ANOVA table for RCBD with rep., omitting interaction
| Source | Df | SS | MS | F value |
|---|---|---|---|---|
| A | \(a-1\) | \(\operatorname{SS}_{\operatorname{A}}\) | \(\operatorname{MS}_{\operatorname{A}}\) | \(F_{\operatorname{A}} = \operatorname{MS}_{\operatorname{A}}/ \operatorname{MS}_{\operatorname{Error}}\) |
| B | \(b-1\) | \(\operatorname{SS}_{\operatorname{B}}\) | \(\operatorname{MS}_{\operatorname{B}}\) | \(F_{\operatorname{B}} = \operatorname{MS}_{\operatorname{B}}/ \operatorname{MS}_{\operatorname{Error}}\) |
| Error | \(\textcolor{red}{abn-a - b + 1}\) | \(\operatorname{SS}_{\operatorname{Error}}\) | \(\operatorname{MS}_{\operatorname{Error}}\) | |
| Total | \(abn-1\) | \(\operatorname{SS}_{\operatorname{Tot}}\) |
- Reject \(H_0\): \(\mu_1 = \dots = \mu_a\) if \(F_{\operatorname{A}} > F_{a-1,\textcolor{red}{abn-a - b + 1},\alpha}\).
- Reject \(H_0\): \(\sigma_B^2 = 0\) if \(F_{\operatorname{B}} > F_{b-1,\textcolor{red}{abn-a - b + 1},\alpha}\).
Expected MS in RCBD with rep., omitting interaction
| Source | Df | Expected mean square |
|---|---|---|
| A | \(a-1\) | \(nb\theta_A^2 + \sigma^2_\varepsilon\) |
| B | \(b-1\) | \(na\sigma_B^2 + \sigma^2_\varepsilon\) |
| Error | \(abn - a - b + 1\) | \(\sigma_\varepsilon^2\) |
In the above \(\theta_A^2 = (a-1)^{-1}\sum_{i=1}^a(\mu_i - \bar\mu_.)^2\).
Cabbage counts
ANOVA table, omitting the interaction:
a <- nlevels(nitr)
b <- nlevels(blk)
n <- 2
y... <- predict(lm(y ~ 1))
yi.. <- predict(lm(y ~ nitr))
y.j. <- predict(lm(y ~ blk))
SST <- sum((y - y...)**2)
SSA <- sum((yi.. - y...)**2)
SSB <- sum((y.j. - y...)**2)
SSE <- sum((y - (yi.. + y.j. - y...))**2)
MSA <- SSA / (a-1)
MSB <- SSB / (b-1)
MSE <- SSE / (a*b*n - a - b + 1)
FA <- MSA / MSE
FB <- MSB / MSE
pA <- 1 - pf(FA,a-1,a*b*n - a - b + 1)
pB <- 1 - pf(FB,b-1,a*b*n - a - b + 1)Cabbage counts
| Source | Df | SS | MS | F value | p value |
|---|---|---|---|---|---|
| A | 4 | 4813.00 | 1203.25 | 23.7160 | 0.0000 |
| B | 1 | 1022.45 | 1022.45 | 20.1525 | 0.0005 |
| Error | 14 | 710.30 | 50.74 | ||
| Total | 19 | 6545.75 |
Cabbage counts
With the interaction term omitting, anova() on lm() output is correct:
MoM for variance components when omitting interaction
Set the mean squares equal to their expectations and solve: \[ \begin{aligned} \operatorname{MS}_{\operatorname{B}}&\overset{\text{set}}{=} na\sigma_B^2 + \sigma_\varepsilon^2\\ \operatorname{MS}_{\operatorname{Error}}&\overset{\text{set}}{=} \sigma_\varepsilon^2 \end{aligned} \] Solve for \(\sigma_B^2\) and \(\sigma_\varepsilon^2\): \[ \begin{aligned} \dot \sigma_B^2 & =\frac{1}{na}(\operatorname{MS}_{\operatorname{B}}- \operatorname{MS}_{\operatorname{Error}})\\ \dot \sigma_\varepsilon^2 &=\operatorname{MS}_{\operatorname{Error}} \end{aligned} \]
Possible to obtain a negative value for \(\dot \sigma_B^2\). In this case use REML.
Cabbage counts
MoM estimates when omitting interaction:
We obtain \(\dot \sigma_B = 9.858\) and \(\dot \sigma_\varepsilon = 7.123\).
Cabbage counts
REML estimates when omitting interaction:
Linear mixed model fit by REML ['lmerModLmerTest']
Formula: y ~ nitr + (1 | blk)
REML criterion at convergence: 111
Random effects:
Groups Name Std.Dev.
blk (Intercept) 9.86
Residual 7.12
Number of obs: 20, groups: blk, 2
Fixed Effects:
(Intercept) nitr50 nitr100 nitr150 nitr200
112.8 32.8 36.3 44.8 36.3 CIs (unadjusted) for diffs in means in RCBD with rep. without interaction
- A \((1-\alpha)100\%\) confidence interval for \(\mu_i - \mu_{i'}\) is \[ \bar Y_{i..} - \bar Y_{i'..} \pm t_{abn - a - b + 1,\alpha/2}\sqrt{\operatorname{MS}_{\operatorname{Error}}}\sqrt{\frac{2}{nb}} \] for \(1 \leq i < i' \leq a\).
Cabbage counts
Compare the means of the nitrogen level \(50\) and \(0\) groups:
Cabbage counts
Can use ls_means() as below:
Least Squares Means table:
Estimate Std. Error df t value lower upper Pr(>|t|)
nitr0 - nitr50 -32.75 5.04 14 -6.50 -43.55 -21.95 1.4e-05 ***
nitr0 - nitr100 -36.25 5.04 14 -7.20 -47.05 -25.45 4.6e-06 ***
nitr0 - nitr150 -44.75 5.04 14 -8.88 -55.55 -33.95 4.0e-07 ***
nitr0 - nitr200 -36.25 5.04 14 -7.20 -47.05 -25.45 4.6e-06 ***
nitr50 - nitr100 -3.50 5.04 14 -0.69 -14.30 7.30 0.498
nitr50 - nitr150 -12.00 5.04 14 -2.38 -22.80 -1.20 0.032 *
nitr50 - nitr200 -3.50 5.04 14 -0.69 -14.30 7.30 0.498
nitr100 - nitr150 -8.50 5.04 14 -1.69 -19.30 2.30 0.114
nitr100 - nitr200 0.00 5.04 14 0.00 -10.80 10.80 1.000
nitr150 - nitr200 8.50 5.04 14 1.69 -2.30 19.30 0.114
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Confidence level: 95%
Degrees of freedom method: Satterthwaite Skin response data, Mohr, Wilson, and Freund (2021)
Galvanic skin responses of five subjects under shock and noise stimuli.
Skin response
resp noise shock subj
1 3 40 0.25 1
2 7 40 0.25 2
3 9 40 0.25 3
4 4 40 0.25 4
5 1 40 0.25 5
6 5 40 0.5 1
7 11 40 0.5 2
8 13 40 0.5 3
9 8 40 0.5 4
10 3 40 0.5 5
11 9 40 0.75 1
12 12 40 0.75 2
13 14 40 0.75 3
14 11 40 0.75 4
15 5 40 0.75 5
16 6 40 1 1
17 11 40 1 2
18 12 40 1 3
19 7 40 1 4
20 4 40 1 5Treatment effects model for two-way factorial RCBD
Assume \[ Y_{ijk} = \mu + \tau_i + \gamma_j + (\tau\gamma)_{ij} + C_k + \varepsilon_{ijk}, \] for \(i = 1,\dots,a\), \(j = 1,\dots,b\), and \(k = 1,\dots,c\), where
- \(Y_{ijk}\) is the response in block \(k\) under treatment combination \(i\times j\).
- \(\mu\) is an overall or baseline mean.
- the \(\tau_i\) are fixed effects for factor A.
- the \(\gamma_j\) are fixed effects for factor B.
- the \((\tau\gamma)_{ij}\) are effects for the AB interaction.
- the \(C_k\) are independent \(\mathcal{N}(0,\sigma_C^2)\) block effects.
- the \(\varepsilon_{ijk}\) are independent \(\mathcal{N}(0,\sigma_\varepsilon^2)\) error terms.
Define the cell means as \[ \mu_{ij} = \mu + \tau_i + \gamma_j + (\tau\gamma)_{ij}, \quad i = 1,\dots,a, \quad j = 1,\dots,b. \]
Sums of squares for the two-way factorial RCBD
| Sum of squares | Symbol | Formula |
|---|---|---|
| Total | \(\operatorname{SS}_{\operatorname{Tot}}\) | \(\sum_{i=1}^a\sum_{j=1}^b\sum_{k=1}^c(Y_{ijk} - \bar Y_{...})^2\) |
| A | \(\operatorname{SS}_{\operatorname{A}}\) | \(bc\sum_{i=1}^a (\bar Y_{i..} - \bar Y_{...})^2\) |
| B | \(\operatorname{SS}_{\operatorname{B}}\) | \(ac\sum_{j=1}^b (\bar Y_{.j.} - \bar Y_{...})^2\) |
| AB | \(\operatorname{SS}_{\operatorname{AB}}\) | \(c\sum_{i=1}^a\sum_{j=1}^b (Y_{ij.} - (\bar Y_{i..} + \bar Y_{.j.} - \bar Y_{...}))^2\) |
| C | \(\operatorname{SS}_{\operatorname{C}}\) | \(ab\sum_{k=1}^c (\bar Y_{..k} - \bar Y_{...})^2\) |
| Error | \(\operatorname{SS}_{\operatorname{Error}}\) | \(\operatorname{SS}_{\operatorname{Tot}}- (\operatorname{SS}_{\operatorname{A}}+ \operatorname{SS}_{\operatorname{B}}+ \operatorname{SS}_{\operatorname{AB}}+ \operatorname{SS}_{\operatorname{C}})\) |
- Then we have \(\operatorname{SS}_{\operatorname{Tot}}= \operatorname{SS}_{\operatorname{A}}+ \operatorname{SS}_{\operatorname{B}}+ \operatorname{SS}_{\operatorname{AB}}+ \operatorname{SS}_{\operatorname{C}}+ \operatorname{SS}_{\operatorname{Error}}\).
- The error \(\operatorname{SS}_{\operatorname{Error}}\) is really the interaction sum of squares \(\operatorname{SS}_{\operatorname{Trt}\times C}\).
- Again, without replication, we cannot estimate this interaction.
- So the interaction serves as the error term.
ANOVA table for two-way factorial RCBD
| Source | Df | SS | MS | F value |
|---|---|---|---|---|
| A | \(a-1\) | \(\operatorname{SS}_{\operatorname{A}}\) | \(\operatorname{MS}_{\operatorname{A}}\) | \(F_{\operatorname{A}}= \operatorname{MS}_{\operatorname{A}}/ \operatorname{MS}_{\operatorname{Error}}\) |
| B | \(b-1\) | \(\operatorname{SS}_{\operatorname{B}}\) | \(\operatorname{MS}_{\operatorname{B}}\) | \(F_{\operatorname{B}}= \operatorname{MS}_{\operatorname{B}}/ \operatorname{MS}_{\operatorname{Error}}\) |
| AB | \((a-1)(b-1)\) | \(\operatorname{SS}_{\operatorname{AB}}\) | \(\operatorname{MS}_{\operatorname{AB}}\) | \(F_{\operatorname{AB}}= \operatorname{MS}_{\operatorname{B}}/ \operatorname{MS}_{\operatorname{Error}}\) |
| C | \(c-1\) | \(\operatorname{SS}_{\operatorname{C}}\) | \(\operatorname{MS}_{\operatorname{C}}\) | \(F_{\operatorname{C}}= \operatorname{MS}_{\operatorname{C}}/ \operatorname{MS}_{\operatorname{Error}}\) |
| Error | \((ab-1)(c-1)\) | \(\operatorname{SS}_{\operatorname{Error}}\) | \(\operatorname{MS}_{\operatorname{Error}}\) | |
| Total | \(abc-1\) | \(\operatorname{SS}_{\operatorname{Tot}}\) |
- Reject \(H_0\): \(\mu_{1.} = \dots = \mu_{a.}\) if \(F_{\operatorname{A}}> F_{a-1,(ab - 1)(c-1),\alpha}\).
- Reject \(H_0\): \(\mu_{.1} = \dots = \mu_{.b}\) if \(F_{\operatorname{B}}> F_{b-1,(ab - 1)(c-1),\alpha}\).
- R. \(H_0\): \(\mu_{ij} = \bar \mu_{i.} + \bar \mu_{.j} - \bar \mu_{..}~ \forall ij\) if \(F_{\operatorname{AB}}> F_{(a-1)(b-1),(ab - 1)(c-1),\alpha}\).
- Reject \(H_0\): \(\sigma_C^2 = 0\) if \(F_{\operatorname{C}}> F_{c-1,(ab - 1)(c-1),\alpha}\).
Skin response
Analysis of Variance Table
Response: resp
Df Sum Sq Mean Sq F value Pr(>F)
noise 1 65 65.0 28.38 1.1e-05 ***
shock 3 219 73.1 31.90 3.6e-09 ***
subj 4 362 90.5 39.48 4.0e-11 ***
noise:shock 3 13 4.2 1.84 0.16
Residuals 28 64 2.3
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Skin response
y <- skin$resp
y... <- predict(lm(resp ~ 1,data = skin))
yi.. <- predict(lm(resp ~ noise,data=skin))
y.j. <- predict(lm(resp ~ shock,data=skin))
y..k <- predict(lm(resp ~ subj,data=skin))
yij. <- predict(lm(resp ~ noise + shock + noise:shock,data=skin))
SSA <- sum((yi.. - y...)^2)
SSB <- sum((y.j. - y...)^2)
SSC <- sum((y..k - y...)^2)
SSAB <- sum((yij. - (yi.. + y.j. - y...))^2)
SST <- sum((y - y...)^2)
SSE <- SST - (SSA + SSB + SSC + SSAB)
a <- 2
b <- 4
c <- 5
MSA <- SSA/(a-1)
MSB <- SSB/(b-1)
MSC <- SSC/(c-1)
MSAB <- SSAB/((a-1)*(b-1))
MSE <- SSE/((a*b-1)*(c-1))
FA <- MSA / MSE
FB <- MSB / MSE
FAB <- MSAB / MSE
FC <- MSC / MSE
pA <- 1 - pf(FA,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))
pC <- 1 - pf(FC,c-1,(a*b-1)*(c-1))Skin response
Expected mean squares in factorial RCBD
| Source | Df | Expected mean square |
|---|---|---|
| A | \(a-1\) | \(bc \theta_A^2 + \sigma_\varepsilon^2\) |
| B | \(b-1\) | \(ac \theta_B^2 + \sigma_\varepsilon^2\) |
| AB | \((a-1)(b-1)\) | \(c \theta_{AB}^2 + \sigma_\varepsilon^2\) |
| C | \(c-1\) | \(ab\sigma_C^2 + \sigma_\varepsilon^2\) |
| Error | \((ab-1)(c-1)\) | \(\sigma_\varepsilon^2\) |
In the above
- \(\theta_A^2 = (a-1)^{-1}\sum_{i=1}^a(\bar \mu_{i.} - \bar \mu_{..})^2\)
- \(\theta_B^2 = (b-1)^{-1}\sum_{j=1}^b(\bar \mu_{.j} - \bar \mu_{..})^2\)
- \(\theta_{AB}^2 = [(a-1)(b-1)]^{-1}\sum_{i=1}^a\sum_{j=1}^b(\mu_{ij} - (\bar \mu_{i.} + \bar \mu_{.j} - \bar \mu_{..}))^2\)
MoMs for variance components in two-way factorial RCBD
- Equating \(\operatorname{MS}_{\operatorname{C}}\) and \(\operatorname{MS}_{\operatorname{Error}}\) with their expectations gives
\[ \dot \sigma_C^2 = \frac{\operatorname{MS}_{\operatorname{C}}- \operatorname{MS}_{\operatorname{Error}}}{ab} \quad \text{ and }\quad \displaystyle \dot \sigma_{\varepsilon}^2 = \operatorname{MS}_{\operatorname{Error}}. \]
- May obtain \(\dot \sigma_C^2 < 0\), so one should use REML estimation.
Skin response
Obtain REML estimators of \(\sigma_C^2\) and \(\sigma_{\varepsilon}^2\) on the skin response data.
Linear mixed model fit by REML ['lmerModLmerTest']
Formula: resp ~ noise + shock + noise:shock + (1 | subj)
Data: skin
REML criterion at convergence: 145
Random effects:
Groups Name Std.Dev.
subj (Intercept) 3.32
Residual 1.51
Number of obs: 40, groups: subj, 5
Fixed Effects:
(Intercept) noise80 shock0.5 shock0.75
4.8 2.0 3.2 5.4
shock1 noise80:shock0.5 noise80:shock0.75 noise80:shock1
3.2 -0.6 2.4 0.4 Skin response
Obtain MoM estimates for \(\sigma_C^2\) and \(\sigma_{\varepsilon}^2\) on the skin response data.
We have \(\dot \sigma_C = 3.320\) and \(\dot \sigma_{\varepsilon}= 1.514\).
Skin response
Skin response
