Rainbow-coloured threads streaming across a dark background

Week: Design of Experiments

Workshop 8A: Difference of Means with t-tests + Permutation Tests

Applying T-tests and ANOVA

Tim Fraser, PhD

Assistant Teaching Professor, Systems Engineering
Cornell University

SYSEN 5300: Systems Engineering & Six Sigma

for Design and Operation of Reliable Systems

Risk Analysis
- Risk assessment and risk characterization
- Failure Modes and Effects Analysis
Modeling Reliability
- Component and system reliability
- Physical acceleration models
Quality Control
- Six sigma and statistical process control
- Catch performance slipping, in time to act
System Improvement
- Design of experiments
- Response surfaces

Today's Class

Recap: Comparing 2 Groups
Comparing 3 Groups
Block Randomization with 3 Groups
Latin Squares

Rows of laboratory glassware on a bench, one flask in focus

Recap: Comparing 2 Groups

Workshop 8A

Designing Experiments

Example 1

Comparing 2 Unpaired Samples with T-tests in R

An industrial process made varying yields from methods A and B.

Does B have a significantly better yield than A?

Code: sysen/code/12_workshop.R

library(dplyr)
library(broom)
data = tibble(
  method = c("A", "A", "A", "A", "A", "A", "A", "A", "A", "A",
             "B", "B", "B", "B", "B", "B", "B", "B", "B", "B"),
  yield = c(89.7, 81.4, 84.5, 84.8, 87.3, 79.7, 85.1, 81.7,
            83.7, 84.5, 84.7, 86.1, 83.2, 91.9, 86.3, 79.3,
            82.6, 89.1, 83.7, 88.5))
# T-test assuming equal variance (unpaired)
t.test(formula = yield ~ method, data = data, var.equal = TRUE)

data:  yield by method
t = -0.88158, df = 18, p-value = 0.3896
alternative hypothesis: true difference in means between group A and group B is not equal to 0

# T-test not assuming equal variance
t.test(formula = yield ~ method, data = data, var.equal = FALSE)

data:  yield by method
t = -0.88158, df = 17.129, p-value = 0.3902
alternative hypothesis: true difference in means between group A and group B is not equal to 0

import pandas as pd
from scipy import stats

data = pd.DataFrame({
  'method': ["A"]*10 + ["B"]*10,
  'yield': [89.7, 81.4, 84.5, 84.8, 87.3, 79.7, 85.1, 81.7,
            83.7, 84.5, 84.7, 86.1, 83.2, 91.9, 86.3, 79.3,
            82.6, 89.1, 83.7, 88.5]})
a = data['yield'][data.method == "A"]
b = data['yield'][data.method == "B"]
# T-test assuming equal variance (unpaired)
r = stats.ttest_ind(a, b, equal_var=True)
print(f"t = {r.statistic:.5f}, df = {r.df:.0f}, "
      f"p-value = {r.pvalue:.4f}")

t = -0.88158, df = 18, p-value = 0.3896

# T-test not assuming equal variance
r = stats.ttest_ind(a, b, equal_var=False)
print(f"t = {r.statistic:.5f}, df = {r.df:.3f}, "
      f"p-value = {r.pvalue:.4f}")

t = -0.88158, df = 17.129, p-value = 0.3902

Example 2

Comparing Variance in 2 Groups with F-tests in R

An industrial process made varying yields from methods A and B.

Is the variance in yield for methods A vs. B significantly different?

library(dplyr)
library(broom)

data = tibble(
  method = c("A", "A", "A", "A", "A", "A", "A", "A", "A", "A",
             "B", "B", "B", "B", "B", "B", "B", "B", "B", "B"),
  yield = c(89.7, 81.4, 84.5, 84.8, 87.3, 79.7, 85.1, 81.7,
            83.7, 84.5, 84.7, 86.1, 83.2, 91.9, 86.3, 79.3,
            82.6, 89.1, 83.7, 88.5))

# F test to compare two variances.
var.test(formula = yield ~ method, data = data)

# Are they significantly different? (No --> p = 0.5049)

