倾向性评分匹配（propensity score matching, PSM）主要是在随机对照试验(Randomized controlled trials,RCT)中用于衡量treat组和control组样本的其他各项特征（如年龄、体重、身高、人种等）的整体均衡性的度量。比如说研究一种药物对疾病的影响，在临床实验中，treat组和control组除了使用药物（安慰剂）不同外，其他的临床特征（如年龄、体重等）都应该基本是相似的，这样treat和control组才有可比性，进而才能验证药物的有效性。

如下图所示，该治疗方法实际上是无效的，但是由于分组中年龄的不平衡导致得出错误的结论。

对于衡量同一特征的组间差异或者距离，我们通常使用标准均值误差（Standardized mean difference，SMD，PMC3144483）。

对于连续性特征变量公式如下：

对于离散型特征变量，公式如下：

下面我们使用tableone包来统计每个变量的标准均值误差SMD,数据来源是右心导管插入(right heart catheterization, RHC)数据，treat组是"RHC"，control是"No RHC"（变量 swang1）

library(tableone) ## PS matching library(Matching) ## Weighted analysislibrary(survey) library(reshape2) library(ggplot2)## Right heart cath dataset rhc <- read.csv("http://biostat.mc.vanderbilt.edu/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") ## Construct a table tabUnmatched <- CreateTableOne(vars = vars, strata = "swang1", data = rhc, test = FALSE) ## Show table with SMD print(tabUnmatched, smd = TRUE)

下图是统计结果，第一列是协变量，第二列是按照有无RHC（treat、contol组）各变量的统计值（mean和SD），最后一列是SMD，可以看出treat和control组的age和sex差异都较小（<10%），income较高（>10%）。

接下来计算通过logit回归计算每个样本的倾向性分数（Propensity Score, PS），也就是被分配为RHC的概率.

## Fit modelpsModel <- glm(formula = 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) ## Predicted probability of being assigned to RHC rhc$pRhc <- predict(psModel, type = "response") ## Predicted probability of being assigned to no RHC rhc$pNoRhc <- 1 - rhc$pRhc ## Predicted probability of being assigned to the treatment actually assigned (either RHC or no RHC) rhc$pAssign <- NA rhc$pAssign[rhc$swang1 == 1] <- rhc$pRhc[rhc$swang1 == 1] rhc$pAssign[rhc$swang1 == 0] <- rhc$pNoRhc[rhc$swang1 == 0] ## Smaller of pRhc vs pNoRhc for matching weight rhc$pMin <- pmin(rhc$pRhc, rhc$pNoRhc)

然后使用Matching包进行匹配一致性样本(1:1匹配)。

