【原】倾向性得分加权分析

倾向性得分（Propensity Score，PS）就是在给定一组协变量的条件下，研究对象被分配到处理组的条件概率，这个概率一般通过logit或者probit模型来获取。倾向性得分本身并不能控制混杂，而是需要通过进一步的匹配、分层、回归校正、加权等方式来达到分析的目的。

既往倾向性得分分析一些笔记：

倾向性得分匹配（PSM）；

倾向性得分匹配续集（R笔记：MatchIt）；

中文网独此一份：多分组资料的倾向性得分匹配（TriMatch）；

SPSS之1:n倾向性得分匹配；

PS Matching安装攻略：SPSS也可以现1:n的倾向性得分分析了；

倾向性得分匹配每次匹配的对象都不一样？提前埋个“种子”试一试吧！；

倾向性得分校正；

还有一篇变量匹配分析：病例对照匹配（CCM）。

此次笔记是演示倾向性得分加权分析。权，然后知轻重，在古代权指的是秤砣，是个非常重要的重量标准。统计学里面的权重不同于一般的比重，它更强调的是指标的相对重要程度，倾向于贡献度或重要性，常见的有抽样权重/概率权重、逆方差权重、频数权重。

Complex Surveys:A Guide to Analysis Using R

理论上在进行IPW后，样本量会膨胀到原数据样本量的2倍，但连续变量和极端权重的截尾处理等情况的存在，大多数情况下加权后的样本量其实是达不到原来的样本量2倍的。样本量的增大往往会带来假阳性，而且对原始样本中的一个个体的加权获得的伪样本是完全相同的，不存在变异，从而使方差估计出现偏差。为获得与原人群相同的标准人群，降低假阳性的发生，常用稳定权重(stabilized weights，SW)来替代IPTW：Wt=Pt/PS，Wc=(1-Pt)/(1-PS)，其中Pt为不考虑协变量的处理组概率[中华流行病学杂志,2010,31(2):223-226]。

中国循证医学杂志.2022,22(3):365-372.

SPSS并未直接提供倾向值权重计算方法，但根据上表中的定义获得相应IPW权重并不难。

【1】获得给定协变量条件下的倾向值

建立logistic回归模型或者Probit回归模型，保存预测概率即可。

Analyze>>Regression>>Binary Logistic…

本例中的分类变量均为0/1型的二分类变量，按连续变量处理结果是一致的。

我们在《倾向性得分匹配（PSM）》中是通过PSM过程来生成倾向值的。

Data>>Propensity Score Matching…

这样我们就可以获得每个个案的倾向值，即下图中的PS：

【2】计算各个观测的IPW权重和SW权重

根据前面表格中的公式生成新的变量即可，当然也可以获取更多种类的权重。

计算Case组的逆概率权重：

Transform>>Computer Variable…

【if】：Include if case satisfies condition：mbsmoke = 1，【Continue】

Target Variable：PSW，Numeric Expression：1/PS

【OK】

计算Control组的逆概率权重：

Transform>>Computer Variable…

【if】：Include if case satisfies condition：mbsmoke = 0，【Continue】

Target Variable：PSW，Numeric Expression：1/(1-PS)

Change existing variable？【OK】

【OK】

计算稳定IPW权重SW。示例有4642例观察对象，其中孕期吸烟864例，Pt=864/4642。

Case组的稳健逆概率权重：

Transform>>Computer Variable…

【if】：Include if case satisfies condition：mbsmoke = 1，【Continue】

Target Variable：SW，Numeric Expression：(864/4642)/PS

【OK】

Control组的稳健逆概率权重：

Transform>>Computer Variable…

【if】：Include if case satisfies condition：mbsmoke = 0，【Continue】

Target Variable：SW，Numeric Expression：((4642-864)/4642)/(1-PS)

Change existing variable？【OK】【OK】

【3】加权样本的组间平衡性及结局变量在组间的差异比较

获得权重后，就可以进行伪样本的组间均衡性考察，以及平衡条件下的组间的差异比较了。

但正如前所言，加权后伪样本的样本量膨胀和加权样本的同质化使得方差估计（相应的标准误、置信区间估计）不能直接使用传统的方法。

