混杂因素的控制？控制变量足够吗？PS(倾向性得分)方法的应用场景？

MingzhaoHuang

最近在看关于倾向性得分相关的文献，由此引申了对混杂因素以及对混杂控制的一些思考
现在很多文章都会使用PS相关方法，来声称 “更好的控制混杂因素”
但是相对于传统的控制变量法，PS是否会更优秀吗？
不知道是否有现成的资料或者参考文献

关于上述这点，我也做了些模拟
正文：

我模拟了一个这样的框架（实在不好意思，不知道如何上传图片）：
A是我们研究的某种干预（二分类）对结局B的影响，其中U1,U2是两个混杂因素，且U1对U2有影响。

rm(list = ls()); gc()

set.seed(123)

n <- 100000

U1 <- rnorm(n, mean = 0, sd = 1)

U2 <- 0.7 * U1 + rnorm(n, mean = 0, sd = 1)

A <- 0.5 * U1 + 0.6 * U2 + rnorm(n, mean = 0, sd = 1) 

A <- ifelse(A >0, 1, 0)

B <- 0.7 * A + 0.3 * U1 + 0.8 * U2 + rnorm(n, mean = 0, sd = 1)

data <- data.frame(U1, U2, A, B)

首先，我们单独分析了干预A对结局B的影响

f0 <- lm(B ~ as.factor(A)) summary(f0)

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -0.677771   0.006181  -109.6   <2e-16 ***
as.factor(A)1  2.063214   0.008731   236.3   <2e-16 ***

A对B影响的回归系数远大于我们假定的0.7，因为有其他的混杂因素的效应被错误的归因到A的效应上。

第一种，最常见的控制混杂因素的方法-多元回归，将U1,U2作为控制变量纳入回归。

lm <- lm(data = data, B ~ as.factor(A) + U1 + U2) summary(lm)

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.0004866  0.0050435   0.096    0.923    
as.factor(A)1 0.7002599  0.0078488  89.218   <2e-16 ***
U1            0.3081652  0.0039986  77.068   <2e-16 ***
U2            0.7903921  0.0033975 232.639   <2e-16 ***

干预A的系数为0.7，与我们模拟数据的真实情况一致，混杂因素的效应被控制住。

第二种，使用倾向性得分加权

weight <- glm(data = data, as.factor(A) ~ U1 + U2, family = binomial(link = "logit"))

ps <- weight[["fitted.values"]]
ps <- unname(ps)
data <- cbind(data, ps)

library(dplyr)
data <- data %>%
    mutate( weight = ifelse(A == 1, 1/ps, 1/(1-ps)))

library(survey)
design <- svydesign(ids = ~1, 
                      weights = ~weight, 
                      data = data)

f1 <- svyglm(B ~ as.factor(A),
               design = design, 
               family = gaussian(link = "identity"))
summary(f1)

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -0.06745    0.01184  -5.696 1.23e-08 ***
as.factor(A)1  0.84162    0.01731  48.607  < 2e-16 ***

A的回归系数为 0.84，与真实值 0.7有一定偏差

第三种，使用双重稳健的倾向性得分加权

f2 <- svyglm(B ~ as.factor(A) + U1 + U2,
               design = design, 
               family = gaussian(link = "identity"))
  summary(f2)

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.004957   0.006603   0.751    0.453    
as.factor(A)1 0.695921   0.009639  72.202   <2e-16 ***
U1            0.310406   0.005776  53.737   <2e-16 ***
U2            0.792583   0.005413 146.427   <2e-16 ***

A的回归系数为 0.69，接近0.7

那么：
①在第一个回归中，U1,U2,A之间存在复杂的相关关系，破坏了回归分析自变量相互独立的假设，为什么不会对估计结果产生影响？为什么不需要纳入U1 U2 的交互项？
②既然将混杂因素纳入回归就可以进行较好的控制，为什么还要作倾向性得分匹配/加权？他们的应用场景到底是什么样的？

还望各位老师不吝赐教！

pyrty-Chen

你写的也太长了，太长的问题就很难引起读者兴趣把。。。不过确实是好问题。

第一个问题，在经典的低维的线性回归模型 $Y = X^{\top} \beta + \epsilon$ ，为了准确的估计真实的 $\beta$ ，其实不需要自变量X分量之间相互独立，相关也是没问题的。低维情况是样本数n远远大于自变量个数p，这里你的 n=100000，X对应的是（1，A，U1，U2)，所以p=4，你的样本量就是非常非常大，大数定理给出的结果就是平均值几乎就是常数了。那确实直接多元回归就可以非常准确的估计了，不管A，U1，U2相关性怎么样。

这里的证明也很简单， $\beta$ 的估计公式是
$\hat{\beta} = \left( \frac{1}{n}\sum_{i=1}^{n} X_i X_i^{\top} \right)^{-1} \frac{1}{n}\sum_{j=1}^{n} X_j Y_j \xrightarrow{p} \left( \Omega^{-1} \Omega \beta \right) =\beta$
其中， $\Omega = \mathbb{E}[X X^{\top}]$ 是X的协方差矩阵。收敛关系 $\xrightarrow{p}$ 表示概率收敛（n充分大的时候，与极限几乎没有差距，比如你这里的100000就是超级充分大了）。分别用大数定理和Slutsky定理就可以了。

MingzhaoHuang

pyrty-Chen
非常感谢您的回复！讲的很清楚！相比而言，我简直是个门外汉，还需要加强学习。
PS：我觉得您说的很对，我准备把我的问题精简一下，确实有点啰嗦了！

pyrty-Chen

第二个问题我再研究一下，但是除非你对A关于U1，U2的关系的估计非常非常自信，那么你应该无脑相信双稳健估计！现在双稳健还可以结合机器学习，深度学习提供非常好的估计。Double robust /AIPW/Double machine learning 都是双稳健的意思，可以看看相关的。简单地说，双稳健估计的偏差是你的 ”倾向性得分估计偏差“ 乘以 ”Y对于A和U的模型的估计偏差“，所以只要你有一个模型偏差接近0，双稳健估计的偏差就接近0，而倾向性得分匹配/加权的偏差就是倾向性得分估计偏差乘一些常数，所以如果你的倾向性得分估计的非常准，那么两者都ok，但是估计的不准，那倾向性得分匹配/加权就很糟糕了。AIPW的A就是Augmented，加强了的IPW（倾向性得分加权）。

pyrty-Chen

对了，这里你的倾向性得分加权的方法偏差比较大就是因为，你估计倾向性得分的时候用的模型是logistic model, 和你真实的模型不一样，所以倾向性得分加权方法在这里的确是有偏差的。

MingzhaoHuang

pyrty-Chen
用一些机器学习算法计算PS是不是会更合适？
另外，这个地方使用 logistic model 为什么会导致这种偏差，老师能展开讲讲吗。