data:  yield by method
F = 0.63193, num df = 9, denom df = 9, p-value = 0.5049
alternative hypothesis: true ratio of variances is not equal to 1

import pandas as pd
from scipy import stats

data = pd.DataFrame({
  'method': ["A"]*10 + ["B"]*10,
  'yield': [89.7, 81.4, 84.5, 84.8, 87.3, 79.7, 85.1, 81.7,
            83.7, 84.5, 84.7, 86.1, 83.2, 91.9, 86.3, 79.3,
            82.6, 89.1, 83.7, 88.5]})
a = data['yield'][data.method == "A"]
b = data['yield'][data.method == "B"]

from scipy.stats import f

# F test to compare two variances.
stat = a.var(ddof=1) / b.var(ddof=1)
p = 2 * min(f.cdf(stat, len(a)-1, len(b)-1),
            f.sf(stat, len(a)-1, len(b)-1))
print(f"F = {stat:.5f}, num df = {len(a)-1}, "
      f"denom df = {len(b)-1}, p-value = {p:.4f}")

# Are they significantly different? (No --> p = 0.5049)

F = 0.63193, num df = 9, denom df = 9, p-value = 0.5049

Example 3

Comparing 2 Paired Samples with T-tests in R

An industrial process made a series of paired products, each using method A and B.

Do the products made by B have a significantly better yield than A?

library(dplyr)
library(broom)

# When data comes in pairs, like this, we use paired t-tests.
data2 = tibble(
  yield_a = c(89.7, 81.4, 84.5, 84.8, 87.3, 79.7, 85.1, 81.7,
              83.7, 84.5),
  yield_b = c(84.7, 86.1, 83.2, 91.9, 86.3, 79.3, 82.6, 89.1,
              83.7, 88.5)
)

# Basic method for paired t-test
t.test(data2$yield_a, data2$yield_b, paired = TRUE, var.equal = TRUE)

data:  data2$yield_a and data2$yield_b
t = -0.9736, df = 9, p-value = 0.3557
alternative hypothesis: true mean difference is not equal to 0

import pandas as pd
from scipy import stats

# When data comes in pairs, we use paired t-tests.
data2 = pd.DataFrame({
  'yield_a': [89.7, 81.4, 84.5, 84.8, 87.3,
              79.7, 85.1, 81.7, 83.7, 84.5],
  'yield_b': [84.7, 86.1, 83.2, 91.9, 86.3,
              79.3, 82.6, 89.1, 83.7, 88.5]})

# Basic method for paired t-test
r = stats.ttest_rel(data2.yield_a, data2.yield_b)
print(f"t = {r.statistic:.4f}, df = {r.df:.0f}, "
      f"p-value = {r.pvalue:.4f}")

t = -0.9736, df = 9, p-value = 0.3557

Example 4

Permutation Test for Difference of Means

Same yields from methods A and B as in Example 1.

Does B have a significantly better yield than A?

Histogram of 1000 shuffled differences of means centered near zero, with the observed difference of -1.3 marked by a red vertical line

library(dplyr)
library(broom)
data = tibble(
  method = c("A", "A", "A", "A", "A", "A", "A", "A", "A", "A",
             "B", "B", "B", "B", "B", "B", "B", "B", "B", "B"),
  yield = c(89.7, 81.4, 84.5, 84.8, 87.3, 79.7, 85.1, 81.7,
            83.7, 84.5, 84.7, 86.1, 83.2, 91.9, 86.3, 79.3,
            82.6, 89.1, 83.7, 88.5))
set.seed(12345)
perms = tibble(rep = 1:1000) %>%
  group_by(rep) %>%
  reframe(data) %>%
  group_by(rep) %>%
  mutate(yield = sample(yield, size = n(), replace = FALSE)) %>%
  group_by(rep) %>%
  summarize(dbar = mean(yield[method == "A"]) -
                   mean(yield[method == "B"]))
obs = data %>%
  summarize(dbar = mean(yield[method == "A"]) -
                   mean(yield[method == "B"]))
mean( abs(perms$dbar) >= abs(obs$dbar) )

