Monte Carlo Simulations

Overview

This vignette validates CausalModel estimators via Monte Carlo simulation. For each estimator, we check:

Simulation 1: Observational Estimators (Binary Treatment)

We use generate_data() with known \(\tau = 10\) and vary sample size.

run_obs_sim <- function(N, n_reps, tau = 10) {
  estimators <- c("OLS", "IPW", "AIPW", "Matching")
  results <- list()
  for (est in estimators) {
    results[[est]] <- matrix(NA, nrow = n_reps, ncol = 2,
                              dimnames = list(NULL, c("ate", "se")))
  }

  for (i in seq_len(n_reps)) {
    data <- generate_data(N = N, k = 2, tau = tau)
    obs <- observational(data$Y, data$Z, data$X)

    r <- est_via_ols(obs)
    results[["OLS"]][i, ] <- c(r$ate, r$se)

    r <- est_via_ipw(obs)
    results[["IPW"]][i, ] <- c(r$ate, r$se)

    r <- est_via_aipw(obs)
    results[["AIPW"]][i, ] <- c(r$ate, r$se)

    r <- est_via_matching(obs, num_matches = 3)
    results[["Matching"]][i, ] <- c(r$ate, r$se)
  }
  results
}

# Simulation settings
tau <- 10
Ns <- c(500, 2000)
n_reps <- 500

sim1_results <- list()
for (N in Ns) {
  cat(sprintf("Running N = %d ...\n", N))
  sim1_results[[as.character(N)]] <- run_obs_sim(N, n_reps, tau)
}

Summary Table

make_summary <- function(results, tau) {
  rows <- list()
  for (est in names(results)) {
    ates <- results[[est]][, "ate"]
    ses <- results[[est]][, "se"]
    bias <- mean(ates) - tau
    rmse <- sqrt(mean((ates - tau)^2))
    mean_se <- mean(ses)
    emp_se <- sd(ates)
    ci_lo <- ates - 1.96 * ses
    ci_hi <- ates + 1.96 * ses
    coverage <- mean(ci_lo <= tau & tau <= ci_hi)
    rows[[est]] <- data.frame(
      Estimator = est, Bias = bias, RMSE = rmse,
      Mean_SE = mean_se, Emp_SE = emp_se, Coverage = coverage
    )
  }
  do.call(rbind, rows)
}

for (N in Ns) {
  cat(sprintf("\n### N = %d\n", N))
  tbl <- make_summary(sim1_results[[as.character(N)]], tau)
  print(knitr::kable(tbl, digits = 3, row.names = FALSE))
}
#> 
#> ### N = 500
#> 
#> 
#> |Estimator |  Bias|  RMSE| Mean_SE| Emp_SE| Coverage|
#> |:---------|-----:|-----:|-------:|------:|--------:|
#> |OLS       | 0.000| 0.068|   0.068|  0.068|    0.952|
#> |IPW       | 0.023| 0.195|   0.630|  0.194|    1.000|
#> |AIPW      | 0.004| 0.123|   0.118|  0.123|    0.952|
#> |Matching  | 0.170| 0.208|   0.117|  0.120|    0.680|
#> 
#> ### N = 2000
#> 
#> 
#> |Estimator |  Bias|  RMSE| Mean_SE| Emp_SE| Coverage|
#> |:---------|-----:|-----:|-------:|------:|--------:|
#> |OLS       | 0.002| 0.033|   0.034|  0.033|    0.960|
#> |IPW       | 0.010| 0.144|   0.323|  0.144|    1.000|
#> |AIPW      | 0.001| 0.066|   0.060|  0.066|    0.950|
#> |Matching  | 0.082| 0.102|   0.060|  0.060|    0.706|

Studentized Statistics

If the estimator and SE are consistent, \(T = (\hat{\tau} - \tau) / \widehat{SE}\) should be approximately standard normal.

