22  Comparing Multiple Groups

So far we have used regression with cardinal explanatory variables. We now turn to factor variables that label which group each observation belongs to: regions, industries, treatment arms. The chapter builds up from the ANOVA decomposition of variation between and within groups to its nonparametric Kruskal-Wallis analog. The next chapter rewrites these comparisons as regressions and develops permutation inference for a factor variable.

Factor Variables

So far we have used cardinal explanatory variables where \(4-3=2-1\); many useful predictors are instead labels of group membership that do not arithmetic this way.

ImportantKey Definition

A factor variable labels which group each observation belongs to. Each level of the factor can be encoded as an indicator variable (or dummy) \(\hat{D}_{c} = \mathbf{1}(\text{Factor} = c)\), which equals \(1\) when the observation is in group \(c\) and \(0\) otherwise.

Factor variables are useful when group identity itself predicts the outcome (region, industry, treatment arm) and we want to compare group means within a regression framework. They split into two flavors:

  • Ordered: refers to ordinal data. The difference between units means something, but not always the same thing. For example, \(4th - 3rd \neq 2nd - 1st\).
  • Unordered: refers to categorical data. The difference between units is meaningless. For example, \(D-C=?\)

One common case is if you have observations of individuals over time periods, then you may have two factor variables. An unordered factor that indicates who an individual is; for example \(\hat{D}_{i}=\mathbf{1}( \text{Individual} = i)\), and an ordered factor that indicates the time period; for example \(\hat{D}_{t}=\mathbf{1}( \text{Time} \in [\text{month}~ t, \text{month}~ t+1) )\). There are many other cases you see factor variables, including spatial ID’s in purely cross sectional data.

Be careful not to handle categorical data as if they were cardinal. E.g., generate city data with Leipzig=1, Lausanne=2, LosAngeles=3, … and then include city as if it were a cardinal number (that’s a big no-no). The same applied to ordinal data; PopulationLeipzig=2, PopulationLausanne=3, PopulationLosAngeles=1.

In this chapter, we use state.region (Northeast, South, North Central, West) as an unordered factor to group the USArrests data by US Census region.

Code 
# Attach the four-level region factor to USArrests; the levels
# are the values the categorical variable can take.
dat <- USArrests
dat$Region <- state.region
levels(dat$Region)
## [1] "Northeast"     "South"         "North Central" "West"

Multiple Groups

With multiple groups, begin with a summary figure (such as a boxplot) and then run a global test of whether at least one group differs.

Code 
# Step 1: Visual summary
boxplot(Murder ~ Region, data = dat,
    col = hcl.colors(4, alpha = .45),
    ylab = 'Murder arrests (per 100k)', xlab = 'Region', main = NA)
title('Murder rates by region', font.main = 1)

22.1 Equal Means (ANOVA)

In Comparing 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 \(Y\) must come from differences between groups rather than within them. We can split each piece out.

ImportantKey Definition

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}_{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{TSS}} &= \underbrace{\sum_{g=1}^{G} n_g (\hat{M}_{g}-\hat{M}_{Y})^2}_{\hat{BSS}} \;+\; \underbrace{\sum_{g=1}^{G}\sum_{i=1}^{n_g}(\hat{Y}_{ig}-\hat{M}_{g})^2}_{\hat{WSS}}, \end{aligned}\] where \(\hat{BSS}\) measures variation in group means around the grand mean and \(\hat{WSS}\) measures variation of observations around their own group mean.