import numpy as np
import pandas as pd

data = pd.DataFrame({
  'method': ["A"]*10 + ["B"]*10,
  'yield': [89.7, 81.4, 84.5, 84.8, 87.3, 79.7, 85.1, 81.7,
            83.7, 84.5, 84.7, 86.1, 83.2, 91.9, 86.3, 79.3,
            82.6, 89.1, 83.7, 88.5]})

# Shuffling is random -- fix the seed so your numbers match.
rng = np.random.default_rng(12345)
y = data['yield'].to_numpy()
is_a = (data.method == "A").to_numpy()
dbar = []
for rep in range(1000):
  z = rng.permutation(y)
  dbar.append(z[is_a].mean() - z[~is_a].mean())
perms = np.array(dbar)
obs = y[is_a].mean() - y[~is_a].mean()
np.mean(np.abs(perms) >= np.abs(obs))

Activity — Practice It (Permutation Tests)

Shuffle the labels yourself: run a permutation test on a new pair of samples.

Open the activity →

Comparing 3 Groups

Workshop 8A

Designing Experiments

Case Study: Blood Coagulation Times

24 animals receiving four different diets A, B, C, D.
Diets randomly assigned to animals
Testing done in random order.

Is there real difference between the mean coagulation times for the 4 diets?

Analysis of Variance (ANOVA) table

data · code

	A	B	C	D
	62⁽²⁰⁾	63⁽¹²⁾	68⁽¹⁶⁾	56⁽²³⁾
	60⁽²⁾	67⁽⁹⁾	66⁽⁷⁾	62⁽³⁾
	63⁽¹¹⁾	71⁽¹⁵⁾	71⁽¹⁾	60⁽⁶⁾
	59⁽¹⁰⁾	64⁽¹⁴⁾	67⁽¹⁷⁾	61⁽¹⁸⁾
	63⁽⁵⁾	65⁽⁴⁾	68⁽¹³⁾	63⁽²²⁾
	59⁽²⁴⁾	66⁽⁸⁾	68⁽²¹⁾	64⁽¹⁹⁾
Treatment average	61	66	68	61
Grand average	64	64	64	64
Difference	−3	+2	+4	−3

Dot plot of the 24 coagulation times, one row per diet with D on top and A at the bottom, each diet's dots on a soft blob in its own colour: A tan, B purple, C blue, D rose, on an axis from 56 to 71

ANOVA: Analysis of Variance

The same dot plot with a dark vertical line at the grand mean of 64 and a coloured tick at each diet's group mean: D 61 red, C 68 light blue, B 66 purple, A 61 yellow

Key Question: How much better are the groups at explaining variation than the grand mean?

y_ti = η + τ_t + ε_ti

y_ti: i-th observation, diet tη: grand meanτ_t: treatment effectε_ti: IID error

	A	B	C	D
	62⁽²⁰⁾	63⁽¹²⁾	68⁽¹⁶⁾	56⁽²³⁾
	60⁽²⁾	67⁽⁹⁾	66⁽⁷⁾	62⁽³⁾
	63⁽¹¹⁾	71⁽¹⁵⁾	71⁽¹⁾	60⁽⁶⁾
	59⁽¹⁰⁾	64⁽¹⁴⁾	67⁽¹⁷⁾	61⁽¹⁸⁾
	63⁽⁵⁾	65⁽⁴⁾	68⁽¹³⁾	63⁽²²⁾
	59⁽²⁴⁾	66⁽⁸⁾	68⁽²¹⁾	64⁽¹⁹⁾
Treatment average	61	66	68	61
Grand average	64	64	64	64
Difference	−3	+2	+4	−3

E(m_T) = Σ τ_t² + σ² E(m_R) = σ² No difference in the four treatments, i.e. H₀: τ₁ = τ₂ = τ₃ = τ₄ = 0 Then both m_T and m_R estimate σ²

Calculating Deviations in ANOVA

