In Youngstown, two first-grade classes in each of eight schools participated in the Children’s Television Workshop experiment. Within each school, one class was randomly assigned to watch The Electric Company. We observe class-average pretest and post-test scores.

We use the same data as Lecture 8 to test Fisher’s sharp null, estimate an average effect, and adjust for pretest scores. The class is the treatment unit; the school defines the pair. Our average-effect target gives each of the 16 classes equal weight.

1. Inspect the eight pairs

pre_C and pre_T are pretest scores; post_C and post_T are observed follow-up outcomes. A pretest score is not the class’s untreated potential outcome at follow-up. All calculations below use base R.

pairs <- data.frame(
  school = 1:8,
  pre_C  = c(12.9, 15.1, 16.8, 15.8, 13.9, 14.5, 17.0, 15.8),
  pre_T  = c(12.0, 12.3, 17.2, 18.9, 15.3, 16.6, 16.0, 20.1),
  post_C = c(54.6, 56.5, 75.2, 75.6, 55.3, 59.3, 87.0, 73.7),
  post_T = c(60.6, 55.5, 84.8, 101.9, 70.6, 78.4, 84.2, 108.6)
)
pairs$D <- pairs$post_T - pairs$post_C
pairs$dX <- pairs$pre_T - pairs$pre_C
J <- nrow(pairs)
pairs

Each row supplies a treatment-minus-control outcome difference \(D_j\) and pretest difference \(\Delta X_j\). Six outcome differences are positive and two are negative. Each compares two different classes, so it need not equal that school’s average causal effect.

2. Test Fisher’s sharp null

Step 1: Enumerate the assignments allowed by pairing

Under \(H_0:Y_{jk}(1)=Y_{jk}(0)\) for every class, the observed outcome supplies both potential outcomes, as in Example 1. Swapping the labels within a pair changes \(D_j\) to \(-D_j\).

Before running the code: How many assignments are possible? Why would shuffling all 16 treatment labels use the wrong design?

# Each row is one reassignment; each column is one school.
# +1 keeps the observed labels; -1 swaps them within that school.
flips <- as.matrix(expand.grid(rep(list(c(-1, 1)), J)))
nrow(flips)  # 2^8 = 256 equally likely assignments
## [1] 256
head(flips)
##      Var1 Var2 Var3 Var4 Var5 Var6 Var7 Var8
## [1,]   -1   -1   -1   -1   -1   -1   -1   -1
## [2,]    1   -1   -1   -1   -1   -1   -1   -1
## [3,]   -1    1   -1   -1   -1   -1   -1   -1
## [4,]    1    1   -1   -1   -1   -1   -1   -1
## [5,]   -1   -1    1   -1   -1   -1   -1   -1
## [6,]    1   -1    1   -1   -1   -1   -1   -1

Step 2: Compare three statistics using the same assignments

For outcome differences, the statistic is \(T=|\bar D|\). The other two statistics replace outcomes by ranks among all 16 classes, or by ranks within each pair. With two classes, the within-pair rank difference is sign(D): +1, -1, or 0 for a tie.

The function below averages the reassigned differences and counts statistics at least as large as the observed statistic. %*% multiplies each row of signs by the vector of differences and sums across pairs.

exact_test <- function(d) {
  T_obs <- abs(mean(d))
  T_rand <- abs(drop(flips %*% d) / J)
  extreme <- T_rand >= T_obs - 1e-10  # Include ties despite rounding.
  c(statistic = T_obs, extreme = sum(extreme), p_value = mean(extreme))
}

global_ranks <- rank(c(pairs$post_C, pairs$post_T))
rank_diff <- global_ranks[J + seq_len(J)] - global_ranks[seq_len(J)]
fisher_results <- rbind(
  Mean_difference = exact_test(pairs$D),
  Global_ranks = exact_test(rank_diff),
  Within_pair_ranks = exact_test(sign(pairs$D))
)
fisher_results
##                   statistic extreme   p_value
## Mean_difference      13.425       8 0.0312500
## Global_ranks          3.750       8 0.0312500
## Within_pair_ranks     0.500      74 0.2890625

The mean and global-rank statistics give \(8/256=0.03125\). Within-pair ranks give \(74/256\approx0.289\): they use the directions of the eight differences but discard their magnitudes. All three tests concern the same sharp null.

Step 3: See where the eight extreme assignments lie

D_rand <- drop(flips %*% pairs$D) / J
hist(D_rand, breaks = 20, col = "grey85", border = "white",
     main = "All 256 assignments under no individual effect",
     xlab = "Mean treatment-minus-control difference")
abline(v = c(-abs(mean(pairs$D)), abs(mean(pairs$D))),
       col = "red", lty = 2)