par(mfrow = c(length(Ns), 4), mar = c(3, 3, 2, 1))
for (N in Ns) {
  res <- sim1_results[[as.character(N)]]
  for (est in names(res)) {
    ates <- res[[est]][, "ate"]
    ses <- res[[est]][, "se"]
    t_stats <- (ates - tau) / ses
    hist(t_stats, breaks = 30, freq = FALSE, col = "lightblue",
         main = sprintf("%s (N=%d)", est, N),
         xlab = "", xlim = c(-4, 4))
    curve(dnorm(x), add = TRUE, col = "red", lwd = 2)
    abline(v = 0, lty = 2)
  }
}

Simulation 2: Continuous Treatment (DML)

run_dml_sim <- function(N, n_reps, tau = 10) {
  results <- matrix(NA, nrow = n_reps, ncol = 2,
                     dimnames = list(NULL, c("ate", "se")))
  for (i in seq_len(n_reps)) {
    data <- generate_data_continuous(N = N, k = 2, tau = tau)
    obs <- observational(data$Y, data$Z, data$X)
    r <- est_via_dml(obs, k_folds = 2)
    results[i, ] <- c(r$ate, r$se)
  }
  results
}

tau_dml <- 10
Ns_dml <- c(1000, 5000)
n_reps_dml <- 500

sim2_results <- list()
for (N in Ns_dml) {
  cat(sprintf("Running DML N = %d ...\n", N))
  sim2_results[[as.character(N)]] <- run_dml_sim(N, n_reps_dml, tau_dml)
}

Summary Table

for (N in Ns_dml) {
  cat(sprintf("\n### DML, N = %d\n", N))
  ates <- sim2_results[[as.character(N)]][, "ate"]
  ses <- sim2_results[[as.character(N)]][, "se"]
  ci_lo <- ates - 1.96 * ses
  ci_hi <- ates + 1.96 * ses
  tbl <- data.frame(
    Bias = mean(ates) - tau_dml,
    RMSE = sqrt(mean((ates - tau_dml)^2)),
    Mean_SE = mean(ses),
    Emp_SE = sd(ates),
    Coverage = mean(ci_lo <= tau_dml & tau_dml <= ci_hi)
  )
  print(knitr::kable(tbl, digits = 3, row.names = FALSE))
}
#> 
#> ### DML, N = 1000
#> 
#> 
#> |   Bias|  RMSE| Mean_SE| Emp_SE| Coverage|
#> |------:|-----:|-------:|------:|--------:|
#> | -0.001| 0.032|   0.032|  0.032|    0.954|
#> 
#> ### DML, N = 5000
#> 
#> 
#> |   Bias|  RMSE| Mean_SE| Emp_SE| Coverage|
#> |------:|-----:|-------:|------:|--------:|
#> | -0.001| 0.014|   0.014|  0.014|    0.954|

Studentized Statistics

par(mfrow = c(1, length(Ns_dml)), mar = c(3, 3, 2, 1))
for (N in Ns_dml) {
  ates <- sim2_results[[as.character(N)]][, "ate"]
  ses <- sim2_results[[as.character(N)]][, "se"]
  t_stats <- (ates - tau_dml) / ses
  hist(t_stats, breaks = 30, freq = FALSE, col = "lightblue",
       main = sprintf("DML (N=%d)", N),
       xlab = "", xlim = c(-4, 4))
  curve(dnorm(x), add = TRUE, col = "red", lwd = 2)
  abline(v = 0, lty = 2)
}

Simulation 3: Experimental Estimators

We compare difference-in-means (DM) and ANCOVA under a completely randomized design. ANCOVA should achieve smaller variance when covariates are predictive.

N_exp <- 1000
tau_exp <- 5
n_reps_exp <- 500

dm_results <- matrix(NA, nrow = n_reps_exp, ncol = 2,
                      dimnames = list(NULL, c("ate", "se")))