Source of Variation	Sum of Squares	Degrees of Freedom	Mean Square
Between treatments	S_T = 228	ν_T = 3	m_T = 76.0	F_3,20 = 13.6
Within treatments	S_R = 112	ν_R = 20	m_R = 5.6
Total about the grand average	S_D = 340	ν_D = 23

Observations y_ti

62	63	68	56
60	67	66	62
63	71	71	60
59	64	67	61
63	65	68	63
59	66	68	64

Y

:

Deviations from grand average of 64

−2	−1	4	−8
−4	3	2	−2
−1	7	7	−4
−5	0	3	−3
−1	1	4	−1
−5	2	4	0

D = Y − 64

340

23

=

Treatment deviations

−3	2	4	−3
−3	2	4	−3
−3	2	4	−3
−3	2	4	−3
−3	2	4	−3
−3	2	4	−3

T

228

3

+

Residuals within-treatment deviations

1	−3	0	−5
−1	1	−2	1
2	5	3	−1
−2	−2	−1	0
2	−1	0	2
−2	0	0	3

R

112

20

	A	B	C	D
	62⁽²⁰⁾	63⁽¹²⁾	68⁽¹⁶⁾	56⁽²³⁾
	60⁽²⁾	67⁽⁹⁾	66⁽⁷⁾	62⁽³⁾
	63⁽¹¹⁾	71⁽¹⁵⁾	71⁽¹⁾	60⁽⁶⁾
	59⁽¹⁰⁾	64⁽¹⁴⁾	67⁽¹⁷⁾	61⁽¹⁸⁾
	63⁽⁵⁾	65⁽⁴⁾	68⁽¹³⁾	63⁽²²⁾
	59⁽²⁴⁾	66⁽⁸⁾	68⁽²¹⁾	64⁽¹⁹⁾
Treatment average	61	66	68	61
Grand average	64	64	64	64
Difference	−3	+2	+4	−3

Plain dot plot of the 24 coagulation times, one row per diet, D on top, on an axis from 56 to 71

Quantities of Interest in ANOVA

Source of Variation	Sum of Squares	Degrees of Freedom	Mean Square
Between treatments	S_T = 228	ν_T = 3	m_T = 76.0	F_3,20 = 13.6
Within treatments	S_R = 112	ν_R = 20	m_R = 5.6
Total about the grand average	S_D = 340	ν_D = 23

Total Sum of Squares measures the total variation in your data, using the grand mean as benchmark.
TSS = Σ (x − x̄_grand)²
Residual Sum of Squares is the unexplained error left within each treatment group.
RSS = Σ (x − x̄_group)²
Explained Sum of Squares is the deviation explained by your group means.
ESS = TSS − RSS

n = sample size (24); k = # of groups (4)
Mean Squares Explained: average variation explained, divided by k − 1 groups.
MSE = (TSS − RSS) / (k − 1)
Mean Squared Error: average variation unexplained, divided by n − k.
MSR = RSS / (n − k)
F-statistic: ratio of variation explained versus unexplained.
F = MSE / MSR

F-statistic in ANOVA

If the error is normally distributed (independence of observation assumption of regression)
Then MSE and MSR are independently distributed.
So, ratio F = MSE / MSR follows an F distribution.
F-statistic: How much better is your model than an intercept model alone? (grand mean)

data · code

Under H₀ that Στ² = 0,

F = m_T / m_R ~ F_3,20

The F distribution with 3 and 20 degrees of freedom: a sharp peak near zero and a long right tail, with arrows at the 5% point (3.10), the 1% point (4.94) and the 0.1% point (8.10), and the observed ratio 13.6 far out in the tail

Example 5

Comparing 3+ group means with F-test in R

Craig, Kim, and Melanie work for a bakery, but Craig’s donuts seem kinda weird.

Does the weight of these donuts vary significantly by the baker?

Calculate an F-statistic and evaluate its significance.

Code: sysen/code/12_workshop.R

library(dplyr); library(readr); library(broom)
url = "https://raw.githubusercontent.com/timothyfraser/sysen/main/"
donuts = read_csv(paste0(url, "workshops/donuts.csv"))

import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
url = "https://raw.githubusercontent.com/timothyfraser/sysen/main/"
donuts = pd.read_csv(url + "workshops/donuts.csv")

