首发，nhanes数据（复杂调查数据）倾向性评分匹配函数（PSM）svypm2发布

写在前面的话，本函数仅适用于nhanes数据（复杂调查数据）的二分类倾向评分匹配，不适用于其他类型，请通读全篇文章，了解倾向评分匹配的优点和不足再决定，电子产品，售出不能退换。
倾向评分匹配（Propensity Score Matching，简称PSM）是一种统计学方法，用于处理观察研究（Observational Study）的数据，在SCI文章中应用非常广泛。在观察研究中，由于种种原因，数据偏差（bias）和混杂变量（confounding variable）较多，倾向评分匹配的方法正是为了减少这些偏差和混杂变量的影响，以便对实验组和对照组进行更合理的比较。
为什么需要做倾向评分匹配？
我们知道RCT的证据力度高，是因为对患者进行了严格的筛选。我们的回顾性研究都是过去的数据，很难像RCT一样进行严格的筛选出两组患者基线相近的基础资料，但我们可以通过倾向评分匹配把回归性的数据进行筛选，把基线资料相近的患者进行匹配，得到近似RCT的效果。
应用场景
 1.基线资料不平
 2.开展病例对照研究病阳性例数较少，如罕见病研究
 3.将众多混杂因素变为一个变量：倾向值
以下为一个实例，没进行匹配前两组患者基线资料相差很大，进行倾向评分匹配后，基线资料近似一致了

目前据我所知，目前尚未有专门的nhanes数据匹配的函数或者R包，应粉丝的要求，开发了svypm2函数，目前只能做2组分类的倾向评分匹配，3组分类的倾向评分匹配的还要等一等。想了解原理的话可以参看这篇文章《R语言3组患者倾向性评分匹配（PSM）》

再开始分析演示之前，我想先引用一位大佬的话

倾向性评分匹配不是万能的，PSM只能部分消除（或者说减弱）暴露因素与各个混杂因素之间的相关性，往往不能做到完全消除相关性。因此，在你协变量很多的时候，往往很难配平衡所有协变量，我们应该理性看待这个方法。给大家多一种选择，仅此而已。

OK，下面咱们正式开始介绍，先导入R包和函数，

library(survey)
library(tableone)
options(survey.lonely.psu = "adjust") 
bc<-read.csv("E:/r/test/subtext.csv",sep=',',header=TRUE)
bc <- na.omit(bc)
str(bc)

我介绍一下数据，SEQN：序列号，RIAGENDR, # 性别， RIDAGEYR, # 年龄，RIDRETH1, # 种族，DMDMARTL, # 婚姻状况，WTINT2YR,WTMEC2YR, # 权重，SDMVPSU, # psu，SDMVSTRA，# strata，LBDGLUSI, #血糖mmol表示，LBDINSI, #胰岛素( pmmol/L)，PHAFSTHR #餐后血糖，LBXGH #糖化血红蛋白，SPXNFEV1, #FEV1：第一秒用力呼气量，SPXNFVC #FVC：用力肺活量，ml（估计肺容量），LBDGLTSI #餐后2小时血糖，factor.FVC是我把肺活量分为了2分类，方便用于测试。

把分类变量转成因子

bc$DMDMARTL<-ifelse(bc$DMDMARTL==1,1,0)
bc$RIAGENDR<-as.factor(bc$RIAGENDR)
bc$RIDRETH1<-as.factor(bc$RIDRETH1)
bc$DMDMARTL<-as.factor(bc$DMDMARTL)
bc$factor.FVC<-as.factor(bc$factor.FVC)

建立抽样调查函数

bcSvy2<- svydesign(ids = ~ SDMVPSU, strata = ~ SDMVSTRA, weights = ~ WTMEC2YR,
                   nest=TRUE,data = bc)

定义一下协变量

cov1<-c("RIAGENDR","RIDAGEYR","LBDGLTSI","RIDRETH1","PHAFSTHR")

进行倾向评分匹配，假设我是以是否结婚这个变量来进行倾向性评分匹配，一句代码就匹配好了

out<-svypm2(data=bc,x="DMDMARTL",y="factor.FVC",covs=cov1,design=bcSvy2)

这里我解释一下参数，data是你的数据，x是你需要进行分组匹配的变量，Y是你研究的结局变量，必须是分二类变量，covs是你要调整的协变量，design这里填入咱们的调查函数

生成的结果再out这里，是一个列表数据文件，fMtchDf_1和fMtchDf_2是分组匹配的结果mbc是合并后的数据，可以看到匹配成功了225例

咱们用tableone包来分析一下匹配前和匹配后的P值和smd，分析过程就不介绍了，直接上代码

allVars <-cov1
fvars<-c("RIAGENDR","RIDRETH1")#分类变量定义为fvars

Svytab2<- svyCreateTableOne(vars = allVars,
                            strata = "DMDMARTL",
                            data =bcSvy2 ,
                            factorVars = fvars)

mbc<-out[["mbc"]]
bcSvy3<- svydesign(ids = ~ SDMVPSU, strata = ~ SDMVSTRA, weights = ~ WTMEC2YR,
                   nest=TRUE,data = mbc)
