A school district conducts a completely randomized experiment with 200 high school students in the same grade. Exactly 100 receive a new mathematics teaching method and 100 receive the standard method. We observe each student’s treatment indicator \(W_i\), pre-treatment achievement level \(X_i\) (Basic, Proficient, or Advanced), and post-test score \(Y_i^{\mathrm{obs}}\).

Our target is the average treatment effect among these 200 students: \[\tau=\frac{1}{n}\sum_{i=1}^n\{Y_i(1)-Y_i(0)\}.\]

We compare three estimates: the difference in means, regression with achievement indicators, and regression with treatment-by-achievement interactions. As in Example 1, distinguish the potential outcomes, the assignment, and the observed data. Here we use large-sample inference for an average effect, rather than an exact randomization test of a sharp null.

1. Generate and inspect a simulated experiment

Step 1: Set up R

We use ggplot2 for plots and sandwich for heteroskedasticity-robust standard errors. If needed, run install.packages(c("ggplot2", "sandwich")) once in the R console.

library(ggplot2)
library(sandwich)
set.seed(123)

Step 2: Generate a fixed table of potential outcomes

Because this is a simulation, we can generate both potential outcomes and know the true average effect. In a real experiment, we would observe only one potential outcome per student.

Students with higher prior achievement have higher baseline scores. The mean treatment effects in the data-generating model are 5, 15, and 3 points for Basic, Proficient, and Advanced students, respectively. Individual effects also vary within each level.

n <- 200
n1 <- 100
Achievement_levels <- c("Basic", "Proficient", "Advanced")
Achievement <- factor(
  sample(Achievement_levels, n, replace = TRUE, prob = c(0.3, 0.5, 0.2)),
  levels = Achievement_levels
)

epsilon <- rnorm(n, mean = 0, sd = 5)
PostTest_0 <- 60 + 10 * (Achievement == "Proficient") +
  20 * (Achievement == "Advanced") + epsilon
Treatment_effect <- 5 + 10 * (Achievement == "Proficient") -
  2 * (Achievement == "Advanced") + rnorm(n, mean = 0, sd = 5)
PostTest_1 <- PostTest_0 + Treatment_effect

true_ATE <- mean(Treatment_effect)
true_ATE
## [1] 9.95407

The realized true_ATE need not equal the population-model average \(0.3(5)+0.5(15)+0.2(3)=9.6\): this particular set of students has random achievement counts and individual effects. Once generated, hold their potential outcomes fixed when reasoning about repeated treatment assignments.

Step 3: Randomize treatment and reveal observed outcomes

Before running the code: Under complete randomization, can the treated group contain 110 students? How would rbinom(n, 1, 0.5) change the assignment design?

sample() shuffles exactly 100 ones and 100 zeros. Independent Bernoulli draws would instead allow the number treated to vary.

Treatment <- sample(c(rep(1, n1), rep(0, n - n1)))
PostTest <- ifelse(Treatment == 1, PostTest_1, PostTest_0)
data <- data.frame(Y = PostTest, W = Treatment, X = Achievement)

head(data)
table(data$W)
## 
##   0   1 
## 100 100
table(Achievement = data$X, Treatment = data$W)
##             Treatment
## Achievement   0  1
##   Basic      29 31
##   Proficient 49 54
##   Advanced   22 15

Complete randomization balances the overall group sizes here; it does not force achievement levels to be balanced between arms. The analysis below uses only the observed columns Y, W, and X. The known true_ATE is a simulation benchmark.

Step 4: Visualize the observed outcomes

Horizontal jitter separates overlapping points without changing the scores. Color identifies a covariate measured before treatment.

ggplot(data, aes(x = W, y = Y, color = X)) +
  geom_point(position = position_jitter(width = 0.1, height = 0, seed = 123),
             alpha = 0.7, size = 2) +
  scale_x_continuous(breaks = c(0, 1), labels = c("Control", "Treatment")) +
  labs(title = "Post-test scores by treatment and prior achievement",
       x = "Assigned teaching method", y = "Post-test score",
       color = "Prior achievement") +
  theme_minimal()

Check your understanding: Does this plot reveal any student’s individual treatment effect? Why might comparing students with similar prior achievement improve precision?

2. Regression without covariate adjustment

Step 1: Connect the regression coefficient to the difference in means

