The HOMEFOOD trial assigned 106 older adults discharged from hospital to nutrition therapy with home-delivered food and supplements, or to standard care. We focus on the follow-up Short Physical Performance Battery (SPPB) score, a measure of physical function. See the original study for the trial background.

This example brings together the methods from Examples 1 and 2. We estimate an average treatment effect using Neyman’s approach, adjust for pre-treatment covariates using post-stratification and regression, and test a sharp null using Fisher’s approach. Unlike Example 2, this is a real dataset: we never observe both potential outcomes or know the true treatment effect.

1. Read and inspect the trial data

Step 1: Set up R

Keep HOMEFOOD.tab in the same folder as this R Markdown file. We use ggplot2 for plots and sandwich for HC2 heteroskedasticity-robust standard errors, as in Example 2. If needed, run install.packages(c("ggplot2", "sandwich")) once in the R console.

library(ggplot2)
library(sandwich)

data_all <- read.table("HOMEFOOD.tab", sep = "\t", header = TRUE)
dim(data_all)
## [1] 106  65
table(data_all$group)
## 
##  0  1 
## 53 53

group == 1 denotes the intervention and group == 0 denotes control. The file contains 53 participants in each arm. The suffixes _0 and _1 indicate baseline and follow-up measurements, respectively: sppb_0 and sppb_1 are not the two potential outcomes. Both are observed measurements taken at different times.

Step 2: Identify the participants included in this analysis

Two participants have missing follow-up SPPB scores, one in each arm. We inspect those rows before restricting the analysis to participants with observed outcomes.

missing_outcome <- is.na(data_all$sppb_1)
data_all[missing_outcome, c("ID", "group", "sex", "BMI_1", "sppb_1")]
# In this file, this selects the same participants as filtering on BMI_1.
stopifnot(identical(missing_outcome, is.na(data_all$BMI_1)))
data <- data_all[!missing_outcome, ]
data <- data[, c("group", "sex", "age", "education_2_cat", "smoker",
                 "alcohol", "nr_of_medicines", "BMI_0", "sppb_0",
                 "ICD10", "MMSE_0", "ISNST_0", "sppb_1")]

stopifnot(!anyNA(data))
n <- nrow(data)
n1 <- sum(data$group == 1)
n0 <- sum(data$group == 0)
c(n = n, treated = n1, control = n0)
##       n treated control 
##     104      52      52
head(data)

The analysis uses 104 participants: 52 treated and 52 controls. We analyze the 104 participants with observed follow-up scores as a randomized comparison, assuming that treatment assignment did not determine which participants had observed follow-up scores. This is an additional assumption: excluding participants after randomization does not automatically preserve randomization, even when equal numbers are excluded from the two arms.

Under that working assumption, our average-effect target is \[\tau=\frac{1}{n}\sum_{i=1}^n\{Y_i(1)-Y_i(0)\},\] where \(Y_i(w)\) is participant \(i\)’s follow-up SPPB score under assignment \(w\). The complete-case calculations alone do not identify the effect for all 106 randomized participants.

Check your understanding: Why is a participant’s baseline score not their untreated potential outcome at follow-up? Why might dropping participants with missing outcomes affect causal interpretation?

2. Inspect pre-treatment covariate balance

Step 1: Compare the sex composition of the two arms

Sex is coded as 1 for female and 2 for male. Create a factor for display and a female indicator whose mean is easy to interpret.

data$sex_label <- factor(data$sex, levels = c(1, 2),
                         labels = c("Female", "Male"))
data$female <- as.numeric(data$sex == 1)
cell_counts <- table(Sex = data$sex_label, Treatment = data$group)
cell_counts
##         Treatment
## Sex       0  1
##   Female 28 38
##   Male   24 14
prop.table(cell_counts, margin = 2)
##         Treatment
## Sex              0         1
##   Female 0.5384615 0.7307692
##   Male   0.4615385 0.2692308

Before running the next code: Does an equal number of treated and control participants guarantee equal proportions of female participants? Which sex is more common in the intervention arm?