在SPSS中有两个地方可以实现加权。一是数据处理中的频数加权（Data>>Weight Cases），另外一个是在一些分析过程中出现的加权选项（比如WLS权重，残差权重）。频数加权其实可以看做复制权重，权重是几样本量就是原来的几倍，很显然这种加权方法几乎适用于所有的SPSS分析过程，如果用IPW进行频数加权，参数值的点估计（效应值估计）结果问题不大，但组内方差的估计有偏甚至是错误的，方差偏小样本量偏大，标准误很明显会偏小，置信区间偏窄，样本量的膨胀也使得更容易获得有显著意义的P值。

在线性回归过程中出现的WLS加权（相当于进行加权最小二乘回归，minimizing sum(w*e^2)），参数值跟采用复制权重后的数据来估算的结果一致，但样本量采用了原始数据的样本量，如果采用IPW作为WLS权重，效应值估计值也没问题，但方差偏小样本量也相对偏小。而且这种方法也无法实现对分类变量的分析。

所以在SPSS里用计算的逆概率作为权重进行分析的结果都是不合适的，虽然效应值问题不大，但是标准误、置信区间都不合适。R语言lm函数对权重参数的释义时也特别指出：The sigma estimate and residual degrees of freedom may be suboptimal; in the case of replication weights, even wrong. Hence, standard errors and analysis of variance tables should be treated with care.

在观察性研究中使用稳定权重理论上可以得到与原始数据样本量相同的伪数据，治疗效应的方差估计可以直接通过传统回归模型进行估计[Value Health. 2010,13(2):273-277]。实际情况是SW加权后的样本量也不绝对与原始数据的样本量完全一样，但非常接近。所以笔者以为在SPSS中可以用SW作为频数加权来获取稳定的逆概率加权近似的分析结果，虽然采用SW作为权重进行WLS加权分析理论上也可以，但因为样本量的原因还是会存在一些偏差。不过两者结果已经非常接近了。

SW作为频数的组间均衡性结果如下。不过比较遗憾的是，经SW加权后年龄和教育水平仍未达到组间均衡。

最终还是把SPSS不能解决IPW加权分析交给R吧，R有包，包治百病。首先说明一下，逆概率加权分析不能直接通过指定lm{stats}或glm{MASS}函数的weights参数来进行，这里的设置weights参数后函数运行的实际上就是加权最小二乘法（weighted least squares），这笔者在前面的SPSS部分也有介绍。

Bwfit.ipw<-lm(bweight~mbsmoke,data=cattaneo,weights=ipw$ipw.weights)

#Bwfit.ipw<-glm(bweight~mbsmoke,family=gaussian,data=cattaneo,weights=ipw$ipw.weights)

#summary(Bwfit.ipw)

anova(Bwfit.ipw)

confint.default(Bwfit.ipw)

#confint(Bwfit.ipw)

#两者在计算置信区间时默认的方法是不一样的：lm进线性回归后，confint()提供了正态近似法（Asymptotic normality）计算的置信区间，而confint.default()可提供轮廓似然法（Profile likelihood）结果，glm进行线性回归后提供的confint()可提供轮廓似然法（Profile likelihood）结果。两者种方法的结果略有差异。SPSS采用的是confint()正态近似法

在R中实现实IPW加权后数据分析的常用包是survey。survey用于基于设计的（而非传统的基于模型）复杂调查样本的抽样的分析。svyglm always returns "model-robust" standard errors; the Horvitz-Thompson-type standard errors used everywhere in the survey package are a generalisation of the model-robust ’sandwich’ estimators. In particular, a quasi-Poisson svyglm will return correct standard errors for relative risk regression models.

##数据导入

library(bruceR)

cattaneo<- import("D:/Temp/cattaneo.sav" )

cattaneo$mmarried<-factor(cattaneo$mmarried,levels=c(0,1),labels=c("No","Yes"))

cattaneo$alcohol<-factor(cattaneo$alcohol,levels=c(0,1),labels=c("No","Yes"))

cattaneo$mbsmoke<-factor(cattaneo$mbsmoke,levels=c(0,1),labels=c("No","Yes"))

cattaneo$mrace<-factor(cattaneo$mrace,levels=c(0,1),labels=c("No","Yes"))

cattaneo$frace<-factor(cattaneo$frace,levels=c(0,1),labels=c("No","Yes"))

cattaneo$fbaby<-factor(cattaneo$fbaby,levels=c(0,1),labels=c("No","Yes"))

cattaneo$prenatal1<-factor(cattaneo$prenatal1,levels=c(0,1),labels=c("No","Yes"))

