匹配法
偏误
我们先突然谈起非参数估计,正如直觉那样,非参数估计就是不对参数做任何假定或者假设。
我们看会估计ATE的例子
ATE=E[Y1−Y0]
如果我们现在存在一个分组X,其会影响到Y,在估计ATE的时候,考虑到X混杂因素的存在,我们会求加权平均数估计ATE。这时,我们用小学的知识就知道,这样估计的ATE会更靠近分组内单位数更多的一组。
后面我们知道,我们可以通过线性回归来控制X,得到处理的系数值,即是我们想要的ATE。但与上面的加权平均数不一样。回归并不是根据样本数量加权的。相反,回归使用的权重与各组中处理变量的方差成正比。而对于连续变量,ATE会往协变量方差更大的方向偏移。
子类估计量 (Subclassification Estimator)
这个方法使用很少,或者基本不使用,在我的理解里就是对不同组进行加权平均数,这个缺陷也很显然,在1D时还好,但随着D的上升,特征空间的稀疏会使得子类估计量的实现出现偏误
如果我们想估计某种因果效应,但由于某些变量 X 的混杂而难以实现,我们需要做的是在 X 相同的小群体内进行处理组与对照组的比较。若具备条件独立性假设Y0,Y1⊥T∣X,则ATE为:
ATE=∫(E[Y∣X,T=1]−E[Y∣X,T=0])dP(x)
此时,会遍历特征X分布的所有空间,计算在Δx的均值差异
对于分类变量,特征X分布的空间里有K个单元,X1,...Xk,此时在离散情形下:
ATE^=k=1∑K(Yˉk1−Yˉk0)∗NNk
匹配估计量(Matching Estimator)
匹配这个概念我们常常听见,即那一个单位去在另一个群中找出一个或几个的类似单元,使得具有可比性,即为每个处理单元匹配相似的对照单元实现可比性。这么做的原因和子类估计量或者线性回归类似,都是为了控制潜在的混杂因素。
设处理组有X个单位,而我们需要知道处理的效应,在控制组中有任意个单位(因为它可以被重复对照),那么最后我们匹配后的处理组应该有≤X个单位。此时我们便可以得到处理组的平均处理效应ATET。
如何匹配
在真实世界里,我们很自然的无法完美匹配,偶尔有几个的正常,如果有很多,那么可能得注意是不是有人偷懒了。
我们常用距离来衡量单元间的相似程度。在日常用到的最多的是欧几里得范数(L2范数) ∣∣Xi−Xj∣∣2。除此之外,还有L1范数,即曼哈顿距离∣∣X∣∣1=∑xi。还有余弦距离等。在使用不同特征尺度的距离时,我们常常会进行缩放处理,如标准化,使处于一个量纲。
在定义了距离之后,我们现在可以将匹配定义为与待匹配样本最接近的邻居。用数学术语表达,匹配估计量可以写作如下形式:
ATE^=N1i=1∑N(2Ti−1)(Yi−Yjm(i))
其中Yjm(i)代表来自另一处理组中与Yi最为相似的样本。我们进行(2Ti−1)双向匹配:既将处理组与对照组匹配,也将对照组与处理组匹配。
在现实应用中,我们会用到大名鼎鼎的,每个机器学习入门级的算法:k-NN,即最近邻算法。
匹配偏误
我们受限考虑处理组的平均处理效应(ATET):
ATET^=N11∑(Yi−Yj(i))
其中N1代表接受处理的个体数量,Yj(i)为处理单元 i 的未处理匹配项。
为了检验偏误,我们想:
N1(ATET^−ATET)→N(0,σ2)
此时我们定义X下的处理组的平均结果为μ0(x)=E[Y∣X=x,T=0],则有
E[N1(ATET^−ATET)]=E[N1(μ0(Xi)−μ0(Xj(i)))]
证明:
∵solve ATET∴Y(0)⊥T∣X∴E[Y(0)∣X=x,T=1]=E[Y(0)∣X=x,T=0]=μ0(x)⇒E[Y(0)∣X=x,T=1]=μ0(x)⇒ATET^−ATET=[Yi(1)−Yj(i)(0)]−[Y(1)−Y(0)]=Y(0)−Yj(i)(0)=μ0(Xi)−μ0(Xj(i))
对于μ0(Xi),即Y0,是处理单元的反事实结果,而μ0(Xj(i))则是匹配的未处理单元的实际结果。而显然我们前面说过,匹配是为了找到最近邻的,而不是完全一样的,所以自然这一差值可以不为0,而我们在匹配背景下,我们通常做的是:
∵Xi≈Xj∴Y0i≈Y0j
当随着样本量N1的增加,可匹配的单元数量增多,即单元i与其匹配单元j之间的差异也会缩小。但显然,收敛成0速度慢,也所以E[N1(μ0(Xi)−μ0(Xj(i)))]可能不会收敛至0,因为N1是不收敛的,所以取决于非收敛函数的发散特征与收敛函数的趋近速度。现在我们知道偏误的具体大小,其来自于匹配的错值想象一下一个东西与另外一个东西不能重合,他们相差的。我们还知道每项观测数据对偏误的贡献为(μ0(Xi)−μ0(Xj(i))),所以在每次匹配中,我们都减去这一数值即可,我们可以通过线性回归的方式获得这个偏误值μ0^(Xj(i))来代替μ0(Xj(i))。
ATET^=N11∑((Yi−Yj(i))−(μ0^(Xi)−μ0^(Xj(i))))
Umm,这看起来复杂了很多,因为我们这次用线性回归不是为了直接获得系数,而是为了校正偏误。按照书上说的,这部分线性回归具有局部性,即它并不试图考察若处理组看起来像未处理组时会如何我觉得这里翻译过来就是不考虑反事实推理,因为此时的样本已经变了,我们通过匹配少了一些未处理组的单元。也就是说,这里的线性回归不进行外推,匹配工作完成了外推这一部分,且其是核心。
当然,校正偏误这里我们所用的局部性线性拟合是可以的,但我想的是,如果偏误的分布是一个曲线型的,又是否可以采用多项式进行拟合?
维度灾难
基于少数特征(如性别)匹配个体很容易,但当我们增加更多特征时,寻找匹配对象就变得越来越困难。更概括地说,特征维度越高,个体与其匹配对象之间的距离就越大。也就意味着填满特征空间所需的数据点数量会随着特征或维度的增加呈指数级增长。例如,若填满 3 个特征空间需要 X 个数据点,那么填满 4 个特征的空间则需要指数级更多的点。
在子类估计器里,当样本少时,其特征空间也自然稀疏,部分单元格将仅有极少数据,甚至可能仅包含处理组或对照组,导致无法在该处估计平均处理效应(ATE)。而在匹配估计器里也是一样,单元间距离的增加到导致匹配对之间的距离增加并引发偏误问题。这时候不妨考虑一下线性回归,因为其通过投影这一巧妙的思想降成了1D空间。
倾向得分
直觉建立
倾向得分的想法是,我们不通过直接控制混杂变量X,即也可以实现条件独立性假设(Y0,Y1)⊥T∣X。
画在DAG上,我们可以更加清晰知道这一想法
那么,这里的e(x)即倾向得分。我们还可以知道,这个得分就是P(T∣X)。
因此,倾向性得分可以让我们不必要对X整体进行条件限制就能够实现潜在结果对处理的独立性,此时的条件独立性假设为(Y0,Y1)⊥T∣e(X)
若我们已知e(x),那么由于阻断效应X⊥T∣e(x),当我们控制e(x)的时候,就等同于控制了X。
对于e(x),我们所能理解的是有多大概率接受处理。想想,当我们从处理组和对照组中各取一个人,且这两个人具有相同的处理接受概率时,他们便具有可比性。这时候便是,两个人都有一样的接受处理的概率,一个人接受,一个人不接受,我们就有理由相信这个不接受处理的人是偶然的,换句话说,这样看起来像是随机分布的一样。
那么,我们现在来证明这一点。严格地说,这个等价关系需要拆成两步:右边这个式子直接证明的是倾向得分的平衡性质,即
X⊥T∣e(X)⟺E[T∣e(X),X]=E[T∣e(X)]
我们已知条件独立性假设(Y0,Y1)⊥T∣X成立,可以得到
(Y0,Y1)⊥T∣e(X)
e(X)是X的一个函数,并且在二元处理变量下有
e(X)=P(T=1∣X)=E[T∣X]
在给定X时,再额外给定e(X)并不会带来新的信息:
E[T∣e(X),X]=E[T∣X]=e(X)
另一方面,根据迭代期望公式:
E[T∣e(X)]=E[E[T∣X]∣e(X)]=E[e(X)∣e(X)]=e(X)
所以可以得到
E[T∣e(X),X]=E[T∣e(X)]
由于T是二元变量,E[T∣⋅]就等于P(T=1∣⋅)。于是上式等价于
P(T=1∣e(X),X)=P(T=1∣e(X))
同时也有
P(T=0∣e(X),X)=P(T=0∣e(X))
这说明在给定倾向得分e(X)以后,T的条件分布不再依赖完整的X,因此
X⊥T∣e(X)
反过来,如果已经知道X⊥T∣e(X),那么T在给定e(X)后的条件分布自然不依赖X,所以也有
E[T∣e(X),X]=E[T∣e(X)]
即
X⊥T∣e(X)⟺E[T∣e(X),X]=E[T∣e(X)]
接下来证明潜在结果是条件独立性。记Y∗=(Y0,Y1),对∀集合A,有
P(Y∗∈A∣T=t,e(X)=e)=∫P(Y∗∈A∣T=t,X=x,e(X)=e)dP(x∣T=t,e(X)=e)=∫P(Y∗∈A∣T=t,X=x)dP(x∣T=t,e(X)=e)=∫P(Y∗∈A∣X=x)dP(x∣T=t,e(X)=e)=∫P(Y∗∈A∣X=x)dP(x∣e(X)=e)=P(Y∗∈A∣e(X)=e)
其中,第三行是因为原始假设(Y0,Y1)⊥T∣X,第四行用到了平衡性质X⊥T∣e(X)。所以最终得到
(Y0,Y1)⊥T∣e(X)
也就是说,倾向得分本身并不是和潜在结果独立性“天然等价”,而是先通过e(X)实现X与T的平衡;在原本控制X即可消除混杂的前提下,控制e(X)才可以替代控制完整的X。
倾向性评分逆概率加权(IPTW)
假设现在我们已经获得了e(x),当我们用倾向性评分加权时,我们获得的ATE为:
E[Y∣X,T=1]−E[Y∣X,T=0]=E[e(x)Y∣X,T=1]P(T)−E[(1−e(x))Y∣X,T=0](1−P(T))
对于接受处理的人,是将所有接受处理的个体按其接受处理的概率的倒数进行加权。换句话说,如果一个人原本不接受处理的概率很高,但它被接受处理了,那么这个人会被赋予更高的权重。这很简单,因为这个人代表了很多接受处理概率不高的人。通过这种方式,我们构造出一个“每个人都接受了处理”的人群,且这个人群的分布和原始人群相同
另外上面的等式可以被进一步简化:
E[Y∣X,T=1]−E[Y∣X,T=0]=E[e(x)Y∣X,T=1]P(T)−E[(1−e(x))Y∣X,T=0](1−P(T))=E[e(x)YTX]−E[(1−e(x))Y(1−T)X]=E[e(x)YT−(1−e(x))Y(1−T)X]=E[Ye(x)(1−e(x))T(1−e(x))−e(x)(1−T)X]=E[Ye(x)(1−e(x))T−e(x)X]
这时候我们对上述等式进行积分,即有:
∫E[Y∣X,T=1]−E[Y∣X,T=0]dX=∫E[Ye(x)(1−e(x))T−e(x)X]dX=E[Ye(x)(1−e(x))T−e(x)]
显然,由于处理概率不可能在(0,1)的范围之外,所以每个人都会有概率接受处理和不接受处理。另外一种表述即是处理组与非处理组需要有重叠区域,即正值性假设。
因果推断中有三大假设:
- 一致性假设:若接受干预T=1,则观测到的结果等于干预状态下的潜在结果Y(1)。
- 可忽略性假设:当我们控制X后,(Y(0),Y(1))⊥T∣X
- 稳定单元干预假设(SUTVA):
- 无干扰:一个个体的处理状态不影响其他个体
- 一致干预:每个个体处理是单一且明确的,不因其他干预而改变
估计e(x)
由于现实情况里处理分配机制未知,我们一般而言,才有e(x)^来近似代替e(x)。常用Logistic回归,但也可以采用梯度提升等机器学习方法。
梯度提升:指的是累加多个弱学习器,比如决策树,来逼近目标函数,进而优化对象(整个模型的函数形式),在每次迭代中都会添加一个新函数。
梯度下降:常用于最小化损失函数,主要用于优化对象,比如权重和偏置,每次更新迭代会修改原有的参数。
当损失函数的SSE=1且步长step=1时,梯度下降更新参数与梯度提升的更新函数等价。
我们会使用scikit-learn库中的一个核心模块进行实现:
from sklearn.linear_model import LogisticRegression
T = 'intervention'
Y = 'outcome'
X = data_with_categ.columns.drop([T,Y])
ps_model = LogisticRegression(C=1e6).fit(data_with_categ[X], data_with_categ[T])
data_ps = data.assign(propensity_score=ps_model.predict_proba(data_with_categ[X])[:, 1])
这样子我们就可以后续直接使用IPTW分析处理与结局的ATE了。
在使用方法之前,我们需要检查处理组与未处理组之间是否存在重叠,观察未处理组和处理组的倾向得分经验分布。
标准误差
如果我们得到了真实的倾向得分,换句话说,我们得到的倾向得分是无偏的,我们可以精确算出来:
σw2=∑i=1nwi∑i=1nwi(yi−μ^)2
但显然,我们很难做到,这时候使用Bootstrap自助法抽样是一个好的选择,来获得ATE的分布。
常见问题
- 倾向得分无需非常精确地预测处理分配,它只需涵盖所有混杂变量。换句话说,我们首先考虑哪些变量与我们的结果有关,我们只纳入哪些就可以了,如果我们还额外纳入预测处理分配的变量,会增加倾向得分估计量的方差。
- 某些情况下,接受处理者的处理概率远高于未处理者,且倾向得分分布的重叠区域非常有限。这时候很有可能会面临没有可对比的情况,这时候我们就不得不进行外推,这时候就会引发偏误。
- 极高或极低倾向得分的实体权重过大,这会增加方差,根据经验法则,当任何权重超过 20 时(例如未处理个体倾向得分达 0.95 或处理个体倾向得分低至 0.05 时),就会面临问题。
倾向得分匹配(PSM)
如前所述,在拥有倾向得分的情况下,我们不再直接控制变量 X,仅控制倾向得分即可。因此,可将倾向得分视为对特征空间的一种降维操作,换句话来说,我们把一个高维的特征空间都压缩到了一个1D变量上,即e(x)。所以,我们可将倾向得分作为其他模型(例如回归模型)的输入特征。
psm_model = smf.ols("achievement_score ~ intervention + propensity_score", data=data_ps)
psm_model.fit()
此外,我们还可基于倾向得分进行匹配。此次,无需寻找所有 X 特征均相似的匹配对象,仅需确保其倾向得分相同即可。
cm = CausalModel(
Y=data["outcomes"].values,
D=data["intervention"].values,
X=data[["propensity_score"]].values
)
cm.est_via_matching(matches=1, bias_adj=True)
print(cm.estimates)
置信区间的求解不能用Bootstrap实现,因为其不适用于匹配的场景。
但在网上搜索时,网上给出的答案是可以使用 Bootstrap 来估计和修正标准误。说明这一场景仍有争议
双重稳健估计
在前面我们知道了,我们可以使用线性回归和倾向得分来估计ATE,但哪个更加准确呢,因为他们俩的优缺点好像刚好互补,于是乎,我们不如两个一起用,这时候便提出了双重稳健估计:
ATE^=N1∑(e^(Xi)Ti(Yi−μ1^(Xi))+μ1^(Xi))−N1∑(1−e^(Xi)(1−Ti)(Yi−μ0^(Xi))+μ0^(Xi))
其中,μ1^(x)是E[Y∣X,T=1]的估计,μ0^(x)是E[Y∣X,T=0]的估计,通过线性回归系数得到。
正如直觉一样,这个等式有两部分组成,第一部分是E[Y1],第二部分是E[Y0]。
我们只看第一部分,因为显然第二部分与第一部分是对称的:
E^[Y1]=N1∑(e^(Xi)Ti(Yi−μ1^(Xi))+μ1^(Xi))
- 假设线性回归是正确的,即μ1^(Xi)是正确的(无偏),那么E[Ti(Yi−μ1^(Xi))]=0,所以μ1^在处理组上的残差均值为零,这最后的结果就变成了E^[Y1]=μ1^(Xi)。
- 接下来我们来考虑倾向性得分是对的,在此之前,我们先变换一下,让结果可以一目了然:
E^[Y1]=N1∑(e^(Xi)Ti(Yi−μ1^(Xi))+μ1^(Xi))=N1∑(e^(Xi)TiYi−e^(Xi)Tiμ1^(Xi)+μ1^(Xi))=N1∑(e^(Xi)TiYi−(e^(Xi)Ti−1)μ1^(Xi))=N1∑(e^(Xi)TiYi−(e^(Xi)Ti−e^(Xi))μ1^(Xi))
如果e^(Xi)是正确的,即它是无偏的,那么E[Ti−P^(Xi)]=0,从而消掉了依赖于μ1^(Xi)的一部分,双重稳健估计量简化为倾向得分加权估计量P^(Xi)TiYi