【原】多分组资料的倾向性得分加权分析，以逆概率加权为例

倾向性得分加权分析的理论概述及两组的倾向性得分加权分析见《倾向性得分加权分析》。文中示例结局变量为连续变量（新生儿体重），研究因素（孕母是否吸烟）为二分类变量，我们借助R包tableone、survey来完成了处理因素为二分类的倾向性得分加权分析，在计算权重时介绍了ipw、PSW、WeighIt、PSweight包。

在真实的临床中，研究因素也可能是连续变量或者多水平的分类变量，结局变量除了连续变量也会存在分类变量的情况。在这种多分组的资料中，协变量如果在组间不均衡，常用多重回归来进行校正，以前也介绍过多分组资料的倾向性得分匹配。此次笔记演示的是如何用倾向性得分逆概率加权来分析这种组间不均衡的多分组资料。

本文中分组变量、处理因素同意，均为主要的研究因素。涉及的几个主要R包：tableone、cobalt、survey、WeighIt、PSweight。

示例：同<<多分组资料的倾向性得分匹配(TriMatch)>>，数据来自TriMatch包。分析目的是考察吸烟程度对年度医疗支出TOTALEXP的影响，考虑的控制因素LASTAGE（调查时年龄；19-94岁）、MALE（性别；0=male，1=female）、RACE3（种族；1=白人white，2=黑black人，3=其他other）、marital（婚姻状况；1=已婚married，2=丧偶widowed，3=离异divorced，4=分居separated，5=未婚never married）、educate（教育水平；1=大学毕业生college graduate，2=大学辍学或在读some college，3=高中毕业生high school graduate，4=其他other）、POVSTALB（贫困状况；1=贫困poor，2=接近贫困near poor，3=低收入low ncome，4=中等收入middle income，5=高收入high income）、SREGION（普查区域；1=Northeast, 2=Mid west, 3=South, 4=West）。

packyears是一种累积吸烟的度量，它结合了受试者自我报告的吸烟程度和持续时间，我们可以根据这个变量定义一个标识吸烟程度的新变量smoke.ord。为了让研究跟具有代表性，我们选择年龄＞=40岁的受试者。

【1】数据导入及预处理

data(nmes,package="TriMatch") #从TriMatch包中载入数据National Medical Expenditure Study数据。需先安装TriMatch包

nmes<-nmes[nmes$age==1, ] #筛选40岁及以上受试者作为研究对象

nmes<-nmes[,c("PIDX","LASTAGE","MALE","RACE3","smoke","packyears","TOTALEXP","lc5","chd5","educate","marital","SREGION","POVSTALB")] #选择要分析的变量

#将分类变量设置成因子

nmes$smoke<-factor(nmes$smoke,levels=c(0,1,2),labels=c("never","smoker","former"))

nmes$MALE<-factor(nmes$MALE,levels=c(0,1),labels=c("male","female"))

nmes$RACE3<-factor(nmes$RACE3,levels=c(1,2,3),labels=c("white","black","other"))

nmes$marital<-factor(nmes$marital,levels=c(1,2,3,4,5),labels=c("married","widowed","divorced","separated","never married"))

nmes$educate<-factor(nmes$educate,levels=c(1,2,3,4),labels=c("college graduate","some college","high school graduate","other"))

nmes$POVSTALB<-factor(nmes$POVSTALB,levels=c(1,2,3,4,5),labels=c("poor","near poor","low income","middle income","high income"))

nmes$SREGION<-factor(nmes$SREGION,levels=c(1,2,3,4),labels=c("Northeast","Mid west","South","West"))

#根据packyears构建表示吸烟程度的三分类变量smoke.ord

subnmes<-nmes[!is.na(nmes$packyears),]

medSF<-median(subnmes[subnmes$smoke !="never",]$packyears) 

medSF

[1] 28

nmes$smoke.ord<-ifelse(nmes$smoke=="never","Never",ifelse(nmes$packyears>=28,"Heavy",ifelse(nmes$packyears>=0, "Moderate",NA))) #根据packyears定义表示吸烟程度的因子变量smoke.ord：0=不吸烟，1=累积吸烟量<packyears中位数，2=累积吸烟量>=packyears中位数

