Comparing Multiple Groups compared means and distributions across several groups with ANOVA and Kruskal-Wallis.
This chapter asks three questions.
None of these methods turns an observational relationship into a causal effect. The word effect in effect coding names a parameterization, not a causal claim.
We continue with the group comparison from Comparing Multiple Groups.
## MAIN EXAMPLE
# 3294 US workers, one row per worker.
# wage is the hourly wage in 1980 dollars.
# school and exper are years of schooling and experience.
data(Wages1, package='Ecdat')
xy_hat <- Wages1
# Three-level schooling factor, as in Comparing Multiple Groups
xy_hat$school_g <- cut(xy_hat$school, breaks=c(0, 11, 12, 16), labels=c('<12', '12', '>12'))
table(xy_hat$school_g)
##
## <12 12 >12
## 1353 1188 753
# The groups differ in size, which matters for effect coding below.ANOVA is a special case of linear regression where the only explanatory variables are group indicators.
The model \(\hat{y}_{ig} = b_{0} + \sum_{g=2}^{G} b_{g} \hat{d}_{ig} + \hat{e}_{ig}\) produces the same \(F\)-statistic as the ANOVA model. After fitting, \(\hat{b}_{0}\) estimates the mean of the reference group and each \(\hat{b}_{g}\) estimates the difference between group \(g\) and the reference. This default representation is called treatment coding.
## MAIN EXAMPLE
# Compare the group means to treatment-coded coefficients.
m_g_hat <- aggregate(wage ~ school_g, data=xy_hat, FUN=mean)[,'wage']
m_g_hat
## [1] 4.990937 5.668354 7.275887
# Mean wages rise with schooling, from about 4.99 to 7.28 dollars.
# The intercept is the reference mean.
m_g_hat - m_g_hat[1]
## [1] 0.0000000 0.6774177 2.2849502
# Each remaining coefficient is a difference from the reference.
reg_lm <- lm(wage ~ school_g, data=xy_hat)
coef(reg_lm)
## (Intercept) school_g12 school_g>12
## 4.9909367 0.6774177 2.2849502
# The reference is the first level, fewer than 12 years of schooling.
# The >12 group earns about 2.28 dollars more per hour than the reference.The ANOVA table for a regression with a single factor is therefore the same as the table from aov().
reg_aov <- aov(wage ~ school_g, data=xy_hat)
anova(reg_aov)
## Analysis of Variance Table
##
## Response: wage
## Df Sum Sq Mean Sq F value Pr(>F)
## school_g 2 2541 1270.27 128.02 < 2.2e-16 ***
## Residuals 3291 32654 9.92
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(reg_lm)
## Analysis of Variance Table
##
## Response: wage
## Df Sum Sq Mean Sq F value Pr(>F)
## school_g 2 2541 1270.27 128.02 < 2.2e-16 ***
## Residuals 3291 32654 9.92
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Both tables report the same F-statistic, about 128.Treatment coding makes every coefficient depend on which group R selects as the reference. Effect coding instead compares each displayed group with the mean across all group means.
Effect coding is useful when no group is a natural reference and the deviations from a common center are easier to discuss than several pairwise differences. When the groups contain different numbers of observations, the unweighted mean of the group means generally differs from the observation-weighted mean across all observations.
## MAIN EXAMPLE
# Display the sum-to-zero contrast matrix used for effect coding.
effect_matrix <- contr.sum(nlevels(xy_hat$school_g))
rownames(effect_matrix) <- levels(xy_hat$school_g)
effect_matrix
## [,1] [,2]
## <12 1 0
## 12 0 1
## >12 -1 -1
# Apply contrasts to this model only, without changing global options.
reg_effect <- lm(
wage ~ school_g,
data=xy_hat,
contrasts=list(school_g='contr.sum')
)
coef(reg_effect)
## (Intercept) school_g1 school_g2
## 5.9783926 -0.9874560 -0.3100382
# Compare each group mean with the unweighted mean of group means.
group_means <- aggregate(wage ~ school_g, data=xy_hat, FUN=mean)
effect_mean <- mean(group_means[,'wage'])
group_effects <- group_means[,'wage'] - effect_mean
data.frame(
school_g=group_means[,'school_g'],
group_mean=group_means[,'wage'],
coded_effect=group_effects
)
## school_g group_mean coded_effect
## 1 <12 4.990937 -0.9874560
## 2 12 5.668354 -0.3100382
## 3 >12 7.275887 1.2974942
# The center is about 5.98 dollars.
# Only workers with more than 12 years of schooling are above it.
# Recover the omitted effect, for more than 12 years of schooling.
# It follows from the sum-to-zero restriction.
-sum(coef(reg_effect)[-1])
## [1] 1.297494Changing the contrasts changes the meanings of the coefficients, but not the fitted values, residuals, \(\hat{R}^2\), or global \(F\)-statistic.
treatment_F <- anova(reg_lm)[['F value']][1]
effect_F <- anova(reg_effect)[['F value']][1]
c(
maximum_fitted_difference=max(abs(fitted(reg_lm) - fitted(reg_effect))),
maximum_residual_difference=max(abs(resid(reg_lm) - resid(reg_effect))),
treatment_F=treatment_F,
effect_F=effect_F
)
## maximum_fitted_difference maximum_residual_difference
## 5.151435e-14 5.134781e-14
## treatment_F effect_F
## 1.280237e+02 1.280237e+02
# The fitted values and residuals differ only by rounding error.
# Both codings give the same F-statistic, about 128.We can extend the group-comparison model by adding continuous predictors alongside the factor. For example, we might ask whether group differences persist after holding a continuous predictor fixed.
## MAIN EXAMPLE
# Add a continuous predictor to the schooling model.
reg_combined <- lm(wage ~ school_g + exper, data=xy_hat)
summary(reg_combined)
##
## Call:
## lm(formula = wage ~ school_g + exper, data = xy_hat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.754 -2.027 -0.497 1.441 35.133
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.15373 0.21102 19.684 < 2e-16 ***
## school_g12 0.64024 0.12520 5.114 3.34e-07 ***
## school_g>12 2.32815 0.14317 16.261 < 2e-16 ***
## exper 0.10452 0.02409 4.339 1.48e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.141 on 3290 degrees of freedom
## Multiple R-squared: 0.07746, Adjusted R-squared: 0.07662
## F-statistic: 92.09 on 3 and 3290 DF, p-value: < 2.2e-16
# Holding experience fixed, the gaps from the <12 group are about 0.64 and 2.33.
# Those are close to the raw gaps of about 0.68 and 2.28.
# Each year of experience goes with about 0.10 dollars more per hour.Holding a predictor fixed changes the conditional association being summarized. It does not by itself remove confounding, reverse causation, measurement error, or other reasons that a regression coefficient may not be causal.
The dummy-variable model treats each group as having its own intercept. This framing is useful when group identity captures features we cannot otherwise measure.
One common case is observations of individuals over several time periods. Then you may have two factor variables. An unordered factor indicates who an individual is, and an ordered factor indicates the time period (Mundlak 1978; Wooldridge 2010). 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. Instead, give each city its own dummy variable, or equivalently its own fixed effect.
Fixed effects are useful when a factor predicts the outcome but the coefficients on its individual levels are not the main object of interest. The coefficient on another predictor then describes its partial relationship with the outcome after removing level differences across the factor.
## MAIN EXAMPLE
# Estimate the experience coefficient with explicit schooling indicators.
reg_dv <- lm(wage ~ -1 + exper + school_g, data=xy_hat)
coef(reg_dv)['exper']
## exper
## 0.1045248
# Estimate the same coefficient by absorbing the schooling group effects.
reg_fe <- fixest::feols(wage ~ exper | school_g, data=xy_hat)
coef(reg_fe)
## exper
## 0.1045248
# Both give about 0.105 dollars per hour for each year of experience.With fixed effects, we can also compute averages for each group: \(\hat{m}_{Yg}=\sum_{i}^{n_{g}}\hat{y}_{ig}/n_{g}\), where group \(g\) has \(n_{g}\) observations denoted \(\hat{y}_{ig}\). We can construct a between estimator from \(\hat{m}_{Yg}=b_{0}+\hat{m}_{Xg}b_{1}\). Alternatively, we can subtract each group average to construct a within estimator: \((\hat{y}_{ig}-\hat{m}_{Yg})=(\hat{x}_{ig}-\hat{m}_{Xg})b_{1}\). Absorbing fixed group differences does not remove confounders that vary within a group.
The ANOVA \(F\)-statistic also compares nested regression models. Under additional parametric assumptions, the statistic follows an \(F\) distribution under the null. We can build a sequence of nested models, each adding one layer of complexity, to see which predictors contribute. The first asks whether group means differ, the second whether a continuous predictor adds explanatory power, and the third whether that predictor’s slope varies across groups (Brambor et al. 2006; Balli and Sørensen 2013).
## MAIN EXAMPLE
# Order the nested models from simplest to most complex.
reg0 <- lm(wage ~ 1, data=xy_hat)
reg1 <- lm(wage ~ school_g, data=xy_hat)
reg2 <- lm(wage ~ school_g + exper, data=xy_hat)
reg3 <- lm(wage ~ school_g * exper, data=xy_hat)
# Test whether the schooling groups add information.
anova(reg0, reg1)
## Analysis of Variance Table
##
## Model 1: wage ~ 1
## Model 2: wage ~ school_g
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 3293 35194
## 2 3291 32654 2 2540.5 128.02 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# F is about 128, and the p-value is essentially zero.
# Test whether experience adds information beyond the schooling groups.
anova(reg1, reg2)
## Analysis of Variance Table
##
## Model 1: wage ~ school_g
## Model 2: wage ~ school_g + exper
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 3291 32654
## 2 3290 32468 1 185.77 18.824 1.476e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# F is about 19, and the p-value is about 0.000015.
# Test whether the experience slope differs by schooling group.
anova(reg2, reg3)
## Analysis of Variance Table
##
## Model 1: wage ~ school_g + exper
## Model 2: wage ~ school_g * exper
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 3290 32468
## 2 3288 32331 2 136.84 6.9581 0.000965 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# F is about 7, and the p-value is about 0.001.
# All three tests reject, although RSS falls by only about 8% in total.# Reproduce the first nested-model F-test by hand.
rss0 <- sum(resid(reg0)^2)
rss1 <- sum(resid(reg1)^2)
df0 <- df.residual(reg0)
df1 <- df.residual(reg1)
F <- ((rss0 - rss1)/(df0 - df1))/(rss1/df1)
# Ask for the upper tail directly, since 1 - pf() would round this p to 0.
p <- pf(F, df0 - df1, df1, lower.tail=FALSE)
cbind(F, p)
## F p
## [1,] 128.0237 2.863426e-54
# F matches the first table above, about 128.
# The p-value, about 3e-54, is below the 2.2e-16 limit that table prints.The theoretical \(F\)-test uses assumptions about independent Normal errors with equal variance. A permutation test instead constructs a null distribution by rearranging the observed factor labels.
Permutation tests are useful when we want a reference distribution that preserves the observed outcomes and the exact number of observations in each group. As in Comparing Multiple Groups, sampling without replacement keeps all the group counts fixed.
The global \(F\)-statistic tests the factor as one object. This is preferable to separately permuting its \(G-1\) coded columns, which would destroy the rule that each observation belongs to exactly one group. It is also coding-invariant, so treatment and effect coding produce the same observed and permuted statistics.
## MAIN EXAMPLE
set.seed(23)
F_stat <- function(y_hat, group) {
reg <- lm(y_hat ~ group)
F_value <- anova(reg)[['F value']][1]
return(F_value)
}
F_obs <- F_stat(y_hat=xy_hat[,'wage'], group=xy_hat[,'school_g'])
n_perm <- 9999
F_perm <- rep(NA, n_perm)
for (b in seq(n_perm)) {
group_b <- sample(xy_hat[,'school_g'], replace=FALSE)
F_perm[b] <- F_stat(y_hat=xy_hat[,'wage'], group=group_b)
}
# Permutation p-value
p_perm <- mean(F_perm >= F_obs)
p_param <- anova(reg_lm)[['Pr(>F)']][1]
c(F_obs=F_obs, p_permutation=p_perm, p_parametric=p_param)
## F_obs p_permutation p_parametric
## 1.280237e+02 0.000000e+00 2.863426e-54
# No shuffle reached the observed F of about 128.
# So p_permutation is 0.
# With 9999 shuffles, that only says the p-value is probably below about 3/9999, or 0.0003.
# A formula has no such limit, so p_parametric can report a far smaller value.
hist(F_perm, breaks=40, border=NA, freq=FALSE, main=NA,
xlim=c(0, F_obs*1.05),
xlab='F-statistic under permuted schooling labels')
abline(v=F_obs, col=rgb(1, 0, 0, .8), lwd=2)
title(paste0('permutation p-value = ',
formatC(p_perm, format='f', digits=4)), font.main=1)
# The observed F sits far to the right of every shuffled F.
# Schooling was not randomly assigned, so this shows association only.
# It does not estimate the effect of moving a worker to another group.An observed \(F\)-statistic in the far-right tail of the permutation distribution provides evidence that the outcome and the factor are associated, under the assumption that the labels are exchangeable under the null.
Repeat the permutation test with murder arrest rates across the census regions.
# All 50 US states, with the four census regions as a factor.
xy2_hat <- USArrests
xy2_hat$Region <- state.region
reg2_lm <- lm(Murder ~ Region, data=xy2_hat)
F2_obs <- F_stat(y_hat=xy2_hat[,'Murder'], group=xy2_hat[,'Region'])
F2_perm <- rep(NA, n_perm)
for (b in seq(n_perm)) {
group2_b <- sample(xy2_hat[,'Region'], replace=FALSE)
F2_perm[b] <- F_stat(y_hat=xy2_hat[,'Murder'], group=group2_b)
}
# Permutation p-value
p2_perm <- mean(F2_perm >= F2_obs)
p2_param <- anova(reg2_lm)[['Pr(>F)']][1]
c(F_obs=F2_obs, p_permutation=p2_perm, p_parametric=p2_param)
## F_obs p_permutation p_parametric
## 1.114389e+01 0.000000e+00 1.282093e-05
hist(F2_perm, breaks=40, border=NA, freq=FALSE, main=NA,
xlim=c(0, F2_obs*1.05),
xlab='F-statistic under permuted region labels')
abline(v=F2_obs, col=rgb(1, 0, 0, .8), lwd=2)
title(paste0('permutation p-value = ',
formatC(p2_perm, format='f', digits=4)), font.main=1)
# The observed F of about 11 lies in the far-right tail of the shuffled F's.
# This is evidence that murder arrest rates and region are associated.
# That reading assumes the labels are exchangeable under the null.
# States may be spatially related, and region was not randomly assigned.
# The p-value does not estimate the effect of moving a state to another region.The simple shuffle above is appropriate for the factor-only model because the null removes the only relationship in the model. In a regression with controls, naively shuffling the factor would also break its relationship with those controls. Adjusted permutation tests therefore require a design-specific restricted permutation or a method that permutes residuals under a reduced model.
Exchangeability is also an assumption about the data, which the shuffle itself cannot check. Observations may be related to one another, and the factor may not have been randomly assigned. A small permutation \(p\)-value does not estimate the effect of moving an observation to another group.
Permutation and subsampling both rerun a statistic on altered data, but they answer different questions. Permutation breaks a relationship to impose a null. Subsampling preserves each row and asks how much an estimate changes when some rows are left out.
Random subsampling is useful for checking whether a regression result is stable across overlapping subsets of the observed data. The table compares it with the bootstrap and jackknife.
| Method | Sample Size per Iteration | Number of Iterations | Resample |
|---|---|---|---|
| Bootstrap | \(n\) | \(B\) | With replacement |
| Jackknife | \(n-1\) | \(n\) | Without replacement |
| Random delete-\(d\) subsample | \(m=n-d\) | \(B\) | Without replacement |
The raw spread of estimates across overlapping subsamples is not the standard error of the full-sample estimator. For the random delete-\(d\) jackknife used here, the variance estimate applies the correction (Politis et al. 1999) \[ \hat{v}^{\text{sub}} = \frac{m}{dB} \sum_{b=1}^{B} \left(\hat{b}^{\text{sub}}_{b}-\hat{b}^{\text{sub}}\right)^2, \] where \(\hat{b}^{\text{sub}}\) is the mean of the \(\hat{b}^{\text{sub}}_{b}\), and \(\hat{\text{SE}}^{\text{sub}}=\sqrt{\hat{v}^{\text{sub}}}\). The example uses a Normal approximation around the full-sample coefficient rather than taking raw percentiles of the overlapping subsample estimates.
## MAIN EXAMPLE
reg <- lm(wage ~ school, data=xy_hat)
b_full <- coef(reg)['school']
# Draw random delete-d samples without replacement.
# Each draw deletes about a fifth of the 3294 workers: d=659 and m=2635.
n <- nrow(xy_hat)
d <- round(n/5)
m <- n - d
n_sub <- 999
rs_coefs <- rep(NA, n_sub)
for (b in seq(n_sub)) {
b_id <- sample(n, m, replace=FALSE)
xy_sub <- xy_hat[b_id,]
reg_b <- lm(wage ~ school, data=xy_sub)
rs_coefs[b] <- coef(reg_b)['school']
}
# Correct the spread for overlapping delete-d samples.
# Here m/d is about 4, so the correction about doubles the raw spread.
rs_mean <- mean(rs_coefs)
rs_var <- (m/d)*mean((rs_coefs - rs_mean)^2)
rs_se <- sqrt(rs_var)
rs_ci <- b_full + qnorm(c(.025, .975))*rs_se
c(b_full=b_full, subsample_SE=rs_se, lower=rs_ci[1], upper=rs_ci[2])
## b_full.school subsample_SE lower upper
## 0.55716169 0.03551907 0.48754560 0.62677779
# Each year of schooling goes with about 0.56 dollars more per hour.
# The subsample SE is about 0.036.
# The interval runs from about 0.49 to 0.63, well above zero.
hist(rs_coefs, breaks=25, border=NA, freq=FALSE, main=NA,
xlim=range(c(rs_coefs, rs_ci)),
xlab='Subsample slope on years of schooling')
abline(v=rs_ci, col=rgb(0, 0, 0, .8), lwd=2)
title(paste0('delete-', d, ' SE = ', round(rs_se, 3)), font.main=1)
# Nearly all raw slopes sit inside the interval, since the subsamples overlap.Repeat the delete-\(d\) standard error for the slope of murder arrest rates on urban population across the states.
# All 50 US states, with percent urban population in the column UrbanPop.
xy2_hat <- USArrests
reg2_full <- lm(Murder ~ UrbanPop, data=xy2_hat)
b2_full <- coef(reg2_full)['UrbanPop']
# Draw random delete-10 samples without replacement, so m=40 states remain.
n2 <- nrow(xy2_hat)
d2 <- 10
m2 <- n2 - d2
n2_sub <- 999
rs2_coefs <- rep(NA, n2_sub)
for (b in seq(n2_sub)) {
b_id <- sample(n2, m2, replace=FALSE)
xy2_sub <- xy2_hat[b_id,]
reg2_b <- lm(Murder ~ UrbanPop, data=xy2_sub)
rs2_coefs[b] <- coef(reg2_b)['UrbanPop']
}
# Correct the spread for overlapping delete-d samples.
rs2_mean <- mean(rs2_coefs)
rs2_var <- (m2/d2)*mean((rs2_coefs - rs2_mean)^2)
rs2_se <- sqrt(rs2_var)
rs2_ci <- b2_full + qnorm(c(.025, .975))*rs2_se
c(b_full=b2_full, subsample_SE=rs2_se, lower=rs2_ci[1], upper=rs2_ci[2])
## b_full.UrbanPop subsample_SE lower upper
## 0.02093466 0.04503122 -0.06732491 0.10919422
# The slope is about 0.02, and the subsample SE is about 0.045.
# The interval runs from about -0.07 to 0.11, so it includes zero.
hist(rs2_coefs, breaks=25, border=NA, freq=FALSE, main=NA,
xlim=range(c(rs2_coefs, rs2_ci)),
xlab='Subsample slope on UrbanPop')
abline(v=rs2_ci, col=rgb(0, 0, 0, .8), lwd=2)
title(paste0('delete-', d2, ' SE = ', round(rs2_se, 3)), font.main=1)
Subsampling also extends to strongly dependent data, including spatial data (Politis and Romano 1994; Sherman 1996; Politis et al. 1999; Choi and Chue 2007; Adamson 2025). The related m-out-of-n bootstrap draws resamples of size \(m<n\) with replacement (Fukuchi 1994; Bickel et al. 1997; Bickel and Sakov 2008; Adamson 2020).
Can we rescale the subsample estimates themselves and then take their percentiles, as in Confidence Intervals? Wu (1986, 1270–71) develops a jackknife percentile method that does this with delete-\(d\) samples. For the sample mean, let \(\hat{m}\) be the full-sample mean and \(\hat{m}^{\text{sub}}_{b}\) the mean of subsample \(b\), which retains \(m=n-d\) observations. The rescaled value is \[ \hat{m}^{\text{scaled}}_{b} = \hat{m} + \sqrt{\frac{m}{d}}\left(\hat{m}^{\text{sub}}_{b} - \hat{m}\right). \] We can take the \(2.5\)th and \(97.5\)th percentiles of these values to construct an approximate \(95\%\) interval for the population mean. For independent observations from the same nonnormal distribution with finite, positive variance, both the number retained and the number deleted must grow with \(n\) for the rescaled mean histogram to approach a Normal distribution. Keeping \(d=1\), or any other fixed number, does not generally satisfy this condition.
For regression coefficients, Wu’s method also uses weights based on the subset design matrices and a scale factor that accounts for the number of coefficients. The mean example above is the simpler special case, so its formula alone does not reproduce Wu’s regression percentile method.
Rescaling the ordinary leave-one-out estimates can match the variance without matching the shape of the sampling distribution. Efron and Tibshirani (1993, 143–44) illustrate an inflated jackknife histogram by multiplying deviations from the average jackknife estimate by \(\sqrt{n-1}\). Their example compares its spread with a bootstrap histogram; it does not establish that its percentiles give valid confidence limits. For the sample mean, their transformation is \[ \hat{m}^{\text{scaled}(-i)} = \hat{m}^{\text{jack}} + \sqrt{n-1}\left(\hat{m}^{(-i)} - \hat{m}^{\text{jack}}\right), \] where \(\hat{m}^{\text{jack}}\) is the average of the leave-one-out means. The variance of these rescaled values, computed with divisor \(n\), equals the usual jackknife variance estimate. However, since \(\hat{m}^{\text{jack}}=\hat{m}\), we can rewrite the transformation as \[\begin{aligned} \hat{m}^{(-i)} &= \frac{n\hat{m} - \hat{x}_{i}}{n-1} \\ \hat{m}^{\text{scaled}(-i)} &= \hat{m} + \frac{\hat{m} - \hat{x}_{i}}{\sqrt{n-1}}. \end{aligned}\] These rescaled values retain the reflected shape of the individual observations, while the sampling distribution of the sample mean approaches a Normal distribution under the assumptions above. Matching the variance therefore does not generally match the percentiles.
Medians and other sample quantiles have an additional difficulty: the ordinary delete-one jackknife can give inconsistent variance estimates (Martin 1990). For smooth statistics, the ordinary jackknife remains useful for estimating standard errors to use in a Normal approximation. Percentile intervals instead require a suitable bootstrap or a justified delete-\(d\) procedure; simply stretching the leave-one-out values is not enough.
A permutation iteration keeps all \(n\) observations but changes which outcome is paired with each factor label. A subsampling iteration keeps the original pairings but uses only \(m<n\) observations. The permutation distribution describes a test statistic under a null relationship. The subsampling distribution describes sensitivity to omitting observations and requires a scaling correction before it estimates full-sample uncertainty.
Treatment coding and effect coding are two coordinate systems for the same fitted group means.
A small permutation \(p\)-value provides evidence against a no-association null under exchangeability. It shows that the observed pairing of outcomes and labels is unusual. It does not show what would happen under an intervention.
The subsampling result shows how a coefficient responds to deleting observations. A stable coefficient can still be confounded, and an unstable coefficient is not necessarily biased.
These limits extend the distinction introduced in Association is not Causation.
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.
Suppose three groups have means \(4\), \(7\), and \(5\). Compute the treatment-coded coefficients using the first group as the reference. Then compute the effect-coded intercept and all three group effects, checking that the effects sum to zero.
Store USArrests with a Region column from state.region as xy2_hat, then construct a permutation test of whether mean Assault arrest rates differ across regions. Use \(B=9999\) permutations, report the global \(F\)-statistic and permutation \(p\)-value, and compare the result with anova(lm(Assault ~ Region, data=xy2_hat)).
Extend the running example by fitting wage ~ school_g + exper + sex and wage ~ school_g * exper + sex. Use anova() to test whether the school_g-exper interaction terms are jointly significant after controlling for sex. Explain why this test does not establish that experience causes higher wages.
Repeat the subsampling example with \(d\) equal to about \(10\%\), \(20\%\), and \(40\%\) of the observations. Apply the delete-\(d\) correction in each case and compare the standard errors. Explain why increasing \(B\) makes the simulation more stable without increasing the information in the original sample.
This chapter rewrote ANOVA as a regression with a factor variable and showed that treatment and effect coding change coefficient meanings without changing the fitted model. The schooling example gave an effect-coded intercept of about \(5.98\), the unweighted mean of the three group means, while its global \(F\)-statistic was identical under both codings. With a fixed effect for each schooling group, a year of experience went with about \(0.10\) dollars more per hour, and the schooling gaps barely moved. Nested-model \(F\)-tests asked whether each added layer adds information, and all three rejected the simpler model. The permutation test shuffled the factor labels to test a no-association null, whereas random delete-\(d\) subsampling preserved each row and measured sensitivity to omitted observations.