#分类变量设置为因子，这几个分类变量都是0/1型二分类变量，不设置因子按连续变量处理也没有问题。后续CreateTableOne和svyCreateTableOne{tableone}的就不需要再指定factorVars了。需要特别说明一下，此处如果采用levels和labels参数设置了标签，标签将代替数字。

##基线考察

library(tableone)

dput(names(cattaneo)) #获得的结果可以方便后续命令的书写

c("id", "bweight", "mmarried", "alcohol", "mage", "medu", "fage", "mbsmoke", "mrace", "frace", "fbaby", "prenatal1", "PS", "PSW","SW")

Sum.vars<-c("mmarried","alcohol","mage","medu","fage","mrace","frace","fbaby","prenatal1","bweight") #指定要分析的变量，去掉id、分组变量等不需要的一些变量

Cat.vars<-c("mmarried","alcohol","mbsmoke","mrace","frace","fbaby","prenatal1") #分类变量

Tab1.stas<-CreateTableOne(vars=Sum.vars,strata="mbsmoke",data=cattaneo,factorVars=Cat.vars,test=TRUE,smd=TRUE,addOverall=TRUE) #创建基线表。var指定所有分析变量，包括分组变量，结局变量可以不用指定，我们在Sum.vars也指定了结局变量bweight，可以顺带看一下单变量分析结果；factorVars指定应该按分类变量处理的那些用数字编码的变量，本例在导入数据后已经将分类变量设置为了因子，此处其实可以无需再指定

print(Tab1.stas,catDigits=2,contDigits=2,pDigits=3,test=TRUE,smd=TRUE) #查看结果。设置分类和连续变量结果的小数点，显示3位数的P值，显示标准化均数差SMD

二分类变量默认只显示其中高水平的结果，多分类变量则显示所有水平。结果中貌似父亲的年龄不成正态。我们假设父亲fage和孕母mage的年龄不满足正态的要求，可以单独把非正态的数据列出。在分类分析中，如果样本量和单元格的期望频数太小，常常需要改为精确检验，这里假设alcohol望频数低于5，需要进行精确检验。

Nonnorm.vars <-c("mage","fage")

Exact.vars<-"alcohol"

Tab1<-print(Tab1.stas,catDigits=2,contDigits=2,pDigits=3,test=TRUE,smd=TRUE,noSpaces=TRUE,exact=Exact.vars,nonnormal=Nonnorm.vars) #设置分类和连续变量结果的小数点，显示3位数的P值，显示smd，去掉空格，对非正态连续变量和样本量比较小的分类变量执行特殊的检验

write.csv(Tab1,file="Table1.csv") #导出可以用excel打开的csv文件

结果显示孕母吸烟组的新生儿的体重明显低于不吸烟组，但各个变量在吸烟者和不吸烟组之间也是不均衡的。为了校正这些干扰因素的影响，多重回归、倾向性得分匹配都是常用的方法。此次笔记笔者演示的是加权分析的过程，要做加权分析首先需要获得权重，我们仍以获得逆概率权重（IPW）及其稳定权重（SW）为例。

##获得倾向值，计算IPW和SW

PSfit<-glm(mbsmoke~mmarried+alcohol+mage+medu+fage+mrace+frace+fbaby+prenatal1,family=binomial(link=logit),data=cattaneo)

PS<-PSfit$fitted.values

iptw<-ifelse(cattaneo$mbsmoke=="Yes",1/PS,1/(1-PS)) #获得IPW权重。特别注意此处ifelse条件是[=="Yes"]而不是[==1]，因为本例在导入数据后已经将分类变量设置为了因子，并且通过levels和labels参数设置了标签。如果用[==1]，则所有个案的值均按1/PS计算。如果在导入数据后没有设置因子或设置为了因子但没有设置标签者则用[==1]

summary(cattaneo$mbsmoke)  #获得sw权重。

0     1

3778 864

sw<-ifelse(cattaneo$mbsmoke=="Yes",(864/4642)/PS,(3778/4642)/(1-PS))

##加权样本的组间均衡性考察

library(tableone)

library(survey)

library(grid);library(survival)

Svy.design<-svydesign(ids=~1,weights=iptw,data=cattaneo) #指定调查设计的信息。ids=~1或~0表示没有标识集聚性（多水平）的变量

