在 R 中复制 Gary King 的仿真研究答案

【问题标题】：Replicating Simulation Study of Gary King in R在 R 中复制 Gary King 的仿真研究
【发布时间】：2021-12-28 14:08:47
【问题描述】：

在this的文章中，作者描述了以下simulation study的Propensity Score Matching，

对于两个协变量中的每一个，我们随机且独立地从 Uniform(0,5) 中抽取 100 个控制单元，从 Uniform(2,6) 中抽取 100 个处理单元。这使得重叠的 [1, 5] x [1, 5] 正方形成为完全随机的实验，并且观察结果超出了增加的不平衡性。我们将结果生成为 $Y = 2T_i + X_{i1} + X_{i2} + \varepsilon$，其中 $ \varepsilon \sim N(0,1)$

但是，这里没有明确说明T 是如何定义的。 $X_{i,k} ~ \forall ~k$ 就论文而言，是一个长度为 100 的向量，而 T 也是一个（自然）向量。

我对此的第一个解释是，

# First Attempt

dgp <- function(n = 100) 
{
  
  
  # Generate Covariates
  x1 <- runif(n, min = 0, max = 5)
  x2 <- runif(n, min = 1, max = 6)
  
  # Generate data with treatment indicators
  # and outcomes
  
  tibble(
    x1 = x1,
    x2 = x2
  ) %>% mutate(
    treatment = as.numeric(x1 >1 | x2 > 1),
    y         = 2*treatment + x1 + x2 + rnorm(n, 0, 1)
  )
}

但是，此尝试将未处理的 runif(n, 0,5) 中的一些 obs 分配给处理组，这不是作者的意图。

我的第二次尝试包括每个协变量的cbind(x1 = runif(100, 0,5), treatment = 0) - 但是，treatment-变量的长度超过了协变量的长度，因此会崩溃。

我如何解释论文中提到的模拟研究，并在 R 中复制它？

编辑：在@whuber 的评论之后完全引用文章。

【问题讨论】：

一个很好的例子，说明为什么期刊应该为所有已发表的研究成果强制开放数据和源代码
我绝对同意这一点。
从This leaves the overlapping [1,5]×[1,5] square as a completely randomized experiment，我猜它对应2T，那么T应该是Uniform( 1,5 )
2. Thus, for unit i (i = 1, . . . , n), denote the treatment variable as T i∈ {0, 1}, where 0 refers to the “control group” and 1 the “treated group”.这也许对你理解也有用。
如果我对论文的理解正确，$T_i$ 是对照组的 0 和治疗组的 1——它只是一个指标。因此，从 U(0, 5) 得出的前 100 个观察值中的每一个都应该具有 T = 0，而从 U(1, 6) 得出的第二组观察值中的每一个都应该具有 T = 1。

标签： r simulation matching

【解决方案1】：

我尝试使用文本中描述的数据生成过程以两种方式复制这项研究，我希望演示的是一个等效的数据生成过程。我的复制是here。

这是我用来生成数据的代码的 sn-p，其中A 是处理：

n <- 200
A <- rep(0:1, each = n/2)
X1 <- ifelse(A==1, runif(n, 1, 6), runif(n, 0, 5))
X2 <- ifelse(A==1, runif(n, 1, 6), runif(n, 0, 5))

Y <- 2*A + X1 + X2 + rnorm(n)

请注意，在论文中，他们首先生成治疗组，然后将协变量值分配给分配的单位。也就是说，治疗不依赖于协变量。相反，治疗是随机的，协变量实际上是中介。

一位同事认为这破坏了他们的整篇论文，我不同意，这就是为什么我着手制定一个等效的数据生成过程，其中治疗是基于协变量创建的，但治疗和协变量具有相同的联合分布在原始配方中。这是执行此操作的代码：

G1 = sample(LETTERS[1:5], n, replace = TRUE, 
            prob = c(5/50, 4/50, 32/50, 4/50, 5/50))
X1 = runif(n, 
           sapply(G1, switch, A = 1, B = 0, C = 1, D = 5, E = 0),
           sapply(G1, switch, A = 6, B = 1, C = 5, D = 6, E = 5))
X2 = runif(n, 
           sapply(G1, switch, A = 5, B = 1, C = 1, D = 1, E = 0),
           sapply(G1, switch, A = 6, B = 5, C = 5, D = 5, E = 1))

#P is propensity score.
P = 1*(X1 > 1)*(X2 > 1) - .5*(X1 > 1)*(X1 < 5)*(X2 > 1)*(X2 < 5)

A <- rep(0, n)
A[P == 1] <- 1
A[P == .5] <- sample(c(rep(1, n/2 - sum(P==1)), 
                       rep(0, n/2 - sum(P==0))), 
                     sum(P == .5), replace = FALSE)

基本上，您创建一个外生辅助分组变量G1，然后根据G1 生成X1 和X2。倾向得分以不同寻常的方式定义；值可以是 0、.5 或 1。通常，对于倾向得分不是 0 或 1 的人，将使用 rbinom() 生成处理，但为了确保处理单元和控制单元的数量相等，我们使用 @987654330 @ 反而。这可确保倾向得分准确且组大小平衡。倾向得分为 0.5 的人对应于“完全随机实验”，倾向得分为 0 或 1 的人对应于那些“增加不平衡”。

无论采用何种方法，都可以验证变量的联合分布是否相同，表明只要变量的联合分布相同，无论使用哪种数据生成过程，论文的消息都适用.这是因为该论文提出了统计论点，而不是因果论点，并且统计数据与人们通常会调用的因果假设无关（即，关于治疗是否导致协变量或协变量导致治疗）。

这与统计估计有关，用于匹配的是 $E[E[Y|X,A=1]] - E[E[Y|X,A =0]]$。这是否对应于因果估计是另一个问题。在混杂的数据生成过程（我描述的第二个）下，因果估计是去混杂的平均治疗效果 $E[Y^1] - E[Y^0] $。在中介数据生成过程（论文中的那个）下，因果估计是治疗的直接影响，而不是通过协变量。作者没有澄清这一点，因为这与他们的论点相切，但对于那些深入调查的人来说，了解实际情况是值得的。

【讨论】：