#nmes$smoke.ord<-factor(nmes$smoke.ord,order=TRUE) #注意此处若设置为有序多分类，后面在采用weightit{WeightIt}或SumStat{PSweight}默认的参数进行权重计算时，计算结果会有所不同。因为两者对有序多分类变量的权重估计采用了不同的方法。此处仅展示并不打算运行

nmes<-na.omit(nmes) #暴力删除有缺失值的记录

【2】基线均衡性考察

library(tableone)

Sum.vars<-c("LASTAGE","MALE","RACE3","educate","marital","SREGION","POVSTALB","TOTALEXP") #需要分析的协变量集合

Cat.vars<-c("MALE","RACE3","educate","marital","SREGION","POVSTALB") #分类变量集合Nonnorm.vars <-"TOTALEXP" #不满足正态分布的变量

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

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

#write.csv(Tab1,file="Table1.csv") #将结果导出到默认目录，保存为名称是Table1的csv文件

虽然不同吸烟程度患者的医疗支出不同（P<0.001）,但大量的协变量在组间同样不均衡（P<0.05,SMD>0.1）。为得到更为严谨的结论，需要对这些干扰因素进行校正。

【3】通过程序包获得倾向值并考察协变量在处理因素各水平间的均衡性

借助程序包不仅可以快速实现各种权重的计算，而且可以对各个协变量在加权前后的均衡性进行考察。处理因素为多分类，推荐考虑weightit{WeightIt}或SumStat{PSweight}。WeightIt包提供了更多的倾向值估计方法，默认的glm处理因素可以是二分类、无序多分类、有序多分类、连续变量及纵向数据，但不能绘制展示均衡性的可视化图形，可借助cobalt包的love.plot()、bal.plot()等函数来弥补。PSweight包倾向值估计方法相对较少，默认的glm处理因素只能是二分类和无序多分类，但可以绘制love图、直方图、密度图等。

PS.fmu<-smoke.ord~LASTAGE+MALE+RACE3+educate+marital+SREGION+POVSTALB 

#为减少后面命令的字符输入，先定义好倾向值的计算模型

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

library(WeightIt)

Wtit<-weightit(formula=PS.fmu,data=nmes,method="glm",link="logit",estimand="ATE",stabilize=TRUE)

#weightit实际上是一个函数调度器，可以调用很多方法来执行均衡权重的估计，method可选"glm"、"gbm"、"cbps"、"npcbps" 、"ebal"、"optweight" 、"super"、"bart" 、"energy"；本例使用了默认的"glm"，也可用"ps"表示，即Propensity score weighting using generalized linear models，该法可用于处理因素为二分类、多分类和连续资料的数据。处理因素为二分类默认链接参数为"logit"（link= "logit")，link可用"logit"或"probit"。处理组为无序多分类，默认调用mlogit包中mlogit函数，link="logit"或"probit"；link="bayes.probit"，将调用MNP包中的mnp函数执行Bayesian probit；link="br.logit"，调用brglm2包的brmultinom函数执行biased-reduced multinomial logistic regression。处理因素如为有序多分类，默认调用MASS包的polr函数，也可通过link="br.logit"调用brglm2包的bracl函数。处理因素为连续变量，则执行线性回归。处理因素为纵向数据，权重为每个时间点的估计权重。其他method可参见WeightIt包的帮助文件

#estimand可选"ATE"、"ATT"、 "ATC"、"ATO"、"ATM"、"ATOS"

#stabilize可指定是否估计稳定权重，等于个体的估计权重乘以其实际所在组的个案数占所有个案数的比例

summary(Wtit)

iptw<-Wtit$weights #因求权重时stabilize=TRUE此处获得的实际上是稳定的IPW

library(cobalt)

bal.tab(Wtit,stats=c("m","v"),continuous="std",binary="std",thresholds=c(m=0.1,v =2),pairwise=TRUE,abs=FALSE,un=TRUE)