Step 2: Summarize the other baseline characteristics

The following code repeats the same calculation for each pre-treatment covariate: compare the two means, estimate the standard error of their difference, and calculate a two-sided normal-approximation p-value.

covariates <- c("female", "age", "education_2_cat", "smoker", "alcohol",
                "nr_of_medicines", "BMI_0", "sppb_0", "ICD10",
                "MMSE_0", "ISNST_0")

balance_results <- t(sapply(covariates, function(v) {
  x_t <- data[[v]][data$group == 1]
  x_c <- data[[v]][data$group == 0]
  se <- sqrt(var(x_t) / length(x_t) + var(x_c) / length(x_c))
  z <- (mean(x_t) - mean(x_c)) / se
  c(treated_mean = mean(x_t), treated_SD = sd(x_t),
    control_mean = mean(x_c), control_SD = sd(x_c),
    z_statistic = z, p_value = 2 * pnorm(-abs(z)))
}))
signif(balance_results, 3)
##                 treated_mean treated_SD control_mean control_SD z_statistic
## female                 0.731      0.448        0.538      0.503      2.0600
## age                   83.300      6.740       81.800      6.030      1.2100
## education_2_cat        0.692      0.466        0.673      0.474      0.2090
## smoker                 1.960      0.194        1.900      0.298      1.1700
## alcohol                1.620      0.491        1.540      0.503      0.7890
## nr_of_medicines       11.800      4.640       12.300      4.230     -0.6410
## BMI_0                 28.500      6.510       27.000      5.330      1.2800
## sppb_0                 2.500      1.840        2.380      1.960      0.3090
## ICD10                 10.300      4.930       10.400      3.800     -0.0891
## MMSE_0                26.100      2.730       25.800      2.880      0.4200
## ISNST_0                5.100      1.730        4.520      1.260      1.9400
##                 p_value
## female           0.0396
## age              0.2260
## education_2_cat  0.8350
## smoker           0.2420
## alcohol          0.4300
## nr_of_medicines  0.5220
## BMI_0            0.2000
## sppb_0           0.7570
## ICD10            0.9290
## MMSE_0           0.6750
## ISNST_0          0.0519

For female, the means are proportions. For the other coded binary variables, the table summarizes their stored numerical codes. These p-values are approximate, separate comparisons; they are not a pass/fail test of whether randomization worked. Chance imbalances are expected, and examining many covariates increases the chance of seeing a small p-value.

We use sex and baseline SPPB to illustrate adjustment.

Check your understanding: Does a small balance p-value establish that treatment was assigned nonrandomly? Does a large one establish that the groups are identical?

3. Estimate the average effect without adjustment

Step 1: Calculate the treated-minus-control difference in means

The unadjusted estimate is \[\widehat\tau=\overline Y_1^{\mathrm{obs}}-\overline Y_0^{\mathrm{obs}}.\] Under the complete-randomization assumptions above, this estimator is unbiased for the study average effect. Individual treatment effects need not be equal.

Y_t <- data$sppb_1[data$group == 1]
Y_c <- data$sppb_1[data$group == 0]
tau_hat <- mean(Y_t) - mean(Y_c)
c(treated_mean = mean(Y_t), control_mean = mean(Y_c), estimate = tau_hat)
## treated_mean control_mean     estimate 
##     3.826923     2.788462     1.038462

Step 2: Use Neyman’s standard error for approximate inference

As in Example 2, estimate the variance using \[\widehat V=\frac{s_1^2}{n_1}+\frac{s_0^2}{n_0}.\] Its square root is a standard error (SE) of the effect estimate, not a standard deviation (SD) of participants’ scores. This variance estimator omits an unidentifiable treatment-effect-heterogeneity term and is conservative in expectation under complete randomization.

For a two-sided test of \(H_0:\tau=0\), use \(z=\widehat\tau/\widehat{\mathrm{SE}}\). We use the same normal approximation for all average-effect confidence intervals and p-values below.

