Propensity score简介

预测causal effect的标准做法是做随机试验,但是这种随机试验通常很难完成,比如要预测吸烟与肺癌的因果效应,那么就必须随机让一部分人强制抽烟(treatment=抽烟),让另一部分人强制不抽烟(treatment=不抽烟),然后这样做是不现实的。此外如果不进行随机试验,那么我们很容易得出一些莫名其妙的结论。比如,你怎么知道是吸烟导致肺癌,而不是因为喜欢吸烟的人本身就容易得肺癌呢?换句话说,这里可能存在一个协变量X(喜欢抽烟),导致了一个人选择treatment=抽烟。

那这个问题怎么解决呢?直觉上,随机试验的本质其实就是将所有的协变量随机打散,分在treatment=抽烟和treatment=不抽烟中,从而使得协变量与是否做干预是独立的treatment⊥Xtreatment \bot XtreatmentX。所以,类似的思路,我们只要找到一种分组方式,使得X与treatment独立就可以了。那么怎么找到这种分组?我们可以定义Propensity score:
e(x) =def Pr⁡(treatment=抽烟∣X=x). {\displaystyle e(x)\ {\stackrel {\mathrm {def} }{=}}\ \Pr(\text{treatment}=抽烟|X=x).} e(x) =def Pr(treatment=X=x).
现在如果有两个人,一个抽样,一个不抽烟,但是他们的Propensity score是一样的(或者相似的),我就将他们分别分到抽烟组和不抽烟组,只要找到足够多这样的pair,这两个分组就会跟协变量独立。更直观来说,两个人的Propensity score一样,意味着他们的X应该是很接近,如果我们能找到两个“一样”的人,他们的生活环境,年龄,性别等等因素(协变量)都是一样的,除了一个人抽烟一个人不抽烟这一点不同,那么,只要有足够多这样的人,我们就能判断吸烟对肺癌的影响。

但是以上分组的过程我们完全采用一种简单matching的技术解决,也就是计算X之间的相似性然后根据相似性进行分组即可,为什么还要使用propensity score呢?原因是,如果X是超级高维的变量,简单的matching就会失效,为了维度诅咒的原因,高维的matching是非常不可靠的,然而对于propensity score而言,不管X有多少维,这个score永远只会是1维

Propensity score theorem

根据上面的分析,利用propensity score的大小就可以帮助我们找到一种分组使得该分组与treatment的效果独立,正式地:
(Y(1),Y(0))⊥T∣X⇒(Y(1),Y(0))⊥T∣e(X) ( Y ( 1 ) , Y ( 0 ) ) \bot T | X \Rightarrow ( Y ( 1 ) , Y ( 0 ) ) \bot T | e ( X ) (Y(1),Y(0))TX(Y(1),Y(0))Te(X)
要理解这个等式,我们可能需要先理解一下,为什么需要

(Y(1),Y(0))⊥T ( Y ( 1 ) , Y ( 0 ) ) \bot T (Y(1),Y(0))T

这一条件也被称为Ignorability/Exchangebility Asummption. 从下图上可以直观理解这一条件。那就是treatment的选择是独立的。
在这里插入图片描述
咦,那为什么,T会跟Y独立,难道不是T导致Y的吗?我们用个例子解释这个问题,假设你在纠结要不要养猫,那么我们用Y(1)Y(1)Y(1)表示得到小猫的幸福程度,用Y(0)Y(0)Y(0)表示没有小猫的幸福程度,显然不管我们选择得到小猫,还是没有小猫,但对于结果Y(1),Y(0)Y(1),Y(0)Y(1),Y(0)它永远是“存在“的,只是你无法“观测”而已,你可以想象这两个结果是上帝决定的,我们选择treatment只是决定观测到其中的哪一个,并不会影响他具体的值。

为什么这一假设如此重要?因为我们只能观测到选择treatment后的Y,而无法知道上帝决定的Y(1),Y(0)Y(1),Y(0)Y(1),Y(0),如果T的选择会影响上帝的结果,那么我们的观测值将不再可信,这意味着存在某些因素会同时影响t和Y。基于该假设,我们就可以预测出average treatment effect。

