因果推断-倾向性评分匹配

所使用的包
创建数据
合并数据框
样本匹配
保存匹配样本

倾向性评分最大的优势是将多个混杂因素的影响用一个综合的值来表示，即倾向性评分值（Propensity Score, PS），从而降低协变量的维度，因此该方法尤其适用于协变量较多的情况。

倾向得分（Propensity Score）的数学原理是基于条件概率和概率密度函数的统计原理。倾向得分是处理组和对照组之间的处理选择概率，用于降低选择偏差。以下是倾向得分的数学原理：

假设我们有一个观察性研究，其中有两组参与者：处理组（接受了某项干预或治疗）和对照组（没有接受干预）。我们希望通过倾向得分来估计每个个体接受处理的概率。

基本原理：

设 Y 表示观察到的结果变量，D 表示处理接受与否的二进制变量（1 表示接受处理，0 表示不接受处理），X 表示观察到的特征或协变量（如年龄、性别、健康状况等）。

倾向得分（Propensity Score，P(X)）是指在给定协变量 X 的条件下，接受处理 D=1 的概率。它可以表示为 P(X) = P(D=1|X)。

估计倾向得分的方法：

通常，倾向得分可以使用逻辑回归模型估计。逻辑回归模型可以表示为：

P(D=1|X) = 1 / (1 + exp(-βX))

其中，β 是模型的系数，X 是协变量的向量。

模型参数 β 可以通过最大似然估计或其他方法来估计。

倾向得分的应用：

一旦估计出倾向得分 P(X)，可以用它来进行如下的应用：

匹配处理组和对照组：通过匹配处理组和对照组的个体，使它们在倾向得分上更相似，从而降低选择偏差。处理效应估计：倾向得分可以用作权重，以估计处理的因果效应，考虑到处理组和对照组的概率分布差异。倾向得分的优点和限制：

优点：倾向得分可以有效降低因选择偏差而引起的观察性研究中的混杂性，同时允许研究人员更准确地估计处理效应。限制：倾向得分方法依赖于模型的假设，如弃治假设，且需要高质量的数据和正确选择的协变量。错误的模型选择可能导致倾向得分的失效。总之，倾向得分的数学原理是基于条件概率和逻辑回归模型的统计原理，它用于处理选择偏差和进行因果推断，从而提高观察性研究的可靠性。

倾向性评分的一般步骤为：

估计 PS 值；
利用 PS 值均衡协变量分布；
均衡性检验及模型评价；
处理效应估计。

中，PS 值的估计是以处理因素作为因变量，其他混杂因素作为自变量，通过建立一个模型（可以是传统的回归模型，也可以是机器学习方法）来估计每个研究对象接受处理因素的可能性。

目前用于估计 PS 值的方法有 logistic 回归，Probit 回归、神经网络、支持向量机、分类与回归数、Boosting 算法、SuperLearner 等。其中 logistic 回归是目前最常用的方法。

所使用的包

pacman::p_load(knitr, wakefield, MatchIt, tableone, captioner)

创建数据

创建一个名为df.patients的数据框，包含250个病人的年龄和性别数据，所有病人的年龄都要在30-78岁之间，并且70%的病人被设定为男性。

set.seed(180717)
library(wakefield)
df.patients <- r_data_frame(n = 250, 
            age(x = 30:78, name = 'Age'), 
            sex(x = c("Male", "Female"), 
            prob = c(0.70, 0.30), name = "Sex"))

## Warning: `tbl_df()` was deprecated in dplyr 1.0.0.
## ℹ Please use `tibble::as_tibble()` instead.
## ℹ The deprecated feature was likely used in the wakefield package.
##   Please report the issue at <https://github.com/trinker/wakefield/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

df.patients$Sample <- as.factor('Patients')

summary(df.patients)

##       Age            Sex           Sample   
##  Min.   :30.00   Male  :156   Patients:250  
##  1st Qu.:41.00   Female: 94                 
##  Median :54.00                              
##  Mean   :53.61                              
##  3rd Qu.:66.00                              
##  Max.   :78.00