se_Neyman <- sqrt(var(Y_t) / n1 + var(Y_c) / n0)
alpha <- 0.05
critical_value <- qnorm(1 - alpha / 2)

# Reuse the same reporting rule for each average-effect estimator.
summarize_effect <- function(estimate, se) {
  c(estimate = estimate, SE = se,
    CI_lower = estimate - critical_value * se,
    CI_upper = estimate + critical_value * se,
    p_value = 2 * pnorm(-abs(estimate / se)))
}

result_unadjusted <- summarize_effect(tau_hat, se_Neyman)
signif(result_unadjusted, 3)
## estimate       SE CI_lower CI_upper  p_value 
##   1.0400   0.4190   0.2170   1.8600   0.0132

The estimated difference is 1.04 SPPB points, with an approximate 95% interval from 0.22 to 1.86. An average-effect null permits positive and negative individual effects that cancel. It is weaker than the sharp null that every individual effect is zero.

Check your understanding: Does a positive estimated average effect imply that every participant benefits?

4. Post-stratify by sex and connect the result to regression

Step 1: Compare treatment and control within each sex

Post-stratification uses a pre-treatment characteristic to define subgroups after assignment. For this analysis, we treat the experiment as stratified by sex, using the observed numbers of treated and control participants in each stratum. Within each sex, we apply the difference in means and Neyman variance calculation from Section 3.

For sex stratum \(k\), calculate its difference in means \(\widehat\tau_k\) and Neyman variance estimate \(\widehat V_k=s_{1k}^2/n_{1k}+s_{0k}^2/n_{0k}\).

strata <- data$sex_label
stratum_levels <- levels(strata)
stopifnot(all(cell_counts > 1))  # Each cell needs a mean and a sample variance.

stratum_results <- t(sapply(stratum_levels, function(k) {
  in_stratum <- strata == k
  y_t <- data$sppb_1[in_stratum & data$group == 1]
  y_c <- data$sppb_1[in_stratum & data$group == 0]
  c(n = sum(in_stratum), n1 = length(y_t), n0 = length(y_c),
    estimate = mean(y_t) - mean(y_c),
    variance = var(y_t) / length(y_t) + var(y_c) / length(y_c))
}))
stratum_results <- data.frame(sex = stratum_levels, stratum_results,
                             row.names = NULL)
stratum_results$SE <- sqrt(stratum_results$variance)
stratum_results

Step 2: Combine the subgroup estimates using sample proportions

Our target averages over participants, so weight each stratum by its share of all 104 analysis participants: \[\widehat\tau_{\mathrm{PS}}=\sum_k\frac{n_k}{n}\widehat\tau_k, \qquad \widehat V_{\mathrm{PS}}=\sum_k\left(\frac{n_k}{n}\right)^2\widehat V_k.\]

Before running the code: Should female and male participants receive equal stratum weights? Why are the weights squared in the variance formula?

weights <- stratum_results$n / n
names(weights) <- stratum_levels
weights
##    Female      Male 
## 0.6346154 0.3653846
tau_PS <- sum(weights * stratum_results$estimate)
se_PS <- sqrt(sum(weights^2 * stratum_results$variance))
result_PS <- summarize_effect(tau_PS, se_PS)
signif(result_PS, 3)
## estimate       SE CI_lower CI_upper  p_value 
##   1.1000   0.4260   0.2640   1.9300   0.0099

The stratum weights are 0.635 for female and 0.365 for male participants. Equal weights would target an average giving the two sex groups equal importance, regardless of their sizes. Precision-based weights would also change the target when subgroup effects differ.

Step 3: Fit the equivalent centered interaction regression

Example 2 centered covariate indicators to make the treatment coefficient represent an overall average fitted effect. Apply the same idea here: center the female indicator using the full analysis sample, then include its interaction with treatment.

