R语言实现基于多重插补的广义模型

数据有缺失，是一件令人不愉快的事情。然而，无论是回顾性研究还是前瞻性的研究，都可能难以避免数据缺失。那么，当你拿到了一份数据，发现其中有个别或一些缺失，你会如何处理呢？目前，可供选择的做法大概有以下两种：1.若缺失较少，且样本量较大，这时选择删除缺失的观测行，也是无伤根本的；2.若样本量不大，缺失比例也不是很小（比如5%-10%的缺失），或者无论样本量有多大，你就是不想删除这5%-10%的缺失，那么可以根据数据情况选择均值/中位数填补、K-mean填补、回归分析填补、决策树填补、随机森林填补和多重插补。上述缺失值补充的方法种的多重插补 (Multivariate Imputation by Chained Equations, MICE)相比其他方法要更加新颖和可靠，想必大家也都听说过。关于MICE的原理，不是本文的重点，这里不多做说明；本文的重点是在多重插补后的数据集上进行广义线性模型、广义相加模型、广义估计方程等广义模型的分析。

在正式进行代码演示之前，先简单介绍一下多重插补及后续分析的主要步骤和本次实践所用的示例数据。如图1（参考文献原图）所示，第一步：在原数据集的基础上进行多次填补，形成多个填补后的数据集；第二步：在每个填补后的数据集上均进行一次分析；第三步：将所有的分析结果进行合并，从而得到一个最终的结果。本次所用数据的情况见表1，是研究H病患者再入院的影响因素，readmission为因变量。

下面开始以代码配合注释的方式，对整个过程进行展示。

####---------- 数据准备阶段 ----------####

#### 载入数据

dat<-read.csv("30_day_readmissions.csv",

sep = ",")

View(dat)

在原始数据集导入后，通过lapply函数和factor函数批量进行定性变量的因子转换，即将数据类型从默认的numeric转换为factor。

#### 将分类变量定义为factor

catvars<-c("ART","Comorbidity")

dat[catvars]<-lapply(dat[catvars], factor)

查看原数据集各变量的缺失情况，以缺失比例和缺失模式图进行展示。从图2可知comorbidity、CD4和Age的缺失比例分别为1.319%、1.237%和0.495%。

#### 查看缺失值

library(VIM)

aggr(dat,plot = T, sortVars=T, only.miss=T)

为减小模拟误差，我们选择生成50个填补后的数据集（m=50），当然你也可以选择更多；根据缺失变量的类型，选择不同的插补方法（defaultMethod）：定量变量为“pmm”，二分类变量为“logreg”，多分类变量为“polyreg”；为保持每次复现时填补的数据集一致，特设置了种子数（seed=316）。

#### 缺失值处理

## 多重插补

library(mice)

dat_mice<-mice(dat, m=50,defaultMethod = c("pmm","pmm","polyreg"),seed = 316)

dat_mice

# 查看补充的数据

dat_mice$imp$CD4

完成数据缺失值的填补不是数据分析的结束，而是分析的开始，下面开始在插补后的50个数据集上进行回归分析。首先是使用with函数，一次性完成50个数据集的分析，接着再利用pool函数进行结果合并。在这里要特别说明一下，之前有朋友来问，可不可以用其中的一个数据集来作分析。我的回答是，多重插补的亮点之一就是生成多个数据集→利用with函数在每个数据集上均进行分析→用pool函数合并分析结果，如果选择其中一个就和多重插补的初衷相悖了。

######---------- 数据分析阶段----------######

###### 广义线性模型 ######

#### 在填补的多个数据集上进行glm

fit_glm<-with(dat_mice,glm(readmission~Age+CD4+ART+Comorbidity,family = binomial))

#### 合并结果

glm_pool<-pool(fit_glm)

#### 输出结果

summary(glm_pool)

###### 广义相加模型 ######

library(mgcv)

### 在填补的多个数据集上进行gam

fit_gam<-with(dat_mice,gam(readmission~Age+CD4+ART+Comorbidity,family = binomial))

#### 合并结果

gam_pool<-pool(fit_gam)

#### 输出结果

summary(gam_pool)

###### 广义估计方程 ######

#### 加载geepack

library(geepack)

#### 在填补的多个数据集上进行gee

fit_gee<-with(dat_mice,geeglm(readmission~Age+CD4+ART+Comorbidity, family=binomial, id=ID, corstr="independence"))

#### 合并结果

gee_pool<-pool(fit_gee)

#### 输出结果

summary(gee_pool)

最后，写一个自定义函数，用于计算OR值及其95%置信区间，并输出β、S.E和P值。

#### 建立自定义函数

poolTable<-function(fit.pool){

 Table<-summary(glm_pool)

 OR<-round(exp(Table$estimate),3)

 LCI<-round(exp(Table$estimate-1.96*Table$std.error),3)

 UCI<-round(exp(Table$estimate+1.96*Table$std.error),3)

 β<-round(Table$estimate,3)

 S.E<-round(Table$std.error,3)

 P<-round(Table$p.value,4)

 library(reshape)

 library(tidyr)

 Table<-rename(Table,c(term="Variable"))

 Table<-cbind(Table[1],β,S.E,OR,LCI,UCI,P)

 Table$k1<-"(";Table$k2<-"-";Table$k3<-")"

 Table<-Table[,c("Variable","β","S.E","OR","k1","LCI","k2","UCI","k3","P")]

 Table<-unite(Table,"OR(95%CI)",c(OR,k1,LCI,k2,UCI,k3),sep = "",remove = T)

 return(Table)

}

## 使用自定义函数

poolTable(gee_pool)