Violin plots of donut weight for three bakers: Craig at the bottom with a long right tail up to about 39, Kim in the middle tightly around 28.5, Melanie at the top with a gap before a few heavy donuts; a large blue dot marks each baker's mean and a vertical blue line marks the grand mean of 30

# use lm()
donuts %>%
  summarize(lm(weight ~ baker) %>% glance()) %>%
  select(sigma, statistic, p.value, df)

sigma	statistic	p.value	df
2.94	6.24	0.00395	2

# or use aov()
donuts %>%
  reframe(aov(weight ~ baker) %>% tidy())

term	df	sumsq	meansq	statistic	p.value
baker	2	108.	54.0	6.24	0.00395
Residuals	47	407.	8.66	NA	NA

# use lm()
m = smf.ols('weight ~ baker', data=donuts).fit()
# sigma, statistic, p.value, df
m.mse_resid ** 0.5, m.fvalue, m.f_pvalue, m.df_model

(2.9424, 6.2398, 0.003951, 2.0)

# or use aov()
sm.stats.anova_lm(m, typ=2)

          sum_sq    df       F    PR(>F)
baker      108.05   2.0  6.2398  0.003951
Residual   406.91  47.0     NaN       NaN

Example 6

Comparing unequal group variances with K² test in R

Craig, Kim, and Melanie work for a bakery, but Craig’s donuts seem kinda weird.

Does the variance in weight vary significantly by baker? Calculate a K² statistic.
If so, do NOT assume equal variance in ANOVA. Calculate your F statistic.

library(dplyr); library(readr); library(broom)
url = "https://raw.githubusercontent.com/timothyfraser/sysen/main/"
donuts = read_csv(paste0(url, "workshops/donuts.csv"))

# Are the variances of my 3+ groups significantly different?
# Homogeneity of Variance - Bartlett's test for K^2
bartlett.test(weight ~ baker, data = donuts)
# If K-squared is significant, the differences
# in variance are significant.

data:  weight by baker
Bartlett's K-squared = 10.64, df = 2, p-value = 0.004892

donuts %>%
  summarize(
    oneway.test(weight ~ baker,
                data = donuts, var.equal = FALSE) %>% tidy())

num.df	den.df	statistic	p.value
2	28.4	6.14	0.00607

import numpy as np
from scipy import stats
# donuts: the same data as Example 5
g = [d['weight'].to_numpy()
     for _, d in donuts.groupby('baker')]
# Bartlett's test for K^2
stats.bartlett(*g)

BartlettResult(statistic=10.64, pvalue=0.004892)

# Welch's ANOVA -- no equal-variance assumption
n = np.array([len(x) for x in g])
mu = np.array([x.mean() for x in g])
w = n / np.array([x.var(ddof=1) for x in g])
mw = (w * mu).sum() / w.sum()
k = len(g)
lam = ((1 - w / w.sum())**2 / (n - 1)).sum()
num = (w * (mu - mw)**2).sum() / (k - 1)
stat = num / (1 + 2 * (k - 2) / (k**2 - 1) * lam)
df2 = 1 / (3 / (k**2 - 1) * lam)
k - 1, df2, stat, stats.f.sf(stat, k - 1, df2)

(2, 28.4, 6.14, 0.00607)

Randomized Block Design with 3+ groups

Workshop 8A

Designing Experiments

Randomized Block Design

To thoroughly control for alternative explanations, we can randomly assign observations to treatments within blocks.
Example: Penicillin
Estimate effect of treatments accounting for blends.

Block	Treatment				Block averages	Block deviations
Block	A	B	C	D	Block averages	Block deviations
Blend 1	89⁽¹⁾	88⁽³⁾	97⁽²⁾	94⁽⁴⁾	92	+6
Blend 2	84⁽⁴⁾	77⁽²⁾	92⁽³⁾	79⁽¹⁾	83	−3
Blend 3	81⁽²⁾	87⁽¹⁾	87⁽⁴⁾	85⁽³⁾	85	−1
Blend 4	87⁽¹⁾	92⁽³⁾	89⁽²⁾	84⁽⁴⁾	88	+2
Blend 5	79⁽³⁾	81⁽⁴⁾	80⁽¹⁾	88⁽²⁾	82	−4
Treatment averages	84	85	89	86	Grand average:
Treatment deviations	−2	−1	+3	0	86

