Suppose we conduct an experiment with complete randomization and obtain the assignment \(W = (1,1,0,1,1,0,0,1)\), and observed outcomes \(Y^{\text{obs}} = (12, 7, 2, 7, 6, 6, 2, 7)\).

1. Test the sharp null of no individual effect

Under complete randomization, exactly five of the eight units receive treatment. We test \[H_0(0): Y_i(1)=Y_i(0) \quad \text{for every unit }i.\]

Step 1: Fill in the missing potential outcomes

Under this null, each unit has the same outcome under treatment and control. Its observed outcome therefore supplies both entries in its potential-outcome table. For example, unit 1 reveals \(Y_1(1)=12\), so we impute \(Y_1(0)=12\). Unit 3 reveals \(Y_3(0)=2\), so we impute \(Y_3(1)=2\).

W <- c(1, 1, 0, 1, 1, 0, 0, 1)
Yobs <- c(12, 7, 2, 7, 6, 6, 2, 7)
n <- length(W)
n1 <- sum(W)

Y0_null <- Yobs
Y1_null <- Yobs
print(data.frame(unit = seq_len(n), W, Yobs, Y0_null, Y1_null))
##   unit W Yobs Y0_null Y1_null
## 1    1 1   12      12      12
## 2    2 1    7       7       7
## 3    3 0    2       2       2
## 4    4 1    7       7       7
## 5    5 1    6       6       6
## 6    6 0    6       6       6
## 7    7 0    2       2       2
## 8    8 1    7       7       7

Step 2: List the assignments allowed by the design

Before running the code: How many ways can we choose five treated units from eight? Would an assignment with four treated units be allowed?

combn(n, n1) lists all possible treated sets. Each column contains the indices of five treated units. We convert each set into a column of zeros and ones, giving \(\binom{8}{5}=56\) equally likely assignments.

treated_all <- combn(n, n1)
permuted_assignments <- apply(treated_all, 2, function(idx) {
  w <- rep(0, n)
  w[idx] <- 1
  w
})

# Rows are units; columns are possible assignment vectors.
permuted_assignments[, 1:3]
##      [,1] [,2] [,3]
## [1,]    1    1    1
## [2,]    1    1    1
## [3,]    1    1    1
## [4,]    1    1    1
## [5,]    1    0    0
## [6,]    0    1    0
## [7,]    0    0    1
## [8,]    0    0    0
ncol(permuted_assignments)
## [1] 56

Step 3: Calculate the statistic for each assignment

Let \(D\) be the treated-minus-control difference in means. Under the zero-effect null, use \(T=|D|\) to measure departures from zero in either direction. The observed difference is \(D_{\mathrm{obs}}=39/5-10/3=67/15\).

For each new assignment w, select outcomes from the imputed table and calculate the same difference. Here Y_new always equals Yobs, but the groups being averaged change.

D_obs <- mean(Yobs[W == 1]) - mean(Yobs[W == 0])
T_obs <- abs(D_obs)

# First examine one hypothetical assignment.
w <- permuted_assignments[, 1]
Y_new <- ifelse(w == 1, Y1_null, Y0_null)
D_new <- mean(Y_new[w == 1]) - mean(Y_new[w == 0])
c(D_new = D_new, T_new = abs(D_new))
## D_new T_new 
##   1.8   1.8
# apply(..., 2, ...) repeats the calculation for every column.
D_rand <- apply(permuted_assignments, 2, function(w) {
  Y_new <- ifelse(w == 1, Y1_null, Y0_null)
  mean(Y_new[w == 1]) - mean(Y_new[w == 0])
})
T_rand <- abs(D_rand)

Step 4: Count results at least as extreme as the observed result

Each allowed assignment has probability \(1/56\), so the exact p-value is the proportion with \(T(w)\geq T_{\mathrm{obs}}\). Include ties. In R, mean() of a logical vector is the proportion of entries that are TRUE.

hist(D_rand, breaks = 10,
     main = "Randomization distribution under no individual effect",
     xlab = "Difference in means")
abline(v = c(-T_obs, T_obs), col = "red", lty = 2)

# A small tolerance prevents numerical rounding from excluding ties.
extreme <- T_rand >= T_obs - 1e-10
sum(extreme)
## [1] 2
p <- mean(extreme)
print(p)  # 2 / 56, approximately 0.0357
## [1] 0.03571429

