R语言| 批量单因素logistic回归
来源：一只勤奋的科研喵
欢迎大家关注我的公众号：一只勤奋的科研喵
本文来源：https://t.1yb.co/k8On
本期介绍使用R语言实现一次性输出所有变量的单因素logistic回归三线表。思路与R语言|6. 批量单因素Cox回归类似。先是构建一套能提取OR、95%CI及P值的函数，然后将该函数至于function(x)循环中，最后将需要进行单因素logistic回归的变量带入循环批量输出结果。 几种函数含义：lrm ()=逻辑回归模型、glm()=广义线性模型、lm ()=线性模型
目 录
1. 常用的logistic回归代码或R包
1-1. glm()手动提取OR、95%CI及P值
glm()的logistic回归需加family=binomial
glm()不能直接数输出95%CI
1-2. rms包的lrm()直接输出95CI及P值
1-3. epiDisplay包直接输出OR、95%CI及P值
2. 批量单因素logistic回归代码详解
1、构建循环函数；
2、将变量带入循环函数；
3、批量生成单因素OR、95%CI及P值三线表；
4、优化三线表
#结果合并需要的包
library(plyr)
#可进行logistic回归的包
library(rms)#可实现逻辑回归模型（lrm）
library(epiDisplay)#快速输出OR、95%CI、P
library(gtsummary)#精美三线表（但，95%CI有误）
#1.清理工作环境
rm(list = ls())
#2.数据放入工作目录，读取
aa<- read.csv('批量单因素逻辑回归.csv')
#3.查看变量名
names(aa)
#4.查看变量性质
str(aa)
#分类变量要变为因子，本期数据源数据不用转变，自己的数据需要转变可用以下函数批量转变。
#只改红字，数字意义是需要进行转换的变量所在源数据的列数
#for (i in names(aa)[c(1,4:12)]){aa[,i] <- as.factor(aa[,i])}
1.png
一、常用的logistic回归代码或R包
01-glm() 手动提取OR、95%CI及P值
glm（）广义线型模型使用注意
glm做logistic回归需加入family=binomial，解释如下
参数 family 规定回归模型的类型
family=gaussian适用于一维连续因变量，服从高斯分布，即正态分布，对应的模型为线性回归
family=mgaussian 说明因变量为服从高斯分布的连续型变量，但是有多个因变量，输入的因变量为一个矩阵，对应的模型为线性回归模型
family=poisson"适用于非负次数因变量，离散型变量，服从泊松分布，对应的模型为泊松回归
family=binomial 适用于二元离散因变量，服从二项分布，对应的模型为逻辑回归
family=multinomial 适用于多元离散因变量，对应的模型为逻辑回归模型
glm1<- glm(status==0~t,
data=aa,
family = binomial)
glm2<- summary(glm1);glm2
2.png
结果：glm（）汇报了p值、回归系数β、标准误SE，根据函数exp(coef(glm1) 和exp(β-1.96*SE)我们可以计算出OR和95%CI
#计算OR并保留两位数
OR<-round(exp(coef(glm1)),2)
#提取SE
SE<-glm2$coefficients[,2]
#计算CI，保留两位数并合并
CI5<-round(exp(coef(glm1)-1.96*SE),2)
CI95<-round(exp(coef(glm1)+1.96*SE),2)
CI<-paste0(CI5,'-',CI95)
#提取P值
P<-round(glm2$coefficients[,4],3)
#OR、95%CI、P值合并
res1<-data.frame(OR,CI,P)[-1,];res1
3.png
CI97<-exp(confint(glm1));CI97
运行结果：
2.5 % 97.5 %
(Intercept) 0.1436722 0.2235527
tT2 1.1585637 1.9212430
tT3 1.1386044 1.8714825
tT4 2.6330820 4.3754725
- 注意2：glm不加family=binomial其结果是错误的
glm1<- glm(status==0~t,
data=aa)
glm2<- summary(glm1);glm2
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.15271 0.01699 8.988 < 2e-16 ***
tT2 0.05854 0.02009 2.914 0.00359 **
tT3 0.05489 0.01963 2.796 0.00520 **
tT4 0.22588 0.02137 10.569 < 2e-16 ***
02. rms包的lrm函数直接输出95%CI及P值
library(rms)
#为后续代码设定环境
bb<-datadist(aa)
options(datadist='bb')
#lrm
lrm1<-lrm(status==0~t,data = aa);lrm1
#结果
lrm2 <-summary(lrm1)
OR<-round(exp(coef(lrm1)),2);OR
lrm（）结果与带有family=binomial的glm（）是一致的，但OR值需要进一步计算。
3-1.png
03-epiDisplay包直接提取OR、95%CI及P值
library(epiDisplay)
glm1<-glm(status==0~t,family = binomial,data = aa)
#提取OR/CI/P值
res<- logistic.display(glm1, simplified=T);re
此方法，最简单且简洁，结果与上述方法一致
运行结果：
OR lower95ci upper95ci Pr(>|Z|)
tT2 1.486011 1.154250 1.913127 2.120391e-03
tT3 1.453608 1.134096 1.863138 3.140919e-03
tT4 3.380248 2.622830 4.356393 4.986736e-21
04. gtsummary包
library(gtsummary)
#glm（）做逻辑回归
glm1<- glm(status==0~t,
family = binomial,
data = aa)
#tbl_regression提取结果
res<-tbl_regression(glm1,
exponentiate=T);res
4.png
小结
上述4种方法均可输出OR、95%CI及P值，但经反复比较后发现，只有glm( )手动计算OR、95%CI运行是报错时最少的，所我选用第一种方法＋function(x)循环进行批量单因素logistic回归制作三线表。
二、批量单因素logistic回归代码详解
说明：带入自己数据时只需要修改代码的3处红色标记即可！！！
e027e1aeea549cfd5a1f891f37666ca.png
7722f9f58856ac10b4884885ede8361.png
08b821f2588d3e128110b2b3206c6a7.png
bf8e912a0dd614b9af96ea453129b6d.png
54396b1ab8e8f5c96adbd4b422d6072.png
