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

Week: Comparing Groups: t-tests & ANOVA

Workshop 12: ANOVA & Latin Squares

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

Recap: Comparing 2 Groups

Workshop 12

Designing Experiments

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?

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))

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]})

Code: sysen/code/12_workshop.R

T-test assuming equal variance

t.test(formula = yield ~ method, data = data, var.equal = TRUE)

a = data['yield'][data.method == "A"]
b = data['yield'][data.method == "B"]

stats.ttest_ind(a, b, equal_var=True)

R console output of a two-sample t-test assuming equal variance

T-test not assuming equal variance (Welch)

t.test(formula = yield ~ method, data = data, var.equal = FALSE)

stats.ttest_ind(a, b, equal_var=False)

R console output of a Welch two-sample t-test

Comparing Variance in 2 Groups with F-tests

Same yields from methods A and B (data from Example 1).
Is the variance in yield for methods A vs. B significantly different?

F test to compare two variances

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

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

from scipy.stats import f

a = data['yield'][data.method == "A"]
b = data['yield'][data.method == "B"]

# F test to compare two variances.
stat = a.var(ddof=1) / b.var(ddof=1)
2 * min(f.cdf(stat, len(a)-1, len(b)-1),
        f.sf(stat, len(a)-1, len(b)-1))

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

R console output of var.test comparing the variance of yield between methods A and B

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))

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 a paired t-test

t.test(data2$yield_a, data2$yield_b,
       paired = TRUE, var.equal = TRUE)

stats.ttest_rel(data2.yield_a, data2.yield_b)

R console output of a paired t-test comparing yield_a and yield_b

Permutation Test for Difference of Means

Same yields from methods A and B (data from Example 1).
Does B have a significantly better yield than A?
Shuffle the labels to see what “no difference” looks like.

Shuffle the labels 1000 times

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"]))

import numpy as np

# Shuffling is random -- your numbers will vary.
y = data['yield'].to_numpy()
is_a = (data.method == "A").to_numpy()

dbar = []
for rep in range(1000):
  z = np.random.permutation(y)
  dbar.append(z[is_a].mean() - z[~is_a].mean())
perms = np.array(dbar)

Compare the observed difference to the null

# Get the observed difference of means
obs = data %>%
  summarize(dbar = mean(yield[method == "A"]) -
                   mean(yield[method == "B"]))

# 2-tailed percentage of random stats
# greater than the observed
mean( abs(perms$dbar) >= abs(obs$dbar) )

# Get the observed difference of means
obs = y[is_a].mean() - y[~is_a].mean()

# 2-tailed percentage of random stats
# greater than the observed (varies by run)
np.mean( np.abs(perms) >= np.abs(obs) )

Histogram of 1000 permuted differences of means with the observed difference marked

Comparing 3 Groups

Workshop 12

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?

Dot plot of coagulation times for diets A, B, C and D with the grand mean of 64 dashed

Diet A	Diet B	Diet C	Diet 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
61	66	68	61
−3	+2	+4	−3

Coagulation times (s). Row 7: diet mean. Row 8: difference from grand mean (64).

Data: sysen/workshops/anova_coagulation.csv · Code: sysen/functions/slides_anova_workshop12.R

ANOVA: Analysis of Variance

Model · one-way ANOVA

y_ti = η + τ_t + ε_ti

y_ti: the i-th observation in diet group t

η: overall (grand) meanτ_t: treatment effectε_ti: error, assumed IID

Expected mean squares

E(m_T) = Σ τ_t² + σ² E(m_R) = σ² H₀: τ₁ = τ₂ = τ₃ = τ₄ = 0

Then both m_T and m_R estimate σ²

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

Calculating Deviations in ANOVA

Worked table of deviations of each observation from the grand mean and from its group mean

Diagram naming the total, explained and residual deviation for a single observation

Three Sums of Squares

Total Sum of Squares measures the total variation in your data, using the grand mean as benchmark.
- TSS = sum( (x - xbbar)^2 )
Residual Sum of Squares measures how much deviation remains within each treatment group. This is the unexplained error in your model.
- RSS = sum( (x - xbar)^2 )
Explained Sum of Squares measures how much deviation was explained by your treatment group means.
- ESS = TSS - RSS

Quantities of Interest in ANOVA

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 by your model.
- F = MSE / MSR

n = sample size (24) k = # of groups (4)

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)

F distribution curve with the upper tail shaded

Comparing 3+ group means with F-test

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.

Load the donuts data

library(dplyr)
library(readr)
library(broom)

donuts = read_csv(
  "https://raw.githubusercontent.com/timothyfraser/sysen/main/workshops/donuts.csv")

import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

donuts = pd.read_csv(
  "https://raw.githubusercontent.com"
  "/timothyfraser/sysen/main"
  "/workshops/donuts.csv")

Route 1 — use lm()

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

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