We reject the sharp null at the 5% level. This provides evidence against the claim that every individual effect is zero; it does not establish that every unit benefits.

Check your understanding: What stays fixed across assignments? Why do we count both positive and negative differences?

2. Test the sharp null that every individual effect is 1

The candidate null is \[H_0(1): Y_i(1)-Y_i(0)=1 \quad \text{for every unit }i.\] This assumption lets us fill in each missing potential outcome. The filled-in values are what the null implies, not additional observations.

Step 1: Start with what we actually observed

The original assignment W tells us which potential outcome was observed:

  • If W[i] == 1, we observed Y1[i]. Keep it and subtract 1 to obtain Y0[i].
  • If W[i] == 0, we observed Y0[i]. Keep it and add 1 to obtain Y1[i].
Unit Original assignment Observed outcome Under the null of effect 1
1 Treated \(Y_1(1)=12\) Impute \(Y_1(0)=12-1=11\)
3 Control \(Y_3(0)=2\) Impute \(Y_3(1)=2+1=3\)

Before running the code: What are the two potential outcomes for unit 2? For unit 6?

tau0 <- 1

# NA marks an outcome we did not observe.
Y0_observed <- ifelse(W == 0, Yobs, NA)
Y1_observed <- ifelse(W == 1, Yobs, NA)

# Keep observed values; fill only the missing values under H0(1).
Y0_null <- Y0_observed
Y1_null <- Y1_observed
Y0_null[W == 1] <- Yobs[W == 1] - tau0
Y1_null[W == 0] <- Yobs[W == 0] + tau0

imputation_table <- data.frame(
  unit = seq_len(n), W_observed = W, Yobs = Yobs,
  Y0_observed, Y1_observed, Y0_null, Y1_null
)
print(imputation_table)
##   unit W_observed Yobs Y0_observed Y1_observed Y0_null Y1_null
## 1    1          1   12          NA          12      11      12
## 2    2          1    7          NA           7       6       7
## 3    3          0    2           2          NA       2       3
## 4    4          1    7          NA           7       6       7
## 5    5          1    6          NA           6       5       6
## 6    6          0    6           6          NA       6       7
## 7    7          0    2           2          NA       2       3
## 8    8          1    7          NA           7       6       7

Check the completed table: every row has Y1_null - Y0_null = 1, and the original assignment recovers Yobs.

stopifnot(all(Y1_null - Y0_null == tau0))
stopifnot(all(ifelse(W == 1, Y1_null, Y0_null) == Yobs))

For any candidate effect \(\tau_0\), the same rule can be written compactly as \[Y_i(0)=Y_i^{\mathrm{obs}}-W_i^{\mathrm{obs}}\tau_0, \qquad Y_i(1)=Y_i(0)+\tau_0.\]

# Equivalent to the explicit missing-value imputation above:
Y0_null <- Yobs - W * tau0
Y1_null <- Y0_null + tau0

Step 2: Hold this table fixed and try a new assignment

Impute once for the candidate null using the original assignment W. Then use each hypothetical assignment w to select an outcome from that fixed table. Do not redo the imputation using w in place of W.

For example, assign treatment to units 1 through 5. Unit 3 switches from control to treatment, so its hypothetical observed outcome becomes 3. Unit 8 switches from treatment to control, so its hypothetical observed outcome becomes 6.

w <- c(1, 1, 1, 1, 1, 0, 0, 0)
Y_new <- ifelse(w == 1, Y1_null, Y0_null)

print(data.frame(
  unit = seq_len(n), W_observed = W, Yobs,
  Y0_null, Y1_null, w_new = w, Y_new
))
##   unit W_observed Yobs Y0_null Y1_null w_new Y_new
## 1    1          1   12      11      12     1    12
## 2    2          1    7       6       7     1     7
## 3    3          0    2       2       3     1     3
## 4    4          1    7       6       7     1     7
## 5    5          1    6       5       6     1     6
## 6    6          0    6       6       7     0     6
## 7    7          0    2       2       3     0     2
## 8    8          1    7       6       7     0     6
D_new <- mean(Y_new[w == 1]) - mean(Y_new[w == 0])
D_new
## [1] 2.333333
D_new - tau0
## [1] 1.333333

Under a nonzero null, a unit’s observed outcome changes when its treatment changes. Under the zero-effect null in Section 1, the two potential outcomes coincide, so the outcomes stay the same under every assignment.

Step 3: Compare deviations from the hypothesized effect