Tabipw.stas<-svyCreateTableOne(vars=Sum.vars,strata="mbsmoke",data=Svy.design,test=TRUE,smd=TRUE,addOverall=TRUE) #本例在导入数据后已经将分类变量设置为了因子，此处无需再指定factorVars参数

Tab.ipw<-print(Tabipw.stas,catDigits=2,contDigits=2,pDigits=3,test=TRUE,smd=TRUE,noSpaces=TRUE,nonnormal=Nonnorm.vars)

write.csv(Tab.ipw, file = "Table.ipw.csv")

结果显示，在进行逆概率加权后，虚拟的样本为9177.82，大部分的干扰因素在组间已达到均衡（P>0.05,标准化均数差(SMD)<0.1），不过孕母的年龄的组间仍未均衡性（P=0.002<0.05，SMD=0.127>0.1）。标准化均数差（SMD）是根据一个或多个变量来衡量两组平均值间距离的度量，在实践中，它经常被用作倾向得分匹配前后个体协变量的均衡性的度量指标。

我们也可以看一下SW的加权结果：

Svy.design<-svydesign(ids=~1,weights=sw,data=cattaneo)

Tabsw.stas<-svyCreateTableOne(vars=Sum.vars,strata="mbsmoke",data=Svy.design,test=TRUE,smd=TRUE,addOverall=TRUE) 

Tab.sw<-print(Tabsw.stas,catDigits=2,contDigits=2,pDigits=3,test=TRUE,smd=TRUE,noSpaces=TRUE,nonnormal=Nonnorm.vars)

年龄的组间均衡性仍然不好，应考虑其他权重，可参见本文最后部分。为演示完整的分析步骤，我们先把结果分析做完。

##加权模型进行结果分析：加权样本的组间效应考察

Smoke.ipw<-svyglm(bweight~mbsmoke,design=Svy.design, family=gaussian())

summary(Smoke.ipw)

几个计算权重的R包

前面我们通过IPW、SW加权，协变量年龄仍未达到组间均衡。这时我们可以尝试其他的一些权重。权重值的计算在R中有很多的程序包。比如ipw、PSW、PSweight、twang、WeightIt、hrIPW等。这些程序包不仅提供了各种各样的权重的计算，而且有些包还可以考察变量的组间均衡性，以及实现处理效应、增强估计等功能。

library(bruceR)

cattaneo<- import("D:/Temp/cattaneo.sav" )

#重新导入数据。前面我们导入的数据中含有带标签的因子，这将使ipwpoint{ipw}无法计算IPTW。psw{PSW}要求所有的分类变量编码为哑变量

##ipwpoint{ipw}:Estimate Inverse Probability Weights

library(ipw)

ipw<-ipwpoint(exposure=mbsmoke, family="binomial", link="logit", denominator=~mmarried+alcohol+mage+medu+fage+mrace+frace+fbaby+prenatal1, data=cattaneo)

ipw$ipw.weights #IPTW权重

##psw{PSW}:Propensity Score Weighting Methods for Dichotomous Treatments

library(PSW)

ipw<-psw(data=cattaneo,form.ps=mbsmoke~mmarried+alcohol+mage+medu+fage+mrace+frace+fbaby+prenatal1,weight="ATE",std.diff=TRUE) #weight参数可指定的权重："ATE"、 "ATT", "ATC"、"MW"、"OVERLAP"、"TRAPEZOIDAL"；std.diff指定是否计算标准化均数差

summary(ipw) #结果会显示倾向值计算模型、估算的倾向值、权重、加权前后的标准化均数差等

ipw$W #权重值

加权后各个变量的标准化均数差大多在-10~10之内，均衡性明显变好。SMD均衡性的界值定义一般为绝对值0.1，但在该程序包中SMD定义如下，是常规定义的100倍：

##weightit{WeightIt}: Weighting for Covariate Balance in Observational Studies

library(WeightIt)

SW<-weightit(mbsmoke~mmarried+alcohol+mage+medu+fage+mrace+frace+fbaby+prenatal1,data=cattaneo,method="glm",estimand="ATE",stabilize=TRUE) #method可选"gbm"、"cbps"、"npcbps" 、"ebal"、"optweight" 、"super"、"bart" 、"energy"；estimand可选"ATE"、"ATT"、 "ATC"、"ATO"、"ATM"、"ATOS"；stabilize可指定是否估计稳定权重，等于个体的估计权重乘以其实际所在组的个案数占所有个案数的比例