ancova_results <- matrix(NA, nrow = n_reps_exp, ncol = 2,
                          dimnames = list(NULL, c("ate", "se")))

for (i in seq_len(n_reps_exp)) {
  X <- matrix(rnorm(N_exp * 2), ncol = 2)
  Z <- sample(c(rep(1, N_exp / 2), rep(0, N_exp / 2)))
  # Covariates strongly predict outcome
  Y <- tau_exp * Z + X %*% c(3, -2) + rnorm(N_exp)
  exp_obj <- experimental(Y, Z, X)

  r <- est_via_dm(exp_obj)
  dm_results[i, ] <- c(r$ate, r$se)

  r <- est_via_ancova(exp_obj)
  ancova_results[i, ] <- c(r$ate, r$se)
}

Summary Table

make_exp_summary <- function(results, tau, name) {
  ates <- results[, "ate"]
  ses <- results[, "se"]
  ci_lo <- ates - 1.96 * ses
  ci_hi <- ates + 1.96 * ses
  data.frame(
    Estimator = name,
    Bias = mean(ates) - tau,
    RMSE = sqrt(mean((ates - tau)^2)),
    Mean_SE = mean(ses),
    Emp_SE = sd(ates),
    Coverage = mean(ci_lo <= tau & tau <= ci_hi)
  )
}

tbl <- rbind(
  make_exp_summary(dm_results, tau_exp, "DM"),
  make_exp_summary(ancova_results, tau_exp, "ANCOVA")
)
knitr::kable(tbl, digits = 3, row.names = FALSE,
             caption = sprintf("Experimental estimators (N=%d, tau=%d)", N_exp, tau_exp))
Experimental estimators (N=1000, tau=5)
Estimator Bias RMSE Mean_SE Emp_SE Coverage
DM 0.008 0.232 0.237 0.232 0.954
ANCOVA 0.004 0.065 0.063 0.065 0.946

ANCOVA achieves a substantially smaller RMSE and empirical SE than DM when covariates are predictive, demonstrating the efficiency gain from covariate adjustment.

Studentized Statistics

par(mfrow = c(1, 2), mar = c(3, 3, 2, 1))

t_dm <- (dm_results[, "ate"] - tau_exp) / dm_results[, "se"]
hist(t_dm, breaks = 30, freq = FALSE, col = "lightblue",
     main = "DM", xlab = "", xlim = c(-4, 4))
curve(dnorm(x), add = TRUE, col = "red", lwd = 2)
abline(v = 0, lty = 2)

t_anc <- (ancova_results[, "ate"] - tau_exp) / ancova_results[, "se"]
hist(t_anc, breaks = 30, freq = FALSE, col = "lightblue",
     main = "ANCOVA", xlab = "", xlim = c(-4, 4))
curve(dnorm(x), add = TRUE, col = "red", lwd = 2)
abline(v = 0, lty = 2)

Simulation 4: Interference (Clustered AIPW)

This simulation validates the clustered AIPW estimator under network interference. We use fewer replications since estimation is more expensive.

n_clusters <- 200
group_struct <- c(2, 3)
tau_int <- 1
gamma_int <- c(0.5, 0.1)
n_reps_int <- 100

# We focus on group 0 and track beta_g at a few g values
# g = (0,0): pure direct effect = tau + 0*0.5 + 0*0.1 = 1.0
# g = (1,0): tau + 1*0.5 + 0*0.1 = 1.5
# g = (0,1): tau + 0*0.5 + 1*0.1 = 1.1
true_beta <- function(g) tau_int + sum(g * gamma_int)

# Track results for group 0 (first group)
g_indices <- list(c(0, 0), c(1, 0), c(0, 1))
g_encoded <- sapply(g_indices, function(g) {
  g[1] + g[2] * (group_struct[1] + 1) + 1  # 1-indexed
})
true_values <- sapply(g_indices, true_beta)