There are three separate quantities to compute for each assignment \(w\):

  1. Hypothetical observed outcomes: \(Y_i^{\mathrm{new}}(w)=w_iY_i(1)+(1-w_i)Y_i(0)\), selected from the imputed table.
  2. Difference in means: \(D(w)=\operatorname{mean}\{Y_i(1):w_i=1\}-\operatorname{mean}\{Y_i(0):w_i=0\}\).
  3. Two-sided test statistic: \(T(w)=|D(w)-\tau_0|\).

Why subtract \(\tau_0\)? Under complete randomization and the candidate null, the randomization expectation of \(D(w)\) is \(\tau_0\). A difference close to the hypothesized effect agrees with the null; a large deviation in either direction is evidence against it. We use the same rule for the original and hypothetical assignments.

For \(\tau_0=1\):

Quantity Original assignment: treat 1, 2, 4, 5, 8 New assignment: treat 1, 2, 3, 4, 5
Treated outcomes 12, 7, 7, 6, 7 12, 7, 3, 7, 6
Control outcomes 2, 6, 2 6, 2, 6
\(D\) \(39/5-10/3=67/15\) \(35/5-14/3=7/3\)
\(D-1\) \(52/15\) \(4/3\)
\(T=|D-1|\) \(52/15\approx3.467\) \(4/3\approx1.333\)

This new assignment is less extreme than the observed assignment. To get the p-value, repeat this comparison for all 56 allowed assignments.

Why not average the individual effects in the imputed table? They all equal 1 by construction, so their average cannot distinguish assignments. The observed difference in means compares different groups of units and varies as the assignment changes.

D_obs <- mean(Yobs[W == 1]) - mean(Yobs[W == 0])
tau_obs <- D_obs - tau0

D_rand <- apply(permuted_assignments, 2, function(w) {
  # Select outcomes from the SAME imputed table for each assignment.
  Y_new <- ifelse(w == 1, Y1_null, Y0_null)
  mean(Y_new[w == 1]) - mean(Y_new[w == 0])
})
tau_rand <- D_rand - tau0
T_obs <- abs(tau_obs)
T_rand <- abs(tau_rand)

hist(tau_rand, breaks = 10,
     main = "Randomization distribution under an effect of 1",
     xlab = "Difference in means minus 1")
abline(v = c(-abs(tau_obs), abs(tau_obs)), col = "red", lty = 2)

# Include ties; tolerance protects against floating-point rounding.
p <- mean(T_rand >= T_obs - 1e-10)
print(p)  # 11 / 56, approximately 0.1964
## [1] 0.1964286

Check your understanding: Why would shuffling treatment labels while keeping Yobs fixed give the wrong reference distribution for this nonzero null?

3. Confidence interval for a common additive treatment effect

Assume \(Y_i(1)-Y_i(0)=\tau\) for every unit. This section constructs a confidence interval for that common additive effect.

Intuition: Try each candidate effect and keep those the data do not reject. Together, these values form the confidence set.

For this bracketing method, assume the exact nonrejected set is one interval, the grid hits it, and rejected grid points bracket it. These conditions ensure that the reported interval encloses the exact set; a finite grid alone does not establish the conditions.

Why does inversion give coverage? If the true common effect is \(\tau^*\), the randomization test rejects it with probability at most 0.05. Thus \[\Pr_W\{\tau^*\in C_{0.95}(W)\} =1-\Pr_W\{p_W(\tau^*)\leq0.05\}\geq0.95.\] Probability here is over repeated treatment assignments with potential outcomes fixed. An interval that encloses this exact set for every assignment has at least the same coverage.

Step 1: Turn the single-candidate test into a function

Sections 1 and 2 tested candidate effects 0 and 1. The function below repeats those same steps for any candidate tau: impute, reassign, and compare test statistics.

Each new candidate requires a new imputed table. Within that candidate’s test, the table stays fixed across assignments.

getPvalue <- function(tau, data, assignments) {
  # Always rebuild from the original observed data for this candidate.
  Y0_null <- data$Yobs - data$W * tau
  Y1_null <- Y0_null + tau
  tau_obs <- mean(data$Yobs[data$W == 1]) -
    mean(data$Yobs[data$W == 0]) - tau

  tau_rand <- apply(assignments, 2, function(w) {
    Y_new <- ifelse(w == 1, Y1_null, Y0_null)
    mean(Y_new[w == 1]) - mean(Y_new[w == 0]) - tau
  })
  mean(abs(tau_rand) >= abs(tau_obs) - 1e-10)
}