summary(SW)

SW$weights 

结果同前面我们手动计算的稳定的IPW（SW）完全一致。

##PSweight{PSweight}/SumStat{PSweight}:Propensity Score Weighting for Causal Inference with Observational Studies and Randomized Trials

library(PSweight)

PSwt<-PSweight(ps.formula=mbsmoke~mmarried+alcohol+mage+medu+fage+mrace+frace+fbaby+prenatal1,yname="bweight",data=cattaneo,weight="IPW",ps.method="glm") #weight参数可指定的权重："IPW"、"treated"、"overlap"、"matching"、"entropy"，默认是重叠加权"overlap"；ps.method指定倾向值的计算法方法，可用"glm"(默认)、"gbm"、"SuperLearner"

PSwt$propensity #阴性结果和阳性结果的倾向值

PSwt[["propensity"]][,2] #阳性结果的倾向值

SumSt<-SumStat(ps.formula=mbsmoke~mmarried+alcohol+mage+medu+fage+mrace+frace+fbaby+prenatal1,data=cattaneo,weight=c("IPW","treated","overlap","matching","entropy"),method="glm") #weight参数可指定的权重："IPW"、"treated"、"overlap"、"matching"、"entropy"，默认是重叠加权"overlap"；ps.method指定倾向值的计算法方法，可用"glm"(默认)、"gbm"、"SuperLearner"

summary(SumSt) #结果中会显示加权前后各组的协变量均值和SMD

SumSt[["propensity"]][,2] #阳性结果的倾向值

SumSt[["ps.weights"]][,3] #IPW权重

tail(SumSt[["ps.weights"]])

SumStat计算的权重明显偏小，应该是经过了某种转化，但这并不影响最终结果的计算。

SumSt[["ess"]] #有效样本量，effective sample sizes。有效样本量定义参照Am J Epidemiol. 2019;188(1): 250–257.

plot(SumSt,type="balance") #可绘制估算的倾向值的直方图（"hist")、估算的倾向值的密度图（"density"）、均衡统计量（("balance"）

其中除了IPW，ATT加权（treated）、重叠加权（overlap）、匹配加权（matching）、熵加权（entropy）都可以使协变量在组间达到均衡。

最后，我们用重叠加权后组间均衡和结果分析作为此次推文的结束吧：

library(tableone)

library(survey)

Sum.vars<-c("mmarried","alcohol","mage","medu","fage","mrace","frace","fbaby","prenatal1","bweight")

Cat.vars<-c("mmarried","alcohol","mbsmoke","mrace","frace","fbaby","prenatal1")

Nonnorm.vars <-c("mage","fage")

Svy.design<-svydesign(ids=~1,weights=SumSt[["ps.weights"]][,5],data=cattaneo)

Tabow.stas<-svyCreateTableOne(vars=Sum.vars,strata="mbsmoke",data=Svy.design,factorVars=Cat.vars,test=TRUE,smd=TRUE,addOverall=TRUE) 

Tab.ow<-print(Tabow.stas,catDigits=2,contDigits=2,pDigits=3,test=TRUE,smd=TRUE,noSpaces=TRUE,nonnormal=Nonnorm.vars)

你要觉得这种虚拟的样本量看起来不舒服，可以根据一开始的表格手动计算OW权重：

PSfit<-glm(mbsmoke~mmarried+alcohol+mage+medu+fage+mrace+frace+fbaby+prenatal1,family=binomial(link=logit),data=cattaneo)

PS<-PSfit$fitted.values

ow<-ifelse(cattaneo$mbsmoke==1,1-PS,PS)

Svy.design<-svydesign(ids=~1,weights=ow,data=cattaneo)

Nonnorm.vars <-c("mage","fage")

Tabow.stas<-svyCreateTableOne(vars=Sum.vars,strata="mbsmoke",data=Svy.design,factorVars=Cat.vars,test=TRUE,smd=TRUE,addOverall=TRUE) 

Tab.ow<-print(Tabow.stas,catDigits=2,contDigits=2,pDigits=3,test=TRUE,smd=TRUE,noSpaces=TRUE,nonnormal=Nonnorm.vars)

这样看起来结果就舒服多了。

最终的分析结果如下：

Smoke.ow<-svyglm(bweight~mbsmoke,design=Svy.design, family=gaussian())

summary(Smoke.ow)

— END —

关注“统计浆糊”

获取更多信息