Comparing Two Groups split a cardinal outcome by a two-level factor and tested whether the two groups differed. This chapter extends that idea to \(G\) groups: how to describe several groups at once, and then two ways to test whether they differ, one built on means and one built on ranks.
So far the variables were all cardinal, where \(4-3=2-1\); many useful variables are instead labels of group membership that do not do arithmetic this way.
Factor variables are useful for splitting a dataset into groups that can be compared: regions, industries, treatment arms. They come in two flavors:
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, and an ordered factor that indicates the time period. 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 applies to ordinal data; PopulationLeipzig=2, PopulationLausanne=3, PopulationLosAngeles=1.
The built-in state.region vector labels the US Census region of each state, in the same order as the rows of USArrests, so it can be paired directly with the murder rates.
# Outcome: murder arrests per 100K, one value per state
y_hat <- USArrests[, 'Murder']
# Group: the four-level region factor, in the same state order as USArrests.
# Its levels are the values the categorical variable can take.
region <- state.region
levels(region)
## [1] "Northeast" "South" "North Central" "West"
table(region)
## region
## Northeast South North Central West
## 9 16 12 13A factor with \(G\) levels splits a cardinal variable into \(G\) conditional distributions, one per group. The group summaries from Bivariate Distributions apply directly, and with more than two or three groups, side-by-side boxplots are usually the clearest figure, ordered here by median murder rate from lowest to highest.
# Order regions by their median murder rate, lowest to highest
region_order <- names(sort(tapply(y_hat, region, median)))
region <- factor(region, levels=region_order)
# Visual summary: one boxplot per region, box width scaled to group size
boxplot(y_hat ~ region,
col=hcl.colors(4, alpha=.45), varwidth=TRUE,
ylab='Murder arrests (per 100,000)', xlab='Region', main=NA)
title('Murder rates by region', font.main=1)
A boxplot shows five numbers per group; varwidth=TRUE also scales each box’s width to \(\sqrt{n_g}\), so a region with more states draws a wider box. 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. The same low-to-high order carries over, so South is last in the legend, and each line’s thickness scales with its group size the same way the boxes do.
# One ECDF per region, on the same axes (region is already ordered by median)
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(y_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(y_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)
The South’s curve sits to the right of the others over its whole range, so its higher center in the boxplot is not driven by a few states.
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.
# Group size, mean, and median of the murder rate by region
summary_fun <- function(y_hat) {
out <- c(n=length(y_hat), mean=mean(y_hat), median=median(y_hat))
return(out)
}
aggregate(y_hat ~ region, FUN=summary_fun)
## region y_hat.n y_hat.mean y_hat.median
## 1 Northeast 9.000000 4.700000 3.400000
## 2 North Central 12.000000 5.700000 5.150000
## 3 West 13.000000 7.030769 6.800000
## 4 South 16.000000 11.706250 12.850000The South has the highest center and the Northeast the lowest. Whether such differences are larger than chance alone would produce is the question the rest of this chapter takes up.
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?
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.
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}_{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{BSS}\) (“Sum Sq” for region), \(\hat{WSS}\) (“Sum Sq” for Residuals).
# Fit ANOVA
fit_aov <- aov(y_hat ~ region)
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 ' ' 1We 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 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.
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=region)
F_obs
## [1] 11.14389To compute a \(p\)-value, we can use permutation or bootstrap methods to build a null distribution, just as in previous chapters.
# Pairs bootstrap: resample group labels to break the Y--group pairing
n <- length(y_hat)
F_nullboot <- vector(length=999)
for( b in seq_along(F_nullboot)){
idx <- sample(n, replace = TRUE)
g_boot <- region[idx]
F_b <- F_fun(y_hat=y_hat, 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 = range(c(0, F_nullboot, F_obs)))
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 = 2, 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).
# Degrees of freedom for the regional data
n <- length(y_hat)
G <- nlevels(region)
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 1.282093e-05 0ANOVA tests whether group means differ; sometimes the entire distributions differ in ways the mean would miss.
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 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~\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:
kw <- kruskal.test(y_hat ~ region)
kw
##
## Kruskal-Wallis rank sum test
##
## data: y_hat by region
## Kruskal-Wallis chi-squared = 19.345, df = 3, p-value = 0.000232A 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. Here, pairwise.wilcox.test() runs a Mann-Whitney test on each pair of regions, and the Holm correction adjusts the \(p\)-values for the number of pairs tested.
pairwise.wilcox.test(y_hat, region, p.adjust.method='holm')
##
## Pairwise comparisons using Wilcoxon rank sum exact test
##
## data: y_hat and region
##
## Northeast North Central West
## North Central 0.6012 - -
## West 0.2048 0.6012 -
## South 0.0011 0.0032 0.0082
##
## P value adjustment method: holmLooking back at the boxplot in Describing Groups, which region stands out?
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{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 = ...)).
This chapter first described several groups with boxplots, ECDFs, and a table of group means, then tested whether the description reflected more than chance. ANOVA’s decomposition of variation between and within groups produced an \(F\)-statistic; the worked three-group example with \(A=\{4,5,6\}, B=\{8,7,9\}, C=\{5,6,7\}\) made the arithmetic concrete, giving \(\hat{BSS}=14\) and \(\hat{WSS}=6\) by hand, which aov(y0_hat ~ g0) reproduces exactly. The Kruskal-Wallis test replaced outcomes with ranks when the equal-variance normal-error assumptions were doubtful.