Fit \(Y_i^{\mathrm{obs}}=\alpha+\beta W_i+e_i\). With a binary treatment and an intercept, the fitted intercept is the control mean and the coefficient of W is exactly the treated-minus-control difference in means.

D_obs <- mean(data$Y[data$W == 1]) - mean(data$Y[data$W == 0])
lm.result1 <- lm(Y ~ W, data = data)
c(difference_in_means = D_obs, regression_coefficient = unname(coef(lm.result1)["W"]))
##    difference_in_means regression_coefficient 
##               9.054639               9.054639

Under complete randomization, this estimator is unbiased for the average effect among the study students. We do not need all individual effects to be equal.

Step 2: Compare conventional and robust standard errors

A standard error describes uncertainty in an estimate across repeated experiments; it is not the standard deviation of students’ scores. The conventional lm() standard error uses a common residual variance. Treatment-effect heterogeneity can produce different variances in the two arms, so we also calculate an HC2 heteroskedasticity-robust standard error.

est <- unname(coef(lm.result1)["W"])
se_conventional <- summary(lm.result1)$coefficients["W", "Std. Error"]
se_HC2 <- sqrt(vcovHC(lm.result1, type = "HC2")["W", "W"])
result1 <- c(estimate = est, SE_conventional = se_conventional, SE_HC2 = se_HC2)
result1
##        estimate SE_conventional          SE_HC2 
##        9.054639        1.488320        1.488320

Here the two standard errors coincide because the treatment arms have equal sizes in this unadjusted regression. That equality does not establish equal outcome variances, and need not hold with unequal arm sizes or covariate adjustment.

For this unadjusted regression, HC2 equals the familiar Neyman standard error: \[\widehat{\mathrm{SE}}_{\mathrm{HC2}} =\sqrt{\frac{s_1^2}{n_1}+\frac{s_0^2}{n_0}},\] where \(s_w^2\) is the sample variance of observed outcomes in arm \(w\).

se_Neyman <- sqrt(var(data$Y[data$W == 1]) / sum(data$W == 1) +
                   var(data$Y[data$W == 0]) / sum(data$W == 0))
c(HC2 = se_HC2, Neyman = se_Neyman)
##     HC2  Neyman 
## 1.48832 1.48832
stopifnot(isTRUE(all.equal(unname(se_HC2), se_Neyman)))

For fixed potential outcomes, the exact randomization variance of the difference in means is \(S_1^2/n_1+S_0^2/n_0-S_\tau^2/n\), where \(S_\tau^2\) is the finite-population variance of individual effects. We cannot estimate that last term from the observed data alone. Omitting it gives a variance estimator that is conservative in expectation; an estimated standard error need not exceed the true one in every realized experiment.

Check your understanding: Does using a robust standard error change the estimated treatment effect? Must a robust standard error be larger than a conventional one?

3. Regression with additive covariate adjustment

Step 1: Understand how R represents achievement

Because X is a factor with Basic as its reference level, R uses two indicators: one for Proficient and one for Advanced. The additive regression is \[Y_i^{\mathrm{obs}}=\alpha+\beta W_i+ \gamma_1\mathbf{1}(X_i=\text{Proficient})+ \gamma_2\mathbf{1}(X_i=\text{Advanced})+e_i.\]

head(model.matrix(~ W + X, data = data))
##   (Intercept) W XProficient XAdvanced
## 1           1 1           1         0
## 2           1 1           0         0
## 3           1 1           1         0
## 4           1 0           0         1
## 5           1 1           0         1
## 6           1 1           1         0
lm.result2 <- lm(Y ~ W + X, data = data)
summary(lm.result2)
## 
## Call:
## lm(formula = Y ~ W + X, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.378  -4.579   0.079   3.702  23.234 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   57.041      1.010  56.470   <2e-16 ***
## W              9.613      0.966   9.952   <2e-16 ***
## XProficient   16.440      1.105  14.880   <2e-16 ***
## XAdvanced     19.724      1.426  13.832   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.803 on 196 degrees of freedom
## Multiple R-squared:  0.6515, Adjusted R-squared:  0.6462 
## F-statistic: 122.2 on 3 and 196 DF,  p-value: < 2.2e-16

This working regression fits a common treatment-control gap across achievement levels. It does not establish that the actual effects are equal. Under complete randomization, the treatment coefficient still consistently estimates the overall average effect under standard large-sample conditions, even when effects differ by achievement. It need not be exactly unbiased in a finite sample.

Step 2: Extract the estimate and standard errors