#stats指定需要输出的指示均衡性的统计量，可使用简写。分类变量可使用"mean.diffs"（简写"m"）、"variance.ratios"（简写"v"）、"ks.statistics"、"ovl.coefficients"（简写"ovl"），连续变量可使用"correlations"（简写"cor"）、"spearman.correlations"（简写"sp"）、"mean.diffs.target"（简写"m"）、"ks.statistics.target"（简写"ks"）。选择使用方差比或KS统计量，所显示的值可能不是来自相同的两两比较，也就是说最大的标准化均值差和最大的方差比可能不是来自相同的比较组

#continuous和binary指定计算统计量时变量标准化("std")还是原始数据 ("raw")。continuous默认"std"，而binar默认是"raw"

#thresholds指定各指标的均衡阈值，比如SMD常用0.1。指定后结果中将额外输出各个协变量是否超过阈值，并汇总统计超过和在阈值内的变量数量等信息

#pairwise处理因素任意两水平之间的均衡性（TRUE，默认）还是每个水平与所有水平组合之间的均衡性（FALSE）

#abs是否显示绝对值

#un是否显示加权前的（未校正）的相应结果

本例均衡性指标不仅选择了SMD，还选择了方差比。结果显示加权前后各协变量在处理因素各水平之间的均衡性，我们可以看到经过IPTW加权后，协变量的组间达到了均衡。默认的是只显示多分组中两两成对比较中差值的最大值（逻辑是最大不均衡如果是可以接受的，那么该变量的所有其他不均衡也是可以容忍的，因此专注于减少最大的不均衡就足以减少总体不均衡），还有另外一种方式是每个组与所有组的比较结果，这个可通过参数pairwise来设置。如果你想查看所有两两比较的结果，可以用以下命令来实现：

bal.tab(Wtit,which.treat=.all,disp="means", un=TRUE)

#which.treat=.all显示治疗因素任意两水平间的比较

#disp显示汇总哪些分布统计信息，可选均值（"means" ）或标准差（"sds"）

#un是否显示未校正的相应结果

我们只展示一部分结果。几个缩写（M/SD.0/1.Un/Adj，Diff.Un//Adj）下释义：M.0.Un表示未校正对照组均值，M.1.Adj则表示校正的处理组均值，Diff.Un表示校正前的标准化均数差SMD。

#绘制可视化均衡性分布图

bal.plot(Wtit,var.name="LASTAGE", which="both",type="density")

#var.name指定绘图变量；

#which显示校正("adjusted")/未校正("unadjusted")/同时显示两者("both")分布均衡性结果

#type指定绘图类型，densities ("density"), histograms ("histogram"), empirical cumulative density function plots ("ecdf")

#绘制love图

love.plot(Wtit,threshold=0.1,agg.fun="range",binary="std") #"mean","max","range"

love.plot(Wtit,thresholds=c(m=0.1),binary="std",which.treat=.all,abs=FALSE)

加权后各协变量的SMD均在-0.1~1之间，均衡性确实得到了非常好的改善。

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

library(PSweight)

SumSt<-SumStat(ps.formula=PS.fmu,data=nmes,weight=c("IPW","treated","overlap","matching","entropy"),method="glm") #weight参数可指定的权重："IPW"、"treated"、"overlap"、"matching"、"entropy"，默认是重叠加权"overlap"；ps.method指定倾向值的计算法方法，可用"glm"(默认)、"gbm"、"SuperLearner"

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

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

SumStat计算的权重明显偏小，应该是对样本量做了归一化的缩放处理（后续svyCreateTableOne{tableone}使用该权重进行分析发现加权基线表中各组的样本量均为1），具体怎么转化的就不纠结了，反正并不影响最终结果的计算。

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

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

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

使用获得的6种权重绘制love图显示协变量的均衡性均得到了较好的提升，其中年龄的均衡性上，匹配加权略微差一点。

【4】加权数据的组间均衡性考察

在第【3】步中其实已经利用计算权重的程序包对协变量在组间的均衡性进行了考察。此处再利用svyCreateTableOne{tableone}进行考察。

