基于 PS (倾向评分)的IPTW 法首先由Rosenbaum作为一种以模型为基础的直接标准化法提出,属于边际结构模型。简单来说,就是把许多协变量和混杂因素打包成一个概率并进行加权,这样的话,我只用计算它的权重就可以了,方便了许多。那么,如何将多个协变量的影响用一个倾向评分值来表示呢? 即如何估计倾向评分值呢? 根据 Rosen-baum 和 Rubin 的定义:倾向评分值为在给定一组协变量(X i )条件下,研究对象 i(i =1,2,…N)被分配到某处理组或接受某暴露因素(Z i =1)的条件概率。理论上,所有可计算得到该条件概率的方法均可用于估计倾向评分值。
采用 logistic 回归模型估计倾向评分值具有模型
简单、容易实现、可直接得到倾向评分值、结果易于解释等显著优势。
我们拿logistic 回归做例子:
logistic 回归是最早提出的估计倾向评分值的方法,由于其原理为人们所熟悉且容易实现,
也是目前最常用的估计方法。logistic 回归模型如下:
假设为二元logistic 回归,右边一系列混杂因素的方程会生成一个0-1之间的目标事件时间发生的概率,概率越大代表事件发生的可能性越大,这样就等于把多个混杂因素做成了一个综合评分来表示。逆概率加权(IPTW)是利用倾向性评分的倒数来处理数据间混杂的一种方法。Robins等给出的加权系数(形)计算方法是:处理组观察单位的权数Wt=1/PS,对照组观察单位的权数Wc=1/(1一PS)。PS为观察单位的倾向评分值。此方法得到的
人群往往与原来人群的数量不同,因此虚拟人群各变量的方差大小可有变化,而且有较低 PS 的处理组对象与有较高 PS 的非处理组对象将会获得很大的权重。由于非常大的权重会诱导不稳定性。Heman等人对计算方法进行调整,将整个研究人群的处理率和非处理率加入公式调整后得到稳定权数(stabilized weights)。具方法是:处理组观察单位的权数Wt=Pt/PS,对照组观察单位的权数Wc=(1–Pt)/(1–PS)。 目前蛮多文章使用稳定权重。
综上所述:倾向评分加权法首先将多个主要混杂变量的信息综合为一个变量倾向评分 ,然后将倾向评分作为需要平衡的混杂因素 ,通过标准化法的原理加权,使各对比组中倾向评分分布一致 ,则达到使各混杂因素在各比较组分布一致的目的。
下图所示,经过处理后,各项数据的分布基本上被配平了,更有利于做时候比较分析,有点起到RCT的效果。
既往我们已经演示了R语言进行倾向评分逆概率加权(IPTW)。后台有粉丝问我,SPSS能不能做倾向评分逆概率加权(IPTW),我研究了一下,也是可以的,步骤和R的差不多。
1.导入数据,整理好数据
2.把需要比较的变量拿出来做结果变量,其他的变量做协变量,组成logistic回归方程,生成预测值
3.通过预测值计算出拟概率权重
4.建立回归方程,加入权重。
今天我们SPSS来演示倾向评分逆概率加权(IPTW),先导入我们的数据,
这是一个关于早产低体重儿的数据(公众号回复:早产数据,可以获得该数据),低于2500g被认为是低体重儿。数据解释如下:low 是否是小于2500g早产儿,age 母亲的年龄,lwt 末次月经体重,race 种族,smoke 孕期抽烟,ptl 早产史(计数),ht 有高血压病史,ui 子宫过敏,ftv 早孕时看医生的次数
bwt 新生儿体重数值。
假设我们想比较有无高血压病史对早产低体重儿的影响,我们先建立没有加权的模型看看。
分析----回归----二元logistic 回归
二元logistic 回归的具体操作我就不演示了,可以看我相关的文章,这里直接上结果,下图是没有加权时候的结果,可以看到ht的系数是1.862,OR值为6.438
步骤的关键就是生成权重W,我们这里以Robins的公式来计算权重,处理组观察单位的权数Wt=1/PS,对照组观察单位的权数Wc=1/(1一PS)。PS为观察单位的倾向评分值。
我们先以ht为结果变量,其他变量为协变量建立二元logistic 回归方程,
在分类中,把分类变量选入分类协变量,变化量选择第一个(就是以第一个为参考的意思)
在保存中选上概率
我们主要就是要这个概率值,其他就不用选了,直接确定,得到概率数据
接下来就是通过公式来计算权重,
点击转换----计算变量
我们把权重命名为W,点击如果
点击:在个案满足条件时包括
双击变量预测概率,然后填入
ht=1
然后点击继续
双击预测概率,输入公式
点击黏贴
进入下图界面,把公式复制下来
然后进行黏贴,注意一下,句号也是要的
我们把第二个公式修改一下,改成这样:
IF (ht=0) W=1 / (1-PRE_1).
这里注意一下,修改公式要在全英文输入法界面修改,括号要在句号前面。
最后点击
运行---------全部
打开数据界面,我们最关键的权重W终于生成了
我们调整一下它的小数位数
我们的权重生成后就可以进行回归方程加权了,不能使用我们的常见的个案加权,这样SPSS会报错(下图)
我们要通过广义线性方程来对它进行加权计算,点击:
分析—广义线性模型—广义线性模型
模型部分选择:二元logistic 回归
响应这里依次选入因变量和因变量设置
权重W填入下面
预测变量这里填入分类和连续变量
选项这里,因为我们是以最低值为参考的,所以这里选择降序
模型这里选入所有变量
最后点击确定,得到最终结果
逆概率加权后,ht的系数为1.941,P还是小于0.05的,证明了加权后对加过影响不大,结果还是有意义的。
为了验证这样计算是否正确,我用R语言重新做了一遍,下面是R的代码:
bc<-read.csv("E:/r/test/zaochan.csv",sep=',',header=TRUE)
names(bc)
#########
summary(bc$race)
bc$race<-ifelse(bc$race=="black",1,ifelse(bc$race=="white",2,3))
bc$smoke<-ifelse(bc$smoke=="nonsmoker",0,1)
bc$low<-factor(bc$low)
bc$race<-factor(bc$race)
bc$ht<-factor(bc$ht)
bc$ui<-factor(bc$ui)
################
fit<-glm(low ~ age + lwt + race + smoke + ptl + ht + ui + ftv,
family = binomial("logit"),
data = bc)
summary(fit)
exp(confint(fit))
exp(coef(fit))
############
pr<- glm(ht ~age + lwt + race + smoke + ptl + ui + ftv, data=bc,
family=binomial(link = "logit"))
pr1<-predict(pr,type = "response")
pr1
W <- (bc$ht==1) * (1/pr1) + (bc$ht==0) * (1)/(1-pr1)
bc$W<-W
##########
fit.IPTW<-glm(low ~ age + lwt + race + smoke + ptl + ht + ui + ftv,
family = binomial("logit"),data = bc,weights = W)
summary(fit.IPTW)
exp(confint(fit.IPTW))
exp(coef(fit.IPTW))
R的结果为:
和我们使用SPSS计算的完全一样,证明我们使用SPSS计算没有问题。
参考文献:
[1]李智文, 刘建蒙, 任爱国. 基于个体的标准化法——倾向评分加权[J]. 中华流行病学杂志, 2010(2):4.
[2] Robins J M , MN Hernán, Brumback B . Marginal Structural Models and Causal Inference in Epidemiology[J]. Epidemiology, 2000, 11.
[3]吴美京, 吴骋, 王睿,等. 倾向性评分法中评分值的估计方法及比较[J]. 中国卫生统计, 2013, 30(3):5.
[4]涂娇, 陈卉, 倪平,等. 倾向评分逆处理概率加权在生存分析中的应用及SAS宏实现[J]. 中国临床药理学与治疗学, 2015, 50(9):9.
[5]冯士雍. 加权回归分析的适用场合,加权回归方程效果的比较及权的估计[J]. 数理统计与管理, 1983(6):8.
[6] 辩证看待倾向性评分法 https://blog.csdn.net/fjsd155