对 MatchIt/vignettes/estimating-effects 中的代码案例感到困惑

Question

关于文档中给出的例子： https://cran.r-project.org/web/packages/MatchIt/vignettes/estimating-effects.html 代码：

gen_X <- function(n) {
  X <- matrix(rnorm(9 * n), nrow = n, ncol = 9)
  X[,5] <- as.numeric(X[,5] < .5)
  X
}

gen_A <- function(X) {
  LP_A <- - 1.2 + log(2)*X[,1] - log(1.5)*X[,2] + log(2)*X[,4] - log(2.4)*X[,5] + log(2)*X[,7] - log(1.5)*X[,8]
  P_A <- plogis(LP_A)
  rbinom(nrow(X), 1, P_A)
}


gen_Y_C <- function(A, X) {
  2*A + 2*X[,1] + 2*X[,2] + 2*X[,3] + 1*X[,4] + 2*X[,5] + 1*X[,6] + rnorm(length(A), 0, 5)
}

gen_Y_B <- function(A, X) {
  LP_B <- -2 + log(2.4)*A + log(2)*X[,1] + log(2)*X[,2] + log(2)*X[,3] + log(1.5)*X[,4] + log(2.4)*X[,5] + log(1.5)*X[,6]
  P_B <- plogis(LP_B)
  rbinom(length(A), 1, P_B)
}

gen_Y_S <- function(A, X) {
  LP_S <- -2 + log(2.4)*A + log(2)*X[,1] + log(2)*X[,2] + log(2)*X[,3] + log(1.5)*X[,4] + log(2.4)*X[,5] + log(1.5)*X[,6]
  sqrt(-log(runif(length(A)))*2e4*exp(-LP_S))
}


set.seed(19599)

n <- 2000
X <- gen_X(n)
A <- gen_A(X)

Y_C <- gen_Y_C(A, X)
Y_B <- gen_Y_B(A, X)
Y_S <- gen_Y_S(A, X)

d <- data.frame(A, X, Y_C, Y_B, Y_S)
#=================================================================
pair_ids <- levels(md$subclass)

est_fun <- function(pairs, i) {
  
#Compute number of times each pair is present
numreps <- table(pairs[i])
  
#For each pair p, copy corresponding md row indices numreps[p] times
ids <- unlist(lapply(pair_ids[pair_ids %in% names(numreps)],
                       function(p) rep(which(md$subclass == p), 
                                              numreps[p])))
 
  
#Subset md with block bootstrapped ids
md_boot <- md[ids,]
  
#Effect estimation
fit_boot <- lm(Y_C ~ A, data = md_boot, weights = weights)
  
#Return the coefficient on treatment
return(coef(fit_boot)["A"])
}

boot_est <- boot(pair_ids, est_fun, R = 499)
boot_est
boot.ci(boot_est, type = "bca")

est_fun有错吗?

运行 pair_ids <- levels(md$subclass) 仅生成 c(1,2,3,4,5,6,7……）（pair_ids 中的每个观察值的大小为 1），因为它只输出 level.

所以运行numreps <- table(pairs[i])，如果i=3，输出[1] "3"； numreps 是 而不是“每对出现的次数”。

运行numreps[p]（这里是numreps[3]），输出NULL，让rep(which(md$subclass == p), numreps[p])显得毫无意义（复制总是 1).

虽然这个案例代码输出了正确的效果预估结果，但也许我们设置的这个函数有一点点无效？

Answer 1

阅读 boot::boot() 的文档。 i不是一个数字，它是由运行

组成的长度为length(pair_ids)的向量

i <- sample(seq_along(pair_ids), length(pair_ids), replace = TRUE)

很多次。因为当 bootstrapping 您正在使用替换采样时，i 有可能（并且几乎可以保证）有重复的条目，在这种情况下 table(pairs[i]) 不会产生全 1。此函数的作用是通过替换对 ID 对进行采样，使用重采样对 ID 中存在的单位创建新数据集，然后估计每个新数据集中的处理效果。这就是集群 bootstrap 应该做的；此实现没有任何无效之处。

自己试试运行：

i <- sample(seq_along(pair_ids), length(pair_ids), replace = TRUE)
table(pair_ids[i])
table(table(pair_ids[i]))

您会看到许多对 ID 出现不止一次；第三行显示每个重复项出现了多少次。

对 MatchIt/vignettes/estimating-effects 中的代码案例感到困惑

Confusion about the code case in MatchIt/vignettes/estimating-effects

r

match

matching

propensity-score-matching