int_results <- array(NA, dim = c(n_reps_int, length(g_indices), 2),
                      dimnames = list(NULL, NULL, c("beta", "se")))

for (i in seq_len(n_reps_int)) {
  cdata <- generate_fixed_cluster(
    clusters = n_clusters, group_struct = group_struct,
    tau = tau_int, gamma = gamma_int
  )
  cl_obj <- clustered(
    Y = cdata$Y, Z = cdata$Z, X = cdata$X,
    cluster_labels = cdata$cluster_labels,
    group_labels = cdata$group_labels,
    ingroup_labels = cdata$ingroup_labels,
    n_matches = 50
  )
  res <- est_via_aipw(cl_obj)

  # Extract group 0 results
  for (j in seq_along(g_indices)) {
    idx <- g_encoded[j]
    int_results[i, j, "beta"] <- res[[1]]$beta_g[idx]
    int_results[i, j, "se"] <- res[[1]]$se[idx]
  }
}
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases
#> Warning in predict.lm(model, newdata = df): prediction from rank-deficient fit;
#> attr(*, "non-estim") has doubtful cases

Summary Table

int_rows <- list()
for (j in seq_along(g_indices)) {
  betas <- int_results[, j, "beta"]
  ses <- int_results[, j, "se"]

  # Remove NaN entries
  valid <- !is.nan(betas) & !is.nan(ses) & ses > 0
  betas <- betas[valid]
  ses <- ses[valid]
  tv <- true_values[j]

  ci_lo <- betas - 1.96 * ses
  ci_hi <- betas + 1.96 * ses

  int_rows[[j]] <- data.frame(
    g = paste0("(", paste(g_indices[[j]], collapse = ","), ")"),
    True_beta = tv,
    Mean_est = mean(betas),
    Bias = mean(betas) - tv,
    Emp_SE = sd(betas),
    Mean_SE = mean(ses),
    Coverage = mean(ci_lo <= tv & tv <= ci_hi),
    N_valid = length(betas)
  )
}
tbl_int <- do.call(rbind, int_rows)
knitr::kable(tbl_int, digits = 3, row.names = FALSE,
             caption = sprintf("Clustered AIPW (group 0, %d clusters, %d reps)",
                               n_clusters, n_reps_int))
Clustered AIPW (group 0, 200 clusters, 100 reps)
g True_beta Mean_est Bias Emp_SE Mean_SE Coverage N_valid
(0,0) 1.0 1.212 0.212 2.752 0.448 0.80 100
(1,0) 1.5 1.430 -0.070 0.491 0.429 0.90 100
(0,1) 1.1 1.134 0.034 0.230 0.243 0.94 100

Studentized Statistics

par(mfrow = c(1, length(g_indices)), mar = c(3, 3, 2, 1))
for (j in seq_along(g_indices)) {
  betas <- int_results[, j, "beta"]
  ses <- int_results[, j, "se"]
  valid <- !is.nan(betas) & !is.nan(ses) & ses > 0
  t_stats <- (betas[valid] - true_values[j]) / ses[valid]
  g_label <- paste0("(", paste(g_indices[[j]], collapse = ","), ")")
  hist(t_stats, breaks = 20, freq = FALSE, col = "lightblue",
       main = sprintf("g = %s", g_label),
       xlab = "", xlim = c(-4, 4))
  curve(dnorm(x), add = TRUE, col = "red", lwd = 2)
  abline(v = 0, lty = 2)
}

Summary

The simulations confirm:

  1. Observational estimators (OLS, IPW, AIPW, matching) are approximately unbiased with well-calibrated standard errors and ~95% coverage at moderate N.
  2. DML correctly estimates the continuous treatment effect with proper coverage.
  3. ANCOVA achieves meaningful efficiency gains over DM when covariates predict the outcome, while maintaining correct coverage.
  4. Clustered AIPW recovers the direct and spillover effects under interference, with coverage approaching nominal levels as cluster counts increase.