R语言| 批量单因素logistic回归

文章来源： https://t.1yb.co/k9eM

本期介绍使用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值

1-4. gtsummary包输出的不是95%而是97.5%CI

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])}

--------------------------------------------------------------------------------------------------------------------------

说明：本期数据5000例，14个变量，结局status=0为死亡，1为生存；13个自变量（一个连续变量，12个分类变量）

一、常用的logistic回归代码或R包

01-glm() 手动提取OR、95%CI及P值

glm（）广义线型模型使用注意

glm不能直接输出OR值及95%CI需根据回归系数和标准误计算；

glm直接用exp(confint())提取到的是97.5%CI；

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

-------------------------------------------------

结果：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

-----------------------------------------------------------

• 注意1：glm直接用exp(confint( ))提取到的是97.5%CI；
• 注意2：glm不加family=binomial其结果是错误的

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值需要进一步计算。

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);res

----------------------------------------------------------------------------

04. gtsummary包输出的不是95%而是97.5%C

虽然结果很美观，但是可以发现其输出的结果与上面提到的注意1：CI97<-exp(confint(glm1));CI97输出的97.5%CI是一致的。因此不建议使用。

library(gtsummary)

#glm（）做逻辑回归

glm1<- glm(status==0~t,

family = binomial,

data = aa)

#tbl_regression提取结果

res<-tbl_regression(glm1,

exponentiate=T);res

----------------------------------------------------------------------

小结

上述4种方法，前三种均可输出OR、95%CI及P值，但经反复比较后发现，只有glm( )手动计算OR、95%CI运行是报错时最少的，所我选用第一种方法＋function(x)循环进行批量单因素logistic回归制作三线表。

===========================================

二、批量单因素logistic回归代码详解

说明：带入自己数据时只需要修改代码的3处红色标记即可！！！

文章来源： https://t.1yb.co/k9ny