data$female_centered <- data$female - mean(data$female)
lm_sex <- lm(sppb_1 ~ group * female_centered, data = data)
coef(lm_sex)
##           (Intercept)                 group       female_centered 
##             2.8176511             1.0978347             0.3035714 
## group:female_centered 
##            -1.2246241
est_sex <- unname(coef(lm_sex)["group"])
se_sex <- sqrt(vcovHC(lm_sex, type = "HC2")["group", "group"])
result_sex <- summarize_effect(est_sex, se_sex)
signif(result_sex, 3)
## estimate       SE CI_lower CI_upper  p_value 
##   1.1000   0.4260   0.2640   1.9300   0.0099

group * female_centered expands to the two main effects and their interaction. The model allows separate treatment-control gaps for female and male participants. Because the centered indicator averages to zero, the coefficient of group equals the sample-proportion-weighted average of those gaps.

With a binary stratum indicator and full interactions, the regression fits the four cell means exactly. Its treatment estimate and HC2 standard error equal the post-stratified calculations.

comparison_sex <- rbind(Post_stratification = result_PS,
                        Centered_interaction = result_sex)
signif(comparison_sex, 3)
##                      estimate    SE CI_lower CI_upper p_value
## Post_stratification       1.1 0.426    0.264     1.93  0.0099
## Centered_interaction      1.1 0.426    0.264     1.93  0.0099

An additive model without the interaction would impose a common fitted gap across sex groups and need not reproduce this estimate. The equality of the two rows above is an algebraic connection between estimators, not two independent pieces of evidence.

Check your understanding: Without centering the female indicator, whose effect would the treatment coefficient estimate? Why use HC2 rather than the conventional standard error from summary(lm_sex) for this comparison?

5. Add baseline SPPB and compare the average-effect estimates

Step 1: Adjust for a second pre-treatment covariate

Baseline physical function may explain substantial variation in follow-up physical function. Center baseline SPPB and add it, together with its treatment interaction, to the sex-adjusted model.

data$sppb_0_centered <- data$sppb_0 - mean(data$sppb_0)
lm_baseline <- lm(sppb_1 ~ group * (female_centered + sppb_0_centered),
                  data = data)
coef(lm_baseline)
##           (Intercept)                 group       female_centered 
##            2.82523361            1.05982742           -0.01903002 
##       sppb_0_centered group:female_centered group:sppb_0_centered 
##            0.66909931           -0.82426396           -0.27133392
est_baseline <- unname(coef(lm_baseline)["group"])
se_baseline <- sqrt(vcovHC(lm_baseline, type = "HC2")["group", "group"])
result_baseline <- summarize_effect(est_baseline, se_baseline)
signif(result_baseline, 3)
## estimate       SE CI_lower CI_upper  p_value 
##  1.06000  0.37200  0.33000  1.79000  0.00443

Both covariates average to zero over the analysis sample, so group again gives the average fitted treatment effect. The interactions allow the fitted effect to vary with sex and baseline SPPB. This is a regression adjustment for the follow-up outcome, not a comparison of each participant’s before-and-after change.

Under standard large-sample conditions, fully interacted covariate adjustment can improve precision under complete randomization. HC2 supports approximate inference; it does not make the regression model exact or guarantee a smaller estimated standard error in every dataset.

Step 2: Compare estimates, standard errors, and intervals

results <- data.frame(
  method = c("No adjustment", "Post-stratification by sex",
             "Regression: sex", "Regression: sex + baseline SPPB"),
  rbind(result_unadjusted, result_PS, result_sex, result_baseline),
  row.names = NULL
)
results

All four rows use the same participants and target the same average effect under the stated assumptions. The first two use Neyman variance estimates; the regression rows use HC2. All intervals and p-values use a normal approximation.

results$method <- factor(results$method, levels = rev(results$method))
ggplot(results, 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 = 0, color = "red", linetype = "dashed") +
  labs(title = "Average-effect estimates and approximate 95% intervals",
       subtitle = "Neyman or HC2 standard errors; dashed line marks no average effect",
       x = "Effect on follow-up SPPB score (points)", y = NULL) +
  theme_minimal()

Including baseline SPPB changes the estimated standard error from 0.43 to 0.37 points. Unlike the simulation in Example 2, we cannot check which estimate is closest to the true effect. The dashed line marks zero, not a known simulation truth.