library(tableone)

library(survey)

library(grid);library(survival)

Svy.design<-svydesign(ids=~1,weights=iptw,data=nmes) #指定调查设计的信息。ids=~1或~0表示没有标识集聚性（多水平）的变量。此处使用的是weighit函数生成的权重

Tabipw.stas<-svyCreateTableOne(vars=Sum.vars,strata="smoke.ord",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") #将结果导出到默认目录，保存为名称是Table1的csv文件

结果显示，在进行逆概率加权后，干扰因素在组间已达到均衡（P>0.05,标准化均数差(SMD)<0.1）。吸烟程度不同确实对医疗支出的影响是不一样的（P<0.001，SMD=0.045）。

【5】加权数据结果分析：加权数据的组间差异比较

加权后的分析结果我们在考察组间均衡性的时候其实一并考察过了。在组间均衡性的基础上，我们也可以单独对此进行分析。结局变量为连续变量，我们可以考虑survey包中的两样本的t检验函数svyttest或者通过广义线性回归模型函数svyglm来实现，svycontrast进行成对比较，但本例的医疗支出TOTALEXP呈明显的偏态分布，以上这些参数检验法是不适用的，应考虑使用函数svyranktest进行非参数检验Kruskal-Wallis H检验。

笔者未找到survey包中正态性的检验方法，可以考虑看下中位值和均值的关系：

svyby(~TOTALEXP,~smoke.ord,Svy.design,svymean, vartype=c("se","var"))

svyby(~TOTALEXP,~smoke.ord,Svy.design,svyquantile, quantiles=0.5)

可以看出均值与中位数偏离非常大。

Kruskal-Wallis H检验如下：

svyranktest(TOTALEXP~smoke.ord,design=Svy.design,test="KruskalWallis") #c("wilcoxon","vanderWaerden","median","KruskalWallis")

结果显示在进行IPW加权后，吸烟程度不同的患者医疗支出的差异具有明显的统计学意义（Chi2=24.13，P<0.001）。

比较遗憾的是未能找到多重比较（Multiple Comparison）的函数，但事后的两两比较在Kruskal-Wallis H检验有统计学意义后也我们是想进一步了解的。我们可以用最原始的办法进行一下变通：分别进行两两比较，对显著水平α值做Bonferroni校正。k组完全进行两两比较的次数=k(k-1)/2，本例三组，两两比较需要比较三次3次，α_adj=0.05/3≈0.0167。

通过survey包中的subset.survey.design选择子集，然后再进行三次的两组比较，显著水平用0.0167而不再是0.05来判断。

NM<-subset(Svy.design, smoke.ord!="Heavy")

NH<-subset(Svy.design, smoke.ord!="Moderate")

MH<-subset(Svy.design, smoke.ord!="Never")

svyranktest(TOTALEXP~smoke.ord,design=NM,test="wilcoxon") #Never vs. Moderate，wilcoxon检验

svyranktest(TOTALEXP~smoke.ord,design=NH,test="vanderWaerde") #Never vs. Heavy，van der Waerden正态计分检验

svyranktest(TOTALEXP~smoke.ord,design=MH,test="median") #Moderate vs. Heavy，median中位数检验

事后两两比较结果显示：中度吸烟和不吸烟（P=0.002<0.0167）、重度吸烟与不吸烟（P<0.001<0.0167）患者的医疗支出差异具有明显的统计学意义，而中度吸烟和重度吸烟的医疗支出差异无显著的统计学意义（P=0.079>0.0167）。

在nmes数据集中，还有两个变量：lc5和chd5。LC包括肺癌、喉癌和慢性阻塞性肺疾病（COPD），吸烟是这些疾病的主要病因；CHD包括冠心病、卒中及一些吸烟为重要诱因的癌症。如果以这两个变量为结局，结局变量就是二分类变量了，在加权平衡协变量后，可以考虑使用svytable函数进行卡方分析，或者svyglm函数进行广义线性回归模型来完成最后的分析，感兴趣的不妨尝试做一下分析。

— END —