library(GAD)
library(dplyr)

Data entry

chemical1 <- c(73, 68, 74, 71, 67)
chemical2 <- c(73, 67, 75, 72, 70)
chemical3 <- c(75, 68, 78, 73, 68)
chemical4 <- c(73, 71, 75, 75, 69)

cloth <- data.frame(
  Chemical = factor(rep(1:4, each = 5)),
  Bolt = factor(rep(1:5, times = 4)),
  Strength = c(
    chemical1,
    chemical2,
    chemical3,
    chemical4
  )
)

cloth
##    Chemical Bolt Strength
## 1         1    1       73
## 2         1    2       68
## 3         1    3       74
## 4         1    4       71
## 5         1    5       67
## 6         2    1       73
## 7         2    2       67
## 8         2    3       75
## 9         2    4       72
## 10        2    5       70
## 11        3    1       75
## 12        3    2       68
## 13        3    3       78
## 14        3    4       73
## 15        3    5       68
## 16        4    1       73
## 17        4    2       71
## 18        4    3       75
## 19        4    4       75
## 20        4    5       69

Descriptive statistics

chemical_summary <- cloth %>%
  group_by(Chemical) %>%
  summarise(
    n = n(),
    Mean = mean(Strength),
    SD = sd(Strength),
    .groups = "drop"
  )

bolt_summary <- cloth %>%
  group_by(Bolt) %>%
  summarise(
    Mean = mean(Strength),
    .groups = "drop"
  )

knitr::kable(
  chemical_summary,
  digits = 3,
  caption = "Tensile Strength by Chemical"
)
Tensile Strength by Chemical
Chemical n Mean SD
1 5 70.6 3.050
2 5 71.4 3.050
3 5 72.4 4.393
4 5 72.6 2.608 
knitr::kable(
  bolt_summary,
  digits = 3,
  caption = "Mean Tensile Strength by Bolt"
)
Mean Tensile Strength by Bolt
Bolt Mean
1 73.50
2 68.50
3 75.50
4 72.75
5 68.50

The chemical means are:

The bolt means range from 68.5 to 75.5. This suggests that there may be substantial bolt-to-bolt variability.

Question 1: Randomized complete block design

Experimental structure

Chemical is a fixed effect because the chemist wants to compare these four specific chemical agents.

Bolt is a random effect because the five bolts represent a sample from a larger population of cloth bolts. Bolt is included as a block to control nuisance variability.

There is one observation for every Chemical–Bolt combination. Therefore, the Chemical-by-Bolt interaction cannot be estimated separately and is included in the error term under the additive RCBD model.

Linear effects model

Let \(Y_{ij}\) be the tensile strength for chemical \(i\) on bolt \(j\). The randomized complete block model is

\[ Y_{ij}=\mu+\alpha_i+\beta_j+\varepsilon_{ij}, \]

where:

  • \(\mu\) is the overall mean;
  • \(\alpha_i\) is the fixed effect of chemical \(i\);
  • \(i=1,2,3,4\);
  • \(\beta_j\) is the random effect of bolt \(j\);
  • \(j=1,2,\ldots,5\);
  • \(\beta_j\overset{iid}{\sim}N(0,\sigma_\beta^2)\);
  • \(\varepsilon_{ij}\overset{iid}{\sim}N(0,\sigma^2)\).

For the fixed chemical effects, use the constraint

\[ \sum_{i=1}^{4}\alpha_i=0. \]

The random bolt effects and experimental errors are assumed to be independent.

Chemical hypotheses

The primary hypotheses are

\[ H_0:\mu_1=\mu_2=\mu_3=\mu_4 \]

against

\[ H_a:\text{at least one chemical mean differs}. \]

Equivalently,

\[ H_0:\alpha_1=\alpha_2=\alpha_3=\alpha_4=0. \]

Bolt hypotheses

Because Bolt is random, the hypotheses concern the block variance:

\[ H_0:\sigma_\beta^2=0 \]

against

\[ H_a:\sigma_\beta^2>0. \]

RCBD analysis using GAD

cloth$Chemical_GAD <- as.fixed(cloth$Chemical)
cloth$Bolt_GAD <- as.random(cloth$Bolt)

rcbd_model <- lm(
  Strength ~ Chemical_GAD + Bolt_GAD,
  data = cloth
)