Step 2: Test a coarse grid and include one extra point at each end

sapply() calls the function once per candidate. Each row reports the test decision for that specific value. Reject when its p-value is at most 0.05.

data <- data.frame(W = W, Yobs = Yobs)
alpha <- 0.05
tau_seq <- seq(-15, 15, by = 1)
pvalue_seq <- sapply(tau_seq, getPvalue, data = data,
                     assignments = permuted_assignments)
print(data.frame(tau = tau_seq, p_value = pvalue_seq,
                 decision_at_5_percent = ifelse(pvalue_seq <= alpha,
                   "Reject", "Do not reject")))
##    tau    p_value decision_at_5_percent
## 1  -15 0.01785714                Reject
## 2  -14 0.01785714                Reject
## 3  -13 0.01785714                Reject
## 4  -12 0.01785714                Reject
## 5  -11 0.01785714                Reject
## 6  -10 0.01785714                Reject
## 7   -9 0.01785714                Reject
## 8   -8 0.01785714                Reject
## 9   -7 0.01785714                Reject
## 10  -6 0.01785714                Reject
## 11  -5 0.01785714                Reject
## 12  -4 0.01785714                Reject
## 13  -3 0.01785714                Reject
## 14  -2 0.01785714                Reject
## 15  -1 0.01785714                Reject
## 16   0 0.03571429                Reject
## 17   1 0.19642857         Do not reject
## 18   2 0.21428571         Do not reject
## 19   3 0.51785714         Do not reject
## 20   4 0.83928571         Do not reject
## 21   5 0.91071429         Do not reject
## 22   6 0.44642857         Do not reject
## 23   7 0.19642857         Do not reject
## 24   8 0.10714286         Do not reject
## 25   9 0.10714286         Do not reject
## 26  10 0.05357143         Do not reject
## 27  11 0.01785714                Reject
## 28  12 0.01785714                Reject
## 29  13 0.01785714                Reject
## 30  14 0.01785714                Reject
## 31  15 0.01785714                Reject

The tested integers 1 through 10 are not rejected. The neighboring grid points 0 and 11 are rejected.

Before running the code: Why might using only the smallest and largest nonrejected grid values cut off part of the confidence set?

Values between 0 and 1 and between 10 and 11 have not yet been tested. Include the neighboring grid point at each end to obtain the coarse conservative approximation \([0,11]\).

# Find the first and last nonrejected grid values.
not_rejected <- which(pvalue_seq > alpha)
stopifnot(length(not_rejected) > 0)
first <- min(not_rejected)
last <- max(not_rejected)

# We need a tested point beyond each end; otherwise expand the search.
stopifnot(first > 1, last < length(tau_seq))
CI_coarse <- c(lower = tau_seq[first - 1], upper = tau_seq[last + 1])
CI_coarse  # [0, 11]
## lower upper 
##     0    11

Zero and 11 are excluded from the exact inverted set, but may be included in an outer approximation to it. The test-decision column describes individual hypotheses; it does not say whether a point can be an endpoint of the wider reported interval.

For example, on a grid \(-5,5,15\), only 5 is not rejected. Reporting \([5,5]\) would ignore the untested values on either side. The extra-point rule instead gives \([-5,15]\), prompting us to refine the grid.

Step 3: Test more values near the boundaries

Refine each boundary region using increments of 0.1. The potential-outcome table is rebuilt for each new candidate, just as in Section 2.

lower_grid <- seq(CI_coarse["lower"], tau_seq[first], by = 0.1)
upper_grid <- seq(tau_seq[last], CI_coarse["upper"], by = 0.1)
p_lower <- sapply(lower_grid, getPvalue, data = data,
                   assignments = permuted_assignments)
p_upper <- sapply(upper_grid, getPvalue, data = data,
                   assignments = permuted_assignments)

print(data.frame(tau = c(lower_grid, upper_grid),
                 p_value = c(p_lower, p_upper),
                 decision_at_5_percent = ifelse(c(p_lower, p_upper) <= alpha,
                   "Reject", "Do not reject")))