It is important to understand that in a randomized block experiment the treatments are randomized within the blocks.

Randomized Block Design (2)

Source of variation	Sum of squares	Degrees of freedom	Mean square	F ratio
Between blocks (blends)	S_B = 264	ν_B = (n − 1) = 4	m_B = 66.0	F_4,12 = 3.51
Between treatments	S_T = 70	ν_T = (k − 1) = 3	m_T = 23.3	F_3,12 = 1.24
Residuals	S_R = 226	ν_R = (n − 1)(k − 1) = 12	m_R = 18.8
Deviations from grand average	S_D = 560	nk − 1 = 19

Observations
y_bt

89	88	97	94
84	77	92	79
81	87	87	85
87	92	89	84
79	81	80	88

Y

Analysis of

observations

:

Deviations from Grand Average
y_bt − ȳ

3	2	11	8
−2	−9	6	−7
−5	1	1	−1
1	6	3	−2
−7	−5	−6	2

D

S_D = 560

ν_D = 19

=

Block Deviations
ȳ_b − ȳ

6	6	6	6
−3	−3	−3	−3
−1	−1	−1	−1
2	2	2	2
−4	−4	−4	−4

B

S_B = 264

ν_B = 4

+

Treatment Deviations
ȳ_t − ȳ

−2	−1	3
−2	−1	3
−2	−1	3
−2	−1	3
−2	−1	3

T

S_T = 70

ν_T = 3

+

Residuals
y_bt − ȳ_b − ȳ_t + ȳ

−1	−3	2	2
3	−5	6	−4
−2	3	−1	0
1	5	−2	−4
−1	0	−5	6

R

S_R = 226

ν_R = 12

Randomized Block Design (3)

Why Use Randomized Blocking?

Of the total sum of squares not associated with treatments or with the mean, almost 50% is just block-to-block variation.
If experiment arranged on a completely randomized basis with no blocks, the error variance would have been much larger.
With randomized block, error much less: of the total of S_D = 560, S_B = 264 removed by blocks.
Randomized block design greatly increased sensitivity of experiment; now able to detect smaller treatment differences.

Formulation

y_ti = η + β_i + τ_t + ε_ti
η_ti = η + β_i + τ_t

Additive: τ₃ = 6 and block 4 adds 5, so both together add 6 + 5 = 11

If not additive: a block × treatment interaction

Randomized Block Design (4)

Diagnostic Tests

Dot plots of the 20 penicillin residuals on an axis from minus 6 to 6, first all together and then split by blend 1 to 5 and by treatment 1 to 4

Figure 4.9 (a) Dot plots of residuals. (b) Residuals identified by block and treatment.

Scatter plot of the penicillin residuals against the predicted values from 80 to 95, scattered evenly around zero

Figure 4.10 Residuals plotted against the predicted values: penicillin experiment. data · code

Latin Squares

Workshop 8A

Designing Experiments

Using 2 Types of Blocks: Air Pollution

How much does adding Chemical A, B, C, or D to gasoline reduce air pollution?

4 treatments tested
4 different drivers
4 different cars
2 block factors: drivers & cars
Latin Squares Design helps eliminate differences in treatment comparisons that are just due to drivers or cars.
Each treatment (A, B, C, or D) randomly appears 1 per row (driver) and 1 per column (car).

Driver	1	2	3	4
I	A 19	B 24	D 23	C 26
II	D 23	C 24	A 19	B 30
III	B 15	D 14	C 15	A 16
IV	C 19	A 18	B 19	D 16

Cars	Drivers	Additives
1: 19	I: 23	A: 18
2: 20	II: 24	B: 22
3: 19	III: 15	C: 21
4: 22	IV: 18	D: 19