Check your understanding: Does the smallest standard error establish that an estimate is closest to the true effect?

6. Test Fisher’s sharp null without post-stratification

Step 1: Impute potential outcomes under the candidate null

We now change the question. Fisher’s sharp null specifies every individual effect: \[H_0(\tau_0):Y_i(1)-Y_i(0)=\tau_0\quad\text{for every participant }i.\] For \(\tau_0=0\), each observed outcome supplies both potential outcomes. More generally, follow Example 1 and impute using the original assignment: \[Y_i(0)=Y_i^{\mathrm{obs}}-W_i^{\mathrm{obs}}\tau_0, \qquad Y_i(1)=Y_i(0)+\tau_0.\]

tau0 <- 0
Y0_null <- data$sppb_1 - data$group * tau0
Y1_null <- Y0_null + tau0
head(data.frame(W_observed = data$group, Yobs = data$sppb_1,
                Y0_null, Y1_null))

The imputed outcomes are what the null implies; they are not extra observations. Keep this completed table fixed as the assignment changes.

Step 2: Draw assignments and compare two-sided deviations

There are too many assignments to enumerate, so draw 10,000 from the complete-randomization model. sample(data$group) preserves 52 treated and 52 control participants. It does not preserve treated counts within sex groups.

For each assignment \(w\), calculate \(D(w)\), the treated-minus-control difference in means from the imputed outcomes, and compare \(|D(w)-\tau_0|\) with the observed value.

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

  tau_rand <- replicate(M, {
    w <- sample(data$group)
    Y_new <- ifelse(w == 1, Y1_null, Y0_null)
    mean(Y_new[w == 1]) - mean(Y_new[w == 0]) - tau
  })

  # Include ties, with a small tolerance for numerical rounding.
  b <- sum(abs(tau_rand) >= abs(tau_obs) - 1e-10)
  list(observed = tau_obs, simulated = tau_rand,
       p_value = (b + 1) / (M + 1))
}

set.seed(0)
fisher_unadjusted <- getPvalueMC(data, tau = tau0, M = 10000)
fisher_unadjusted$p_value
## [1] 0.01829817

As in Example 1, report \((b+1)/(M+1)\) when \(b\) of the \(M\) simulated assignments are at least as extreme as the observed assignment. This includes the observed assignment in the comparison and prevents reporting a zero Monte Carlo p-value.

Step 3: Inspect the randomization distribution

hist(fisher_unadjusted$simulated, breaks = 20,
     main = "Randomization distribution under no individual effect",
     xlab = "Difference in means", col = "grey85", border = "white")
abline(v = c(-abs(fisher_unadjusted$observed),
              abs(fisher_unadjusted$observed)), col = "red", lty = 2)

The dashed lines mark equally extreme deviations in both directions. The Monte Carlo p-value is 0.0183. It approximates the full randomization p-value; increasing the number of draws improves simulation accuracy, not the size of the trial.

7. Test Fisher’s sharp null after post-stratification

Both tests below shuffle treatment labels within sex; they differ in the statistic used to compare assignments.

7.1. Weighted within-sex difference in means

Step 1: Generate assignments using stratified randomization

After post-stratification by sex, we treat the experiment as stratified. When generating new treatment assignments for Fisher’s test, we therefore use a stratified randomization mechanism: shuffle treatment labels separately within each sex, preserving its observed numbers of treated and control participants.

We use the weighted within-sex difference as the test statistic. In every simulated assignment, the 66 female participants include 38 treated and 28 controls; the 38 male participants include 14 treated and 24 controls.

set.seed(0)
w <- data$group
for (k in stratum_levels) {
  idx <- which(strata == k)
  w[idx] <- sample(data$group[idx])
}
table(Sex = strata, Simulated_treatment = w)
##         Simulated_treatment
## Sex       0  1
##   Female 28 38
##   Male   24 14

Before running the test: What stays fixed in this simulation that was allowed to vary in Section 6? Do the stratum weights change across assignments?