# Check the model structure and expected mean squares
estimates(rcbd_model)
## $tm
##              Chemical_GAD Bolt_GAD   n
## Chemical_GAD            0        5 2.5
## Bolt_GAD                4        1 2.5
## Res                     1        1 1.0
## 
## $mse
##              Mean square estimates         
## Chemical_GAD "1*Residuals + 5*Chemical_GAD"
## Bolt_GAD     "1*Residuals + 4*Bolt_GAD"    
## Residuals    "1*Residuals"                 
## 
## $f.versus
##                 Numerator Denominator
## Chemical_GAD Chemical_GAD   Residuals
## Bolt_GAD         Bolt_GAD   Residuals
# Produce the ANOVA table
rcbd_table <- gad(rcbd_model)
rcbd_table
## $anova
## Analysis of Variance Table
## 
## Response: Strength
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## Chemical_GAD  3  12.95   4.317  2.3761    0.1211    
## Bolt_GAD      4 157.00  39.250 21.6055 2.059e-05 ***
## Residuals    12  21.80   1.817                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

RCBD ANOVA results

The RCBD ANOVA results are approximately:

Source df SS MS F p-value
Chemical 3 12.950 4.317 2.376 0.1211
Bolt 4 157.000 39.250 21.606 0.0000206
Error 12 21.800 1.817
Total 19 191.750

For Chemical, the test statistic is

\[ F_0= \frac{MS_{\text{Chemical}}}{MS_E} = \frac{4.3167}{1.8167} = 2.376. \]

The p-value is

\[ p=0.1211. \]

Because

\[ 0.1211<0.15, \]

we reject \(H_0\) at the specified significance level.

Chemical conclusion: At \(\alpha=0.15\), there is sufficient evidence that the four chemical agents do not all have the same mean tensile strength.

This result would not be significant at the more commonly used \(\alpha=0.05\) level because \(0.1211>0.05\).

The sample means indicate that Chemicals 3 and 4 produced the highest average tensile strengths.

Test for the Bolt effect

For Bolt,

\[ F_0= \frac{MS_{\text{Bolt}}}{MS_E} = \frac{39.25}{1.8167} = 21.606. \]

The p-value is

\[ p=0.0000206. \]

Because

\[ 0.0000206<0.15, \]

we reject

\[ H_0:\sigma_\beta^2=0. \]

Bolt conclusion: Bolt contributes a statistically significant amount of variability to tensile strength.

Variance components

MS_bolt <- 39.25
MS_error_rcbd <- 1.8167
number_chemicals <- 4

sigma2_bolt <- (
  MS_bolt - MS_error_rcbd
) / number_chemicals

sigma2_error <- MS_error_rcbd

bolt_proportion <- sigma2_bolt /
  (sigma2_bolt + sigma2_error)

sigma2_bolt
## [1] 9.358325
sigma2_error
## [1] 1.8167
bolt_proportion
## [1] 0.8374321

The estimated bolt variance is

\[ \widehat{\sigma}_\beta^2 = \frac{MS_{\text{Bolt}}-MS_E}{4}. \]

Therefore,

\[ \widehat{\sigma}_\beta^2 = \frac{39.25-1.8167}{4} = 9.3583. \]

The estimated error variance is

\[ \widehat{\sigma}^2 = MS_E = 1.8167. \]

The proportion of random variability associated with bolts is

\[ \frac{9.3583}{9.3583+1.8167} \times100\% = 83.7\%. \]

Approximately 83.7% of the estimated random variability is associated with differences among bolts.

RCBD model diagnostics

par(mfrow = c(2, 2))
plot(rcbd_model)

par(mfrow = c(1, 1))

Test the residuals for normality:

shapiro.test(residuals(rcbd_model))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(rcbd_model)
## W = 0.8996, p-value = 0.04054

The Shapiro–Wilk test gives approximately

\[ W=0.8996 \]

and

\[ p=0.0405. \]

Because

\[ 0.0405<0.15, \]

the normality assumption is questionable at the assignment’s significance level.

The residual plots should be examined for unequal spread, unusual observations, and patterns. Because this experiment contains only 20 observations, the diagnostic results should be interpreted cautiously.

Question 2: Completely randomized design

If the experiment uses unrelated random pieces of cloth, Bolt is no longer included as a blocking factor. Chemical remains a fixed effect.

Linear effects model

Let \(Y_{ij}\) be the tensile strength of replication \(j\) receiving chemical \(i\). The completely randomized design model is

\[ Y_{ij}=\mu+\alpha_i+\varepsilon_{ij}, \]

where

\[ \varepsilon_{ij}\overset{iid}{\sim}N(0,\sigma^2). \]

Hypotheses

The hypotheses are

\[ H_0:\mu_1=\mu_2=\mu_3=\mu_4 \]

against

\[ H_a:\text{at least one chemical mean differs}. \]