E[Y(1)]−E[Y(0)]=E[Y(1)∣T=1]−E[Y(0)∣T=0]=E[Y∣T=1]−E[Y∣T=0] \left. \begin{array}{l}{ E [ Y ( 1 ) ] - E [ Y ( 0 ) ] = E [ Y ( 1 ) | T = 1 ] - E [ Y ( 0 ) | T = 0 ] }\\{ = E [ Y | T = 1 ] - E [ Y | T = 0 ] }\end{array} \right. E[Y(1)]E[Y(0)]=E[Y(1)T=1]E[Y(0)T=0]=E[YT=1]E[YT=0]

该公式表明,我们可以从观测数据中学习出average treatment effect,然而现实中,T往往与Y(0),Y(1)不独立的,所以我们往往会做一个弱一点的假设,也就是Conditional Exchangeability/Unconfoundendess Assumption:

(Y(1),Y(0))⊥T∣X ( Y ( 1 ) , Y ( 0 ) ) \bot T|X (Y(1),Y(0))TX

这意味着,当给定协变量X的时候,他们就是独立的,如下图所示

在这里插入图片描述
类似的,我们就可以推导出:

E[Y(1)−Y(0)∣X]=E[Y(1)∣X]−E[Y(0)∣X]=E[Y(1)∣T=1,X]−E[Y(0)∣T=0,X]=E[Y∣T=1,X]−E[Y∣T=0,X] \left. \begin{array}{l}{ E [ Y ( 1 ) - Y ( 0 ) | X ] = E [ Y ( 1 ) | X ] - E [ Y ( 0 ) | X ] }\\{ = E [ Y ( 1 ) | T = 1 , X ] - E [ Y ( 0 ) | T = 0 , X ] }\\{ = E [ Y | T = 1 , X ] - E [ Y | T = 0 , X ] }\end{array} \right. E[Y(1)Y(0)X]=E[Y(1)X]E[Y(0)X]=E[Y(1)T=1,X]E[Y(0)T=0,X]=E[YT=1,X]E[YT=0,X]

根据iterated expectations: EX(X)=EY[EX(X∣y)]E_X(X)=E_Y[E_X(X|y)]EX(X)=EY[EX(Xy)]

EX(X)=∫xp(x)dx=∫∫xp(x,y)dxdy=∫∫xp(x∣y)dxp(y)dy=∫EX[X∣y]p(y)dy=EY[EX[X∣y]] \begin{aligned} E_{X}( X) & =\int xp( x) dx\\ & =\int \int xp( x,y) dxdy\\ & =\int \int xp( x|y) dxp( y) dy\\ & =\int E_{X}[ X|y] p( y) dy\\ & =E_{Y}[ E_{X}[ X|y]] \end{aligned} EX(X)=xp(x)dx=xp(x,y)dxdy=xp(xy)dxp(y)dy=EX[Xy]p(y)dy=EY[EX[Xy]]

我们可以进一步得到:

E[Y(1)−Y(0)]=EX[E[Y(1)−Y(0)∣X]]=EX[E[Y∣T=1,X]−E[Y∣T=0,X]] \begin{aligned} E[Y(1)-Y(0)] & =E_{X}[ E[Y(1)-Y(0)|X]]\\ & =E_{X}[ E[Y|T=1,X]-E[Y|T=0,X]] \end{aligned} E[Y(1)Y(0)]=EX[E[Y(1)Y(0)X]]=EX[E[YT=1,X]E[YT=0,X]]

于是,我们就成功得到了average treatment effect. 接下来回到正题,对于propensity score的作用,

(Y(1),Y(0))⊥T∣X⇒(Y(1),Y(0))⊥T∣e(X) ( Y ( 1 ) , Y ( 0 ) ) \bot T | X \Rightarrow ( Y ( 1 ) , Y ( 0 ) ) \bot T | e ( X ) (Y(1),Y(0))TX(Y(1),Y(0))Te(X)

显然,e(X)e(X)e(X)起到了替代X的作用!也就是说,用一个一维的东西代替了一个高维的变量。那么从直觉上,这个东西为什么成立呢?

在这里插入图片描述
看上图,上面的W是我们的X,因为我们已知给定可以block掉confounder,因为e(W)是W的函数,可以认为e(W)也block掉了路径。

该正式证明过程如下:

首先因为T是二值变量,所以概率可以写成期望