创建另一个名为df.population的数据框。我们希望这个数据集的数据和患者df.patients有些不同，因此正常人群的年龄区间被设定为18-80岁，并且男女各占一半。

set.seed(180717)

df.population <- r_data_frame(n = 1000, age(x = 18:80, name = 'Age'), 
sex(x = c("Male", "Female"),
prob = c(0.50, 0.50), 
name = "Sex"))

df.population$Sample <- as.factor('Population')
 
#下方表格显示样本平均年龄为50.1岁，男女比例也大致相等。

summary(df.population)

##       Age            Sex             Sample    
##  Min.   :18.00   Male  :489   Population:1000  
##  1st Qu.:33.00   Female:511                    
##  Median :48.00                                 
##  Mean   :48.81                                 
##  3rd Qu.:65.00                                 
##  Max.   :80.00

合并数据框

在匹配样本之前，我们需要把两个数据框合并。首先，生成一个新变量Group来代表研究对象来自哪个全体（逻辑型变量），再添加另一个变量Distress来表示个体的痛苦程度。Distress变量是利用wakefield包中的age函数创建的，可以发现，女性承受的痛苦级别更高。

mydata <- rbind(df.patients, df.population)
mydata$Group <- as.logical(mydata$Sample == 'Patients')
mydata$Distress <- ifelse(mydata$Sex == 'Male', 
                age(nrow(mydata), x = 0:42, name = 'Distress'),
                age(nrow(mydata), x = 15:42, name = 'Distress'))

mydata

## # A tibble: 1,250 × 5
##      Age Sex    Sample   Group Distress
##    <int> <fct>  <fct>    <lgl>    <int>
##  1    57 Female Patients TRUE        24
##  2    77 Male   Patients TRUE        30
##  3    54 Male   Patients TRUE        24
##  4    72 Female Patients TRUE        31
##  5    44 Male   Patients TRUE        14
##  6    74 Female Patients TRUE        41
##  7    62 Male   Patients TRUE         6
##  8    46 Female Patients TRUE        25
##  9    44 Male   Patients TRUE        22
## 10    55 Female Patients TRUE        28
## # ℹ 1,240 more rows

当我们比较两组样本的年龄和性别分布时，我们可以发现明显的区别：

library(pacman)
pacman::p_load(tableone)
table1 <- CreateTableOne(vars = c('Age', 'Sex', 'Distress'), 
                         data = mydata, 
                         factorVars = 'Sex', 
                         strata = 'Sample')
table1 <- print(table1, 
                printToggle = FALSE, 
                noSpaces = TRUE)

library(knitr)
kable(table1[,1:3],  
      align = 'c', 
      caption = 'Table 1: Comparison of unmatched samples')

Table 1: Comparison of unmatched samples
	Patients	Population	p
n	250	1000
Age (mean (SD))	53.61 (14.19)	48.81 (18.39)	<0.001
Sex = Female (%)	94 (37.6)	511 (51.1)	<0.001
Distress (mean (SD))	24.35 (11.17)	24.42 (11.23)	0.928

样本匹配

现在，我们已经完成了全部的准备工作，可以开始使用MatchIT包中的matchit函数来匹配两组样本了。

library(MatchIt)
set.seed(180717)

match.it <- matchit(Group ~ Age + Sex, 
data = mydata, 
method="nearest", 
ratio=1)

a <- summary(match.it)

a