Cars 1–4 across, drivers I–IV down. Grand average: 20

Using 2 Types of Blocks (2): Air Pollution

No convincing evidence for differences between the treatments
Latin square design effective in eliminating variation due to drivers

Source of Variation	Sum of Squares	df	Mean Square	Ratio of Mean Squares	P
Cars (columns)	S_C = 24	3	m_C = 8.00	F_3,6 = m_C/m_R = 1.5	0.31
Drivers (rows)	S_D = 216	3	m_D = 72.00	F_3,6 = m_D/m_R = 13.5	<0.01
Treatments (additives)	S_T = 40	3	m_T = 13.33	F_3,6 = m_T/m_R = 2.5	0.16
Residuals	S_R = 32	6	m_R = 5.33
Total	S_V = 312	15

Observations (cars 1–4)

A 19	B 24	D 23	C 26
D 23	C 24	A 19	B 30
B 15	D 14	C 15	A 16
C 19	A 18	B 19	D 16

Y

:

Deviations from grand average (20)

−1	4	3	6
3	4	−1	10
−5	−6	−5	−4
−1	−2	−1	−4

V

312

15

=

Columns (cars)

−1	−1	2
−1	−1	2
−1	−1	2
−1	−1	2

C

24

3

+

Rows (drivers)

3	3	3	3
4	4	4	4
−5	−5	−5	−5
−2	−2	−2	−2

D

216

3

+

Treatments (additives)

−2	2	−1	1
−1	1	−2	2
2	−1	1	−2
1	−2	2	−1

T

40

3

+

Residuals

−1	−1	2	0
1	−1	−2	2
−1	0	0	1
1	2	0	−3

R

32

6

Using 3 Types of Blocks: Graeco-Latin Squares

A Graeco-Latin square is a k by k pattern
Allows study of k treatment simultaneously, with 3 blocking variables, each at k levels
Example: 1 extra blocking variable in car emissions
Eg. day of test

Driver	Car 1	Car 2	Car 3	Car 4
I	Aα	Bβ	Cγ	Dδ
II	Bδ	Aγ	Dβ	Cα
III	Cβ	Dα	Aδ	Bγ
IV	Dγ	Cδ	Bα	Aβ

Additives: A, B, C, D · Days: α, β, γ, δ

Using 3+ Types of Blocks: Hyper-Graeco Squares

Machine used to test wearing quality of types of cloth. Records weight loss from test piece in 1 machine cycle (rubbed against a std. grade of emery paper)
Treatments: 4 types of cloth A, B, C, D
Block Type 1: Mounted in 4 specimen holder 1, 2, 3, 4.
Block Type 2: Each holder can be in 1 of 4 positions P1, P2, P3, P4.
Block Type 3: Each emery paper sheet alpha, beta, gamma, delta was cut into 4 quarters.
Objective: (1) Make accurate comparison of the treatments; (2) Understand variability caused by various factors-holders, positions, emery papers, and cycles

A smiling corgi looking out of the round door of a stainless-steel washing machine

Martindale Wear Tester: the design

Cycles	P₁	P₂	P₃	P₄	Legend
C₁	αA1 / 320	βB2 / 297	γC3 / 299	δD4 / 313	Replicate I Cycles: C₁–C₄ Treatments: A, B, C, D Holders: 1, 2, 3, 4 Emery paper sheets: α, β, γ, δ
C₂	βC4 / 266	αD3 / 227	δA2 / 260	γB1 / 240
C₃	γD2 / 221	δC1 / 240	αB4 / 267	βA3 / 252
C₄	δB3 / 301	γA4 / 238	βD1 / 243	αC2 / 290

Table 4.12. Cell = emery paper, treatment, holder / wear.

Cycles	P₁	P₂	P₃	P₄	Legend
C₅	εA1 / 285	ξB2 / 280	θC3 / 331	κD4 / 311	Replicate II Cycles: C₅–C₈ Treatments: A, B, C, D Holders: 1, 2, 3, 4 Emery paper sheets: ε, ξ, θ, κ
C₆	ξC4 / 268	εD3 / 233	κA2 / 291	θB1 / 280
C₇	θD2 / 265	κC1 / 273	εB4 / 234	ξA3 / 243
C₈	κB3 / 306	θA4 / 271	ξD1 / 270	εC2 / 272