P(T=1,∣Y(t),e(X))=E[T∣Y(t),e(X)] P(T=1,|Y(t),e(X))=E[T|Y(t),e(X)] P(T=1,Y(t),e(X))=E[TY(t),e(X)]

根据iterated expectations,可以额外增加一个X\displaystyle XX作为条件

E[T∣Y(t),e(X)]=EX[E[T∣Y(t),e(X),X]∣Y(t),e(X)] E[T|Y(t),e(X)]=E_{X}[ E[T|Y(t),e(X),X]|Y(t),e(X)] E[TY(t),e(X)]=EX[E[TY(t),e(X),X]Y(t),e(X)]

因为我们给定了X\displaystyle XX,所以e(X)\displaystyle e(X)e(X)只是一个X的函数,所以是可以去掉的:

=EX[E[T∣Y(t),X]∣Y(t),e(X)] =E_{X}[ E[T|Y(t),X]|Y(t),e(X)] =EX[E[TY(t),X]Y(t),e(X)]

根据(Y(1),Y(0))⊥T∣X\displaystyle (Y(1),Y(0))\bot T|X(Y(1),Y(0))TX

=EX[E[T∣X]∣Y(t),e(X)] =E_{X}[ E[T|X]|Y(t),e(X)] =EX[E[TX]Y(t),e(X)]

又因为T是二值的,所以期望可以写成概率

=EX[p(T=1∣X)∣Y(t),e(X)]=EX[e(X)∣Y(t),e(X)]=e(X) \begin{aligned} & =E_{X}[ p( T=1|X) |Y(t),e(X)]\\ & =E_{X}[ e( X) |Y(t),e(X)]\\ & =e( X) \end{aligned} =EX[p(T=1X)Y(t),e(X)]=EX[e(X)Y(t),e(X)]=e(X)

最后,因为条件集出现e(X)\displaystyle e(X)e(X),所以这个期望最终的就是e(X)\displaystyle e(X)e(X),从而与T\displaystyle TT是独立的。

Inverse Probability Weighting (IPW)

从上面可以看到,阻碍我们正确预测Y的就是因为存在confounder,即
在这里插入图片描述
然后,我们希望的是没有confounder的情况,也就是,
在这里插入图片描述
这两幅图对应分布的区别是什么呢?那就是T和W的独立性,换句话说,我们希望,

P(T,Y,W)=P(Y∣T,W)P(T)P(W) P(T,Y,W)=P(Y|T,W)P(T)P(W) P(T,Y,W)=P(YT,W)P(T)P(W)

然后实际上,有confounder的情况下他的分布应该是:
P^(T,Y,W)=P(Y∣T,W)P(T∣W)P(W) \hat{P}(T,Y,W)=P(Y|T,W)P(T|W)P(W) P^(T,Y,W)=P(YT,W)P(TW)P(W)

那有没有可能找到一种方法使得他们相等?解决这个问题的一般方法是reweight,也就是说:

P(T,Y,W)=P(T)P(T∣W)P^(T,Y,W)=βP^(T,Y,W) P(T,Y,W)=\frac{P(T)}{P(T|W)}\hat{P}(T,Y,W)=\beta \hat{P}(T,Y,W) P(T,Y,W)=P(TW)P(T)P^(T,Y,W)=βP^(T,Y,W)

通过乘以一个权重,这两个分布就可以做一个转换了,实际上,更方便的方法是对T作一个“干预操作“,这意味着T的取值一定等于某个确定的数,于是此时, 不管W取什么值,T等于某个值的概率为1,这意味着,P(T∣W)=P(T)=1P(T|W)=P(T)=1P(TW)=P(T)=1。因此我们只要将P^(T,Y,W)/p(T∣W)=P(T,Y,W)\hat{P}(T,Y,W)/p(T|W)=P(T,Y,W)P^(T,Y,W)/p(TW)=P(T,Y,W)就可以变成分布PPP!也就是说,对于均值,我们有以下关系E[I(T=t)y]⟹E[I(T=t)Yp(t∣w)]E[I(T=t)y]\Longrightarrow E[\frac{I(T=t)Y}{p(t|w)}]E[I(T=t)y]E[p(tw)I(T=t)Y], 这里的I(T=t)I(T=t)I(T=t)其实就可以理解为一种干预操作。. 并且我们可以证明IPW的目标跟standardization treatment effect(另一种说法是back-door adjustment, 一种经过debiased的treatment effect计算方法)是等价的:

∑wE[Y∣T=t,W=w]P(W=w) \sum _{w} E[ Y|T=t,W=w] P( W=w) wE[YT=t,W=w]P(W=w)

于是

=∑wE[Y∣T=t,W=w]P(W=w)p(t∣w)p(t∣w)=∑w∑yyP(y∣T=t,W=w)P(W=w)p(t∣w)p(t∣w)=∑w∑yyP(y,T=t,W=w)1p(t∣w)=∑w∑y∑tI(T=t,W=w)yP(y,T=t,W=w)1p(t∣w)=∑wEt,y[I(T=t,W=w)y]1p(t∣w)=Et,y[I(T=t)yp(t∣w)] \begin{aligned} & =\sum _{w} E[ Y|T=t,W=w] P( W=w)\frac{p( t|w)}{p( t|w)}\\ & =\sum _{w}\sum _{y} yP( y|T=t,W=w) P( W=w)\frac{p( t|w)}{p( t|w)}\\ & =\sum _{w}\sum _{y} yP( y,T=t,W=w)\frac{1}{p( t|w)}\\ & =\sum _{w}\sum _{y}\sum _{t} I( T=t,W=w) yP( y,T=t,W=w)\frac{1}{p( t|w)}\\ & =\sum _{w} E_{t,y}[ I( T=t,W=w) y]\frac{1}{p( t|w)}\\ & =E_{t,y}\left[\frac{I( T=t) y}{p( t|w)}\right] \end{aligned} =wE[YT=t,W=w]P(W=w)p(tw)p(tw)=wyyP(yT=t,W=w)P(W=w)p(tw)p(tw)=wyyP(y,T=t,W=w)p(tw)1=wytI(T=t,W=w)yP(y,T=t,W=w)p(tw)1=wEt,y[I(T=t,W=w)y]p(tw)1=Et,y[p(tw)I(T=t)y]

最后一条等式考虑了p(t∣w)\displaystyle p( t|w)p(tw)中的条件w\displaystyle ww,默认了是W=w\displaystyle W=wW=w,所以为了方便起见,就把I\displaystyle II中的W\displaystyle WW忽略掉,其实是需要加上去的,不过好像大家都是这么写的,就忽略掉了。

最后我们的计算公式就是如下:

τ≜E[Y(1)−Y(0)]=E[1(T=1)Ye(W)]−E[1(T=0)Y1−e(W)]\tau \triangleq E [ Y ( 1 ) - Y ( 0 ) ] = E [ \frac { 1 ( T = 1 ) Y } { e ( W ) } ] - E [ \frac { 1 ( T = 0 ) Y } { 1 - e ( W ) } ]τE[Y(1)Y(0)]=E[e(W)1(T=1)Y]E[1e(W)1(T=0)Y]

其估计值:

τ^=1n∑i(1(ti=1)yie^(wi)−1(ti=0)yi1−e^(wi))=1n1∑i:ti=1yie^(wi)−1n0∑i:ti=0yi1−e^(wi) \left. \begin{array}{l}{ \hat { \tau } = \frac { 1 } { n } \sum _ { i } ( \frac { 1 ( t _ { i } = 1 ) y _ { i } } { \hat { e } ( w _ { i } ) } - \frac { 1 ( t _ { i } = 0 ) y _ { i } } { 1 - \hat { e } ( w _ { i } ) } ) }\\{ = \frac { 1 } { n _ { 1 } } \sum _ { i : t _ { i } = 1 } \frac { y _ { i } } { \hat { e } ( w _ { i } ) } - \frac { 1 } { n _ { 0 } } \sum _ { i : t _ { i } = 0 } \frac { y _ { i } } { 1 - \hat { e } ( w _ { i } ) } }\end{array} \right. τ^=n1i(e^(wi)1(ti=1)yi1e^(wi)1(ti=0)yi)=n11i:ti=1e^(wi)yin01i:ti=01e^(wi)yi

参考资料

https://www.bradyneal.com/causal-inference-course

https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/

评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值