关于批量单因素logistic回归结果的整合并优化输出

1 GLM模型简介

参数 family 规定回归模型的类型

链接函数	Y的分布	回归类型
binomial(link = "logit")	二元离散 二项分布	逻辑回归
gaussian(link = "identity")	连续 正态分布	线性回归
Gamma(link = "inverse")	非负 右偏误差分布	Gamma回归
inverse.gaussian(link = "1/mu^2")	连续的正值 逆高斯分布	逆高斯回归
poisson(link = "log")	非负次数 离散型 泊松分布	泊松回归
quasi(link = "identity", variance = "constant")	连续性 非负性 常数方差	平方根误差回归
quasibinomial(link = "logit")	二元离散 非负性	逻辑回归
quasipoisson(link = "log")	非负次数 离散型	广义线性模型

2 拟合glm模型

我们拟合常规的logistic回归

# 1 转化为因子
dat1$sex <- factor(dat1$sex)
# 2 拟合模型
mod <- glm(insomnia~sex,family = binomial,data=dat1)
# 3 总结模型
summary(mod)

可以知道，拟合的logistic回归模型，相对于Sex赋值为1的人群，Sex赋值为2的人群，OR=EXP(0.68023),P值<0.001

3 抓取回归的结果

接下来我们对其进行组合

为了方便接下来的工作，我们加载dplyr包

# 加载包
library("dplyr")
mod_dat <- mod %>% summary %>% coef %>% as.data.frame() 

此时的结果长这个样子，已经有

Estimate：β回归系数

Std. Error：S.E (标准误)

z value： Z值

Pr(>|z|)：P值

4 OR的计算

很明显，这个结果还不够，下面我们需要计算OR值和95%CI

 mod_dat <- mod %>% summary %>% coef %>% as.data.frame() %>% 
  mutate(Beta=sprintf("%0.2f",Estimate)) %>% 
  mutate(S.E=sprintf("%0.2f",`Std. Error`)) %>% 
  mutate(Z=sprintf("%0.2f",`z value`)) %>% 
  mutate(P=ifelse(`Pr(>|z|)`<0.001,'<.001',sprintf("%0.3f",`Pr(>|z|)`))) %>% 
  mutate(OR=sprintf("%0.2f",exp(Estimate))) %>% 
  mutate(lower=sprintf("%0.2f",exp(Estimate - qnorm(0.975)*`Std. Error`))) %>% 
  mutate(upper=sprintf("%0.2f",exp(Estimate + qnorm(0.975)*`Std. Error`)))
mod_dat

到这里，其实结果就告一段落了，此时的结果为

5 OR (95%CI)的计算

Beta、S.E、Z、P、OR、lower和upper都有了，各位到这里就可以结束了

细心的同学可以发现，我在计算95%CI的时候用的是：

 exp(Estimate - qnorm(0.975)*`Std. Error`)

这里为什么要这样去用呢，当计算的CI很贴近于1，你用1.96*S.E去计算CI，会得到Lower=1.00，更甚至会得到虽然Lower<1,Upper>1,但是P值是显著的，此时你会怀疑自己，或者怀疑软件了，本质就是因为1.96的精度不够。因此，我们采用qnorm(0.975)去计算

qnorm (0.975)= 1.959964，这个和1.96还是有差别的，用1.96去计算会得到更宽的区间，这也就解释了Lower<1,Upper>1,但是明明P值是显著的，这个结果了。

然后，我们想做的更卷一点

想想，如果，OR=1.000001,四舍五入是不是就变成了1.00,然后你又怀疑自己了，OR=1.00，但是P却显著，因此，我们继续做一个转化

mutate(OR=exp(Estimate)) %>% 
  mutate(OR=ifelse(OR<1.01 & OR>1 & `Pr(>|z|)`<0.05,1.01,#进行特殊情况的转化
                   ifelse(OR>0.99 & OR<1 & `Pr(>|z|)`<0.05,0.99,OR)))

