Introduction

A chemist wants to test four chemical agents on the strength of a cloth. Because the bolts of cloth can be different from each other, she first runs the experiment as a randomized block design using the bolts as blocks. Then we look at what happens if the same numbers come from a completely randomized design with no blocking.

All the analysis uses the GAD package and alpha = 0.15.

library(GAD)

Question 1: Randomized block design

Data

The data is entered in tidy format: one column for the chemical, one for the bolt and one for the strength.

strength <- c(73, 68, 74, 71, 67,
              73, 67, 75, 72, 70,
              75, 68, 78, 73, 68,
              73, 71, 75, 75, 69)

chemical <- as.fixed(factor(rep(1:4, each = 5)))
bolt <- as.random(factor(rep(1:5, times = 4)))

dat <- data.frame(chemical, bolt, strength)

head(dat)
##   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
str(dat)
## 'data.frame':    20 obs. of  3 variables:
##  $ chemical: Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 2 2 2 2 2 ...
##  $ bolt    : Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5 1 2 3 4 5 ...
##  $ strength: num  73 68 74 71 67 73 67 75 72 70 ...

as.fixed() and as.random() are from GAD. They tell R which factor is fixed and which one is random.

The chemical is fixed because the chemist picked those four agents on purpose and wants to compare those four. The bolt is random because the five bolts were just taken from a bigger supply. We do not care about these five bolts, we care about how much the strength changes from one bolt to another in general.

Mean strength of each chemical and of each bolt:

tapply(dat$strength, dat$chemical, mean)
##    1    2    3    4 
## 70.6 71.4 72.4 72.6
tapply(dat$strength, dat$bolt, mean)
##     1     2     3     4     5 
## 73.50 68.50 75.50 72.75 68.50

Linear effects model

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

where

  • \(y_{ij}\) is the strength of chemical \(i\) on bolt \(j\)
  • \(\mu\) is the overall mean
  • \(\alpha_i\) is the effect of chemical \(i\), a fixed effect, with \(i = 1, 2, 3, 4\)
  • \(\beta_j\) is the effect of bolt \(j\), a random effect, with \(j = 1, 2, 3, 4, 5\) and \(\beta_j \sim N(0, \sigma^2_\beta)\)
  • \(\varepsilon_{ij}\) is the error, with \(\varepsilon_{ij} \sim N(0, \sigma^2)\)

Hypotheses

For the chemical, which is fixed, we test if the effects are zero:

\[H_0: \alpha_1 = \alpha_2 = \alpha_3 = \alpha_4 = 0 \qquad H_a: \text{at least one } \alpha_i \neq 0\]

For the bolt, which is random, we test if the variance of the block effect is zero:

\[H_0: \sigma^2_\beta = 0 \qquad H_a: \sigma^2_\beta > 0\]

Analysis

model1 <- lm(strength ~ chemical + bolt, data = dat)

gad(model1)
## $anova
## Analysis of Variance Table
## 
## Response: strength
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## chemical   3  12.95   4.317  2.3761    0.1211    
## bolt       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

GAD also shows which mean square is used in the bottom of each F ratio:

estimates(model1)
## $tm
##          chemical bolt   n
## chemical        0    5 2.5
## bolt            4    1 2.5
## Res             1    1 1.0
## 
## $mse
##           Mean square estimates     
## chemical  "1*Residuals + 5*chemical"
## bolt      "1*Residuals + 4*bolt"    
## Residuals "1*Residuals"             
## 
## $f.versus
##          Numerator Denominator
## chemical  chemical   Residuals
## bolt          bolt   Residuals

Both F ratios use the residual mean square at the bottom. That is why making the bolt random does not change the numbers in this model, but it does change what the bolt hypothesis means: with a random block we are testing a variance and not four separate effects.

Model check

qqnorm(residuals(model1), main = "Normal Q-Q Plot of the Residuals")
qqline(residuals(model1))

plot(fitted(model1), residuals(model1),
     main = "Residuals vs Fitted Values",
     xlab = "Fitted Value",
     ylab = "Residual",
     pch = 19)
abline(h = 0, lty = 2)

The Q-Q plot is close to the line. There is a flat part in the middle because several residuals have almost the same value, but nothing is far away from the line. In the residuals vs fitted plot the points are spread about the same from left to right, with no funnel. The model looks fine.

Conclusions

For the chemical, F = 2.3761 and the p-value is 0.1211. The p-value is smaller than 0.15, so we reject \(H_0\). There is evidence that the chemical agent changes the strength of the cloth.

For the bolt, F = 21.6055 and the p-value is 0.0000206, much smaller than 0.15, so we reject \(H_0\) there too. The variance between bolts is not zero, so blocking on the bolt was a good idea.

Question 2: Completely randomized design

Now we assume the same 20 numbers came from random pieces of cloth, so there are no blocks. The columns are just replications.

chemical2 <- as.fixed(factor(rep(1:4, each = 5)))

