Comparing Two Groups split a cardinal outcome by a two-level factor and tested whether the two groups differed. We now extend that idea to \(G\) groups.
This chapter asks three questions.
- How does \(\hat{y}\) differ across \(G\) groups?
- Group ECDFs and conditional means: these use the ECDF and the conditional mean from earlier chapters, with one for each group.
- How sure are we?
- Notched boxplot: this is new, and it draws a bootstrap confidence interval for each group’s median. It uses the same recipe as the difference between two medians.
- Are the differences larger than we would expect by chance?
- ANOVA \(F\)-statistic and permutation test: the statistic is new, and it compares variation between groups to variation within groups. The test uses the same recipe as two groups.
- Kruskal-Wallis test: this is a new test that compares average ranks, so it only requires ordinal data.
- Categorical ANOVA: this is a new statistic for a factor outcome, which measures variation with one minus the Herfindahl index.
Describing Groups
A cardinal variable can be grouped into a factor with \(G\) levels. We pair these groups with a cardinal outcome.
Code
## MAIN EXAMPLE
# Load the data: 3294 US workers
data(Wages1, package='Ecdat')
# Outcome: hourly wage in 1980 dollars, one value per worker
y_hat <- Wages1[, 'wage']
# Group: years of schooling, cut into three levels (G=3).
# Its levels are the values the categorical variable can take.
# Here they are fewer than 12 years, exactly 12 years, and more than 12 years.
school_g <- cut(Wages1[, 'school'], breaks=c(0, 11, 12, 16),
labels=c('<12', '12', '>12'))
levels(school_g)
## [1] "<12" "12" ">12"
table(school_g)
## school_g
## <12 12 >12
## 1353 1188 753
Group Summaries
A factor with \(G\) levels splits a cardinal variable into \(G\) conditional distributions, one per group. The group summaries from Bivariate Distributions apply directly. The ECDF shows every observation, and drawing one per group on the same axes compares the whole distributions at once: a curve further to the right has higher values at every quantile. Each line’s thickness scales with its group size.
Code
## MAIN EXAMPLE
# One ECDF per schooling group, on the same axes
cols <- c(rgb(.8, 0, 0, .5), rgb(0, .8, 0, .5), rgb(0, 0, .8, .5))
n_g <- table(school_g)
lwd_g <- 1.5 + 2.5 * (n_g - min(n_g)) / (max(n_g) - min(n_g))
# The axis stops just past 20 dollars, so the 16 largest wages are not drawn
plot(NULL, xlim=c(0, 20), ylim=c(0, 1),
xlab='Hourly wage (1980 dollars)', ylab='Proportion of workers', main=NA)
for (g in seq_along(levels(school_g))) {
F_g_hat <- ecdf(y_hat[school_g == levels(school_g)[g]])
plot(F_g_hat, add=TRUE, col=cols[g], lwd=lwd_g[g],
verticals=TRUE, do.points=FALSE, col.01line=NA)
}
legend('bottomright', bty='n',
legend=paste0(levels(school_g), ' years (n=', n_g, ')'), col=cols, lty=1, lwd=lwd_g)
Code
# The curve for more than 12 years sits to the right of the other two.
# That holds over the whole range shown.
# So its higher wages are not driven by a few workers.
A table of group sizes and centers puts numbers on the figure. The conditional mean \(\hat{m}_{Y}(g)\) is the average of \(\hat{y}_{i}\) among the \(n_{g}\) observations in group \(g\), and aggregate() computes it for every group at once.
Code
## MAIN EXAMPLE
# Group size, mean, and median of the hourly wage by schooling group
summary_fun <- function(y_hat) {
out <- c(n=length(y_hat), mean=mean(y_hat), median=median(y_hat))
return(out)
}
aggregate(y_hat ~ school_g, FUN=summary_fun)
## school_g y_hat.n y_hat.mean y_hat.median
## 1 <12 1353.000000 4.990937 4.593337
## 2 12 1188.000000 5.668354 5.166288
## 3 >12 753.000000 7.275887 6.736894
# Workers with more than 12 years of schooling have the highest center.
# Those with fewer than 12 years have the lowest.
Are the differences between group centers larger than we would expect by chance?
Here is another example, with murder arrest rates across the four census regions, an unordered factor with \(G=4\) levels. To make the bootstrap intervals in the later callouts meaningful, pretend we observed only a random sample of \(35\) states, as we did with workers in Sampling.
Code
# USArrests covers all 50 states, so it is the whole population of states.
# Draw one random sample of 35 of the 50 states, and use only it from here on
set.seed(111)
state_id <- sample(nrow(USArrests), 35)
# Outcome: murder arrests per 100K, one value per sampled state
murder_hat <- USArrests[state_id, 'Murder']
# Group: the four-level region factor, for the same sampled states.
# Its levels are the values the categorical variable can take.
# The built-in state.region lists them in the same state order as USArrests.
# So it can be paired directly with the murder rates.
region <- state.region[state_id]
levels(region)
## [1] "Northeast" "South" "North Central" "West"
table(region)
## region
## Northeast South North Central West
## 5 14 8 8
The regions are ordered by median murder rate from lowest to highest, and each line’s thickness scales with its group size.
Code
# Order regions by their median murder rate, lowest to highest
region_order <- names(sort(tapply(murder_hat, region, median)))
region <- factor(region, levels=region_order)
# One ECDF per region, on the same axes
line_cols <- hcl.colors(4)
n_g <- table(region)
lwd_g <- 1.5 + 2.5 * (n_g - min(n_g)) / (max(n_g) - min(n_g))
par(mar=c(4, 4, 4.5, 1))
plot(NULL, xlim=range(murder_hat), ylim=c(0, 1),
xlab='Murder arrests (per 100,000)', ylab='Proportion of states', main=NA)
for (g in seq_along(levels(region))) {
F_g_hat <- ecdf(murder_hat[region == levels(region)[g]])
plot(F_g_hat, add=TRUE, col=line_cols[g], lwd=lwd_g[g],
verticals=TRUE, do.points=FALSE, col.01line=NA)
}
legend('topleft', inset=c(0, -.2), xpd=TRUE, ncol=2, bty='n',
legend=paste0(levels(region), ' (n=', n_g, ')'), col=line_cols, lty=1, lwd=lwd_g)
Code
# The South has the highest median, so it is last in the legend.
# Its curve sits to the right of the others over its whole range.
# So its higher values are not driven by a few states.
A table of group sizes and centers puts numbers on the figure.
Code
# Group size, mean, and median of the murder rate by region
aggregate(murder_hat ~ region, FUN=summary_fun)
## region murder_hat.n murder_hat.mean murder_hat.median
## 1 Northeast 5.00000 3.08000 3.30000
## 2 North Central 8.00000 5.48750 4.95000
## 3 West 8.00000 7.15000 7.35000
## 4 South 14.00000 11.58571 12.85000
# The South has the highest center and the Northeast the lowest.
Equal Means (ANOVA)
In Comparing Two Groups, we tested whether two groups had the same mean. We now extend this to \(G\) groups. The Analysis of Variance (ANOVA) approach asks: are the observed differences in group means larger than we would expect from random variation alone?
Decomposition
If group means differ, some of the spread in \(\hat{y}\) must come from differences between groups rather than within them. We can split each piece out.
Label each observation \(\hat{y}_{ig}\) for individual \(i\) in group \(g\), where group \(g=1,\dots,G\) has \(n_{g}\) observations with group mean \(\hat{m}_{Y}(g)\), and the grand mean across all \(n=\sum_{g} n_{g}\) observations is \(\hat{m}_{Y}\). The ANOVA decomposition splits the total sum of squares into between-group and within-group pieces: \[\begin{aligned}
\underbrace{\sum_{g=1}^{G}\sum_{i=1}^{n_{g}}(\hat{y}_{ig}-\hat{m}_{Y})^2}_{\hat{\text{TSS}}}
&=
\underbrace{\sum_{g=1}^{G} n_{g} (\hat{m}_{Y}(g)-\hat{m}_{Y})^2}_{\hat{\text{BSS}}}
\;+\;
\underbrace{\sum_{g=1}^{G}\sum_{i=1}^{n_{g}}(\hat{y}_{ig}-\hat{m}_{Y}(g))^2}_{\hat{\text{WSS}}},
\end{aligned}\] where \(\hat{\text{BSS}}\) measures variation in group means around the grand mean and \(\hat{\text{WSS}}\) measures variation of observations around their own group mean.
The decomposition is useful for testing whether group means differ. If they do not, \(\hat{\text{BSS}}\) should be small relative to \(\hat{\text{WSS}}\). The equality holds exactly because the cross-term \(\sum_{g,i}(\hat{y}_{ig}-\hat{m}_{Y}(g))(\hat{m}_{Y}(g)-\hat{m}_{Y})\) vanishes, since the within-group deviations sum to zero around their own group mean.
In R, aov() fits the ANOVA model and summary() returns the standard ANOVA table. The ANOVA table reports \(\hat{\text{BSS}}\) (“Sum Sq” for the group factor) and \(\hat{\text{WSS}}\) (“Sum Sq” for Residuals).
Code
# Fit ANOVA
reg_aov <- aov(y_hat ~ school_g)
summary(reg_aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## school_g 2 2541 1270.3 128 <2e-16 ***
## Residuals 3291 32654 9.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
For example, suppose we have \(G=3\) groups with the following data:
Data for Three Groups
| 4 |
8 |
5 |
| 5 |
7 |
6 |
| 6 |
9 |
7 |
Group means: \(\hat{m}_{Y}(A) = 5\), \(\hat{m}_{Y}(B) = 8\), \(\hat{m}_{Y}(C) = 6\). Grand mean: \(\hat{m}_{Y} = (5+8+6)/3 = 19/3 \approx 6.33\). Each group has \(n_{g} = 3\).
Between-group sum of squares: \[\begin{aligned}
\hat{\text{BSS}} &= 3(5-19/3)^2 + 3(8-19/3)^2 + 3(6-19/3)^2 \\
&= 3(16/9) + 3(25/9) + 3(1/9) = 48/9 + 75/9 + 3/9 = 14
\end{aligned}\] Within-group sum of squares: \[\begin{aligned}
\hat{\text{WSS}} &= [(4-5)^2+(5-5)^2+(6-5)^2] + [(8-8)^2+(7-8)^2+(9-8)^2] \\
& + [(5-6)^2+(6-6)^2+(7-6)^2] = 2+2+2 = 6
\end{aligned}\]
Code
# Data
y0_hat <- c(4, 5, 6, 8, 7, 9, 5, 6, 7)
g0 <- factor(c('A', 'A', 'A', 'B', 'B', 'B', 'C', 'C', 'C'))
# Verify by hand
m0_hat <- mean(y0_hat)
m0_g_hat <- aggregate(y0_hat, list(g0), mean)$x
n0_g <- aggregate(y0_hat, list(g0), length)$x
BSS <- sum(n0_g * (m0_g_hat - m0_hat)^2)
WSS <- sum(aggregate(y0_hat, list(g0), function(y_hat) sum((y_hat - mean(y_hat))^2))$x)
cbind(BSS, WSS)
## BSS WSS
## [1,] 14 6
The F-statistic
We have the decomposition; now we need a single number that tells us when the between-group piece is “large” relative to the within-group piece.
The ANOVA \(F\)-statistic is the ratio of average between-group variation to average within-group variation: \[\hat{F} = \frac{\hat{\text{BSS}}/(G-1)}{\hat{\text{WSS}}/(n-G)}.\] Writing \(Y_{ig}\) for the random outcome of individual \(i\) in group \(g\), the null is \(H_{0}: \mathbb{E}[Y_{i1}]=\mathbb{E}[Y_{i2}]=\dots=\mathbb{E}[Y_{iG}]\).
The ANOVA F-statistic is useful as the global “do any means differ?” test: one number and one \(p\)-value, regardless of how many groups you have. A large \(\hat{F}\) means the group means are more spread out than individual-level noise would predict.
Code
## MAIN EXAMPLE
F_fun <- function(y_hat,g){
m_y_hat <- mean(y_hat)
m_g_hat <- aggregate(y_hat, list(g), mean)$x
n_g <- aggregate(y_hat, list(g), length)$x
BSS <- sum(n_g * (m_g_hat - m_y_hat)^2)
WSS <- sum(aggregate(y_hat, list(g), function(y_hat) sum((y_hat - mean(y_hat))^2))$x)
n <- length(y_hat)
G <- nlevels(g)
F_obs <- (BSS/(G-1)) / (WSS/(n-G))
return(F_obs)
}
F_obs <- F_fun(y_hat=y_hat, g=school_g)
F_obs
## [1] 128.0237
To compute a \(p\)-value, we use a permutation test: shuffle the group labels, which keeps each group’s number of observations fixed, and recompute \(\hat{F}\) each time. Only large values of \(\hat{F}\) count against the null. As with KS and CVM in Comparing Two Groups, we reject at the \(5\%\) level when \(\hat{F}\) exceeds the \(95\)th percentile of the permutation distribution.
Code
## MAIN EXAMPLE
# Null distribution: shuffle schooling labels to break the Y--group pairing
F_perm <- vector(length=999)
for( b in seq_along(F_perm)){
g_perm <- sample(school_g, replace=FALSE)
F_perm[b] <- F_fun(y_hat=y_hat, g=g_perm)
}
# Permutation p-value
p_perm <- mean(F_perm >= F_obs)
p_perm
## [1] 0
# One-sided critical value: 95th percentile of the permutation distribution
F_crit <- quantile(F_perm, probs=.95)
hist(F_perm, breaks = 40,
border = NA,
freq = FALSE, main = NA, xlab = 'Null Permutation F-statistic',
xlim = range(c(0, F_perm, F_obs)))
title('Null Permutation F-test', font.main = 1)
abline(v = F_crit, col = rgb(0, 0, 0, .8), lwd = 2)
abline(v = F_obs, col = rgb(1, 0, 0, .8), lwd = 2)
legend('top', legend = c('Observed F', '95th percentile'),
col = c(rgb(1, 0, 0, .8), rgb(0, 0, 0, .8)), lwd = 2, bty='n')
Code
# No shuffle comes close to the observed F of about 128.
# So the permutation p-value is zero.
Pairing the observed outcomes y_hat with a shuffled vector of group labels breaks the link between outcome and group: under the null, group labels are exchangeable, so this approximates the distribution of \(\hat{F}\) when no group effect exists.
Under the null and additional parametric assumptions, \(\hat{F}\) follows an \(F_{G-1,\;n-G}\) distribution. The parametric \(F_{G-1,\;n-G}\) null distribution requires three assumptions: independent observations, normal errors within each group, and equal variances across groups (homoscedasticity). R reports \(p\)-values from this theoretical distribution by default.
Here is the same test for the sampled states from Describing Groups.
Code
# Fit ANOVA
summary(aov(murder_hat ~ region))
## Df Sum Sq Mean Sq F value Pr(>F)
## region 3 357.3 119.09 10.07 8.64e-05 ***
## Residuals 31 366.5 11.82
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# F-statistic
F_obs <- F_fun(y_hat=murder_hat, g=region)
F_obs
## [1] 10.0743
# Null distribution: shuffle region labels to break the Y--group pairing
F_perm <- vector(length=999)
for( b in seq_along(F_perm)){
g_perm <- sample(region, replace=FALSE)
F_perm[b] <- F_fun(y_hat=murder_hat, g=g_perm)
}
# Permutation p-value
p_perm <- mean(F_perm >= F_obs)
p_perm
## [1] 0
# One-sided critical value: 95th percentile of the permutation distribution
F_crit <- quantile(F_perm, probs=.95)
hist(F_perm, breaks = 40,
border = NA,
freq = FALSE, main = NA, xlab = 'Null Permutation F-statistic',
xlim = range(c(0, F_perm, F_obs)))
title('Null Permutation F-test', font.main = 1)
abline(v = F_crit, col = rgb(0, 0, 0, .8), lwd = 2)
abline(v = F_obs, col = rgb(1, 0, 0, .8), lwd = 2)
legend('top', legend = c('Observed F', '95th percentile'),
col = c(rgb(1, 0, 0, .8), rgb(0, 0, 0, .8)), lwd = 2, bty='n')
The figure below draws the theoretical \(F_{G-1,\;n-G}\) density over the permutation distribution for the sampled states, from the callout above.
Code
# Degrees of freedom for the regional data
n <- length(murder_hat)
G <- nlevels(region)
p_param <- 1 - pf(F_obs, G-1, n-G)
# Compare parametric and permutation p-values
cbind(F_obs, p_param, p_perm)
## F_obs p_param p_perm
## [1,] 10.0743 8.639652e-05 0
# Overlay: permutation histogram vs. theoretical F density
hist(F_perm, breaks = 40,
border = NA,
freq = FALSE, main = NA, xlab = 'Null F-statistic',
xlim = range(c(0, F_perm, F_obs)))
title('Permutation vs. Theoretical Null Distribution', font.main = 1)
F_grid <- seq(0, max(c(F_perm, F_obs)), length = 200)
lines(F_grid, df(F_grid, G-1, n-G), col = rgb(0, 0, 1, .8), lwd = 1, lty = 2)
abline(v = F_obs, col = rgb(1, 0, 0, .8), lwd = 2)
legend('top', legend = c('Null permutation', paste0('F(', G-1, ', ', n-G, ')'), 'Observed F'),
fill = c(grey(.83), NA, NA), border = NA, # grey(.83) matches the default histogram fill
lwd = c(NA, 1, 2), lty = c(NA, 2, 1), col = c(NA, rgb(0, 0, 1, .8), rgb(1, 0, 0, .8)), bty = 'n')
If these assumptions are doubtful, the Kruskal-Wallis test below is a safer alternative. As a quick check, compare side-by-side boxplots: if the spread differs across groups or the distributions show strong skewness or kurtosis, be cautious with the parametric test.
For the three-group example above, \(\hat{F}=(14/2)/(6/6)=7\). With \(G-1=2\) and \(n-G=6\) degrees of freedom, we compare it to an \(F_{2,6}\) distribution: \(p = 1 - F_{2,6}(7) \approx 0.027\).
Code
# Data (the three-group example from the decomposition above)
y0_hat <- c(4, 5, 6, 8, 7, 9, 5, 6, 7)
g0 <- factor(c('A', 'A', 'A', 'B', 'B', 'B', 'C', 'C', 'C'))
# Computation
summary(aov(y0_hat ~ g0))
## Df Sum Sq Mean Sq F value Pr(>F)
## g0 2 14 7 7 0.027 *
## Residuals 6 6 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Verify by hand
n0 <- length(y0_hat)
G0 <- nlevels(g0)
m0_hat <- mean(y0_hat)
m0_g_hat <- aggregate(y0_hat, list(g0), mean)$x
n0_g <- aggregate(y0_hat, list(g0), length)$x
BSS <- sum(n0_g * (m0_g_hat - m0_hat)^2)
WSS <- sum(aggregate(y0_hat, list(g0), function(y_hat) sum((y_hat - mean(y_hat))^2))$x)
Fstat <- (BSS/(G0-1)) / (WSS/(n0-G0))
pval <- 1 - pf(Fstat, G0-1, n0-G0)
cbind(BSS, WSS, Fstat, pval)
## BSS WSS Fstat pval
## [1,] 14 6 7 0.027
Equal Distributions
ANOVA tests whether group means differ; sometimes the entire distributions differ in ways the mean would miss.
Kruskal-Wallis Test
When data are skewed, heavy-tailed, or only ordinal, a test built on means may overlook differences that show up in the rank order.
The Kruskal-Wallis test (Kruskal and Wallis 1952) is a nonparametric test for equal distributions across \(G\) groups, testing \(H_{0}: F_{1} = F_{2} = \dots = F_{G}\) against \(H_{A}:\) at least one \(F_{g}\) differs (where \(F_{g}\) is the continuous distribution of group \(g=1,...G\)). It uses the average ranks within each group rather than the raw observations, so it only requires the data to be ordinal. With two groups it reduces to the Mann-Whitney U (equivalently Wilcoxon rank-sum) test (Wilcoxon 1945; Mann and Whitney 1947).
The Kruskal-Wallis test is useful when the equal-variance normal-error assumption is implausible, as with income, response times, or other heavy-tailed outcomes. Rank-based tests need only the orderings of the data. A significant result tells us that some distribution differs, but not which one; the global rejection has to be followed up with pairwise tests.
To conduct the test, first denote individuals \(i=1,...n\) with overall ranks \(\hat{r}_{1},....\hat{r}_{n}\). Each individual belongs to group \(g=1,...G\), and each group \(g\) has \(n_{g}\) individuals with average rank \(\bar{r}_{g} = \sum_{i} \hat{r}_{i} /n_{g}\). The Kruskal Wallis statistic is \[\begin{aligned}
\hat{\text{KW}} &= (n-1) \frac{\sum_{g=1}^{G} n_{g}( \bar{r}_{g} - \bar{r} )^2 }{\sum_{i=1}^{n} ( \hat{r}_{i} - \bar{r} )^2},
\end{aligned}\] where \(\bar{r} = \frac{n+1}{2}\) is the grand mean rank.
In the special case with only two groups, the Kruskal Wallis test reduces to the Mann–Whitney U test (also known as the Wilcoxon rank-sum test). In this case, we can write the hypotheses in terms of individual outcomes in each group, \(Y_{i}\) in one group and \(Y_{j}\) in the other: \(H_{0}: \text{Prob}(Y_{i} > Y_{j})=\text{Prob}(Y_{j} > Y_{i})\) versus \(H_{A}: \text{Prob}(Y_{i} > Y_{j}) \neq \text{Prob}(Y_{j} > Y_{i})\). The corresponding test statistic is \[\begin{aligned}
\hat{U} &= \min(\hat{U}_{1}, \hat{U}_{2}) \\
\hat{U}_{g} &= \sum_{i~\text{in}~g}\sum_{j~\text{in}~-g}
\Bigl[\mathbf 1( \hat{y}_{i} > \hat{y}_{j}) + \tfrac12\mathbf 1(\hat{y}_{i} = \hat{y}_{j})\Bigr].
\end{aligned}\]
Interpretation template:
- If the Kruskal-Wallis p-value is large, keep the global null and stop.
- If the p-value is small, report that at least one group differs.
- Report both statistical and practical importance (distribution plots + effect size summaries), not only p-values.
Code
## MAIN EXAMPLE
kw <- kruskal.test(y_hat ~ school_g)
kw
##
## Kruskal-Wallis rank sum test
##
## data: y_hat by school_g
## Kruskal-Wallis chi-squared = 259.3, df = 2, p-value < 2.2e-16
The test above does not tell us which groups differ. To answer which, use a post-hoc pairwise test together with a multiple-comparisons correction. Here, pairwise.wilcox.test() runs a Mann-Whitney test on each pair of groups, and the Holm correction (Holm 1979) adjusts the \(p\)-values for the number of pairs tested.
Code
pairwise.wilcox.test(y_hat, school_g, p.adjust.method='holm')
##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: y_hat and school_g
##
## <12 12
## 12 2.6e-10 -
## >12 < 2e-16 < 2e-16
##
## P value adjustment method: holm
# All three adjusted p-values are tiny, so every pair of groups differs.
Looking back at the boxplot in Equal Medians, which two groups are closest?
Here are the same tests for the sampled states from Describing Groups.
Code
kw <- kruskal.test(murder_hat ~ region)
kw
##
## Kruskal-Wallis rank sum test
##
## data: murder_hat by region
## Kruskal-Wallis chi-squared = 16.394, df = 3, p-value = 0.0009415
pairwise.wilcox.test(murder_hat, region, p.adjust.method='holm')
##
## Pairwise comparisons using Wilcoxon rank sum exact test
##
## data: murder_hat and region
##
## Northeast North Central West
## North Central 0.582 - -
## West 0.075 0.582 -
## South 0.001 0.031 0.075
##
## P value adjustment method: holm
Looking back at the boxplot in the Equal Medians callout, which region stands out?
Worked example: \(G=3\) groups with values \(A=\{1,2,3\}\), \(B=\{4,5,6\}\), \(C=\{7,8,9\}\).
Pooling and ranking gives the ranks \(1,\ldots,9\) in order. Group rank averages: \(\bar{r}_{A}=2\), \(\bar{r}_{B}=5\), \(\bar{r}_{C}=8\). Grand mean rank: \(\bar{r}=(n+1)/2=5\). Plug into \(\hat{\text{KW}}\) with \(n=9\), \(n_{g}=3\): \[
\hat{\text{KW}} = (n-1) \cdot \frac{3(2-5)^2 + 3(5-5)^2 + 3(8-5)^2}{\sum_{i=1}^{9}(\hat{r}_{i} - 5)^2}
= 8 \cdot \frac{27 + 0 + 27}{60} = 7.2.
\]
Code
y0_hat <- c(1,2,3, 4,5,6, 7,8,9)
g0 <- factor(rep(c('A','B','C'), each=3))
kruskal.test(y0_hat ~ g0)
##
## Kruskal-Wallis rank sum test
##
## data: y0_hat by g0
## Kruskal-Wallis chi-squared = 7.2, df = 2, p-value = 0.02732
For the two-group case, the Mann-Whitney statistic counts pairwise wins. Compare \(A=\{1,3,5\}\) to \(B=\{2,4,6\}\): out of \(3\times 3 = 9\) pairs, \(A\) wins three (\((3,2), (5,2), (5,4)\)), so \(\hat{U}_{A}=3\), \(\hat{U}_{B}=6\), giving \(\hat{U}=\min(3,6)=3\).
Code
a0_hat <- c(1,3,5); b0_hat <- c(2,4,6)
wilcox.test(a0_hat, b0_hat)
##
## Wilcoxon rank sum exact test
##
## data: a0_hat and b0_hat
## W = 3, p-value = 0.7
## alternative hypothesis: true location shift is not equal to 0
The Mann-Whitney test compares medians when the two groups have the same shape, so that one distribution is a shifted copy of the other. Otherwise, it tests whether a draw from one group is as likely to exceed a draw from the other as the reverse, and it can reject even when the two medians are equal.
Same shape. Suppose hourly wages in two labor markets have the same right-skewed shape, but every wage in the second market is \(2\) dollars higher. The population medians are \(15\) and \(17\). Consider the probability that a worker drawn from the second market earns more than a worker drawn from the first. Its sample analog is the share of pairs that the second group wins, \(\hat{U}_{2}/(n_{1}n_{2})\) in the notation above. With conf.int=TRUE, wilcox.test also reports a shift estimate, the median of all differences \(\hat{y}_{i}-\hat{y}_{j}\) between a worker \(i\) in the second group and a worker \(j\) in the first. Because the shapes match, the shift estimate and the difference in sample medians both estimate the \(2\) dollar gap. So here, a Mann-Whitney rejection is evidence that the medians differ.
Code
# Same right-skewed shape, with every wage in market 2 shifted up by 2 dollars
n_g <- 2000
y1_sim <- 15 + 3 * (rexp(n_g) - log(2)) # population median 15
y2_sim <- 17 + 3 * (rexp(n_g) - log(2)) # population median 17
# Difference in sample medians
median(y2_sim) - median(y1_sim)
## [1] 1.934255
# Share of pairs that market 2 wins
mean(outer(y2_sim, y1_sim, '>'))
## [1] 0.7364205
# A market 2 worker earns more than a market 1 worker about 74% of the time.
# Mann-Whitney test and its shift estimate
mw <- wilcox.test(y2_sim, y1_sim, conf.int=TRUE)
c(shift=unname(mw[['estimate']]), p_value=mw[['p.value']])
## shift p_value
## 1.930319 0.000000
# The market 2 ECDF is the market 1 ECDF moved right by 2
cols <- c(rgb(.8, 0, 0, .5), rgb(0, 0, .8, .5))
F1_hat <- ecdf(y1_sim)
F2_hat <- ecdf(y2_sim)
plot(F1_hat, col=cols[1], lwd=2, verticals=TRUE, do.points=FALSE, col.01line=NA,
xlim=c(5, 35), main=NA,
xlab='Hourly wage', ylab='Proportion of workers')
plot(F2_hat, add=TRUE, col=cols[2], lwd=2, verticals=TRUE, do.points=FALSE, col.01line=NA)
abline(h=.5, col=grey(.5), lty=3) # each curve crosses this line at its median
legend('bottomright', legend=c('Market 1', 'Market 2'),
lty=1, lwd=2, col=cols, bty='n')
Different shapes. Now suppose wages in the first market are symmetric around \(15\), while the second market keeps the right-skewed shape from above, also with median \(15\). The medians are equal, but the second market has no wages below about \(12.92\), while the first market has a long left tail. When both workers earn less than \(15\), the worker from the second market earns more about \(70\%\) of the time. Overall, the second market wins more than half of all pairs. So the Mann-Whitney test can reject, with a shift estimate above zero, even though the medians are equal.
Code
# Same population median of 15, different shapes
y1_sim <- 15 + 3 * rnorm(n_g) # symmetric
y2_sim <- 15 + 3 * (rexp(n_g) - log(2)) # right-skewed, no wages below 12.92
# Difference in sample medians is close to zero
median(y2_sim) - median(y1_sim)
## [1] -0.1661474
# Yet market 2 wins more than half of the pairs
mean(outer(y2_sim, y1_sim, '>'))
## [1] 0.5432795
# The sample share is about 54%, and the population share is about 56%.
# So the Mann-Whitney test rejects
mw <- wilcox.test(y2_sim, y1_sim, conf.int=TRUE)
c(shift=unname(mw[['estimate']]), p_value=mw[['p.value']])
## shift p_value
## 4.178188e-01 2.132699e-06
# The shift estimate is about 0.42, and the population shift is about 0.55.
# Both are well above zero, even though the medians are equal.
# The ECDFs cross at the median but differ below it
F1_hat <- ecdf(y1_sim)
F2_hat <- ecdf(y2_sim)
plot(F1_hat, col=cols[1], lwd=2, verticals=TRUE, do.points=FALSE, col.01line=NA,
xlim=c(5, 35), main=NA,
xlab='Hourly wage', ylab='Proportion of workers')
plot(F2_hat, add=TRUE, col=cols[2], lwd=2, verticals=TRUE, do.points=FALSE, col.01line=NA)
abline(h=.5, col=grey(.5), lty=3) # both curves cross this line at 15
legend('bottomright', legend=c('Market 1', 'Market 2'),
lty=1, lwd=2, col=cols, bty='n')
Categorical ANOVA
When the outcome is a factor, we can still ask whether the groups spread their observations across the categories in the same proportions. Let \(\hat{p}(y)\) be the proportion of all \(n\) observations in category \(y\), and \(\hat{p}(y \mid g)\) the proportion of the \(n_{g}\) observations in group \(g\). The categorical ANOVA measures variation with one minus the Herfindahl index, and splits it into between-group and within-group pieces: \[\begin{aligned}
\underbrace{\frac{n}{2}\left(1-\hat{H}\right)}_{\hat{\text{TSS}}}
&=
\underbrace{\frac{1}{2}\sum_{g=1}^{G} n_{g} \sum_{y}\left(\hat{p}(y \mid g)-\hat{p}(y)\right)^2}_{\hat{\text{BSS}}}
\;+\;
\underbrace{\sum_{g=1}^{G} \frac{n_{g}}{2}\left(1-\hat{H}(g)\right)}_{\hat{\text{WSS}}},
\end{aligned}\] where \(\hat{H}=\sum_{y}\hat{p}(y)^2\) is the Herfindahl index of the pooled data and \(\hat{H}(g)=\sum_{y}\hat{p}(y \mid g)^2\) is the Herfindahl index within group \(g\).
The categorical ANOVA is useful when the outcome is an unordered factor, such as occupation, brand, or vote choice, where there are no means to compare and no ranks to sort. Recall that \(1-\hat{H}\) is zero when every observation falls in one category and grows as the categories become equally common, so it plays the role of the variance. The pieces are also ordinary ANOVA sums of squares. Make one 0/1 indicator for each category and compute its \(\hat{\text{TSS}}\), \(\hat{\text{BSS}}\), and \(\hat{\text{WSS}}\) as in the decomposition above. Adding each piece up across the categories and halving it gives the three pieces here. With a binary outcome, the two indicators have the same sums of squares, so the categorical ANOVA is exactly the ANOVA of a single 0/1 variable. The share \(\hat{\text{BSS}}/\hat{\text{TSS}}\) reports the fraction of the variation in \(\hat{y}\) that lies between groups.
For example, cut a cardinal outcome into brackets and sort each group’s observations into them. The brackets are an ordered factor, but the categorical ANOVA ignores that order, so it treats them as it would any unordered factor.
Code
## MAIN EXAMPLE
# Outcome: hourly wage in four brackets, an ordered factor
# The brackets are up to 4 dollars an hour, 4 to 6, 6 to 8, and above 8.
wage_bracket <- cut(y_hat, breaks=c(0, 4, 6, 8, Inf), ordered_result=TRUE)
tab <- table(school_g, wage_bracket)
tab
## wage_bracket
## school_g (0,4] (4,6] (6,8] (8,Inf]
## <12 542 430 236 145
## 12 365 375 251 197
## >12 126 177 194 256
# Proportions within each schooling group: each bar sums to 1
p_bracket_school <- prop.table(tab, 1)
cols <- hcl.colors(4, alpha=.45)
par(mar=c(4, 4, 3.5, 1))
barplot(t(p_bracket_school), col=cols, border=NA,
xlab='Years of schooling', ylab='Proportion of workers', main=NA)
legend('top', inset=c(0, -.15), xpd=TRUE, ncol=4, bty='n',
fill=cols, border=NA, legend=levels(wage_bracket))
Code
# Workers with more schooling are less often in the lowest bracket.
# They are more often in the highest.
The function below computes the three pieces and the between-group share.
Code
## MAIN EXAMPLE
# Categorical ANOVA: variation is measured by one minus the Herfindahl index
cat_anova <- function(y_hat, g) {
# Total variation, from the pooled proportions
n <- length(y_hat)
p_hat <- table(y_hat) / n
TSS <- n / 2 * (1 - sum(p_hat^2))
# Within-group variation, from each group's own proportions
WSS_g <- tapply(y_hat, g, function(y_g) {
p_g_hat <- table(y_g) / length(y_g)
return(length(y_g) / 2 * (1 - sum(p_g_hat^2)))
})
WSS <- sum(WSS_g)
# Between-group variation, from the decomposition
BSS <- TSS - WSS
return(c(BSS=BSS, WSS=WSS, TSS=TSS, share=BSS / TSS))
}
cat_anova(wage_bracket, school_g)
## BSS WSS TSS share
## 3.017363e+01 1.183800e+03 1.213974e+03 2.485526e-02
# About 2.5% of the variation in wage brackets lies between schooling groups.
Suppose we record how \(n=8\) workers commute, \(n_{g}=4\) from each of two neighborhoods. Neighborhood \(A\) has car, car, car, bus, and neighborhood \(B\) has car, bus, bike, bike.
Pooled, the proportions are \(\hat{p}(\text{car})=4/8\), \(\hat{p}(\text{bus})=2/8\), and \(\hat{p}(\text{bike})=2/8\), so \(\hat{H}=1/4+1/16+1/16=3/8\) and \(\hat{\text{TSS}}=\frac{8}{2}(1-3/8)=5/2\). Within \(A\), \(\hat{H}(A)=(3/4)^2+(1/4)^2=5/8\), and within \(B\), \(\hat{H}(B)=(1/4)^2+(1/4)^2+(2/4)^2=3/8\). So \(\hat{\text{WSS}}=\frac{4}{2}(1-5/8)+\frac{4}{2}(1-3/8)=3/4+5/4=2\). For \(\hat{\text{BSS}}\), each neighborhood’s proportions differ from the pooled ones by \(1/4\) in two categories and by \(0\) in the third, so \(\hat{\text{BSS}}=\frac{1}{2}[4(1/16+0+1/16)+4(1/16+0+1/16)]=1/2\). The pieces add up, \(5/2=1/2+2\), and the between-group share is \(\hat{\text{BSS}}/\hat{\text{TSS}}=1/5\).
Code
# Commute mode of 8 workers, 4 from each neighborhood
y0_hat <- factor(c('car', 'car', 'car', 'bus', 'car', 'bus', 'bike', 'bike'))
g0 <- factor(c('A', 'A', 'A', 'A', 'B', 'B', 'B', 'B'))
table(g0, y0_hat)
## y0_hat
## g0 bike bus car
## A 0 1 3
## B 2 1 1
cat_anova(y0_hat, g0)
## BSS WSS TSS share
## 0.5 2.0 2.5 0.2
As with the ANOVA \(F\)-statistic, we get a \(p\)-value by shuffling the group labels and recomputing the share each time.
Code
# Observed between-group share
share_obs <- cat_anova(wage_bracket, school_g)[['share']]
# Null distribution: shuffle schooling labels to break the outcome--group pairing
share_perm <- rep(NA, 999)
for (b in seq_along(share_perm)) {
g_perm <- sample(school_g, replace=FALSE)
share_perm[b] <- cat_anova(wage_bracket, g_perm)[['share']]
}
# Permutation p-value
p_perm <- mean(share_perm >= share_obs)
p_perm
## [1] 0
# One-sided critical value: 95th percentile of the permutation distribution
share_crit <- quantile(share_perm, probs=.95)
hist(share_perm, breaks=40,
border=NA,
freq=FALSE, main=NA, xlab='Null Permutation Between-Group Share',
xlim=range(c(0, share_perm, share_obs)))
title('Null Permutation Categorical ANOVA', font.main=1)
abline(v=share_crit, col=rgb(0, 0, 0, .8), lwd=2)
abline(v=share_obs, col=rgb(1, 0, 0, .8), lwd=2)
legend('top', legend=c('Observed share', '95th percentile'),
col=c(rgb(1, 0, 0, .8), rgb(0, 0, 0, .8)), lwd=2, bty='n')
Because the brackets are ordered, the Kruskal-Wallis test also applies.
Code
# Kruskal-Wallis ranks the brackets, so it uses their order
kruskal.test(wage_bracket ~ school_g)
##
## Kruskal-Wallis rank sum test
##
## data: wage_bracket by school_g
## Kruskal-Wallis chi-squared = 241.42, df = 2, p-value < 2.2e-16
# Both this test and the categorical ANOVA reject.
# The schooling groups differ along one direction, from low to high brackets.
Ranks concentrate on differences along one direction, from low to high brackets. The categorical ANOVA spreads its attention over every way the proportions could differ, so a rank test is usually more powerful when the categories are ordered. The categorical ANOVA can still detect differences that ranks miss, such as one group concentrated in the middle brackets and another split between the two extremes.
Here is the same analysis for the sampled states from Describing Groups. The urban population bands from Ordered Factors sort each region’s states into mostly rural, mixed, and mostly urban.
Code
# Outcome: urban population in three bands, an ordered factor, for the sampled states
urban_band <- cut(USArrests[state_id, 'UrbanPop'], breaks=c(0, 50, 75, 100),
labels=c('mostly rural', 'mixed', 'mostly urban'), ordered_result=TRUE)
tab <- table(region, urban_band)
tab
## urban_band
## region mostly rural mixed mostly urban
## Northeast 1 1 3
## North Central 1 7 0
## West 0 3 5
## South 5 7 2
# Proportions within each region: each bar sums to 1
p_band_region <- prop.table(tab, 1)
cols <- c(rgb(.8, 0, 0, .5), rgb(0, .8, 0, .5), rgb(0, 0, .8, .5))
par(mar=c(4, 4, 3.5, 1))
barplot(t(p_band_region), col=cols, border=NA,
xlab='Region', ylab='Proportion of states', main=NA)
legend('top', inset=c(0, -.15), xpd=TRUE, ncol=3, bty='n',
fill=cols, border=NA, legend=levels(urban_band))
Code
# The South and North Central lean rural.
# The Northeast and West lean urban.
The function from above computes the three pieces and the between-group share.
Code
cat_anova(urban_band, region)
## BSS WSS TSS share
## 2.3785714 8.3642857 10.7428571 0.2214096
# About 22% of the variation in urban bands lies between regions.
As above, we get a \(p\)-value by shuffling the group labels and recomputing the share each time.
Code
# Observed between-group share
share_obs <- cat_anova(urban_band, region)[['share']]
# Null distribution: shuffle region labels to break the outcome--group pairing
share_perm <- rep(NA, 999)
for (b in seq_along(share_perm)) {
g_perm <- sample(region, replace=FALSE)
share_perm[b] <- cat_anova(urban_band, g_perm)[['share']]
}
# Permutation p-value
p_perm <- mean(share_perm >= share_obs)
p_perm
## [1] 0.01801802
# One-sided critical value: 95th percentile of the permutation distribution
share_crit <- quantile(share_perm, probs=.95)
hist(share_perm, breaks=40,
border=NA,
freq=FALSE, main=NA, xlab='Null Permutation Between-Group Share',
xlim=range(c(0, share_perm, share_obs)))
title('Null Permutation Categorical ANOVA', font.main=1)
abline(v=share_crit, col=rgb(0, 0, 0, .8), lwd=2)
abline(v=share_obs, col=rgb(1, 0, 0, .8), lwd=2)
legend('topright', legend=c('Observed share', '95th percentile'),
col=c(rgb(1, 0, 0, .8), rgb(0, 0, 0, .8)), lwd=2, bty='n')
Because the bands are ordered, the Kruskal-Wallis test also applies.
Code
# Kruskal-Wallis ranks the bands, so it uses their order
kruskal.test(urban_band ~ region)
##
## Kruskal-Wallis rank sum test
##
## data: urban_band by region
## Kruskal-Wallis chi-squared = 9.284, df = 3, p-value = 0.02574
# The categorical ANOVA gives the smaller p-value, about 0.018 against 0.026.
# The regions differ partly along one direction, from rural to urban.
# Ranks concentrate on that direction.
# In this sample, 7 of North Central's 8 states are in the middle band.
# That pattern contributes more to the categorical ANOVA than any other region.
# It contributes the least to the Kruskal-Wallis statistic.
For the sampled states in the callout above, make a binary outcome that marks the mostly urban states. Check that the categorical ANOVA of this outcome gives the same \(\hat{\text{BSS}}\) and \(\hat{\text{WSS}}\) as the ordinary ANOVA of the 0/1 variable.
Code
# 1 for mostly urban states, 0 otherwise
d_hat <- as.numeric(urban_band == 'mostly urban')
cat_anova(factor(d_hat), region)
## BSS WSS TSS share
## 2.353571 4.789286 7.142857 0.329500
summary(aov(d_hat ~ region))
## Df Sum Sq Mean Sq F value Pr(>F)
## region 3 2.354 0.7845 5.078 0.00563 **
## Residuals 31 4.789 0.1545
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Exercises
Comment the script you wrote for this chapter, then restart R and check that it runs from a clean session, then check the script with AI as explained in Working with AI. Write three sentences from memory on the main statistical idea of this chapter, and ask the assistant what is wrong, vague, or missing. Finish with your own questions about whatever you found hardest.
Load Wages1 from the Ecdat package. Draw a boxplot of wage by sex, then add an ECDF of wage for each sex on the same axes. Compute the mean and median wage for each sex with aggregate(). Does one sex’s ECDF sit to the right of the other’s over the whole range, or do the curves cross?
An ANOVA \(F\)-test rejects the null that all group means are equal. Does this tell you which group is different, or only that at least one differs?
Using the USArrests dataset with state.region as the grouping variable, compute the between-group sum of squares (\(\hat{\text{BSS}}\)) and within-group sum of squares (\(\hat{\text{WSS}}\)) for Assault by hand (using aggregate()). Verify your \(\hat{F}\) statistic matches the output of summary(aov(Assault ~ state.region, data = ...)).
Recall
This chapter described several groups with ECDFs and a table of group means. We compared their medians with bootstrap notches on boxplots. We then tested whether the differences reflected more than chance. ANOVA splits variation into between-group and within-group pieces and combines them into an \(F\)-statistic. A permutation test got its \(p\)-value by shuffling the group labels. In the three-group example with \(A=\{4,5,6\}, B=\{8,7,9\}, C=\{5,6,7\}\), we found \(\hat{\text{BSS}}=14\) and \(\hat{\text{WSS}}=6\) by hand. The function aov(y0_hat ~ g0) gives the same numbers. The Kruskal-Wallis test replaced outcomes with ranks when the equal-variance normal-error assumptions were doubtful. The categorical ANOVA measured the variation of a factor outcome with one minus the Herfindahl index; for the commutes in two neighborhoods it gave \(\hat{\text{TSS}}=5/2\), \(\hat{\text{BSS}}=1/2\), and \(\hat{\text{WSS}}=2\) by hand.
Holm, Sture. 1979.
“A Simple Sequentially Rejective Multiple Test Procedure.” Scandinavian Journal of Statistics 6: 65–70.
https://www.jstor.org/stable/4615733.
Kruskal, William H., and W. Allen Wallis. 1952.
“Use of Ranks in One-Criterion Variance Analysis.” Journal of the American Statistical Association 47 (260): 583–621.
https://doi.org/10.1080/01621459.1952.10483441.
Mann, H. B., and D. R. Whitney. 1947.
“On a Test of Whether One of Two Random Variables Is Stochastically Larger Than the Other.” The Annals of Mathematical Statistics 18 (1): 50–60.
https://doi.org/10.1214/aoms/1177730491.
McGill, Robert, John W. Tukey, and Wayne A. Larsen. 1978.
“Variations of Box Plots.” The American Statistician 32 (1): 12–16.
https://doi.org/10.1080/00031305.1978.10479236.
Tukey, John Wilder. 1977. Exploratory Data Analysis. Vol. 2. Springer.
Wilcoxon, Frank. 1945.
“Individual Comparisons by Ranking Methods.” Biometrics Bulletin 1 (6): 80–83.
https://doi.org/10.2307/3001968.