The dashed lines mark the observed absolute difference, 13.425, in both directions. Count assignments on or beyond either line. These are exact p-values, so no Monte Carlo correction is needed.

3. Estimate the average effect using Neyman’s approach

Step 1: Average the pair differences and estimate uncertainty

Now the target is an average effect, without assuming all individual effects are equal. Use \[\widehat\tau=\bar D,\qquad \widehat V=\frac{s_D^2}{J},\qquad \widehat{\mathrm{SE}}=\frac{s_D}{\sqrt J}.\] For fixed pairs, this variance estimator is conservative in expectation when pair-average effects differ.

tau_hat <- mean(pairs$D)
SE <- sd(pairs$D) / sqrt(J)
c(estimate = tau_hat, SE = SE)
##  estimate        SE 
## 13.425000  4.636337
intervals <- rbind(
  Normal = tau_hat + qnorm(c(0.025, 0.975)) * SE,
  Paired_t = tau_hat + qt(c(0.025, 0.975), df = J - 1) * SE
)
colnames(intervals) <- c("lower", "upper")
round(intervals, 2)
##          lower upper
## Normal    4.34 22.51
## Paired_t  2.46 24.39

The estimate is 13.43 points with SE 4.64. With only eight pairs, the normal approximation may be rough. The paired \(t\) interval uses a larger cutoff; it is exact under i.i.d. Normal pair differences, not from randomization alone.

Step 2: Compare with the standard error that ignores pairing

SE_ignore_pairs <- sqrt(var(pairs$post_T) / J + var(pairs$post_C) / J)
c(paired = SE, ignoring_pairs = SE_ignore_pairs)
##         paired ignoring_pairs 
##       4.636337       7.848152

Ignoring pairs leaves the point estimate unchanged but gives SE 7.85 here. This two-sample formula is shown for comparison; it does not use the actual paired assignment mechanism.

Check your understanding: Why can removing shared school variation reduce the SE? Do similar control potential outcomes within pairs imply equal treatment effects across pairs?

4. Adjust for pretest differences

Step 1: Fit one observation per school

With a common conditional treatment effect, fit \(D_j=\tau+\beta\Delta X_j+e_j\). The intercept estimates the effect at zero pretest imbalance. We also fit the full interaction model from Lecture 8 by adding the centered pair mean \(\bar X_j-\bar X\).

Before running the code: Treated classes started with higher pretest scores on average. If pretest scores predict post-test scores positively, would adjustment increase or decrease the estimated effect?

pairs$mean_X_centered <- (pairs$pre_T + pairs$pre_C) / 2 -
  mean(c(pairs$pre_T, pairs$pre_C))
fit0 <- lm(D ~ 1, data = pairs)
fit1 <- lm(D ~ dX, data = pairs)
fit2 <- lm(D ~ dX + mean_X_centered, data = pairs)

regression_row <- function(fit) {
  row <- coef(summary(fit))["(Intercept)", ]
  c(estimate = unname(row["Estimate"]),
    model_SE = unname(row["Std. Error"]), df = df.residual(fit))
}
round(rbind(Unadjusted = regression_row(fit0),
            Pretest_difference = regression_row(fit1),
            Difference_and_pair_mean = regression_row(fit2)), 3)
##                          estimate model_SE df
## Unadjusted                 13.425    4.636  7
## Pretest_difference          8.994    1.410  6
## Difference_and_pair_mean    8.522    1.475  5

The simple adjusted estimate is \(\bar D-\widehat\beta\,\overline{\Delta X}\approx8.99\). Adjustment removes the estimated contribution of the pretest imbalance. These conventional model SEs assume independent pair errors with common variance. The full model uses another degree of freedom; its estimated SE need not be smaller.

Step 2: Connect to regression on all 16 classes

We can instead include a separate intercept for each school, with slopes shared across schools. For the simple adjustment, the treatment coefficient equals the intercept from fit1.

classes <- data.frame(
  school = factor(rep(pairs$school, 2)), W = rep(c(0, 1), each = J),
  X = c(pairs$pre_C, pairs$pre_T), Y = c(pairs$post_C, pairs$post_T)
)
fit_units <- lm(Y ~ school + W + X, data = classes)
c(pair_regression = unname(coef(fit1)["(Intercept)"]),
  unit_regression = unname(coef(fit_units)["W"]))
## pair_regression unit_regression 
##        8.994303        8.994303

For fit2, the equivalent unit-level model adds W:I(X - mean(X)). Its treatment coefficient averages fitted effects over the observed covariates. For PATE inference, the full interaction model also requires accounting for sampling variability in the covariate mean; its conditional model SE alone does not do that.

Check your understanding: Why can we estimate common slopes using all pairs even though each individual pair contains only two classes?