R console output showing sigma, statistic, p.value and df from glance

Route 2 — or use aov()

donuts %>%
  reframe(aov(weight ~ baker) %>% tidy())

m = smf.ols('weight ~ baker', data=donuts).fit()
sm.stats.anova_lm(m, typ=2)

R console output of a tidied analysis of variance table for weight by baker

Comparing unequal group variances with K² test

Same donuts as Example 5 (donuts data).
Does the variance in donut weight differ significantly by baker? Calculate a K² statistic.
If so, do NOT assume equal variance in ANOVA. Calculate your F statistic.

Homogeneity of Variance — Bartlett's test for K²

# Are the variances of my 3+ groups significantly different?
bartlett.test(weight ~ baker, data = donuts)

# If K-squared is significant,
# the differences in variance are significant.

from scipy import stats

# Are the variances of my 3+ groups different?
g = [d['weight'] for _, d in donuts.groupby('baker')]
stats.bartlett(*g)

# If K-squared is significant,
# the differences in variance are significant.

R console output of Bartlett's test of homogeneity of variances

ANOVA without the equal-variance assumption

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

# Welch's ANOVA -- no equal-variance assumption
import numpy as np
g = [d['weight'].to_numpy()
     for _, d in donuts.groupby('baker')]
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)
stat, stats.f.sf(stat, k - 1, df2)

R console output of a tidied one-way test that does not assume equal variance

Randomized Block Design with 3+ groups

Workshop 12

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.
N = 20 animals, five blends of 4.

Diagram of 20 animals split into five blends of four, each blend receiving treatments A, B, C and D

Randomized Block Design (2)

R = D - B - T; the vectors R, B, and T are mutually orthogonal.

Vector decomposition of the penicillin data into block, treatment and residual components

Analysis of variance table for the randomized block penicillin experiment

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 the experiment were arranged on a completely randomized basis with no blocks, the error variance would have been much larger.
With randomized blocks, error is much less: of the total of SD = 560, SB = 264 is removed by blocks.
Randomized block design greatly increased sensitivity of the experiment; now able to detect smaller treatment differences.

The Randomized Block Model

Observed response

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

grand level + block i + treatment t + error

Expected response

η_ti = η + β_i + τ_t

Additive: effects simply add.

Worked · additive

τ₃ = 6, block 4 adds 5 6 + 5 = 11

If not additive: a block × treatment interaction

Randomized Block Design (4) — Diagnostic Tests

Diagnostic plots for the randomized block penicillin model

Residuals from the randomized block model plotted by blend and by treatment

Latin Squares

Workshop 12

Designing Experiments

Example: Air Pollution — the question

How much does adding Chemical A, B, C, or D to gasoline reduce air pollution?
Hypothesized relationship: chemical type → air pollution.
Confounding variables: driver differences, car differences.

The design rule

4 treatments, 4 drivers, 4 cars: 2 block factors, drivers & cars.
Latin squares remove driver and car differences from the comparison.
Each treatment (A, B, C, or D) randomly appears once per row (driver) and once per column (car).

Air Pollution — the 4 × 4 Latin square data

Driver	Car 1	Car 2	Car 3	Car 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

Emissions by driver (row) and car (column); the italic letter is the additive.

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

Averages. Grand average: 20

Data: sysen/workshops/anova_latin_square.csv · Code: sysen/functions/slides_anova_workshop12.R

Air Pollution — analysis of variance

Source	Sum of Squares	df	Mean Square	F	p
Cars (columns)	24	3	8.00	1.5	0.31
Drivers (rows)	216	3	72.00	13.5	0.004
Treatments (additives)	40	3	13.33	2.5	0.16
Residuals	32	6	5.33
Total	312	15

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

Code: sysen/functions/slides_anova_workshop12.R

Using 3 Types of Blocks — Graeco-Latin Squares

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

Graeco-Latin square layout adding Greek-letter blocks to the Latin square of drivers and cars

Using 3+ Types of Blocks — Example: Martindale Wear Tester

Machine used to test wearing quality of types of cloth.
Records weight loss from a 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 holders 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 holders, positions, emery papers, and cycles.

Hyper-Graeco Latin Squares — Example: Martindale Wear Tester

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.

Hyper-Graeco-Latin square layout for the Martindale wear tester experiment

Hyper-Graeco Latin Squares — what the analysis showed

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

Analysis of variance table for the Martindale wear tester experiment

Effect estimates for cloth, holder, position, emery paper and cycle

Hyper-Graeco Latin Squares — graphical analysis

In graphical analysis, position P2 gives much less wear than the others, indicating a need of improvement.

Dot plot of wear by holder position showing position P2 far below the others

Balanced Incomplete Block Design

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

Doubly Balanced Incomplete Block Designs — Youden Squares

Compare 7 treatments in seven blocks of size 4 (e.g. test 7 types of cloth A–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.