Svytab3<- svyCreateTableOne(vars = allVars,
                            strata = "DMDMARTL",
                            data =bcSvy3 ,
                            factorVars = fvars)
print(Svytab2,smd=T)
print(Svytab3,smd=T)

我们可以看到，匹配效果并不好，还不如匹配前

这是因为咱们还没有设置卡钳，默认的卡钳是0.5，我们可以把它设置得小一点，我设置成0.2

out<-svypm2(data=bc,x="DMDMARTL",y="factor.FVC",covs=cov1,design=bcSvy2,CALIP=0.2) 

allVars <-cov1
fvars<-c("RIAGENDR","RIDRETH1")#分类变量定义为fvars

Svytab2<- svyCreateTableOne(vars = allVars,
                            strata = "DMDMARTL",
                            data =bcSvy2 ,
                            factorVars = fvars)

mbc<-out[["mbc"]]
bcSvy3<- svydesign(ids = ~ SDMVPSU, strata = ~ SDMVSTRA, weights = ~ WTMEC2YR,
                   nest=TRUE,data = mbc)
Svytab3<- svyCreateTableOne(vars = allVars,
                            strata = "DMDMARTL",
                            data =bcSvy3 ,
                            factorVars = fvars)
print(Svytab2,smd=T)
print(Svytab3,smd=T)

咱们可以看到（下图），除了性别，其他得P都大于0.05了

好像这个例子说服力差点，我换个性别来匹配

定义协变量

out<-svypm2(data=bc,x="DMDMARTL",y="factor.FVC",covs=cov1,design=bcSvy2,CALIP=0.1)
生成匹配数据
out<-svypm2(data=bc,x="RIAGENDR",y="factor.FVC",covs=cov1,design=bcSvy2,CALIP=0.2)

查看匹配前和匹配后的基线表

allVars <-cov1
fvars<-c("RIAGENDR","RIDRETH1")#分类变量定义为fvars

Svytab2<- svyCreateTableOne(vars = allVars,
                            strata = "RIAGENDR",
                            data =bcSvy2 ,
                            factorVars = fvars)


mbc<-out[["mbc"]]
bcSvy3<- svydesign(ids = ~ SDMVPSU, strata = ~ SDMVSTRA, weights = ~ WTMEC2YR,
                   nest=TRUE,data = mbc)

Svytab3<- svyCreateTableOne(vars = allVars,
                            strata = "RIAGENDR",
                            data =bcSvy3 ,
                            factorVars = fvars)
print(Svytab2,smd=T)
print(Svytab3,smd=T)

最后咱们来个粉丝的数据，这样有点代表性，特多协变量的，试一下
先导入数据

library(survey)
setwd("E:/r/fensi/20240725P交互")
data <- read.csv("Clean_SUA_SHBG.csv",sep=',',header=TRUE)
data<-as.data.frame(data)

str(data)

把分类变量转成因子

data$gender <- factor(data$gender)
data$race <- factor(data$race)
data$education <- factor(data$education)
data$moderate_activity <- factor(data$moderate_activity)
data$BMI_Class <- factor(data$BMI_Class)

定义协变量

cov1<-c("age","gender","race","education","PIR","BMI","Cholesterol_mg_dL","total_protein_g_dL",
        "Triglycerides_mg_dL","ALT_U_L","albumin_g_dL","ALP_IU_L","AST_U_L","urea_nitrogen_mg_dL",
        "Ca_mg_dL","r_IU_L","LDH_IU_L","P_mg_dL","cotinine_ng_mL")

进行倾向评分匹配，这里协变量非常多，咱们尽量把卡钳设置小点，只有15例能匹配成功

out<-svypm2(data=data,x="gender",y="moderate_activity",covs=cov1,design=nhs,CALIP=0.005)

比较一下匹配后和匹配前的基线表，这里就直接上代码了

allVars <-cov1
fvars<-c("gender","race","education")#分类变量定义为fvars

Svytab2<- svyCreateTableOne(vars = allVars,
                            strata = "gender",
                            data =nhs,
                            factorVars = fvars)

mbc<-out[["mbc"]]

bcSvy3<- svydesign(
  data = mbc,
  ids =  ~ psu,
  strata =  ~ strata,
  weights =  ~ wtmec2yr,
  nest = TRUE,
  survey.lonely.psu = "adjust")

Svytab3<- svyCreateTableOne(vars = allVars,
                            strata = "gender",
                            data =bcSvy3,
                            factorVars = fvars)

我们可以看到大部分变量都被配平了，也算尽力了把。
最后我总结一下，所谓匹配，就是找到相似的数据，因此协变量越多的话你的卡钳应该尽量的小，因此匹配得到的数据也会比卡钳大的数据小。所以要是匹配的效果不理想，增加数据量也是一个途径。

如果你还不会，可以看下面视频操作

首发，nhanes数据倾向评分匹配

需要获得nhanes数据（复杂调查数据）倾向性评分匹配函数（PSM）svypm2函数代码的，请看这篇文章

首发，nhanes数据（复杂调查数据）倾向性评分匹配函数（PSM）svypm2发布