Step 2: Repeat the test using the weighted statistic

For every candidate effect, rebuild the imputed table from the original data. Then, for every simulated assignment, calculate the within-stratum differences and combine them using the fixed sample proportions.

getPvalueMCStrata <- function(data, tau, strata, M = 10000) {
  strata <- factor(strata)
  stratum_levels <- levels(strata)
  weights <- prop.table(table(strata))
  stopifnot(all(table(strata, data$group) > 0))

  Y0_null <- data$sppb_1 - data$group * tau
  Y1_null <- Y0_null + tau

  # The same statistic is used for observed and simulated assignments.
  weighted_difference <- function(Y, w) {
    differences <- sapply(stratum_levels, function(k) {
      idx <- strata == k
      mean(Y[idx & w == 1]) - mean(Y[idx & w == 0])
    })
    sum(weights * differences)
  }
  tau_obs <- weighted_difference(data$sppb_1, data$group) - tau

  tau_rand <- replicate(M, {
    w <- data$group
    for (k in stratum_levels) {
      idx <- which(strata == k)
      w[idx] <- sample(data$group[idx])
    }
    Y_new <- ifelse(w == 1, Y1_null, Y0_null)
    weighted_difference(Y_new, w) - tau
  })

  b <- sum(abs(tau_rand) >= abs(tau_obs) - 1e-10)
  list(observed = tau_obs, simulated = tau_rand,
       p_value = (b + 1) / (M + 1))
}

set.seed(0)
fisher_PS <- getPvalueMCStrata(data, tau = tau0, strata = strata, M = 10000)
fisher_PS$p_value
## [1] 0.0109989
stopifnot(isTRUE(all.equal(fisher_PS$observed, tau_PS - tau0)))
hist(fisher_PS$simulated, breaks = 20,
     main = "Stratified randomization under no individual effect",
     xlab = "Weighted within-sex difference in means",
     col = "grey85", border = "white")
abline(v = c(-abs(fisher_PS$observed), abs(fisher_PS$observed)),
       col = "red", lty = 2)

7.2. Regression coefficient

Step 1: Use the estimated treatment coefficient

Fisher’s test can use the estimated treatment coefficient in place of the difference in means. Use the same model as in Section 5: regress follow-up SPPB on treatment, centered sex and baseline SPPB, and their treatment interactions.

For a hypothetical assignment \(w\), let \(\widehat\beta(w)\) be the fitted treatment coefficient. Our two-sided statistic is

\[T(w)=|\widehat\beta(w)-\tau_0|.\]

The coefficient measures the adjusted gap in SPPB points. Under the zero-effect null, we compare the absolute values of the observed and simulated coefficients.

The function below returns the signed deviation from the hypothesized effect. We take absolute values when counting results at least as extreme as the observed one.

covariates_reg <- data[c("female_centered", "sppb_0_centered")]

regression_statistic <- function(Y, w, covariates, tau = 0) {
  reg_data <- data.frame(Y = Y, group = w, covariates)
  fit <- lm(Y ~ group * (female_centered + sppb_0_centered),
            data = reg_data)
  estimate <- unname(coef(fit)["group"])
  estimate - tau
}

reg_observed <- regression_statistic(data$sppb_1, data$group,
                                     covariates_reg, tau = tau0)
reg_observed
## [1] 1.059827

For tau0 = 0, the observed coefficient is 1.06, the same estimate as in Section 5. We now compare it with coefficients obtained from new treatment assignments.

Step 2: Refit the regression for every simulated assignment

Use the same stratified randomization mechanism as in Section 7.1: shuffle treatment labels separately within each sex, preserving the observed treated and control counts. Hold sex, baseline SPPB, and their full-sample centering fixed. For each new assignment:

  1. Select hypothetical observed outcomes from the potential-outcome table imputed under the null.
  2. Replace the treatment indicator with the simulated assignment and refit the entire regression, including the interaction terms.
  3. Extract the treatment coefficient and subtract the hypothesized effect.