## 
## Call:
## matchit(formula = Group ~ Age + Sex, data = mydata, method = "nearest", 
##     ratio = 1)
## 
## Summary of Balance for All Data:
##           Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance         0.2174        0.1957          0.3921     0.7863    0.1125
## Age             53.6120       48.8080          0.3385     0.5958    0.0826
## SexMale          0.6240        0.4890          0.2787          .    0.1350
## SexFemale        0.3760        0.5110         -0.2787          .    0.1350
##           eCDF Max
## distance     0.210
## Age          0.198
## SexMale      0.135
## SexFemale    0.135
## 
## Summary of Balance for Matched Data:
##           Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance         0.2174        0.2173          0.0004     1.0027    0.0003
## Age             53.6120       54.1760         -0.0397     0.9712    0.0101
## SexMale          0.6240        0.6080          0.0330          .    0.0160
## SexFemale        0.3760        0.3920         -0.0330          .    0.0160
##           eCDF Max Std. Pair Dist.
## distance     0.012          0.0012
## Age          0.020          0.0617
## SexMale      0.016          0.0495
## SexFemale    0.016          0.0495
## 
## Sample Sizes:
##           Control Treated
## All          1000     250
## Matched       250     250
## Unmatched     750       0
## Discarded       0       0

函数中method=’nearest’的设定指明了使用近邻法
其他的方法包括:method = c(“exact”, “subclass”, “optimal”, ““genetic”, “full”)）。
ratio=1意味着这是一一配对。在我们的例子里，对于患者组中的每一个病例，将恰好匹配正常人群样本中的一个人。
同时也请注意Group变量需要是逻辑型变量（TRUE与FALSE）。

为了后续工作的便利，我们将summary函数的输出赋值给名为a的变量。在匹配完样本后，正常人群样本量减少到了和患者样本一致

kable(a$nn, digits = 2, align = 'c', caption = 'Table 2: Sample sizes')

Table 2: Sample sizes
	Control	Treated
All (ESS)	1000	250
All	1000	250
Matched (ESS)	250	250
Matched	250	250
Unmatched	750	0
Discarded	0	0

根据输出结果，匹配后的年龄和性别分布基本一致了。

kable(a$sum.matched, digits = 2, align = 'c', caption = 'Table 3: Summary of balance for matched data')

Table 3: Summary of balance for matched data
	Means Treated	Means Control	Std. Mean Diff.	Var. Ratio	eCDF Mean	eCDF Max	Std. Pair Dist.
distance	0.22	0.22	0.00	1.00	0.00	0.01	0.00
Age	53.61	54.18	-0.04	0.97	0.01	0.02	0.06
SexMale	0.62	0.61	0.03	NA	0.02	0.02	0.05
SexFemale	0.38	0.39	-0.03	NA	0.02	0.02	0.05

倾向性评分的分布可以使用MatchIt包中的plot函数进行绘制。

plot(match.it, type = 'jitter', interactive = FALSE)

## Warning in check_font_path(italic, "italic"): 'italic' should be a length-one
## vector, using the first element

保存匹配样本

最后，匹配的样本将保存到名为df.match的新数据框中。

# 获得匹配数据集(含distance和weights列)的第1~ncol(mydata)列

df.match <- match.data(match.it)[1:ncol(mydata)]

# rm()函数：删除变量函数
rm(df.patients, df.population)

# 我们可以对比两类人群间痛苦程度的差异是否依旧显著
pacman::p_load(tableone)
table4 <- CreateTableOne(vars = c('Age', 'Sex', 'Distress'), 
                         data = df.match, 
                         factorVars = 'Sex', 
                         strata = 'Sample')

# printToggle：是否打印输出。如果为FALSE，则不创建输出，并返回一个矩阵。
# noSpaces：是否删除为对齐添加的空格。如果您希望自己在其他软件中对齐数字，请使用此选项。

table4 <- print(table4, 
                printToggle = FALSE, 
                noSpaces = TRUE)

# align/列对齐: 'l' (left), 'c' (center) and/or 'r' (right). 默认则align = NULL, 数值列右对齐, 其他列左对齐.
# kable函数：Create Tables In LaTeX, HTML, Markdown And ReStructuredText

kable(table4[,1:3],  
      align = 'c', 
      caption = 'Table 4: Comparison of matched samples')

Table 4: Comparison of matched samples
	Patients	Population	p
n	250	250
Age (mean (SD))	53.61 (14.19)	54.18 (14.40)	0.659
Sex = Female (%)	94 (37.6)	98 (39.2)	0.783
Distress (mean (SD))	24.35 (11.17)	23.94 (11.50)	0.685