The computation is the same as in Section 2. Select the coefficient by its name, "W", rather than its position in the coefficient table.

est <- unname(coef(lm.result2)["W"])
se_conventional <- summary(lm.result2)$coefficients["W", "Std. Error"]
se_HC2 <- sqrt(vcovHC(lm.result2, type = "HC2")["W", "W"])
result2 <- c(estimate = est, SE_conventional = se_conventional, SE_HC2 = se_HC2)
result2
##        estimate SE_conventional          SE_HC2 
##       9.6133536       0.9659801       0.9578819

A prognostic pre-treatment covariate can reduce uncertainty by explaining outcome variation. A smaller estimated standard error in one simulated dataset is not a guarantee of improved precision in every experiment, nor proof that the point estimate is closer to the truth.

Before moving on: Our simulation gives Proficient students a larger mean effect. Which term is missing from the additive regression if we want to fit different gaps by achievement?

4. Regression with treatment-by-achievement interactions

Step 1: See what an uncentered interaction coefficient means

W * X expands to W + X + W:X. With Basic as the reference level, the coefficient of W estimates the gap for Basic students, not the overall average effect.

lm.uncentered <- lm(Y ~ W * X, data = data)
coef(lm.uncentered)
##   (Intercept)             W   XProficient     XAdvanced W:XProficient 
##     59.769393      4.333236     10.618134     20.290344     11.180770 
##   W:XAdvanced 
##     -2.845305

Step 2: Center the achievement indicators using all students

To make the coefficient of W estimate the overall average effect, subtract each indicator’s full-sample mean. Do not center separately within treatment arms, or use the arbitrary numerical codes of the factor.

data$X1 <- as.numeric(data$X == "Proficient")
data$X2 <- as.numeric(data$X == "Advanced")
data$X1 <- data$X1 - mean(data$X1)
data$X2 <- data$X2 - mean(data$X2)
colMeans(data[c("X1", "X2")])  # Zero, up to numerical rounding.
##            X1            X2 
## -1.998401e-17  2.220446e-18

Fit \[Y_i^{\mathrm{obs}}=\alpha+\beta W_i+\gamma_1X_{1i}+\gamma_2X_{2i} +\delta_1W_iX_{1i}+\delta_2W_iX_{2i}+e_i.\] The fitted treatment effect for student \(i\) is \(\widehat\beta+\widehat\delta_1X_{1i}+\widehat\delta_2X_{2i}\). Because both centered covariates average to zero, the average fitted effect is exactly \(\widehat\beta\).

Centering changes the interpretation of the coefficients, not the fitted outcomes. It makes W the estimated average fitted effect; it does not make it equal to the unknown true average effect in every sample.

lm.result3 <- lm(Y ~ W * (X1 + X2), data = data)
summary(lm.result3)
## 
## Call:
## lm(formula = Y ~ W * (X1 + X2), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.8257  -3.4943  -0.3868   3.3854  20.4270 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  68.9914     0.6101 113.088  < 2e-16 ***
## W             9.5650     0.8633  11.080  < 2e-16 ***
## X1           10.6181     1.4242   7.455 2.90e-12 ***
## X2           20.2903     1.7187  11.806  < 2e-16 ***
## W:X1         11.1808     1.9760   5.658 5.43e-08 ***
## W:X2         -2.8453     2.5709  -1.107     0.27    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.079 on 194 degrees of freedom
## Multiple R-squared:  0.7246, Adjusted R-squared:  0.7175 
## F-statistic: 102.1 on 5 and 194 DF,  p-value: < 2.2e-16
est <- unname(coef(lm.result3)["W"])
se_conventional <- summary(lm.result3)$coefficients["W", "Std. Error"]
se_HC2 <- sqrt(vcovHC(lm.result3, type = "HC2")["W", "W"])
result3 <- c(estimate = est, SE_conventional = se_conventional, SE_HC2 = se_HC2)
result3
##        estimate SE_conventional          SE_HC2 
##       9.5649511       0.8632937       0.8461876

Step 3: Verify the coefficient using achievement-specific comparisons

With only this categorical covariate and full interactions, the fitted gap for each achievement level is its observed difference in means. The overall estimate weights these gaps by the achievement proportions among all 200 students.

Before running the code: Should the three achievement-specific gaps receive equal weights? What would equal weighting target instead?