怎么理解呢？

当OR<1.01并且OR>1,且P值显著时，那就将OR手动转为1.01,反之

当OR>0.99并且OR<1,且P值显著时，那就将OR手动转为0.99

虽然不太合适，但这样是不是可以解除误会了？

6 结果的组合

同样，我们将Lower和Upper也进行转化，我们上一套完整的代码

 mod_dat <- mod %>% summary %>% coef %>% as.data.frame() %>% 
  mutate(Beta=sprintf("%0.2f",Estimate)) %>% 
  mutate(S.E=sprintf("%0.2f",`Std. Error`)) %>% 
  mutate(Z=sprintf("%0.2f",`z value`)) %>% 
  mutate(P=ifelse(`Pr(>|z|)`<0.001,'<.001',sprintf("%0.3f",`Pr(>|z|)`))) %>% 
  mutate(OR=exp(Estimate)) %>% 
  mutate(OR=ifelse(OR<1.01 & OR>1 & `Pr(>|z|)`<0.05,1.01,#进行特殊情况的转化
                   ifelse(OR>0.99 & OR<1 & `Pr(>|z|)`<0.05,0.99,OR))) %>% 
  mutate(lower=exp(Estimate - qnorm(0.975)*`Std. Error`)) %>% 
  mutate(lower=ifelse(lower>1 & lower<1.01 & `Pr(>|z|)`<0.05,1.01,#进行特殊情况的转化
                      ifelse(lower<1 & lower>0.99 & `Pr(>|z|)`<0.05,0.99,lower))) %>% 
  mutate(upper=exp(Estimate + qnorm(0.975)*`Std. Error`)) %>% 
  mutate(upper=ifelse(upper>1 & upper<1.01 & `Pr(>|z|)`<0.05,1.01,#进行特殊情况的转化
                      ifelse(upper<1 & upper>0.99 & `Pr(>|z|)`<0.05,0.99,upper))) %>% 
  mutate(ORCI=paste0(sprintf("%0.2f",OR)," (",
                     sprintf("%0.2f",lower)," - ",
                     sprintf("%0.2f",upper),")")) %>% 
  dplyr::select(Beta,S.E,Z,P,ORCI,OR,lower,upper) #仅保留的生成的新变量

这个结果就很清爽了

7 结果的优化

虽然清爽，但是我们还想更卷一些，看看哪里还可以卷，我们刚才的解释为：相对于Sex赋值为1的人群，Sex赋值为2的人群，OR=1.97, P值<0.001，很明显，这个结果，我们看不到参照的水平，那就继续优化

1. 生成一行，变量名为sex，但是水平是1，并添加Reference字样的值
2. 表格合并

7. 1参照水平

怎么获取参照水平？R语言默认以第一个水平为参照，因此只需要

 levels(dat1[,"sex"])[1]

生成一行，并且将ORCI的值赋值为“1.00 (Reference)”

 data.frame(Variables= c("sex",paste0("  ",levels(dat1[,"sex"])[1])),
           Beta=c("",""),S.E=c("",""),Z=c("",""),P=c("",""),ORCI=c("","1.00 (Reference)"),
           OR=c(NA,NA),lower=c(NA,NA),upper=c(NA,NA))

结果为

7.2 表格合并

用这个结果，再去和mod_dat(去除第一行的常数项)进行纵向合并就好了

 rest_rowt <- rbind(data.frame(Variables= c("sex",paste0("  ",levels(dat1[,"sex"])[1])),
                              Beta=c("",""),S.E=c("",""),Z=c("",""),P=c("",""),ORCI=c("","1.00 (Reference)"),
                              OR=c(NA,NA),lower=c(NA,NA),upper=c(NA,NA)),
                   data.frame(Variables= paste0("  ",levels(dat1[,"sex"])[-1]),mod_dat[-1,])) %>% as_tibble()
