The previous chapter compared means and distributions across several groups. This chapter rewrites those group comparisons as regressions, first with treatment coding and then with effect coding. It then develops two resampling methods with different purposes: permutation tests impose a null relationship by shuffling a factor variable, while subsampling keeps each observation intact and checks how estimates vary when some observations are omitted. These methods sharpen what a regression can establish, but none turns an observational relationship into a causal effect. The word effect in effect coding names a parameterization, not a causal claim.
We continue with USArrests and the four-level state.region factor from Comparing Multiple Groups.
# Recreate the regional data so this chapter runs independently.
dat <- USArrests
dat$Region <- state.region
table(dat$Region)
##
## Northeast South North Central West
## 9 16 12 13ANOVA 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 in the previous chapter. 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.
# Compare the regional means to treatment-coded coefficients.
M_g <- aggregate(Murder ~ Region, data=dat, FUN=mean)[,'Murder']
M_g
## [1] 4.700000 11.706250 5.700000 7.030769
# The intercept is the reference mean.
M_g - M_g[1]
## [1] 0.000000 7.006250 1.000000 2.330769
# Each remaining coefficient is a difference from the reference.
fit_lm <- lm(Murder ~ Region, data=dat)
coef(fit_lm)
## (Intercept) RegionSouth RegionNorth Central RegionWest
## 4.700000 7.006250 1.000000 2.330769The ANOVA table for a regression with a single factor is therefore the same as the table from aov().
fit_aov <- aov(Murder ~ Region, data=dat)
anova(fit_aov)
## Analysis of Variance Table
##
## Response: Murder
## Df Sum Sq Mean Sq F value Pr(>F)
## Region 3 391.24 130.412 11.144 1.282e-05 ***
## Residuals 46 538.32 11.703
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(fit_lm)
## Analysis of Variance Table
##
## Response: Murder
## Df Sum Sq Mean Sq F value Pr(>F)
## Region 3 391.24 130.412 11.144 1.282e-05 ***
## Residuals 46 538.32 11.703
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Treatment 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. Because the four regions contain different numbers of states, the unweighted mean of their four means is not the observation-weighted mean across all 50 states.
# Display the sum-to-zero contrast matrix used for effect coding.
effect_matrix <- contr.sum(nlevels(dat$Region))
rownames(effect_matrix) <- levels(dat$Region)
effect_matrix
## [,1] [,2] [,3]
## Northeast 1 0 0
## South 0 1 0
## North Central 0 0 1
## West -1 -1 -1
# Apply contrasts to this model only, without changing global options.
fit_effect <- lm(
Murder ~ Region,
data=dat,
contrasts=list(Region='contr.sum')
)
coef(fit_effect)
## (Intercept) Region1 Region2 Region3
## 7.284255 -2.584255 4.421995 -1.584255
# Compare each regional mean with the unweighted mean of group means.
region_means <- aggregate(Murder ~ Region, data=dat, FUN=mean)
effect_mean <- mean(region_means[,'Murder'])
region_effects <- region_means[,'Murder'] - effect_mean
data.frame(
region=region_means[,'Region'],
group_mean=region_means[,'Murder'],
coded_effect=region_effects
)
## region group_mean coded_effect
## 1 Northeast 4.700000 -2.5842548
## 2 South 11.706250 4.4219952
## 3 North Central 5.700000 -1.5842548
## 4 West 7.030769 -0.2534856
# Recover the omitted West effect from the sum-to-zero restriction.
-sum(coef(fit_effect)[-1])
## [1] -0.2534856Changing the contrasts changes the meanings of the coefficients, but not the fitted values, residuals, \(\hat{R}^2\), or global \(F\)-statistic.
treatment_F <- anova(fit_lm)[['F value']][1]
effect_F <- anova(fit_effect)[['F value']][1]
c(
maximum_fitted_difference=max(abs(fitted(fit_lm) - fitted(fit_effect))),
maximum_residual_difference=max(abs(resid(fit_lm) - resid(fit_effect))),
treatment_F=treatment_F,
effect_F=effect_F
)
## maximum_fitted_difference maximum_residual_difference
## 3.552714e-15 2.442491e-15
## treatment_F effect_F
## 1.114389e+01 1.114389e+01We can extend the group-comparison model by adding continuous predictors alongside the factor. For example, we might ask whether regional differences in murder arrest rates persist after holding Assault fixed.
# Add a continuous predictor to the regional model.
reg_combined <- lm(Murder ~ Region + Assault, data=dat)
summary(reg_combined)
##
## Call:
## lm(formula = Murder ~ Region + Assault, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.4434 -1.2104 0.0058 1.6280 6.0123
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.216516 0.936455 0.231 0.818201
## RegionSouth 3.702630 1.021206 3.626 0.000731 ***
## RegionNorth Central 1.224174 0.986690 1.241 0.221151
## RegionWest 0.187047 1.007021 0.186 0.853481
## Assault 0.035396 0.004474 7.912 4.63e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.237 on 45 degrees of freedom
## Multiple R-squared: 0.7578, Adjusted R-squared: 0.7363
## F-statistic: 35.2 on 4 and 45 DF, p-value: 2.511e-13Holding Assault 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.
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. The dummy-variable approach and a dedicated fixed-effects estimator are algebraically equal, giving the same coefficient and residuals.
# Estimate the Assault coefficient with explicit region indicators.
reg_dv <- lm(Murder ~ -1 + Assault + Region, data=dat)
coef(reg_dv)['Assault']
## Assault
## 0.03539592
# Estimate the same coefficient by absorbing the region effects.
reg_fe <- fixest::feols(Murder ~ Assault | Region, data=dat)
coef(reg_fe)
## Assault
## 0.03539592With 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.
# Order the nested models from simplest to most complex.
reg0 <- lm(Murder ~ 1, data=dat)
reg1 <- lm(Murder ~ Region, data=dat)
reg2 <- lm(Murder ~ Region + Assault, data=dat)
reg3 <- lm(Murder ~ Region * Assault, data=dat)
# Test whether Region adds information.
anova(reg0, reg1)
## Analysis of Variance Table
##
## Model 1: Murder ~ 1
## Model 2: Murder ~ Region
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 49 929.55
## 2 46 538.32 3 391.24 11.144 1.282e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Test whether Assault adds information beyond Region.
anova(reg1, reg2)
## Analysis of Variance Table
##
## Model 1: Murder ~ Region
## Model 2: Murder ~ Region + Assault
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 46 538.32
## 2 45 225.12 1 313.19 62.605 4.626e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Test whether the Assault slope differs by Region.
anova(reg2, reg3)
## Analysis of Variance Table
##
## Model 1: Murder ~ Region + Assault
## Model 2: Murder ~ Region * Assault
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 45 225.12
## 2 42 212.05 3 13.072 0.863 0.4678# 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)
p <- 1 - pf(F, df0 - df1, df1)
cbind(F, p)
## F p
## [1,] 11.14389 1.282093e-05The 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. They differ from the label bootstrap in the previous chapter because sampling without replacement keeps all four regional 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.
set.seed(23)
F_stat <- function(Y, group) {
fit <- lm(Y ~ group)
F_value <- anova(fit)[['F value']][1]
return(F_value)
}
F_obs <- F_stat(Y=dat[,'Murder'], group=dat[,'Region'])
B <- 9999
F_perm <- rep(NA, B)
for (b in seq(B)) {
group_b <- sample(dat[,'Region'], replace=FALSE)
F_perm[b] <- F_stat(Y=dat[,'Murder'], group=group_b)
}
# Add one simulated result to avoid reporting a zero Monte Carlo p-value.
p_perm <- (sum(F_perm >= F_obs) + 1)/(B + 1)
p_param <- anova(fit_lm)[['Pr(>F)']][1]
c(F_obs=F_obs, p_permutation=p_perm, p_parametric=p_param)
## F_obs p_permutation p_parametric
## 1.114389e+01 1.000000e-04 1.282093e-05
hist(F_perm, breaks=40, border=NA, freq=FALSE, main=NA,
xlim=c(0, F_obs*1.05),
xlab='F-statistic under permuted region 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 regional \(F\)-statistic lies in the far-right tail of the permutation distribution. The test therefore provides evidence that murder arrest rates and region are associated, under the assumption that the labels are exchangeable under the null.
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 Region 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 a substantive assumption rather than a computational fact. States may be spatially related, and region was not randomly assigned. The small permutation \(p\)-value does not estimate the effect of moving a state to another region.
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}-\bar{\hat{B}}^{\text{sub}}\right)^2, \] and \(\hat{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.
xy <- USArrests[,c('Murder', 'UrbanPop')]
colnames(xy) <- c('y', 'x')
reg <- lm(y ~ x, data=xy)
b_full <- coef(reg)['x']
# Draw random delete-10 samples without replacement.
n <- nrow(xy)
d <- 10
m <- n - d
B_sub <- 999
rs_coefs <- rep(NA, B_sub)
for (b in seq(B_sub)) {
b_id <- sample(n, m, replace=FALSE)
xy_b <- xy[b_id,]
reg_b <- lm(y ~ x, data=xy_b)
rs_coefs[b] <- coef(reg_b)['x']
}
# Correct the spread for overlapping delete-d samples.
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.x subsample_SE lower upper
## 0.02093466 0.04493402 -0.06713441 0.10900373
hist(rs_coefs, breaks=25, border=NA, freq=FALSE, main=NA,
xlim=range(c(rs_coefs, rs_ci)),
xlab='Subsample slope on UrbanPop')
abline(v=b_full, lwd=2, col=rgb(1, 0, 0, .8))
abline(v=rs_ci, lty=2, col=rgb(0, 0, 1, .8))
title(paste0('delete-', d, ' SE = ', round(rs_se, 3)), font.main=1)
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. The word effect does not change the statistical model or supply a causal interpretation.
The permutation test provides evidence against a no-association null under exchangeability. It identifies a surprising arrangement of outcomes and labels, not 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 and Prediction Are Not Causation. The later chapters on Observational Data and Experimental Data develop the designs and assumptions required for causal claims.
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.
Using USArrests and state.region, 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=dat)).
Extend the running example by fitting Murder ~ Region + Assault + UrbanPop and Murder ~ Region * Assault + UrbanPop. Use anova() to test whether the Region-Assault interaction terms are jointly significant after controlling for UrbanPop. Explain why this test does not establish that assault arrests cause murder arrests.
Repeat the subsampling example with \(d=5\), \(d=10\), and \(d=20\). 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 regional example gave an effect-coded intercept of about \(7.28\), the unweighted mean of the four regional means, while its global \(F\)-statistic was identical under both codings. 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. None of these operations made the regional regression causal. The next chapter audits fitted regressions for outliers, leverage, collinearity, and misspecification.