cell_counts <- table(data$X, data$W)
stopifnot(all(cell_counts > 0))  # Each level needs both treated and control students.
cell_means <- with(data, tapply(Y, list(X, W), mean))
subgroup_effects <- cell_means[, "1"] - cell_means[, "0"]
weights <- as.numeric(prop.table(table(data$X)))

print(data.frame(
  achievement = levels(data$X), proportion = weights,
  control_mean = cell_means[, "0"], treated_mean = cell_means[, "1"],
  estimated_effect = subgroup_effects
))
##            achievement proportion control_mean treated_mean estimated_effect
## Basic            Basic      0.300     59.76939     64.10263         4.333236
## Proficient  Proficient      0.515     70.38753     85.90153        15.514006
## Advanced      Advanced      0.185     80.05974     81.54767         1.487931
weighted_effect <- sum(weights * subgroup_effects)
c(weighted_subgroup_effect = weighted_effect,
  centered_W_coefficient = unname(coef(lm.result3)["W"]))
## weighted_subgroup_effect   centered_W_coefficient 
##                 9.564951                 9.564951
stopifnot(isTRUE(all.equal(weighted_effect, unname(coef(lm.result3)["W"]))))

Under standard large-sample conditions for complete randomization, this fully interacted adjustment has asymptotic precision at least as good as the unadjusted difference in means. That is not a guarantee for every small sample. Sparse achievement-by-treatment cells can make the fitted gaps unstable; an empty cell prevents estimating its gap this way. HC2 supports large-sample, generally conservative inference for the study average effect under these conditions; it does not make the regression model exact.

Check your understanding: Why do we center before interpreting the coefficient of W as an overall effect? Would centering change the predictions from the interaction model?

5. Compare estimates, confidence intervals, and p-values

Step 1: Use HC2 standard errors for approximate inference

For each method, form an approximate 95% confidence interval \[\widehat\tau\ \pm\ z_{0.975}\widehat{\mathrm{SE}}_{\mathrm{HC2}}.\] For the two-sided null \(H_0:\tau=0\), use \(z=\widehat\tau/\widehat{\mathrm{SE}}_{\mathrm{HC2}}\) and the standard normal approximation. This tests an average-effect null, which permits positive and negative individual effects that average to zero. It is different from Example 1’s sharp null that every individual effect is zero.

result <- data.frame(
  method = c("No adjustment", "Additive adjustment", "Centered interactions"),
  rbind(result1, result2, result3), row.names = NULL
)
alpha <- 0.05
critical_value <- qnorm(1 - alpha / 2)
result$CI_lower <- result$estimate - critical_value * result$SE_HC2
result$CI_upper <- result$estimate + critical_value * result$SE_HC2
result$z_statistic <- result$estimate / result$SE_HC2
result$p_value <- 2 * pnorm(-abs(result$z_statistic))
result

The conventional summary(lm(...)) output above uses conventional standard errors and t-based inference. The comparison table here consistently uses HC2 standard errors and a normal approximation. These intervals are not finite-sample exact randomization intervals.

Step 2: Plot the intervals against the simulation truth

result$method <- factor(result$method,
  levels = c("Centered interactions", "Additive adjustment", "No adjustment"))

ggplot(result, aes(y = method, x = estimate)) +
  geom_segment(aes(x = CI_lower, xend = CI_upper, yend = method),
               linewidth = 0.8) +
  geom_point(size = 3) +
  geom_vline(xintercept = true_ATE, color = "red", linetype = "dashed",
             linewidth = 0.8) +
  labs(title = "Average-effect estimates and approximate 95% intervals",
       subtitle = "HC2 standard errors; dashed line is the true study ATE",
       x = "Effect on mathematics score (points)", y = NULL) +
  theme_minimal()

Compare the estimates and interval widths. All three methods target the same study average effect, but they use prior achievement differently. We know the dashed line only because we generated both potential outcomes. In a real experiment, it would be unavailable.

In this run, the true study ATE is 9.95 points. The HC2 standard errors are 1.49, 0.96, 0.85 for no adjustment, additive adjustment, and centered interactions, respectively. Both adjusted methods give narrower intervals here.

A 95% confidence procedure describes coverage across repeated experiments; a single interval either contains the fixed true effect or does not. One simulation cannot establish a method’s coverage or relative precision.

Check your understanding: Which method has the smallest HC2 standard error in this dataset? Is it also closest to the true effect? Would either observation necessarily hold after a new random assignment?