This chapter compares two groups: how the distribution of one variable changes across the categories of another. We begin with a discrete data and then turn to continuous data where, in both cases, a cardinal outcome \(\hat{y}_{i}\) is split into groups by a factor \(\hat{x}_{i}\). In each case we summarize the groups, compute a statistic for the difference between them, and test that difference using the resampling methods from Hypothesis Testing.
When both variables are factors, a natural way to summarize their association is to look at how the average of one (treated as cardinal) changes across categories of the other. When its categories are ordered, we can also look at how its quantiles change.
Comparing \(\hat{m}_{Y}(x)\) across values of \(x\) is useful as the simplest summary of association between a cardinal outcome and a factor: a single number per subgroup shows whether \(\hat{y}\) tends to be higher in one category. The calculation breaks into two steps: row-by-row, compute the conditional distribution \(\hat{p}(y\mid x)\), then take its weighted average using the values of \(y\).
For example, consider years of schooling (treated as cardinal) across sex in the Wages1 data from the Ecdat package. We first split schooling into one group per sex, and then average within each group.
library('Ecdat')
# Years of schooling within each sex
D_female <- Wages1[, 'sex'] == 'female'
y1_hat <- Wages1[D_female, 'school']
y2_hat <- Wages1[!D_female, 'school']
# Conditional mean of schooling within each sex
m1_hat <- mean(y1_hat)
m2_hat <- mean(y2_hat)
c(m1_hat=m1_hat, m2_hat=m2_hat)
## m1_hat m2_hat
## 11.83748 11.44232Women in this sample average about \(11.8\) years of schooling and men about \(11.4\) years.
Alternatively, aggregate computes the conditional mean for every group at once. To see the two steps in the definition, we first compute the conditional distribution of schooling within each sex, so that each column sums to 1. We then take the weighted average of the schooling values within each column.
# Mean years of school by sex
aggregate(school ~ sex, Wages1, mean)
## sex school
## 1 female 11.83748
## 2 male 11.44232
## Breakdown
# Counts of years of schooling (rows) by sex (columns)
school_sex <- table(Wages1[, 'school'], Wages1[, 'sex'])
# Conditional distribution p(y|x): each column sums to 1
p_school_sex <- prop.table(school_sex, 2)
round(p_school_sex, 3)
##
## female male
## 3 0.000 0.001
## 4 0.000 0.001
## 5 0.001 0.002
## 6 0.002 0.006
## 7 0.005 0.009
## 8 0.014 0.037
## 9 0.038 0.059
## 10 0.106 0.135
## 11 0.196 0.205
## 12 0.379 0.344
## 13 0.120 0.087
## 14 0.089 0.072
## 15 0.043 0.038
## 16 0.007 0.003
# Conditional mean of schooling within each sex
school_values <- as.numeric(rownames(p_school_sex))
m_school_sex <- colSums(school_values * p_school_sex)
m_school_sex
## female male
## 11.83748 11.44232Men are more common at \(10\) years or fewer, while women are more common at \(12\) years or more.
The conditional mean summarizes each subgroup with one number, but we can also ask where the lower, middle, and upper parts of each subgroup lie.
Conditional quantiles are useful when \(\hat{y}\) is ordered but not cardinal: unlike the conditional mean, they only use the ordering of the categories, as with ordered factors. In R, quantile with type=1 inverts the ECDF within each group, so it returns an observed number of years.
# Conditional medians within each sex
med1_hat <- median(y1_hat)
med2_hat <- median(y2_hat)
c(med1_hat=med1_hat, med2_hat=med2_hat)
## med1_hat med2_hat
## 12 12
# Conditional quartiles within each sex
q1_hat <- quantile(y1_hat, probs=c(.25, .5, .75), type=1)
q2_hat <- quantile(y2_hat, probs=c(.25, .5, .75), type=1)
rbind(q1_hat, q2_hat)
## 25% 50% 75%
## q1_hat 11 12 13
## q2_hat 10 12 12The medians are both \(12\) years, because over a third of each group stopped at exactly \(12\) years of schooling. The quartiles still show a difference: the middle half of women have between \(11\) and \(13\) years, while the middle half of men have between \(10\) and \(12\) years. With lumpy discrete data, a single median can hide a difference that several quantiles reveal.
Alternatively, aggregate computes the conditional quantiles for every group at once. To see where these quantiles come from, cumsum down each column of p_school_sex gives \(\hat{F}_{Y}(y \mid x)\) for each sex. Each quartile is the first row where its column reaches \(p\).
# Conditional Quantiles
aggregate(school ~ sex, Wages1, median)
## sex school
## 1 female 12
## 2 male 12
aggregate(school ~ sex, Wages1, quantile, probs=c(.25, .5, .75), type=1)
## sex school.25% school.50% school.75%
## 1 female 11 12 13
## 2 male 10 12 12
## Breakdown
# Conditional ECDF within each sex: cumulative proportions down each column
F_school_sex <- apply(p_school_sex, 2, cumsum)
round(F_school_sex, 3)
##
## female male
## 3 0.000 0.001
## 4 0.000 0.002
## 5 0.001 0.004
## 6 0.003 0.010
## 7 0.008 0.020
## 8 0.022 0.057
## 9 0.060 0.115
## 10 0.166 0.250
## 11 0.361 0.456
## 12 0.740 0.800
## 13 0.860 0.887
## 14 0.950 0.959
## 15 0.993 0.997
## 16 1.000 1.000Consider a mixed dataset where \(\hat{y}_{i}\) is a continuous cardinal variable and \(\hat{x}_{i}\) is a factor variable (typically unordered). For such “grouped data”, we analyze associations via group comparisons. The basic idea is best seen in a comparison of two samples, which corresponds to an \(\hat{x}_{i}\) with two categories. For example, we might compare the heights of men and women in Canada or the homicide rates in two different American states.
We will consider the wages for people with and without completing a degree, again from Wages1. To compare the two groups, we start by making a histogram for each.
# Workers with 12 or 16 years of schooling, and which of them are high school only
xy_hat <- Wages1[Wages1[, 'school'] %in% c(12, 16), c('school', 'wage')]
D_hs <- xy_hat[, 'school'] == 12
# Hourly wages within each level of schooling
y1_hat <- xy_hat[D_hs, 'wage'] # high school only
y2_hat <- xy_hat[!D_hs, 'wage'] # university graduates
# Initial summary figure: two overlaid histograms on a common grid
bks <- seq(0, 34, by=1) # shared bin edges
dlim <- c(0, .2) # shared y limit
cols <- c(rgb(.8, 0, 0, .5), rgb(0, 0, .8, .5)) # semi-transparent fill
hist(y1_hat, breaks=bks, ylim=dlim,
col=cols[1], xlab='Wages',
freq=FALSE, border=NA, main=NA)
rug(y1_hat, col=cols[1])
hist(y2_hat, breaks=bks, ylim=dlim,
col=cols[2],
freq=FALSE, border=NA, add=TRUE)
rug(y2_hat, col=cols[2])
legend('topright',
col=cols, pch=15,
legend=c('12 Years (High School)', '16 Years (University)'),
title='School Completed')
The histogram may show several differences between groups or none at all. In general, we test whether the difference in some statistic \(\hat{g}\) equals \(0\). As in LOTUS and Data Transformations, \(g\) denotes a function, but here \(\hat{g}\) is a function of a group’s whole sample rather than of a single observation. Often, the first statistic we investigate for hypothesis testing is the mean, \(\hat{m}\).
The mean difference is useful as the workhorse single-number summary of “how different are these two groups, on average”, because it is on the same scale as the outcome itself and the natural null hypothesis \(\hat{d}=0\) corresponds to “no difference in population means”. We can test that hypothesis by inverting a confidence interval. The code below computes \(\hat{m}_1\) and \(\hat{m}_2\) from the wage data and reports \(\hat{d}\): the university group has the higher mean wage, so \(\hat{d}<0\) here, and the next step is asking whether the gap is bigger than sampling noise alone would produce. To answer that, we resample each group separately and recompute \(\hat{d}\) each time.
# Difference between means
m1_hat <- mean(y1_hat)
m2_hat <- mean(y2_hat)
d_hat <- m1_hat - m2_hat
round(c(m1_hat=m1_hat, m2_hat=m2_hat, d_hat=d_hat), 2)
## m1_hat m2_hat d_hat
## 5.67 7.54 -1.87
# Bootstrap Distribution
bootstrap_diff <- rep(NA, 9999)
for(b in seq(bootstrap_diff) ){
# keep each worker's group fixed and resample wages within each group
y1_boot <- sample(xy_hat[D_hs, 'wage'], replace=TRUE)
y2_boot <- sample(xy_hat[!D_hs, 'wage'], replace=TRUE)
m1_b_hat <- mean(y1_boot)
m2_b_hat <- mean(y2_boot)
d_b <- m1_b_hat - m2_b_hat
bootstrap_diff[b] <- d_b
}
hist(bootstrap_diff, breaks=50,
border=NA, freq=FALSE,
main=NA, xlab='Difference in Means')
title('Bootstrapped Difference in Means', font.main=1)
# 2-Sided Test via Confidence Interval
boot_ci <- quantile(bootstrap_diff, probs=c(.025, .975))
abline(v=boot_ci, col=rgb(0, 0, 0, .8), lwd=2)
abline(v=0, col=rgb(0, 0, 1, .8), lty=2)
abline(v=d_hat, col=rgb(1, 0, 0, .8), lwd=2)
Such tests make purely statistical statements that only speak to how frequent differences are generated by random chance. They do not say why there are any differences or how large they are. Even on purely statistical grounds, however, we would want to be cautious about making data-driven decisions when
For such reasons, applied statisticians consider many statistics.
The above procedure generalizes from differences in means to differences in medians. To start, we plot the ECDF’s for both groups.
# Quantile Comparison
## Distribution 1
F1_hat <- ecdf(y1_hat)
plot(F1_hat, col=cols[1],
pch=16, xlab='Wages', xlim=c(0, 24),
main=NA, bty='n')
title('Comparing Medians', font.main=1)
## Median 1
med1_hat <- median(y1_hat) # equivalently, quantile(F1_hat, probs=0.5)
segments(med1_hat, 0, med1_hat, 0.5, col=cols[1], lty=2)
abline(h=0.5, lty=2)
## Distribution 2
F2_hat <- ecdf(y2_hat)
plot(F2_hat, add=TRUE, col=cols[2], pch=16)
## Median 2
med2_hat <- median(y2_hat) # equivalently, quantile(F2_hat, probs=0.5)
segments(med2_hat, 0, med2_hat, 0.5, col=cols[2], lty=2)
## Legend
legend('bottomright',
col=cols, pch=15,
legend=c('12 Years (High School)', '16 Years (University)'),
title='School Completed')
# Bootstrap Distribution Function
boot_quant <- function(y_hat, D_group, n_boot=9999, probs=0.5, ...){
bootstrap_diff <- rep(NA, n_boot)
for(b in seq(bootstrap_diff)){
y1_boot <- sample(y_hat[D_group], replace=TRUE)
y2_boot <- sample(y_hat[!D_group], replace=TRUE)
q1_b <- quantile(y1_boot, probs=probs, ...)
q2_b <- quantile(y2_boot, probs=probs, ...)
d_b <- q1_b - q2_b
bootstrap_diff[b] <- d_b
}
return(bootstrap_diff)
}
# 2-Sided Test for Median Differences
d_med <- med1_hat - med2_hat
boot_d <- boot_quant(xy_hat[, 'wage'], D_hs, n_boot=9999, probs=0.5)
hist(boot_d, breaks=50, border=NA, freq=FALSE,
main=NA, xlab='Difference in Medians')
title('Bootstrap Difference in Medians', font.main=1)
abline(v=quantile(boot_d, probs=c(.025, .975)), col=rgb(0, 0, 0, .8), lwd=2)
abline(v=0, col=rgb(0, 0, 1, .8), lty=2)
abline(v=d_med, col=rgb(1, 0, 0, .8), lwd=2)
Changing probs in boot_quant tests other quantiles, such as the quartiles. Note that bootstrap tests can perform poorly with highly unequal variances or skewed data. To see this yourself, make a simulation with skewed data and unequal variances.
In principle, we can also examine whether there are differences in spread (sd or IQR) or shape (skew or kurtosis). We use the same hypothesis testing procedure as above.
We can also examine whether there are any differences between the entire distributions. We typically start by plotting the two ECDF’s, and then calculate a statistic for hypothesis testing. A quantile-quantile plot is another way to compare the two distributions visually.
To go from a picture of two distributions to a single-number summary of how different they are, we measure the gap between their ECDFs.
\(\hat{\text{KS}}\) is useful when we expect distributions to differ by a shift: the gap piles up at one location and the max catches it. The cost is that \(\hat{\text{KS}}\) can miss differences where small positive and small negative gaps cancel along the support. \(\hat{\text{CVM}}\) is useful for accumulated small differences spread across the support, since every pooled sample point contributes to the sum.
# KS and CVM statistics for any two-group split
ks_cvm <- function(y1_hat, y2_hat){
y_values <- sort(c(y1_hat, y2_hat)) # pooled sample points
F1_vals <- ecdf(y1_hat)(y_values)
F2_vals <- ecdf(y2_hat)(y_values)
gap <- abs(F1_vals - F2_vals)
# KS: largest vertical gap (p=1); CVM: summed squared gaps (p=2)
return(c(KS=max(gap), CVM=sum(gap^2)))
}
stat_obs <- ks_cvm(y1_hat, y2_hat)
round(stat_obs, 2)
## KS CVM
## 0.42 76.16
Suppose there are two groups of summer workers, with data on the earnings of each individual. In thousands of Canadian Dollars, specifically, Group \(A\) has earnings \(\{18,20,25,30\}\) and group \(B\) has earnings \(\{19,21,24,35\}\).
We can compute the CVM and KS statistics for further statistical testing. Start by drawing the ECDF’s.
# Earnings of each group, in thousands of Canadian Dollars
yA0_hat <- c(18, 20, 25, 30)
yB0_hat <- c(19, 21, 24, 35)
FA_hat <- ecdf(yA0_hat)
FB_hat <- ecdf(yB0_hat)
# ECDFs as step functions; shape and line type also mark the group
plot(FA_hat, col=cols[1], pch=16, lwd=2,
xlab='Earnings', ylab='ECDF',
main=NA, bty='n')
plot(FB_hat, add=TRUE, col=cols[2], pch=17, lty=2, lwd=2)
# Vertical gaps at the 8 pooled values, one per row of the table below
x_values <- sort(c(yA0_hat, yB0_hat))
segments(x_values, FA_hat(x_values), x_values, FB_hat(x_values),
lwd=3, col=grey(0, .75), lty=3)
cvm0_hat <- sum((FA_hat(x_values) - FB_hat(x_values))^2)
title(paste0('CVM: ', cvm0_hat), adj=1, font.main=1)
legend('topleft', bty='n', inset=c(0, .1),
col=cols, pch=c(16, 17), lty=c(1, 2), lwd=2,
legend=c('Group A', 'Group B'))
Each dotted segment is a gap of \(1/4\), and at the other four pooled values the two ECDFs meet. Then compute \(\hat{F}_A(x) - \hat{F}_B(x)\) and its square at each of the 8 pooled sorted values:
| \(x\) | \(\hat{F}_A(x)\) | \(\hat{F}_B(x)\) | \(\hat{F}_A-\hat{F}_B\) | \((\hat{F}_A-\hat{F}_B)^2\) |
|---|---|---|---|---|
| 18 | 1/4 | 0 | 1/4 | 1/16 |
| 19 | 1/4 | 1/4 | 0 | 0 |
| 20 | 2/4 | 1/4 | 1/4 | 1/16 |
| 21 | 2/4 | 2/4 | 0 | 0 |
| 24 | 2/4 | 3/4 | -1/4 | 1/16 |
| 25 | 3/4 | 3/4 | 0 | 0 |
| 30 | 4/4 | 3/4 | 1/4 | 1/16 |
| 35 | 4/4 | 4/4 | 0 | 0 |
\(\text{CVM} = 1/16 + 1/16 + 1/16 + 1/16 = 1/4\).
The KS statistic is the largest absolute value in the \(\hat{F}_A-\hat{F}_B\) column. \(\text{KS} = \max\{1/4, 0, 1/4, 0, |-1/4|, 0, 1/4, 0\} = 1/4\).
Another useful visualization for two groups is to plot the quantiles against one another: a quantile-quantile plot. I.e., the first data point on the bottom left shows the first quantile for both distributions.
# Percentiles of each group's wages
quants <- seq(0, 1, length.out=101)
Q1 <- quantile(y1_hat, probs=quants)
Q2 <- quantile(y2_hat, probs=quants)
# Compare Quantiles
ry <- range(c(y1_hat, y2_hat))
plot(Q1, Q2, xlim=ry, ylim=ry,
xlab=expression(Q[1]),
ylab=expression(Q[2]),
main=NA,
pch=16, col=grey(0, .25))
title('Quantile-Quantile Plot', font.main=1)
abline(a=0, b=1, lty=2)
Return to the summer workers from the Must Know above, with group \(A\) earning \(\{18,20,25,30\}\) and group \(B\) earning \(\{19,21,24,35\}\). Compare the two distributions by drawing a QQ plot.
We sort both samples and plot Group \(A\) quantiles (x-axis) against Group \(B\) quantiles (y-axis): \((18,19),\ (20,21),\ (25,24),\ (30,35)\).
# Sorting pairs the k-th smallest earnings in group A with the k-th smallest in group B
ry0 <- range(c(yA0_hat, yB0_hat))
plot(sort(yA0_hat), sort(yB0_hat), xlim=ry0, ylim=ry0,
xlab='Group A Earnings', ylab='Group B Earnings',
main=NA, pch=16, col=grey(0, .5))
title('Quantile-Quantile Plot', font.main=1)
abline(a=0, b=1, lty=2) # 45 degree line: identical distributions
The 45° reference line indicates identical distributions: point 1 lies above the line (18,18) and so do points 2 and 4.
The mean and median tests above invert a bootstrap confidence interval. We can also test by “imposing the null”, where it is more common to use \(p\)-values instead of confidence intervals.
The permutation test (also called a randomization test) is useful whenever the null hypothesis is “no association”, because reshuffling labels imposes that null exactly by construction, without needing to assume a distributional form. Every reshuffle uses the same pooled data, so the only variability comes from re-splitting the labels. If using a hard decision rule, it is most common to reject the null at the \(5\%\) level when \(p<0.05\) and fail to reject when \(p>0.05\). Alternatively, we can reject when the observed statistic falls outside the middle \(95\%\) of the permutation distribution, as with the null bootstrap in Hypothesis Testing.
Consider again the two groups of summer workers, with earnings in thousands of Canadian Dollars. Group \(A\) has earnings \(\{18,20,25,30\}\) and group \(B\) has earnings \(\{19,21,24,35\}\). The observed mean difference is \(\hat{d}=23.25-24.75=-1.5\). A permutation keeps the eight pooled values but reassigns which four belong to group \(A\).
# Pooled earnings and group labels (1=A, 2=B)
y0_hat <- c(yA0_hat, yB0_hat) # group A's four earnings, then group B's
g0 <- c(1, 1, 1, 1, 2, 2, 2, 2)
d0_hat <- mean(y0_hat[g0==1]) - mean(y0_hat[g0==2])
d0_hat
## [1] -1.5
# A few reshuffled labels, written out by hand
g0_perm1 <- c(2, 1, 1, 2, 1, 2, 1, 2)
g0_perm2 <- c(1, 2, 2, 1, 1, 2, 1, 2)
g0_perm3 <- c(2, 2, 1, 1, 1, 2, 2, 1)
mean(y0_hat[g0_perm1==1]) - mean(y0_hat[g0_perm1==2])
## [1] -4
mean(y0_hat[g0_perm2==1]) - mean(y0_hat[g0_perm2==2])
## [1] -2.5
mean(y0_hat[g0_perm3==1]) - mean(y0_hat[g0_perm3==2])
## [1] 6.5
# With 8 values there are only choose(8, 4)=70 ways to pick group A
id_A <- combn(8, 4)
d0_perm <- apply(id_A, 2, function(id){
yA0_perm <- y0_hat[id]
yB0_perm <- y0_hat[-id]
mean(yA0_perm) - mean(yB0_perm)
})
hist(d0_perm, breaks=25, border=NA, freq=FALSE,
main=NA, xlab='Difference in Means')
# 2-Sided Test via the middle 95% of the relabelings
perm0_ci <- quantile(d0_perm, probs=c(.025, .975))
perm0_ci
## 2.5% 97.5%
## -7.1375 7.1375
abline(v=perm0_ci, col=rgb(0, 0, 0, .8), lwd=2)
abline(v=d0_hat, col=rgb(1, 0, 0, .8), lwd=2)
# Two-sided p-value: share of relabelings at least as extreme as observed
mean(abs(d0_perm) >= abs(d0_hat))
## [1] 0.7714286About \(77\%\) of the \(70\) relabelings produce a gap at least as large as \(1.5\), so this small difference is easy to explain by chance. Likewise, \(\hat{d}=-1.5\) sits well inside the middle \(95\%\) of relabelings, from about \(-7.1\) to \(7.1\).
For the wage data, there are far too many relabelings to list, so we draw many of them at random.
# Permutation Distribution: re-split the pooled wages of the two groups
permute_diff <- rep(NA, 9999)
for(b in seq(permute_diff) ){
#randomize data
D_hs0 <- sample(D_hs)
y1_perm <- xy_hat[D_hs0, 'wage']
y2_perm <- xy_hat[!D_hs0, 'wage']
# recompute statistic
m1_perm_hat <- mean(y1_perm)
m2_perm_hat <- mean(y2_perm)
d_b <- m1_perm_hat - m2_perm_hat
permute_diff[b] <- d_b
}
hist(permute_diff, breaks=50,
border=NA, freq=FALSE,
main=NA, xlab='Difference in Means')
title('Permuted Difference in Means', font.main=1)
# 2-Sided Test via the middle 95% of the permutation distribution
perm_ci <- quantile(permute_diff, probs=c(.025, .975))
abline(v=perm_ci, col=rgb(0, 0, 0, .8), lwd=2)
abline(v=d_hat, col=rgb(1, 0, 0, .8), lwd=2)
# Two-sided p-value: share of permutations at least as extreme as observed
p <- mean(abs(permute_diff) >= abs(d_hat))
p
## [1] 0.02060206The same reshuffling works for medians: only the statistic that we recompute on each permutation changes. We pool the wages of the two groups, randomly pick which of them belong to the high-school group, and compute the difference in medians.
# Observed difference in medians
d_med <- med1_hat - med2_hat
# Permutation Distribution: same re-split as for means, recomputing medians
permute_med_diff <- rep(NA, 9999)
for(b in seq(permute_med_diff) ){
D_hs0 <- sample(D_hs)
y1_perm <- xy_hat[D_hs0, 'wage']
y2_perm <- xy_hat[!D_hs0, 'wage']
permute_med_diff[b] <- median(y1_perm) - median(y2_perm)
}
hist(permute_med_diff, breaks=50,
border=NA, freq=FALSE,
main=NA, xlab='Difference in Medians')
title('Permuted Difference in Medians', font.main=1)
# 2-Sided Test via the middle 95% of the permutation distribution
perm_med_ci <- quantile(permute_med_diff, probs=c(.025, .975))
abline(v=perm_med_ci, col=rgb(0, 0, 0, .8), lwd=2)
abline(v=d_med, col=rgb(1, 0, 0, .8), lwd=2)
# Two-sided p-value: share of permutations at least as extreme as observed
p_med <- mean(abs(permute_med_diff) >= abs(d_med))
p_med
## [1] 0.00350035The \(p\)-value is the share of permutations with a median gap at least as large as the observed one. A small \(p\)-value means the observed gap in medians is hard to explain by chance.
The same reshuffling tests whether the \(\hat{\text{KS}}\) and \(\hat{\text{CVM}}\) statistics from Distributional Comparisons are bigger than chance alone would produce. To impose the null hypothesis of identical distributions, we permute the group labels: each resample keeps the pooled data but randomly re-splits it into two groups of the original sizes. Each permutation that produces a statistic at least as large as the observed one is a case where random label-shuffling alone created a gap this big, so the \(p\)-value is the fraction of such permutations. A small \(p\)-value means the observed gap is hard to explain by chance. Both statistics are never negative, and only large values count against the null, so we use a one-sided interval: we reject at the \(5\%\) level when the observed statistic exceeds the \(95\)th percentile of its permutation distribution.
# ks_cvm() and the observed statistics stat_obs come from Distributional Comparisons
# Permutation null distribution: reshuffle group labels
perm_stats <- matrix(NA, nrow=999, ncol=2,
dimnames=list(NULL, c('KS', 'CVM')))
for(b in seq(nrow(perm_stats))){
D_hs0 <- sample(D_hs)
y1_perm <- xy_hat[D_hs0, 'wage']
y2_perm <- xy_hat[!D_hs0, 'wage']
perm_stats[b, ] <- ks_cvm(y1_perm, y2_perm)
}
# Permutation p-values
p_ks <- mean(perm_stats[, 'KS'] >= stat_obs['KS'])
p_cvm <- mean(perm_stats[, 'CVM'] >= stat_obs['CVM'])
# One-sided critical values: 95th percentile of each permutation distribution
perm_crit <- apply(perm_stats, 2, quantile, probs=.95)
round(c(KS_stat=stat_obs[['KS']], KS_crit=perm_crit[['KS']], KS_pval=p_ks,
CVM_stat=stat_obs[['CVM']], CVM_crit=perm_crit[['CVM']], CVM_pval=p_cvm), 3)
## KS_stat KS_crit KS_pval CVM_stat CVM_crit CVM_pval
## 0.418 0.330 0.007 76.164 34.893 0.003Comment 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.
A permutation test shuffles the group labels to build a null distribution, while a bootstrap test resamples with replacement from each group separately. Explain why permutation is appropriate for testing “no difference between groups” and when you might prefer the bootstrap instead.
Using the Wages1 data from the Ecdat package, compute the difference in mean wages \(\hat{d} = \hat{m}_1 - \hat{m}_2\) between men and women (the sex variable). Then compute the Kolmogorov-Smirnov statistic \(\hat{\text{KS}}\) between the two wage distributions by hand (i.e., find the maximum absolute difference between the two ECDF’s evaluated at the pooled sorted data).
Using the mtcars dataset, split cars into two groups based on transmission type (am: 0 = automatic, 1 = manual). Write a permutation test with \(B = 9999\) iterations to test whether the difference in mean mpg between automatic and manual cars is statistically different from zero. Report the \(p\)-value.
This chapter compared two groups in three ways: mean differences, quantile differences, and distributional comparisons via \(\hat{\text{KS}}\) and \(\hat{\text{CVM}}\). The worked example throughout used wages from the Wages1 dataset split by schooling (12 vs 16 years), with permutation tests for \(p\)-values; the small worked example of summer earnings \(A=\{18,20,25,30\}\) and \(B=\{19,21,24,35\}\) walked through a hand permutation of the mean difference (\(\hat{d}=-1.5\), \(p\approx 0.77\)) and the cell-by-cell \(\hat{\text{CVM}}=1/4\) calculation.