The design was effective both in removing sources of extraneous variation and in indicating their relative importance.
Because of the elimination of these disturbances, the residual variance was reduced by a factor of 8
We could detect much smaller differences in treatment

data · code

Martindale Wear Tester: the analysis

Treatments	A: 270.0	B: 275.6	C: 279.9	D: 260.4
Holders	1: 268.9	2: 272.0	3: 274.0	4: 271.0
Positions	P₁: 279.0	P₂: 257.4	P₃: 274.4	P₄: 275.1
Emery papers	α: 276.0	β: 264.5	γ: 249.5	δ: 278.5	ε: 256.0	ξ: 265.3	θ: 286.8	κ: 295.3
Cycles	C₁: 307.3	C₂: 248.3	C₃: 245.0	C₄: 268.0	C₅: 301.8	C₆: 268.0	C₇: 253.8	C₈: 279.8
Replicates	I: 267.1	II: 275.8	Grand average = 271.5

Source	df	Sum of Squares	Mean Squares	Ratio of Mean Sq.
Σ(y − ȳ)²	31	26,463.97
Replications	1	603.78	m_Rep = 603.78	m_Rep/m_R = 5.73
Cycles	6	14,770.44	m_C = 2,461.74	m_C/m_R = 23.35
Positions	3	2,217.34	m_P = 739.11	m_P/m_R = 7.01
Emery papers	6	6,108.94	m_E = 1,018.16	m_E/m_R = 9.66
Holders	3	109.09	m_H = 36.36	m_H/m_R = 0.34
Treatments	3	1,705.34	m_T = 568.45	m_T/m_R = 5.39
Residuals	9	949.03	m_R = 105.45

The F-stat=5.39 with 3, 9 DOF, significant at 2% level.
By using a design which makes it possible to remove the effects of many larger disturbing factors, differences between treatment were made detectable.

data · code

Martindale Wear Tester: graphical analysis

Graphical ANOVA: dots for the replicates, cycles, emery papers, positions, holders, treatments and residuals of the Martindale wear example; position P2 sits far to the left

Graphical ANOVA for the Martindale wear example (Box, Hunter & Hunter, Figure 4.14). data · code

In graphical analysis, position P2 gives much less wear than others, indicating a need of improvement
The analysis identified the large contributions to the total variation due to cycles and to emery papers.
This suggested improvements, which later led to changes in the design of the machine.

Balanced Incomplete Block Design

Suppose the Martindale wear tester allowed only three samples to be included in each cycle, but you had 4 treatment A, B, C, and D to compare
You have 4 treatments but a block size of 3, too small to accommodate all the treatments simultaneously → balanced incomplete block design.
Property: every pair of treatment occurs together in a block the same number of times.

Table 4.14. A Balanced Incomplete Block Design, t = 4 Treatments in b = 4 Blocks of Size k = 3

Block (cycle) of 1000 revolutions	Treatments
1	A	B	C
2	A	B	D
3	A	C	D
4	B	C	D

Block (cycle) of 1000 revolutions	A	B	C	D
1	x	x	x
2	x	x		x
3	x		x	x
4		x	x	x

Doubly Balanced Designs: Youden Squares

Compare 7 treatment in seven blocks of size 4 (e.g., test 7 types of cloth A, B, C, D, E, F, and G)
But only 4 test pieces could be compared simultaneously in a single machine cycle)
Also had the opportunity to eliminate a second source of block variation, machine positions

Block (cycle)	A	B	C	D	E	F	G
1		α 627		β 248		γ 563	δ 252
2	α 344		β 233			δ 442	γ 226
3			α 251	γ 211	δ 160		β 297
4	β 337	δ 537			γ 195		α 300
5		γ 520	δ 278		β 199	α 595
6	γ 369			δ 196	α 185	β 606
7	δ 396	β 602	γ 240	α 273