##     tau    p_value decision_at_5_percent
## 1   0.0 0.03571429                Reject
## 2   0.1 0.03571429                Reject
## 3   0.2 0.03571429                Reject
## 4   0.3 0.03571429                Reject
## 5   0.4 0.03571429                Reject
## 6   0.5 0.03571429                Reject
## 7   0.6 0.03571429                Reject
## 8   0.7 0.03571429                Reject
## 9   0.8 0.03571429                Reject
## 10  0.9 0.03571429                Reject
## 11  1.0 0.19642857         Do not reject
## 12 10.0 0.05357143         Do not reject
## 13 10.1 0.01785714                Reject
## 14 10.2 0.01785714                Reject
## 15 10.3 0.01785714                Reject
## 16 10.4 0.01785714                Reject
## 17 10.5 0.01785714                Reject
## 18 10.6 0.01785714                Reject
## 19 10.7 0.01785714                Reject
## 20 10.8 0.01785714                Reject
## 21 10.9 0.01785714                Reject
## 22 11.0 0.01785714                Reject

At the lower boundary, 0.9 is rejected and 1 is not rejected. At the upper boundary, 10 is not rejected and 10.1 is rejected. Again include one extra tested point on each side:

first_lower <- min(which(p_lower > alpha))
last_upper <- max(which(p_upper > alpha))
stopifnot(first_lower > 1, last_upper < length(upper_grid))
CI_fine <- c(lower = lower_grid[first_lower - 1],
             upper = upper_grid[last_upper + 1])
CI_fine  # [0.9, 10.1]
## lower upper 
##   0.9  10.1

Using the 0.1 grid and including one extra tested point at each end, we report \([0.9,10.1]\), the same conservative confidence interval shown on Lecture 3, slide 31. For a finer interval, test more candidate effects near the boundaries.

Check your understanding: Does a rejected endpoint contradict the definition of exact inversion? Explain the difference between the exact set and its grid approximation.

4. Monte Carlo when enumeration is expensive

For a large experiment, draw assignments from the design instead of enumerating them all. The imputation is unchanged: construct one table for the candidate effect, then reuse it for every simulated assignment.

Step 1: Draw one assignment from the design

sample(data$W) randomly rearranges the original zeros and ones. Every draw has exactly five treated units, and each of the 56 assignment vectors is equally likely. Separate draws can repeat an assignment.

Before running the code: Could this command produce an assignment with six treated units?

set.seed(0)
w <- sample(data$W)
w
## [1] 0 1 1 0 1 1 0 1
sum(w)
## [1] 5

Step 2: Repeat the same imputation and statistic calculation

Compare this function with getPvalue() in Section 3. The imputed table and statistic are the same. replicate(M, ...) draws and evaluates M random assignments in place of apply() over all assignments.

The imputation stays outside the simulation loop: the candidate effect stays fixed throughout this test.

getPvalueMC <- function(data, tau, M = 10000) {
  Y0_null <- data$Yobs - data$W * tau
  Y1_null <- Y0_null + tau
  tau_obs <- mean(data$Yobs[data$W == 1]) -
    mean(data$Yobs[data$W == 0]) - tau

  tau_rand <- replicate(M, {
    # Complete randomization: preserve the number treated.
    w <- sample(data$W)
    Y_new <- ifelse(w == 1, Y1_null, Y0_null)
    mean(Y_new[w == 1]) - mean(Y_new[w == 0]) - tau
  })
  # Count simulated statistics at least as extreme as the observed one.
  b <- sum(abs(tau_rand) >= abs(tau_obs) - 1e-10)
  (b + 1) / (M + 1)
}

Step 3: Compare simulated and exact p-values

If b of the M simulated assignments are at least as extreme as the observed assignment, report \((b+1)/(M+1)\). The plus-one correction includes the observed assignment in the comparison and prevents a zero Monte Carlo p-value. Here enumeration is still feasible, so we can check the simulation against the exact answers.

set.seed(0)
print(data.frame(
  candidate_effect = c(0, 1),
  exact = sapply(c(0, 1), getPvalue, data = data,
                 assignments = permuted_assignments),
  Monte_Carlo = sapply(c(0, 1), function(tau) {
    getPvalueMC(data, tau = tau, M = 10000)
  })
))
##   candidate_effect      exact Monte_Carlo
## 1                0 0.03571429  0.03329667
## 2                1 0.19642857  0.19408059

Here sample(data$W) implements complete randomization. An independent Bernoulli design would require drawing each treatment indicator with its specified probability instead. Increasing M improves simulation accuracy; it does not increase the number of study units.

Check your understanding: If you increase M from 1,000 to 10,000, what improves? Does this add study units or change the exact p-value?