CRD analysis using GAD

cloth$Chemical_CRD <- as.fixed(cloth$Chemical)

crd_model <- lm(
  Strength ~ Chemical_CRD,
  data = cloth
)

# Check the model structure
estimates(crd_model)
## $tm
##              Chemical_CRD n
## Chemical_CRD            0 5
## Res                     1 1
## 
## $mse
##              Mean square estimates         
## Chemical_CRD "1*Residuals + 5*Chemical_CRD"
## Residuals    "1*Residuals"                 
## 
## $f.versus
##                 Numerator Denominator
## Chemical_CRD Chemical_CRD   Residuals
# Produce the ANOVA table
crd_table <- gad(crd_model)
crd_table
## $anova
## Analysis of Variance Table
## 
## Response: Strength
##              Df Sum Sq Mean Sq F value Pr(>F)
## Chemical_CRD  3  12.95  4.3167  0.3863 0.7644
## Residuals    16 178.80 11.1750

CRD ANOVA results

The CRD ANOVA results are approximately:

Source df SS MS F p-value
Chemical 3 12.950 4.317 0.386 0.7644
Error 16 178.800 11.175
Total 19 191.750

The test statistic is

\[ F_0= \frac{MS_{\text{Chemical}}}{MS_E} = \frac{4.3167}{11.175} = 0.386. \]

The p-value is

\[ p=0.7644. \]

Because

\[ 0.7644>0.15, \]

we fail to reject \(H_0\).

CRD conclusion: At \(\alpha=0.15\), there is insufficient evidence that the four chemical agents have different mean tensile strengths.

CRD model diagnostics

par(mfrow = c(2, 2))
plot(crd_model)

par(mfrow = c(1, 1))
shapiro.test(residuals(crd_model))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(crd_model)
## W = 0.93129, p-value = 0.1635

The Shapiro–Wilk test gives approximately

\[ p=0.1635. \]

Because

\[ 0.1635>0.15, \]

there is insufficient evidence that the CRD residuals are nonnormal. However, the result is close to the specified significance level and should be interpreted cautiously.

Question 3: Comparison of the two designs

The RCBD and CRD analyses produce different conclusions.

Design Chemical F p-value Decision at \(\alpha=0.15\)
RCBD 2.376 0.1211 Reject \(H_0\)
CRD 0.386 0.7644 Fail to reject \(H_0\)

In the RCBD, the chemical effect is statistically significant at the specified \(\alpha=0.15\) level.

In the CRD, the chemical effect is not statistically significant.

Why the findings differ

The chemical sum of squares is the same in both analyses:

\[ SS_{\text{Chemical}}=12.95. \]

The chemical sum of squares does not change because both analyses contain the same chemical observations and chemical means.

The difference is in the error term.

For the RCBD:

\[ SS_E=21.80 \]

and

\[ MSE_{\text{RCBD}}=1.8167. \]

For the CRD:

\[ SS_E=178.80 \]

and

\[ MSE_{\text{CRD}}=11.175. \]

The RCBD separates the variation caused by Bolt from the experimental error:

\[ SS_{\text{Bolt}}=157.00. \]

The reduction in error mean square is

\[ \frac{MSE_{\text{CRD}}} {MSE_{\text{RCBD}}} = \frac{11.175}{1.8167} = 6.15. \]

Thus, the CRD error mean square is approximately 6.15 times the RCBD error mean square.

Does Bolt represent significant nuisance variability?

Yes.

The test for Bolt gives

\[ F_{\text{Bolt}}=21.606 \]

and

\[ p=0.0000206<0.15. \]

Therefore, Bolt represents a significant amount of nuisance variability.

Blocking on Bolt was effective because it:

  1. separated bolt-to-bolt variability from experimental error;
  2. reduced the error mean square from 11.175 to 1.8167; and
  3. increased the precision of the chemical comparison.

Final answers

  1. For the RCBD, the chemical effect gives

\[ F=2.376,\qquad p=0.1211. \]

Because \(p<0.15\), reject the hypothesis of equal chemical means at the specified significance level.

  1. The random Bolt effect gives

\[ F=21.606,\qquad p=0.0000206. \]

Therefore, Bolt contributes a significant amount of variability.

  1. For the CRD, the chemical effect gives

\[ F=0.386,\qquad p=0.7644. \]

Fail to reject the hypothesis of equal chemical means.

  1. Blocking reduced the error mean square from

\[ 11.175 \]

to

\[ 1.8167. \]

Therefore, Bolt represents substantial nuisance variability, and the RCBD gives a much more precise comparison of the chemical agents.