listMatch <- Match(Tr = (rhc$swang1 ==1), # Need to be in 0,1 ## logit of PS,i.e., log(PS/(1-PS)) as matching scale X = log(rhc$pRhc / rhc$pNoRhc), ## 1:1 matching M = 1, ## caliper = 0.2 * SD(logit(PS)) caliper = 0.2, replace = FALSE, ties = TRUE, version = "fast") # determining if balance exists in any unmatched dataset and in matched datasetsmb <- MatchBalance(psModel$formula, data=rhc, match.out=listMatch, nboots=50)

将匹配的样本提取出来：

rhcMatched <- rhc[unlist(listMatch[c("index.treated","index.control")]), ]

再看下现在匹配后的SMD，现在所有变量的SMD都小于10%了。

## Construct a tabletabMatched <- CreateTableOne(vars = vars, strata ="swang1", data = rhcMatched, test = FALSE) ## Show table with SMD print(tabMatched, smd = TRUE)

然后给样本进行加权，使得各组中的倾向性评分基本一致，进而消除混杂因素，作为标准平衡数据参考。一般有两种加权方法：逆概率处理加权法（the inverse probability of treatment weighting，IPTW）和标准化死亡比加权法（the standardized mortality ratio weighting，SMRW）,本次我们是有IPTW的进阶版（PMID:26238958）:

## Matching weight rhc$mw <- rhc$pMin / rhc$pAssign # IPTW: rhc$mw1=ifelse(rhc$swang1==1,1/(rhc$pRhc),1/(1-rhc$pRhc)) ## Weighted data rhcSvy <- svydesign(ids = ~ 1, data = rhc, weights = ~ mw) ## Construct a table (This is a bit slow.) tabWeighted <- svyCreateTableOne(vars = vars, strata = "swang1", data = rhcSvy, test =FALSE) ## Show table with SMD print(tabWeighted, smd = TRUE)

加权后变量组间差异（很小）：

进行作图比较 Unmatched、Mathced和Weighted结果：

library(data.table) ## Construct a data frame containing variable name and SMD from all methodsdataPlot <- data.table(variable = rownames(ExtractSmd(tabUnmatched)), Unmatched = ExtractSmd(tabUnmatched), Matched = ExtractSmd(tabMatched), Weighted = ExtractSmd(tabWeighted)) colnames(dataPlot) <- c("variable","Unmatched","Matched","Weighted") ## Create long-format data for ggplot2 dataPlotMelt <- melt(data = dataPlot, id.vars = c("variable"), variable.name = "Method", value.name = "SMD") ## Order variable names by magnitude of SMD varNames <- as.character(dataPlot$variable)[order(dataPlot$Unmatched)]## Order factor levels in the same orderdataPlotMelt$variable <- factor(dataPlotMelt$variable, levels = varNames)## Plot using ggplot2ggplot(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())

可以看出加权后的"标准数据"和我们PSM后的结果基本是一致的。最后还看看右心导管插不插对于生存是否有影响，使用ShowRegTable函数计算风险比（hazard ratio,HR）[95% CI)]和pvalue。

## Unmatched model (unadjsuted)glmUnmatched <- glm(formula = (death =="Yes") ~ swang1, family = binomial(link = "logit"), data = rhc)## Matched model glmMatched <- glm(formula = (death == "Yes") ~ swang1, family = binomial(link ="logit"), data = rhcMatched) ## Weighted model glmWeighted <- svyglm(formula = (death == "Yes") ~ swang1, family = binomial(link ="logit"), design = rhcSvy) ## Show results togetherresTogether <-list(Unmatched = ShowRegTable(glmUnmatched, printToggle = FALSE), Matched = ShowRegTable(glmMatched, printToggle =FALSE), Weighted = ShowRegTable(glmWeighted, printToggle = FALSE)) print(resTogether, quote =FALSE)

。

本文转载自生信杂谈，请支持原创

小结

感谢大家耐心看完，自己的文章都写的很细，代码都在原文中，希望大家都可以自己做一做，也欢迎大家的意见和建议。

如果你是一个大学本科生或研究生，如果你正在因为你的统计作业、数据分析、论文、报告、考试等发愁，如果你在使用SPSS,R，Python，Mplus, Excel中遇到任何问题，都可以联系我。因为我可以给您提供最好的，最详细和耐心的数据分析服务。

如果你对Z检验，t检验，方差分析，多元方差分析，回归，卡方检验，相关，多水平模型，结构方程模型，中介调节等等统计技巧有任何问题，请私信我，获取最详细和耐心的指导。

If you are a student and you are worried about you statistical #Assignments, #Data #Analysis, #Thesis, #reports, #composing, #Quizzes, Exams.. And if you are facing problem in #SPSS, #R-Programming, #Excel, Mplus, then contact me. Because I could provide you the best services for your Data Analysis.

Are you confused with statistical Techniques like z-test, t-test, ANOVA, MANOVA, Regression, Logistic Regression, Chi-Square, Correlation, Association, SEM, multilevel model, mediation and moderation etc. for your Data Analysis...??

Then Contact Me. I will solve your Problem...

加油吧，打工人！

往期内容

R数据分析：用R语言做meta分析

R数据分析：贝叶斯定理的R语言模拟

R数据分析：著名的“三门问题”的R语言模拟

R数据分析：用R语言做潜类别分析LCA

。。。。

本文转载自生信杂谈，请支持原创