rest_rowt 

到这里，我们就卷完了

7.3 变量的批量处理

处理连续变量，我们就不需要额外的去生成一个新的参照水平的数据，直接就可以用了，最后我们把所有的变量进行合并输出就行了，见全部的代码

var=c("education","sex","huji","smoking","drinking")
cat_var=c("education","sex","huji","smoking","drinking")
df <- dat1
####因子转化####
for (fv in cat_var){
  df[,fv] <- factor(df[,fv])
}
rest_row <- list()
for (i in 1:length(var)){
  print(paste0("正在分析",var[i]),quote=F)
  formula=paste0("insomnia~",var[i])
  mod <- glm(formula,family = binomial,data=df)
  mod_dat <- mod %>% summary %>% coef %>% as.data.frame() %>% 
    mutate(Beta=sprintf("%0.2f",Estimate)) %>% 
    mutate(S.E=sprintf("%0.2f",`Std. Error`)) %>% 
    mutate(Z=sprintf("%0.2f",`z value`)) %>% 
    mutate(P=ifelse(`Pr(>|z|)`<0.001,'<.001',sprintf("%0.3f",`Pr(>|z|)`))) %>% 
    mutate(OR=exp(Estimate)) %>% 
    mutate(OR=ifelse(OR<1.01 & OR>1 & `Pr(>|z|)`<0.05,1.01,#进行特殊情况的转化
                     ifelse(OR>0.99 & OR<1 & `Pr(>|z|)`<0.05,0.99,OR))) %>% 
    mutate(lower=exp(Estimate - qnorm(0.975)*`Std. Error`)) %>% 
    mutate(lower=ifelse(lower>1 & lower<1.01 & `Pr(>|z|)`<0.05,1.01,#进行特殊情况的转化
                        ifelse(lower<1 & lower>0.99 & `Pr(>|z|)`<0.05,0.99,lower))) %>% 
    mutate(upper=exp(Estimate + qnorm(0.975)*`Std. Error`)) %>% 
    mutate(upper=ifelse(upper>1 & upper<1.01 & `Pr(>|z|)`<0.05,1.01,#进行特殊情况的转化
                        ifelse(upper<1 & upper>0.99 & `Pr(>|z|)`<0.05,0.99,upper))) %>% 
    mutate(ORCI=paste0(sprintf("%0.2f",OR)," (",
                       sprintf("%0.2f",lower)," - ",
                       sprintf("%0.2f",upper),")")) %>% 
    dplyr::select(Beta,S.E,Z,P,ORCI,OR,lower,upper)
  rownames (mod_dat)<-NULL
  mod_dat <- mod_dat[-1,]
  if (df[,var[i]] %>% is.factor ==T){
    rest_row[i] <- rbind(data.frame(Variables= c(var[i],paste0("  ",levels(df[,var[i]])[1])),
                                  Beta=c("",""),S.E=c("",""),Z=c("",""),P=c("",""),ORCI=c("","1.00 (Reference)"),
                                  OR=c(NA,NA),lower=c(NA,NA),upper=c(NA,NA)),
                       data.frame(Variables= paste0("  ",levels(df[,var[i]])[-1]),mod_dat)) %>% list()
  }
  else if (df[,var[i]] %>% is.factor ==F){
    rest_row[i] <- data.frame(Variables= var[i],mod_dat)%>% list()
  }
}
final_single <- do.call(rbind,rest_row) %>% as_tibble()

结果见下图

7.4 结果的输出

结果导出word，结果

final_single <- do.call(rbind,rest_row) %>% 
  select(- OR,-lower,-upper)%>% as_tibble()
library("flextable")
c <- flextable(
  final_single,
  cwidth = 1,
  cheight = 0.1,
  defaults = list(),
  theme_fun = theme_booktabs
)
save_as_docx(c,path="sing.DOCX")

 更多的精彩，下回分享。