dat2 <- data.frame(chemical2, strength)

head(dat2)
##   chemical2 strength
## 1         1       73
## 2         1       68
## 3         1       74
## 4         1       71
## 5         1       67
## 6         2       73

Linear effects model

\[y_{ij} = \mu + \alpha_i + \varepsilon_{ij}\]

where

  • \(y_{ij}\) is replication \(j\) of chemical \(i\)
  • \(\mu\) is the overall mean
  • \(\alpha_i\) is the effect of chemical \(i\), a fixed effect, with \(i = 1, 2, 3, 4\)
  • \(\varepsilon_{ij}\) is the error, with \(\varepsilon_{ij} \sim N(0, \sigma^2)\) and \(j = 1, 2, 3, 4, 5\)

There is no \(\beta_j\) now because there are no blocks.

Hypotheses

\[H_0: \alpha_1 = \alpha_2 = \alpha_3 = \alpha_4 = 0 \qquad H_a: \text{at least one } \alpha_i \neq 0\]

Analysis

model2 <- lm(strength ~ chemical2, data = dat2)

gad(model2)
## $anova
## Analysis of Variance Table
## 
## Response: strength
##           Df Sum Sq Mean Sq F value Pr(>F)
## chemical2  3  12.95  4.3167  0.3863 0.7644
## Residuals 16 178.80 11.1750

Model check

qqnorm(residuals(model2), main = "Normal Q-Q Plot of the Residuals (CRD)")
qqline(residuals(model2))

boxplot(strength ~ chemical2, data = dat2,
        col = "lightblue",
        main = "Strength by Chemical Agent",
        xlab = "Chemical Agent",
        ylab = "Tensile Strength")

The residuals look normal and the four boxes are about the same height, so the model is fine. But the boxes are very tall and they overlap almost completely, which already shows that this design will not find a difference.

Conclusions

F = 0.3863 and the p-value is 0.7644. The p-value is much bigger than 0.15, so we fail to reject \(H_0\). With this design there is no evidence that the chemical agent changes the strength of the cloth.

Question 3: Comparison

The two analyses use exactly the same 20 numbers and give opposite answers.

Randomized block Completely randomized
SS chemical 12.95 12.95
SS bolt 157.00 not separated
SS error 21.80 178.80
df error 12 16
MS error 1.82 11.18
F chemical 2.3761 0.3863
p-value 0.1211 0.7644
Decision at 0.15 Reject H0 Fail to reject H0

The sum of squares of the chemical is 12.95 in both cases, because the chemical means do not change. What changes is the error. In the block design the bolt takes 157.00 out of the total, and only 21.80 is left as error. In the completely randomized design there is no block to take that part, so it goes into the error: 21.80 + 157.00 = 178.80.

That makes the mean square error jump from 1.82 to 11.18, about six times bigger. The F ratio has the mean square error at the bottom, so the F for the chemical falls from 2.3761 to 0.3863 and the p-value goes from 0.1211 to 0.7644. The real effect of the chemical is still there, but it is hidden inside the noise coming from the bolts.

So yes, the bolt of cloth is a significant source of nuisance variability. Its sum of squares (157.00) is more than twelve times the sum of squares of the chemical (12.95), its own test gives a p-value of 0.0000206, and ignoring it is what makes the chemical effect disappear. Blocking on the bolt is what lets the experiment see the difference between the chemical agents.

Complete code

library(GAD)

# ---- Question 1: randomized block design ----
strength <- c(73, 68, 74, 71, 67,
              73, 67, 75, 72, 70,
              75, 68, 78, 73, 68,
              73, 71, 75, 75, 69)

chemical <- as.fixed(factor(rep(1:4, each = 5)))
bolt <- as.random(factor(rep(1:5, times = 4)))

dat <- data.frame(chemical, bolt, strength)

head(dat)
str(dat)

tapply(dat$strength, dat$chemical, mean)
tapply(dat$strength, dat$bolt, mean)

model1 <- lm(strength ~ chemical + bolt, data = dat)

gad(model1)

estimates(model1)

qqnorm(residuals(model1), main = "Normal Q-Q Plot of the Residuals")
qqline(residuals(model1))

plot(fitted(model1), residuals(model1),
     main = "Residuals vs Fitted Values",
     xlab = "Fitted Value",
     ylab = "Residual",
     pch = 19)
abline(h = 0, lty = 2)

# ---- Question 2: completely randomized design ----
chemical2 <- as.fixed(factor(rep(1:4, each = 5)))

dat2 <- data.frame(chemical2, strength)

head(dat2)

model2 <- lm(strength ~ chemical2, data = dat2)

gad(model2)

qqnorm(residuals(model2), main = "Normal Q-Q Plot of the Residuals (CRD)")
qqline(residuals(model2))

boxplot(strength ~ chemical2, data = dat2,
        col = "lightblue",
        main = "Strength by Chemical Agent",
        xlab = "Chemical Agent",
        ylab = "Tensile Strength")