Under the zero-effect null, the outcomes stay fixed, but the treatment indicator and interaction terms change. Refit the regression for every assignment, rather than reusing the observed coefficient.

getPvalueMCRegression <- function(data, tau, strata, M = 10000) {
  strata <- factor(strata)
  stratum_levels <- levels(strata)
  covariates <- data[c("female_centered", "sppb_0_centered")]
  Y0_null <- data$sppb_1 - data$group * tau
  Y1_null <- Y0_null + tau
  tau_obs <- regression_statistic(data$sppb_1, data$group,
                                  covariates, tau)

  tau_rand <- replicate(M, {
    w <- data$group
    for (k in stratum_levels) {
      idx <- which(strata == k)
      w[idx] <- sample(data$group[idx])
    }
    Y_new <- ifelse(w == 1, Y1_null, Y0_null)
    regression_statistic(Y_new, w, covariates, tau)
  })

  b <- sum(abs(tau_rand) >= abs(tau_obs) - 1e-10)
  list(observed = tau_obs, simulated = tau_rand,
       p_value = (b + 1) / (M + 1))
}

set.seed(0)
fisher_regression <- getPvalueMCRegression(data, tau = tau0,
                                          strata = strata, M = 10000)
data.frame(
  statistic = "Treatment coefficient",
  observed = fisher_regression$observed,
  p_value = fisher_regression$p_value
)

Resetting the seed to 0 gives the same 10,000 within-sex assignments as the test in Section 7.1. The regression statistic additionally adjusts for baseline SPPB, so the p-values can differ even though the assignment mechanism is the same.

Step 3: Inspect the randomization distribution of the coefficient

hist(fisher_regression$simulated, breaks = 20,
     main = "Regression coefficient under no individual effect",
     xlab = "Treatment coefficient (SPPB points)",
     col = "grey85", border = "white")
abline(v = c(-abs(fisher_regression$observed),
              abs(fisher_regression$observed)), col = "red", lty = 2)

The coefficient-based Fisher p-value using stratified randomization is 0.0061.

Here regression supplies a covariate-adjusted statistic for testing the sharp null. Its randomization justification does not require the linear outcome model to be correct. For background on regression statistics in Fisher tests, see Zhao and Ding (2021).

Check your understanding: Why must we refit the regression for every simulated assignment?

8. Compare the p-values

test_comparison <- data.frame(
  method = c("Neyman: unadjusted", "Neyman: post-stratified",
             "Regression: sex + baseline SPPB (HC2 normal approximation)",
             "Fisher: unadjusted", "Fisher: post-stratified by sex",
             "Fisher: regression coefficient (sex + baseline SPPB)"),
  null = c(rep("Average effect = 0", 3),
           rep("Every individual effect = 0", 3)),
  p_value = c(result_unadjusted["p_value"], result_PS["p_value"],
              result_baseline["p_value"],
              fisher_unadjusted$p_value, fisher_PS$p_value,
              fisher_regression$p_value)
)
knitr::kable(test_comparison, digits = 4, row.names = FALSE,
             col.names = c("Method", "Null hypothesis", "p-value"))
Method Null hypothesis p-value
Neyman: unadjusted Average effect = 0 0.0132
Neyman: post-stratified Average effect = 0 0.0099
Regression: sex + baseline SPPB (HC2 normal approximation) Average effect = 0 0.0044
Fisher: unadjusted Every individual effect = 0 0.0183
Fisher: post-stratified by sex Every individual effect = 0 0.0110
Fisher: regression coefficient (sex + baseline SPPB) Every individual effect = 0 0.0061

Both regression rows use treatment, the centered sex indicator, centered baseline SPPB (sppb_0), and treatment interactions with both covariates.

All six p-values are below 0.05, so the analyses agree in rejecting their respective null hypotheses. In all three comparisons here, Neyman/HC2 gives a similar but slightly smaller p-value than Fisher: without adjustment (0.0132 versus 0.0183), after post-stratification by sex (0.0099 versus 0.011), and with regression adjustment for sex and baseline SPPB (0.0044 versus 0.0061).