The decomposition is useful as the foundation for any test of “do group means differ?”: if they do not, \(\hat{BSS}\) should be small relative to \(\hat{WSS}\). The equality holds exactly because the cross-term \(\sum_{g,i}(\hat{Y}_{ig}-\hat{M}_{g})(\hat{M}_{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{BSS}\) (“Sum Sq” for Region), \(\hat{WSS}\) (“Sum Sq” for Residuals).

Code 
# Fit ANOVA
fit_aov <- aov(Murder ~ Region, data = dat)
fit_aov
## Call:
##    aov(formula = Murder ~ Region, data = dat)
## 
## Terms:
##                   Region Residuals
## Sum of Squares  391.2357  538.3171
## Deg. of Freedom        3        46
## 
## Residual standard error: 3.420898
## Estimated effects may be unbalanced

For example, suppose we have \(G=3\) groups with the following data:

Group A Group B Group C
4 8 5
5 7 6
6 9 7

Group means: \(\hat{M}_A = 5\), \(\hat{M}_B = 8\), \(\hat{M}_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{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{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
Y <- c(4, 5, 6, 8, 7, 9, 5, 6, 7)
g <- factor(c('A', 'A', 'A', 'B', 'B', 'B', 'C', 'C', 'C'))

# Verify by hand
n <- length(Y)
G <- nlevels(g)

M_Y <- mean(Y)
M_g <- aggregate(Y, list(g), mean)$x
n_g <- aggregate(Y, list(g), length)$x

BSS <- sum(n_g * (M_g - M_Y)^2)
WSS <- sum(aggregate(Y, list(g), function(y) sum((y - mean(y))^2))$x)

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.

ImportantKey Definition

The ANOVA \(F\)-statistic is the ratio of average between-group variation to average within-group variation: \[\hat{F} = \frac{\hat{BSS}/(G-1)}{\hat{WSS}/(n-G)}.\] Under the null \(H_0: \mu_1=\mu_2=\dots=\mu_G\) (with the usual parametric assumptions) it follows an \(F_{G-1, n-G}\) distribution; we can also get a \(p\)-value nonparametrically by bootstrapping the labels.

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. The numerator is the average between-group variation per degree of freedom, the denominator is the average within-group variation, and a large \(\hat{F}\) means group means are more spread out than individual-level noise would predict.

Code 
F_fun <- function(Y,g){
    M_Y <- mean(Y)
    M_g <- aggregate(Y, list(g), mean)$x
    n_g <- aggregate(Y, list(g), length)$x

    BSS <- sum(n_g * (M_g - M_Y)^2)
    WSS <- sum(aggregate(Y, list(g), function(y) sum((y - mean(y))^2))$x)

    n <- length(Y)
    G <- nlevels(g)
    F_obs <- (BSS/(G-1)) / (WSS/(n-G))
}
F_obs <- F_fun(Y=dat$Murder, g=dat$Region)
F_obs
## [1] 11.14389

To compute a \(p\)-value, we can use permutation or bootstrap methods to build a null distribution, just as in previous chapters. The Permutation Tests section in the next chapter instead shuffles the observed labels without replacement, preserving the number of states in each region.

Pairing the observed Murder values with a resampled Region vector 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.

Code 
# Pairs bootstrap: resample group labels to break the Y--group pairing
n <- nrow(dat)

F_nullboot <- vector(length=999)
for( b in seq_along(F_nullboot)){
    idx <- sample(n, replace = TRUE)
    g_boot <- dat[idx, 'Region']
    F_b <- F_fun(Y=dat$Murder, g=g_boot)
    F_nullboot[b] <- F_b
}

# Bootstrap p-value
p_boot <- mean(F_nullboot >= F_obs)
p_boot
## [1] 0

hist(F_nullboot, breaks = 40,
    col = grey(0.5, 0.5), border = NA,
    freq = FALSE, main = NA, xlab = 'Null Bootstrap F-statistic',
    xlim = c(0,14))
title('Null Bootstrap F-test', font.main = 1)
abline(v = F_obs, col = rgb(1, 0, 0, .8), lwd = 2)
legend('topright', legend = 'Observed F',
    col = rgb(1, 0, 0, .8), lwd =0, bty='n')

Under additional parametric assumptions (independent observations, equal variances, normal errors), \(\hat{F}\) (under the null) follows an \(F_{G-1,\;n-G}\) distribution. From that theoretical null distribution, we can then compute p-values (which is the default report method in R).

Code 
summary(fit_aov)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Region       3  391.2   130.4   11.14 1.28e-05 ***
## Residuals   46  538.3    11.7                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p_param <- 1 - pf(F_obs, G-1, n-G)
 
 # Compare parametric and bootstrap p-values
cbind(F_obs, p_param, p_boot)
##         F_obs      p_param p_boot
## [1,] 11.14389 0.0001093623      0

The parametric F-test assumes (i) independent observations, (ii) normal errors within each group, and (iii) equal variances across groups (homoscedasticity). If these assumptions are doubtful, then the nonparametric Kruskal-Wallis test (below) is a safer alternative. As a quick diagnostic, visually compare side-by-side boxplots: if the spread differs across groups or the distribution has excessive skew/kurtosis, be cautious with the parametric test.

With \(G-1=2\) and \(n-G=6\) degrees of freedom, we compare to an \(F_{2,6}\) distribution: \(p = 1 - F_{2,6}(7) \approx 0.027\).

Code 
# Computation
summary(aov(Y~g))
##             Df Sum Sq Mean Sq F value Pr(>F)  
## g            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 (from Before)
Fstat <- (BSS/(G-1)) / (WSS/(n-G))
pval  <- 1 - pf(Fstat, G-1, n-G)
cbind(BSS, WSS, Fstat, pval)
##      BSS WSS    Fstat         pval
## [1,]  14   6 54.83333 5.155876e-13

22.2 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.

ImportantKey Definition

The Kruskal-Wallis test 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.

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{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 \(Y_j\) in the other; \(H_0: Prob(Y_i > Y_j)=Prob(Y_j > Y_i)\) versus \(H_A: Prob(Y_i > Y_j) \neq 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\in g}\sum_{j\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 
kw <- kruskal.test(Murder ~ Region, data = dat)
kw
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Murder by Region
## Kruskal-Wallis chi-squared = 19.345, df = 3, p-value = 0.000232

A significant Kruskal-Wallis result tells us that some group differs, not which. To answer which, use a post-hoc pairwise test together with a multiple-comparisons correction (see Multiple Testing in a later chapter).

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{KW}\) with \(n=9\), \(n_g=3\): \[ \hat{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 
Y <- c(1,2,3, 4,5,6, 7,8,9)
g <- factor(rep(c('A','B','C'), each=3))
kruskal.test(Y ~ g)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Y by g
## 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 
A <- c(1,3,5); B <- c(2,4,6)
wilcox.test(A, B)
## 
##  Wilcoxon rank sum exact test
## 
## data:  A and B
## W = 3, p-value = 0.7
## alternative hypothesis: true location shift is not equal to 0

The Kruskal-Wallis test rejected the null on the regional data, so at least one region has a different distribution of murder rates. A post-hoc check pinpoints which pair differs.

Code 
pairwise.wilcox.test(dat$Murder, dat$Region, p.adjust.method='holm')
## 
##  Pairwise comparisons using Wilcoxon rank sum exact test 
## 
## data:  dat$Murder and dat$Region 
## 
##               Northeast South  North Central
## South         0.0011    -      -            
## North Central 0.6012    0.0032 -            
## West          0.2048    0.0082 0.6012       
## 
## P value adjustment method: holm

Looking back at the boxplot, which region stands out?

22.3 Exercises

  1. 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.

  2. 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?

  3. Using the USArrests dataset with state.region as the grouping variable, compute the between-group sum of squares (\(\hat{BSS}\)) and within-group sum of squares (\(\hat{WSS}\)) for Assault by hand (using aggregate()). Verify your \(\hat{F}\) statistic matches the output of summary(aov(Assault ~ state.region, data = ...)).

Further Reading

Recall

This chapter introduced ANOVA: the decomposition of variation between and within groups, the \(F\)-statistic that compares the two, and the nonparametric Kruskal-Wallis analog. The worked three-group example with \(A=\{4,5,6\}, B=\{8,7,9\}, C=\{5,6,7\}\) made the arithmetic concrete: by hand we got \(\hat{BSS}=14\) and \(\hat{WSS}=6\), which aov(Y ~ g) reproduces exactly. The Kruskal-Wallis test replaced outcomes with ranks when the equal-variance normal-error assumptions were doubtful. The next chapter rewrites group comparisons as regressions, compares treatment and effect coding, and builds a permutation test for the factor as a whole.