匹配与倾向评分，地表最强实操
本章思维导图：
前言
上一章我们系统介绍了PS的来龙去脉与匹配的血缘关系。纸上得来终觉浅得知此事要躬行。本章节，笔者进行SPSS与R全方位实操。SPSS目的是为了同学们学习成本最低解决80%PSM问题，R目的是为了解决100%遇到的PS各类问题。我们采用SPSS 26.0和R4.1一步步带大家实战并配以讲解，让大家依葫芦画瓢的同时加深对前文理论的理解。本章节提供SPSS 的PS分析（主要是PSM）操作以及R语言的配套操作。对于经典的PSM采用SPSS足够方便，R则是更全面但是稍微需要理解代码含义。
3.1 SPSS的3个PS模块
SPSS每年升级，以往补丁形式的Essentials都被集成中心中更加方便。在SPSS26.0实现PSM分析有三种方法：.SPS宏、python插件、R插件。
方法① .SPS宏[25]，优点：扩展性，可结合任意数据挖掘等model新方法生成Pre。这个宏是基于贪婪算法的1:1匹配。如果想进行1：N 匹配则在抽取完首批对照后，通过横向合并变量的方法从原始文件中去除已经配比的对照，重复宏命令即可，此外，要把工作文件夹中中的结果和过程文件移走，以便产生新的配比文件，也可以用程序第1行重新定义一个工作夹。因为手工作业太多，而且结果与方法②③不一致，所以笔者不推荐。
方法② python插件。SPSS安装时候勾选python，后期扩展-扩展中心-PSM单独可以安装，模块安装好后出现在“数据-倾向评分匹配和个案控制匹配”，本插件最大的优势可以进行个体匹配也可以PSM，δ是绝对值，例如0.02或0.03。运行时候有个细节需要注意，治疗组必须赋值是1对照组是0，治疗组不能是2或0，也就是X=1,0（治疗1，对照0）不可以X01,X21,X12。python插件如图13。
图13 python插件有两种
方法③ R插件。登录网站https://sourceforge.net/projects/psmspss/files/ ，下载 psmatching3.04.spe，安装扩展-本地扩展束。这里要说一下SPSS26.0需要安装R3.4，而旧版需要安装Ressentials及对应的R版本，以及所需的安装R包'SparseM', 'xtable', 'svd', 'abind', 'Rcpp', 'RItools', 'digest'。如果高版本遇到扩展束出错，可以在低版本SPSS中打开扩展束，复制对话框构建程序里面的语法模板内容到新版SPSS中更新.spe文件即可。这个插件需要是covx缺失值必须是填补过的，有missing数据跑出来会错误。此插件提供最邻近匹配、完全匹配、最佳匹配三种算法（后面两种有bugs），pre产生模型也有logistic和GAM两种，以及相对卡钳值匹配的选择。
3.2 SPSS实战
笔者这里引用https://biostat.app.vumc.org/wiki/Main/DataSets → rhc.csv数据集。rhc.csv 数据集[26] 采集1996年就使用了倾向评分方法证明右心导管手术在重症患者初始护理中的有效性，论文PMID: 8782638，rhc论文中采用PSM匹配因素CovX由7名重症监护医师判定，包括age、sex、race…amihx等等。主要结局变量为death，手术方式变量为swang1。rhc进入PSM过程分析变量精简如下表2。Python和R插件均不支持字符型，19个字符型变量（sex race income ninsclas cat1 dnr1 ca resp card neuro gastr renal meta hema seps trauma ortho swang1 death）需要转换为分类变量。
表2 rhc.csv数据进行PSM分析变量一览
3.2.1 SPSS倾向评分匹配
论文原文中描述，2184例rhc患者，对照3551例，进行pre绝对值δ=0.03的最邻近匹配，匹配后结果是1008对。方法②：我们尝试python插件（绝对值δ=0.03）与R插件（δ=0.2S）同时测试，如下图。Python插件fuzzy match匹配1750对，治疗组匹配率81.4%, 仔细看输出文档会发现python插件没有考虑哑变量问题，这在无序分类变量比例高的时候会导致一定局限性，Pre计算错误。
图14 python插件运行1:1 PSM。注：绝对δ
图15 python插件模糊匹配
方法③：先将swang1_01转换为连续变量。与python需要输入绝对值 δ=|pre|不同，R插件中卡钳值是隐藏的，默认δ=0.2*S(logitpre) ，pre评分进行logitP变换后的0.2倍标准差，是不是感觉像隐藏关卡？在插件中只需要输入0.2系数SPSS自动会计算δ值。为了更好的让大家理解，我们运行语法COMPUTE logitp=LG10((1-ps)/ps)，发现原始数据logitp的S=0.609，绝对值δ=0.2S=0.122，。matched_cases_in_control=1为匹配成功。R插件匹配1687对，与python匹配不一致的原因，一是δ不一样，二是因为pre产生模型R插件更优秀考虑了哑变量形式，所以不能简单的将相对系数在两个插件中直接换算。由于python插件没有考虑哑变量形式，我们可以发现python_ps与R得到的ps不一致，所以只推荐用来做个体匹配，PSM需要用R插件才能得出正确pre值！
图16 R插件运行1:1~1:N参数。注：match ratio栏设置比例。
图17 R插件1:1PSM结果
3.2.2 PS匹配效果评价：效应量
R插件匹配还能输出均衡性检验，如下图11。标准化差值法直观反映匹配前后的组间差异。这个方法分别计算各个变量匹配前和匹配后的标准化均数差(standardised mean differences,SMD或cohen’s d)，来判断是否满足均衡的要求。其实cohen’s d是effect size效应量的一种，effect size大约有15种，感兴趣的小伙伴可以自行google相关知识。“效应量”用于解决“大样本数据即使是微不足道的关联因素也可能得出P<0.05”这样的质疑。通常cohen’s d<0.2代表差异微小“trivial”，cohen’s d 0.2~0.5 代表有差异很小但是还算有意义（meaningful but small），cohen’s d 0.5-0.8代表差异中等（medium），cohen’s d >0.8代表差异很大，公式如下图18。一般来说PS分析中cohen’s d 最严格时取0.1算匹配后均衡，SMD 超过 0.25 是均衡性不佳。
图18 SMD (cohen’d) 公式
也有学者认为PS还要看方差比VR， ，接近1才能表示covx组间均衡，VR<0.5或者>2认为不均衡。SPSS可输出PS匹配前后密度图（图19a）、点图（图19b）、线图（图19c）。显示两组有较大的共同支持域，线图显示匹配后cohen’s d差异减小；还可输出jitter抖动图；输出“overalll balance test”的imbalance χ2 test[27] 和“Relative multivariate imbalance L1”中，L1统计量理论上介于0~1之间，匹配前和匹配后相比L1 measure统计量越小。overalll balance test P>0.999；本例L1=1，此指标失效原因是当covx过多时，L1会变成1。所以统计量仅仅是个参考，看图为准。点图中匹配后没有发现变量|d| > 0.25，提示匹配效果越好。采用R插件实施1:1~1:N PSM极其便捷，但遗憾的是尚无分层、加权PS以及其他方法。还需要注意，尽管SPSS对中文比较友好，但是PS插件需要英文，所以大家处理的时候尽量都使用英文变量名。
图19 a.匹配前后pre分布密度 b.校正多因素后cohen’s d值匹配前后分布点图 c.校正多因素后cohen’s d值升降趋势
3.2.3 拓展：SPSS个体匹配
若我们拓展一下，假设此队列第二研究终点是rhc术后急性肾损伤（AKI），想要证明接受rhc手术的糖尿病患者比非糖尿病患者可能更容易发生AKI。我们前期收集了一些数据和大部分血生化标本，但是由于纳入排除标准和采样缺陷导致并非全部病例及对照都采样。此时，我们可以采用post-hoc individual match后期弥补成nested case-control研究，将最终纳入分析的case进行1:1匹配。例如设置匹配条件是患者年龄相差δ为±5岁，eGFR± 5 ml/min/1.73 m2，操作参数设置如图20，连续变量age δ=5，eGFR δ=5，如果有分类变量δ=0，随机seed设置1234是为了能重现结果。假如设置age更严格相差3年，则匹配上的样本量会变少。个体匹配完成后变量hege会出现“.00”“. ”，就是未匹配上的变量，hege≥1才是匹配的。如果想要完成1：N多个对照匹配，需要删减首次hege≥1的对象后再跑一次。
图20 python插件运行个体匹配
3.3 R实战
话说，R最牛B的地方就在于可以踩在巨人肩膀上进行开源二次开发，笔者踩在大佬们的肩膀耕耘2个月，也希望诸位踩在笔者肩膀，用起本程序的时候不忘多来公众号看看笔者。
下面程序需求一定R语言基础。有了前面SPSS抛砖引玉，笔者对rhc.csv数据进行PS进一步大刀阔斧地从各种角度验证分析，PSM主要采用R MatchIt包实现，结果展示依据tableone包。MatchIt包中函数matchit，与SPSS插件不同其caliper采用的是Spooled，见SMD公式 。MatchIt包匹配method：有最相近匹配法(method="nearest")，其他mating方法还有method="exact", method="full", method="optimal" (optimal matching), method="subclass"(subclassification), method="genetic" (genetic matching), method="cem" (coarsened exact matching)。distance：默认采用的计算倾向值采用的是logistic回归(distance="logit")，如想采用其他链接函数如Probit回归来计算倾向值，可添加distance="probit"。当然PSM实现还有许多包，nonrandom(1729对)，Matching(1563对)，结果迥异，感兴趣同学可以自行尝试。倾向评分加权PSW，主要采用自编程序实现。因为目前CRAN中加权的包有很多，其中PSweight包，涵盖了所有加权类型，而且可以做多水平治疗组的PSW，也可以做PSW方法的修剪trimming，但是缺点是没有中间步骤，只有协变量的SMD输出；twang包可以做gbm等机器学习的PSW，但是没有OW等新方法，所以考虑后期进行加权后logistic以及图表输出，还是自编程序更方便。
表3 CRAN上已有的PSW相关包
对PS来说，方法学开发远不止如此，比如机器学习计算PS评分再进行后续匹配、PS如何应用在多次治疗、剂量治疗、结构方程、时序性治疗、多水平数据分析上。顿时大家是不是觉得PS简直是个无底洞？不用怕，相信大家经过本节学习已经打通PS任督二脉，加上google搜索一切都不是问题。对于匹配来说，匹不匹配，匹配参数多少，匹配方法选什么，其实都是可以很灵活可调节的，大家可以在论文中做各种敏感性分析。所以针对匹配的数据可以计算一个结果，针对不匹配的数据也可以计算一个结果，针对不同匹配方法的结果进行比较。归根结底，我们的目的是通过这些方法，发现真实的规律，只有敏感性结果一致才好下稳定结论。不同PS方法笔者根据自己经验给出了推荐等级。
最后，就是全文的压轴部分，奉上笔者呕心沥血总结的R代码，大家可以follow一步步重现，有助于理解前面的概念与原理。
#安装并加载包
install.packages("pacman")
pacman::p_load(dplyr,MatchIt,cobalt,gmodels,ggplot2,survival,survMisc,survminer,knitr,kableExtra, tidyverse, tableone,survey,reshape2)
select <- dplyr::select
#library(haven)
#rhc <- read_sav("D:/temp/rhc.sav")
#读者可自行前去网页下载其他格式
rhc <- read.csv("https://biostat.app.vumc.org/wiki/pub/Main/DataSets/rhc.csv")
vars <- c("age","sex","race","edu","income","ninsclas","cat1","das2d3pc","dnr1",
"ca","surv2md1","aps1","scoma1","wtkilo1","temp1","meanbp1","resp1",
"hrt1","pafi1","paco21","ph1","wblc1","hema1","sod1","pot1","crea1",
"bili1","alb1","resp","card","neuro","gastr","renal","meta","hema",
"seps","trauma","ortho","cardiohx","chfhx","dementhx","psychhx",
"chrpulhx","renalhx","liverhx","gibledhx","malighx","immunhx",
"transhx","amihx")
#加入X变量，swang1
Varsx <- c("age","sex","race","edu","income","ninsclas","cat1","das2d3pc","dnr1",
"ca","surv2md1","aps1","scoma1","wtkilo1","temp1","meanbp1","resp1",
"hrt1","pafi1","paco21","ph1","wblc1","hema1","sod1","pot1","crea1",
"bili1","alb1","resp","card","neuro","gastr","renal","meta","hema",
"seps","trauma","ortho","cardiohx","chfhx","dementhx","psychhx",
"chrpulhx","renalhx","liverhx","gibledhx","malighx","immunhx",
"transhx","amihx","swang1")
library(tableone)
#连续变量
conVars <- c("age","edu","das2d3pc","surv2md1","aps1","scoma1","wtkilo1","temp1","meanbp1",
"resp1","hrt1","pafi1","paco21","ph1","wblc1","hema1","sod1","pot1","crea1","bili1",
"alb1","cardiohx","chfhx","dementhx","psychhx","chrpulhx","renalhx","liverhx","gibledhx",
"malighx","immunhx","transhx","amihx")
#分类变量→因子
catVars <- c("sex","race","income","ninsclas","cat1","dnr1","ca","resp","card",
"neuro","gastr","renal","meta","hema","seps","trauma","ortho","swang1","death")
rhc[catVars] <- lapply(rhc[catVars], factor)
#匹配前一览表，匹配前标准化均数差异
tabjixian <- CreateTableOne(vars = Varsx, data = rhc, test = FALSE)
summary(tabjixian)
print(tabjixian, showAllLevels = TRUE, formatOptions = list(big.mark = ","))
#tableone使用正态分布方法分析资料，因此会出现“(mean (sd))”的描述，两组计量资料用oneway.test函数，t检验分析，
#分类资料用 chisq.test函数，卡方检验分析，默认有矫正卡方，精确性检验用fisher.test函数Fisher检验
#正态性检验,按照X治疗分组，rhc数据集,normtest.p>0.05正态
library(pastecs)
by(rhc[conVars],rhc$swang1,function(x)stat.desc(x,norm=TRUE))
#conVars清一色偏态分布，基线total描述中位数±IQR
tabnonnor <- conVars
print(tabjixian, nonnormal= tabnonnor, formatOptions = list(big.mark = ","))
#table1 基线汇总，分group查看分布，中位数±IQR，精确卡方p
tabUnmatched <- CreateTableOne(vars = Varsx, strata = "swang1" , data = rhc,factorVars = catVars)
print(tabUnmatched, nonnormal = tabnonnor, exact = catVars,smd = TRUE,formatOptions = list(big.mark = ","))
图21 tableone命令获取基线表。注：可见age等基线变量极度不均衡。
GetConfInt <- function(obj) {
logitsticModel <- FALSE
if (identical(class(obj), c("glm", "lm")) == TRUE) {
mat <- coef(summary(obj))
logitsticModel <- TRUE }
else if (identical(class(obj), c("geeglm", "gee", "glm")) == TRUE) {
mat <- coef(summary(obj)) }
else if (identical(class(obj), c("coeftest")) == TRUE) {
mat <- obj }
else if (identical(class(obj), c("matrix")) == TRUE) {
mat <- obj }
else {
stop("Not a supported object")
}
#OR点估计，1.96 * SE， LL UL 95CI%，p-value
matRes <- mat[, 1, drop = F]
matSe <- mat[, 2, drop = F] * qnorm(0.975)
matRes <- cbind(matRes, (matRes - matSe), (matRes + matSe))
colnames(matRes) <- c("OR","lower","upper")
matRes <- exp(matRes)
matRes <- cbind(matRes, mat[, 3:4, drop = F])
if (logitsticModel == TRUE) {
matRes[, c("lower","upper")] <- exp(suppressMessages(confint(obj)))
}
matRes
}
#单因素logistic+ OR95%CI, death ~ rhc，手术影响死亡
glmCrude <- glm(formula=death~swang1,family= binomial(link = "logit"),data=rhc)
GetConfInt(glmCrude)
#多因素full logistic， death ~ CovX +rhc
glmFull <- glm(formula=death ~ swang1+ age + sex + race + edu + income + ninsclas + cat1 + das2d3pc + dnr1 + ca + surv2md1 + aps1 + scoma1 + wtkilo1 + temp1 + meanbp1 + resp1 + hrt1 + pafi1 + paco21 + ph1 + wblc1 + hema1 + sod1 + pot1 + crea1 + bili1 + alb1 + resp + card + neuro + gastr + renal + meta + hema + seps + trauma + ortho + cardiohx + chfhx + dementhx + psychhx + chrpulhx + renalhx + liverhx + gibledhx + malighx + immunhx + transhx + amihx,family= binomial(link = "logit"),data=rhc)
library(MASS)
stepAIC(glmFull, direction="backward")
#② 常规多因素stepwise后退法 logistic， death ~ CovX +rhc
glmstepwise <-glm(formula=death ~ swang1 + age + sex + income + ninsclas + cat1 + das2d3pc + dnr1 + ca + surv2md1 + aps1 + wtkilo1 + temp1 + hrt1 + pafi1 + wblc1 + hema1 + bili1 + resp + card + neuro + gastr + hema + seps + cardiohx + chfhx + dementhx + chrpulhx + renalhx + liverhx + immunhx, family = binomial(link = "logit"), data = rhc)
GetConfInt(glmstepwise)
图22 常规多因素分析。注：校正多因素后，手术影响死亡P <0.001。
#-------------logistic rhc~ covx 生成pre值, 前提需要哑变量，psform1与上面formula等同
match_vars1 <- colnames(select(rhc, vars, -swang1))
psform1 <- f.build("swang1", match_vars1)
#必须考虑了dummy哑变量，才pre正确
psmodel1 <- glm(psform1, data = rhc, family = binomial())
summary(psmodel1)
#预测pre值
rhc$p1 <- predict(psmodel1, newdata =rhc, type = "response")
#实际原理是logitPS的δ值匹配
rhc$ps1 <- log((1 - rhc$p1) / (rhc$p1))
## PS ROC达到0.80,PSM 效果好
library(pROC)
rocPsModel <- roc(swang1 ~ p1, data = rhc)
plot(rocPsModel,print.auc=TRUE,plot=TRUE,print.thres=TRUE)
图23 rhc~50个covx的ROC。注：AUC>0.65 PSM效果才好。
#linear PS adjustment , rhc P<0.001
glmPsAdjLinear <- glm(formula = death ~ swang1 + p1, family = binomial(link = "logit"), data= rhc)
GetConfInt(glmPsAdjLinear)
图24 PS直接校正logistic
#10分位probs = seq(from = 0.1, to = 0.9, by = 0.1))
quintilePoints <- quantile(x = rhc$p1, probs = seq(from = 0.2, to = 0.8, by = 0.2))
quintilePoints <- c(0,quintilePoints,1)
rhc$psQunintile <- cut(x = rhc$p1, breaks = quintilePoints, labels = paste0("q",1:5))
glmPsAdjQuintile <-glm(formula=death~swang1+psQunintile,family=binomial(link="logit"),data= rhc)
GetConfInt(glmPsAdjQuintile)
#增加交互项后，P=0.2，所以不要随便加交互项
glmPsAdjQuintileInt<-glm(formula=death~swang1+psQunintile+swang1:psQunintile, family=binomial(link="logit"),data= rhc)
GetConfInt(glmPsAdjQuintileInt)
#测试5分位层间rhc的OR值
library(plyr)
logistic.stratified <- dlply(.data = rhc, .variables = "psQunintile",
.fun = function(DF) {
glm(death ~ swang1, data = DF, family = binomial)
})
res.strata.logit <- lapply(logistic.stratified[1:5], function(X){
GetConfInt(X)
})
print( res.strata.logit)
#print(do.call(rbind, res.strata.logit))
图25 按照PS5分位分层logistic。注：P20层手术
xtabs.stratified <- dlply(.data =rhc, .variables="psQunintile",
.fun=function(DF) {
xtabs(~ swang1 + death, data = DF)[2:1, 2:1]
})
xtabs.stratified[1:5]
#死亡分层分布
library(PSAgraphics)
cat.psa(categorical=rhc$death, treatment=rhc$swang1,strata=rhc$psQunintile)
#其他变量，如年龄分层分布
#box.psa(continuous=rhc$age, treatment=rhc$swang1,strata=rhc$psQunintile)
图26 按照PS5分位分层logistic各层频数构成比
#-定义function,平衡性测试，各PS方法比较
get_bal <- function(out)
cobalt::love.plot(out,binary = "std",stats = c("mean.diffs"),threshold = c(.1), var.order = "unadjusted", line = TRUE)
#--------方法1 Nearest neighbor within caliper .2*SD
#默认最邻近匹配法、干预组和对照组1:1进行匹配、卡钳值相对系数设为0.03
set.seed(1234)
PSM11_logit.out <-matchit(psform1, data =rhc, method = "nearest",distance=rhc$ps1, m.order= "random",ratio=1,caliper=0.03,replace=FALSE)
summary(PSM11_logit.out)
#caliper,ATT,绘制平衡图
PSM11_logit.data <- MatchIt::match.data(PSM11_logit.out)
get_bal(PSM11_logit.out)
#logistic1:1匹配后数据,PSM11_tabMatched，后面整和需要改 test = FALSE
PSM11_tabMatched <- CreateTableOne(vars = Varsx, strata = "swang1" , data = PSM11_logit.data, factorVars = catVars)
print(PSM11_tabMatched, nonnormal = tabnonnor, exact = Varsx, smd = TRUE,formatOptions = list(big.mark = ","))
图27 logistic1:1 PSM 协变量平衡图
图28 logistic1:1 PSM匹配后SMD结果。注意：SMD均<0.1，均衡。
match_vars1m <- colnames(select(rhc, vars, ps1, -group))
#match_vars1m <- colnames(select(rhc, vars, -group))
psform1m <- f.build("swang1", match_vars1)
set.seed(1234)
PSM11_maha.out <- matchit(psform1m, data = rhc,distance = "mahalanobis")
summary(PSM11_maha.out)
PSM11_maha.data <- match.data(PSM11_maha.out)
#考虑pre评分后Mahalanobis匹配后数据,PSM11_mahatabMatched
PSM11_mahatabMatched <- CreateTableOne(vars = Varsx, strata = "swang1" , data = PSM11_maha.data, factorVars = catVars)
print(PSM11_mahatabMatched, nonnormal = tabnonnor, exact = Varsx,
smd = TRUE,formatOptions = list(big.mark = ","))
#MachIt包的其他method方法如下
#PSM11_opti.out <- matchit(psform1, data = rhc,method = "optimal",ratio=1)
#summary(PSM11_opti.out)
#PSM11_full.out<-matchit(psform1, data=rhc,method= "full",min.controls=1, max.controls=10,discard="both")
#summary(PSM11_full.out)
#设置治疗权重，采用的都是pre=p1
rhc$ptreat <- rhc$p1
rhc$pNotreat <- 1 - rhc$ptreat
#----①IPTW接受治疗A的患者以1/pre加权，而接受对照B的患者以1/(1 -pre)加权
rhc$iptw <- NA
#swang是字符RHC，一般是1/0; treat因素有3水平及以上同理，1/level
rhc$iptw <- NA
rhc$iptw[rhc$swang1 == "RHC"] <- 1/rhc$p1[rhc$swang1 == "RHC"]
rhc$iptw[rhc$swang1 != "RHC"] <- (1/(1 - rhc$p1))[rhc$swang1 != "RHC"]
#IPSW加权，均衡性检验
library(survey)
Svy_IPTW <- svydesign(ids = ~ 1, data = rhc, weights = ~ iptw)
tabWeightedipsw <- svyCreateTableOne(vars = Varsx, strata = "swang1", data = Svy_IPTW, factorVars = catVars)
print(tabWeightedipsw, smd = TRUE)
#IPTW加权分析，加权logistic分析匹配后rhc
glmiptw <- glm(formula=death~swang1,family=binomial(link="logit"),data=rhc,weights=iptw)
GetConfInt(glmiptw)
#同理cox分析，time为随访时间
#library(survival)
#library(reportReg)
#coxiptw<- coxph(Surv(time,death)~swang1,data=rhc,weights=iptw)
#reportReg(coxiptw)
图29 IPTW加权后特征SMD与加权logistic评估手术效果
rhc$smrw <- NA
rhc$smrw[rhc$swang1 == "RHC"] <- 1
rhc$smrw[rhc$swang1 != "RHC"] <- (rhc$p1/(1 - rhc$p1))[rhc$swang1 != "RHC"]
#smrw均衡性检验
library(survey)
Svy_smrw <- svydesign(ids = ~ 1, data = rhc, weights = ~ smrw)
tabWeightedsmrw <- svyCreateTableOne(vars = Varsx, strata = "swang1", data = Svy_smrw, factorVars = catVars)
print(tabWeightedsmrw, smd = TRUE)
##smrw加权分析，加权logistic分析匹配后rhc
glmsmrw <- glm(formula=death~swang1,family=binomial(link="logit"),data=rhc,weights=smrw)
GetConfInt(glmsmrw)
图30 SMRW加权后特征SMD与加权logistic评估手术效果
#④ matching weights IPTW,#-----设置治疗权重，采用的都是pre=p1
rhc$matchWeightNumerator <- pmin(rhc$p1, 1 - rhc$p1)
rhc$matchWeight <- rhc$matchWeightNumerator*rhc$iptw
#加权
library(survey)
Svy_matchWeight <- svydesign(ids = ~ 1, data = rhc, weights = ~ matchWeight)
tabWeightedmatchWeight <- svyCreateTableOne(vars = Varsx, strata = "swang1", data = Svy_matchWeight, factorVars = catVars)
print(tabWeightedmatchWeight, smd = TRUE)
##mating加权分析，加权logistic分析匹配后rhc
glmmatchWeight <- glm(formula=death~swang1,family=binomial(link="logit"),data=rhc,weights=matchWeight)
GetConfInt(glmmatchWeight)
图31 MW加权后特征SMD与加权logistic评估手术效果。
rhc$ptreat <- rhc$p1
rhc$pNotreat <- 1 - rhc$ptreat
rhc$ow[rhc$swang1 == "RHC"]<- rhc$pNotreat[rhc$swang1 == "RHC"]
rhc$ow[rhc$swang1 == "No RHC"]<- rhc$ptreat[rhc$swang1 == "No RHC"]
#over Weighting 最强加权， tabWeightedOw
Svy_Ow <- svydesign(ids = ~ 1, data = rhc, weights = ~ ow)
tabWeightedOw <- svyCreateTableOne(vars = Varsx, strata = "swang1", data = Svy_Ow, factorVars = catVars)
print(tabWeightedOw, smd = TRUE)
addmargins(table(ExtractSmd(tabWeightedOw) > 0.1))
#ow加权分析，加权logistic分析匹配后rhc
glmowWeight <- glm(formula=death~swang1,family=binomial(link="logit"),data=rhc,weights=ow)
GetConfInt(glmowWeight)
图 32 OW加权后特征SMD与加权logistic评估手术效果。
#"ASD" refers to the pairwise absolute standardized difference
#https://github.com/thuizhou/PSweight, trtgrp=2 两水平01
library(PSweight)
msstat<-SumStat(psform1, data=rhc, weight=c("IPW","overlap","treated","entropy", "matching"))
#输出overlap 重叠直方图
plot(msstat, type="hist")
#loveplot图
plot(msstat, type="balance", weighted.var=TRUE, threshold=0.1, metric="ASD")
summary(msstat)
图33手术组与对照组pre值的overlap重叠直方图。注：poor overlap用OW分析效果好。
tabUnmatched <- CreateTableOne(vars = Varsx, strata = "swang1" , data = rhc)
#match要有命令test = FALSE，不然匹配不上的组不会显示
PSM11_tabMatched <- CreateTableOne(vars = Varsx, strata = "swang1" , data = PSM11_logit.data)
PSM11_mahatabMatched <- CreateTableOne(vars = Varsx, strata = "swang1" , data = PSM11_maha.data)
tabWeightediptw <- svyCreateTableOne(vars = Varsx, strata = "swang1", data = Svy_IPTW)
tabWeightedsmrw<- svyCreateTableOne(vars = Varsx, strata = "swang1", data = Svy_smrw)
tabWeightedmatchWeight<- svyCreateTableOne(vars = Varsx, strata = "swang1", data = Svy_matchWeight)
tabWeightedOw <- svyCreateTableOne(vars = Varsx, strata = "swang1", data = Svy_Ow)
#匹配后table汇总, 方便统计分析
resCombo<-NA
resCombo <- cbind(print(tabUnmatched,printToggle= FALSE),
print(PSM11_tabMatched,printToggle= FALSE),
print(PSM11_mahatabMatched,printToggle= FALSE),
print(tabWeightediptw,printToggle= FALSE),
print(tabWeightedsmrw,printToggle= FALSE),
print(tabWeightedmatchWeight,printToggle= FALSE),
print(tabWeightedOw,printToggle= FALSE))
#重复N次，N方法
resCombo <- rbind(Group = rep(c("No treat","treat"),1), resCombo)
#中间加空格“”
colnames(resCombo) <- c("Unmatch","","P","","logit","","P","","Maha","","P","","IPTW","","P","","SMRW","","P","","MW","","P","","OW","","P","")
print(resCombo, quote = FALSE)
write.csv(resCombo, "resCombo.csv")
图34 7种方法合集。注：resCombo.csv编辑而来
#①举例测试logistic模型，绘制单个匹配前后差异SMD其他方法相同，带连线
get_bal(PSM11_logit.out)
#不带连线黑白图loveplot
love.plot(bal.tab(PSM11_logit.out, m.threshold=0.1),stat = "mean.diffs", var.names =rhc, abs = F)
plot(summary(PSM11_logit.out))
#按照sex分组，Density function密度图、直方图
bal.plot(PSM11_logit.out,var.name = 'sex',which = 'both')
#金字塔图
library(Hmisc)
histbackback(split(rhc$ps1,rhc$swang1), main= "Propensityscorebeforematching",xlab=c("control","treatment"))
#Quantile-quantile (QQ) plot，匹配前后平衡性的改善程度。
plot(PSM11_logit.out,type = "qq", interactive = FALSE,which.xs = c("age","sex"))
#绘制直方图
plot(PSM11_logit.out,type="hist",interactive = FALSE)
#绘制抖动图
plot(PSM11_logit.out,type="jitter",interactive = FALSE)
#绘制loveplot合成图，多种方法差异
dataPlot <- data.frame(variable = rownames(ExtractSmd(tabUnmatched)), Unmatched = as.numeric(ExtractSmd(tabUnmatched)), logitMatched = as.numeric(ExtractSmd(PSM11_tabMatched)), mahaMatched = as.numeric(ExtractSmd(PSM11_mahatabMatched)), Weightediptw = as.numeric(ExtractSmd(tabWeightediptw)), Weightedsmrw = as.numeric(ExtractSmd(tabWeightedsmrw)), Weightedmw = as.numeric(ExtractSmd(tabWeightedmatchWeight)),
WeightedOw = as.numeric(ExtractSmd(tabWeightedOw)))
dataPlotMelt <- melt(data = dataPlot,id.vars=c("variable"),variable.name ="Method",value.name="SMD")
varNames <- as.character(dataPlot$variable)[order(dataPlot$Unmatched)]
dataPlotMelt$variable <- factor(dataPlotMelt$variable,levels = varNames)
ggplot(data = dataPlotMelt, mapping = aes(x = variable, y = SMD, group = Method, color = Method)) + geom_line() +geom_point() + geom_hline(yintercept = 0.1, color = "black", size = 0.1) +coord_flip() + theme_bw() + theme(legend.key = element_blank())
图35 PS展示图合集
#单因素，Unmatched model
glmUnmatched <- glm(formula = death~ swang1,family = binomial(link = "logit"),data= rhc)
#1:1匹配后，条件/非条件都可logistic，logistic Matched model → PSM11_logit.data
glmMatched <- glm(formula = death~ swang1,family = binomial(link = "logit"),data= PSM11_logit.data)
#maha logistic，logistic Matched model → PSM11_logit.data
glmmaha <- glm(formula = death~ swang1,family = binomial(link = "logit"),data= PSM11_maha.data)
#logitMatched<-ShowRegTable(glmmaha)
#汇总7种方法OR值
GetConfInt(glmUnmatched)
GetConfInt(glmMatched)
GetConfInt(glmmaha)
GetConfInt(glmiptw)
GetConfInt(glmsmrw)
GetConfInt(glmmatchWeight)
GetConfInt(glmowWeight)
图36 7种主流PS方法手术效果。各种方法OR稍微不同，P方向一致。
表4 SPSS插件与R包运用PSM匹配结果
小结
最近感慨，统计简直就是个黑洞，越学越不懂。
谢谢你耐心看完本章。
那么问题来了，你懂了吗？
参考文献：
Carignan A, Valiquette L, Grenier C, Musonera JB, Nkengurutse D, Marcil-Heguy A, Vettese K, Marcoux D, Valiquette C, Xiong WT et al: Anosmia and dysgeusia associated with SARS-CoV-2 infection: an age-matched case-control study. CMAJ 2020, 192(26):E702-E707. Cosma S, Carosso AR, Cusato J, Borella F, Carosso M, Bovetti M, Filippini C, D'Avolio A, Ghisetti V, Di Perri G et al: Coronavirus disease 2019 and first-trimester spontaneous abortion: a case-control study of 225 pregnant patients. Am J Obstet Gynecol 2021, 224(4):391 e391-391 e397. Viera AJ: Odds ratios and risk ratios: what's the difference and why does it matter? South Med J 2008, 101(7):730-734. SUNY Downstate Medical Center. (2004). The Evidence Pyramid. Retrieved 1 June. 2021 from https://guides.downstate.edu/c.php?g=856794&p=6831536. Grimes DA: "Case-control" confusion: mislabeled reports in obstetrics and gynecology journals. Obstet Gynecol 2009, 114(6):1284-1286. Ramos A, Mendoza LC, Rabasa F, Bolibar I, Puig T, Corcoy R: Case-control studies in diabetes. Do they really use a case-control design? Acta Diabetol 2017, 54(7):631-634. Schwartz S, Campbell UB, Gatto NM, Gordon K: Toward a clarification of the taxonomy of "bias" in epidemiology textbooks. Epidemiology 2015, 26(2):216-222. Grool AM, Aglipay M, Momoli F, Meehan WP, 3rd, Freedman SB, Yeates KO, Gravel J, Gagnon I, Boutis K, Meeuwisse W et al: Association Between Early Participation in Physical Activity Following Acute Concussion and Persistent Postconcussive Symptoms in Children and Adolescents. JAMA 2016, 316(23):2504-2514. Zeng C, Dubreuil M, LaRochelle MR, Lu N, Wei J, Choi HK, Lei G, Zhang Y: Association of Tramadol With All-Cause Mortality Among Patients With Osteoarthritis. JAMA 2019, 321(10):969-982. Mehta N, Kalra A, Nowacki AS, Anjewierden S, Han Z, Bhat P, Carmona-Rubio AE, Jacob M, Procop GW, Harrington S et al: Association of Use of Angiotensin-Converting Enzyme Inhibitors and Angiotensin II Receptor Blockers With Testing Positive for Coronavirus Disease 2019 (COVID-19). JAMA Cardiol 2020, 5(9):1020-1026. Azoulay L, Yin H, Filion KB, Assayag J, Majdan A, Pollak MN, Suissa S: The use of pioglitazone and the risk of bladder cancer in people with type 2 diabetes: nested case-control study. BMJ 2012, 344:e3645. Lee CC, Lee MT, Chen YS, Lee SH, Chen YS, Chen SC, Chang SC: Risk of Aortic Dissection and Aortic Aneurysm in Patients Taking Oral Fluoroquinolone. JAMA Intern Med 2015, 175(11):1839-1847. Chen HH, Huang N, Chen YM, Chen TJ, Chou P, Lee YL, Chou YJ, Lan JL, Lai KL, Lin CH et al: Association between a history of periodontitis and the risk of rheumatoid arthritis: a nationwide, population-based, case-control study. Ann Rheum Dis 2013, 72(7):1206-1211. Fortner RT, Eliassen AH, Spiegelman D, Willett WC, Barbieri RL, Hankinson SE: Premenopausal endogenous steroid hormones and breast cancer risk: results from the Nurses' Health Study II. Breast Cancer Res 2013, 15(2):R19. Chakraborty S, Sun CL, Francisco L, Sabado M, Li L, Chang KL, Forman S, Bhatia S, Bhatia R: Accelerated telomere shortening precedes development of therapy-related myelodysplasia or acute myelogenous leukemia after autologous transplantation for lymphoma. J Clin Oncol 2009, 27(5):791-798. Pearce N: Analysis of matched case-control studies. BMJ 2016, 352:i969. Delcoigne B, Stoer NC, Reilly M: Valid and efficient subgroup analyses using nested case-control data. Int J Epidemiol 2018, 47(3):841-849. Noah MA, Peek GJ, Finney SJ, Griffiths MJ, Harrison DA, Grieve R, Sadique MZ, Sekhon JS, McAuley DF, Firmin RK et al: Referral to an extracorporeal membrane oxygenation center and mortality among patients with severe 2009 influenza A(H1N1). JAMA 2011, 306(15):1659-1668. Sadique MZ, Grieve R, Harrison DA, Cuthbertson BH, Rowan KM: Is Drotrecogin alfa (activated) for adults with severe sepsis, cost-effective in routine clinical practice? Crit Care 2011, 15(5):R228. Austin PC: An Introduction to Propensity Score Methods for Reducing the Effects of Confounding in Observational Studies. Multivariate Behav Res 2011, 46(3):399-424. Thomas L, Li F, Pencina M: Using Propensity Score Methods to Create Target Populations in Observational Clinical Research. JAMA 2020, 323(5):466-467. Thomas LE, Li F, Pencina MJ: Overlap Weighting: A Propensity Score Method That Mimics Attributes of a Randomized Clinical Trial. JAMA 2020, 323(23):2417-2418. Hemkens LG, Contopoulos-Ioannidis DG, Ioannidis JP: Agreement of treatment effects for mortality from routinely collected data and subsequent randomized trials: meta-epidemiological survey. BMJ 2016, 352:i493. McMurry TL, Hu Y, Blackstone EH, Kozower BD: Propensity scores: Methods, considerations, and applications in the Journal of Thoracic and Cardiovascular Surgery. J Thorac Cardiovasc Surg 2015, 150(1):14-19. 李智文, 李宏田, 张乐,等. 用SPSS宏程序实现观察对象的倾向评分配比[J]. 中国卫生统计, 2011, 28(001):89-90. Connors AF, Jr., Speroff T, Dawson NV, Thomas C, Harrell FE, Jr., Wagner D, Desbiens N, Goldman L, Wu AW, Califf RM et al: The effectiveness of right heart catheterization in the initial care of critically ill patients. SUPPORT Investigators. JAMA 1996, 276(11):889-897. Hansen, B. B., & Bowers, J. (2008). Covariate balance in simple, stratified and clustered comparative studies. Statistical Science, 23, 219-236.
原创不易，转载请说明来自本公众号。