只需3分钟！输出一篇logistic回归报告，准确度比肩R语言
在医学科研、特别是观察性研究领域，无论是现况调查、病例对照研究、还是队列研究，经常遇到分类的健康结局，包括二分类（如：生存与死亡、阳性与阴性、发病与未发病）或者一些可进行分类的生理生化指标等（如：血压值、血镁值、血脂和胆固醇等）时，线性回归分析往往无法进行，此时可以考虑Logistic回归模型。
logistic回归分析报告的统计模块主要包括3部分内容：统计描述、差异性分析、logistic回归。完成这三步，基本就可以形成一份统计报告了！
实际中，许多人习惯性使用SPSS进行回归分析，但是SPSS无法进行批量单因素分析，还需要手动绘制三线表，费时又费力。下面通过一份实操数据为大家展示R语言操作过程！
一、案例介绍
数据集来源于R自带MASS数据集birthwt，这是一份于1986年在在马萨诸塞州收集的与婴儿出生体重低相关的危险因素的数据。根据婴儿出生体重是否小于2.5kg，分为低出生组（59例）与正常组（130例）。研究的暴露因素见下表。
| 变量名称 | 变量类型 | 变量赋值 |
|---|---|---|
| low | 二分类变量 | 0：出生体重不小于2.5kg，1：出生体重小于2.5kg |
| age | 连续型变量 | 母亲的年龄，取值范围：14~45 |
| lwt | 连续型变量 | 母亲最后一次月经时的体重（磅），取值范围：80~250 |
| race | 无序多分类变量 | 1：白种人，2：黑种人，3：其他种族 |
| smoke | 二分类变量 | 0：怀孕期间不吸烟，1：怀孕期间吸烟 |
| ptl | 多分类变量 | 早产次数，取值0，1，2，3 |
| ht | 二分类变量 | 0：没有高血压病史，1：有高血压病史 |
| ui | 二分类变量 | 0：不存在子宫易怒，1：存在子宫易怒 |
二、 R语言软件复现
1.安装并加载R包
tableone 包，简化了”表 1”的构建，可以将生物医学研究论文中常见的患者数据汇总成一个表中混合的连续变量和分类变量。分类变量可以概括为计数和/或百分比。连续变量可以用”正态”（平均值和标准差）或”非正态”（中位数和四分位距）来描述。
compareGroups包，可以按组对多个变量进行描述。根据这些变量的性质，酌情计算不同的检验（t检验、方差分析、Kruskall-Wallis、Fisher、秩和…），从而可以很轻松的制作出SCI论文基线资料表或单因素分析表，也能做出SCI论文中多个模型比较的多因素分析表，甚至是线性趋势（P for trend）。今天的文章主要通过这个R包来进行差异性分析，是文章中的表二。
autoReg包，一款功能十分强大的R包，不仅可以快捷完成基线表的制作，还可以直接一行代码输出回归分析（支持线性模型、广义线性模型和比例风险模型）的表格。
install.packages("tableone")
install.packages("compareGroups")
install.packages("autoReg")
library(tableone)
library(compareGroups)
library(autoReg)
2.读取示例数据
library(MASS)
data(birthwt)
3.描述性统计
采用tableone包
①对定量数据进行正态性检验，本案例的样本量较小（n=189），采用W检验
shapiro.test(birthwt$lwt) #结果为偏态
##
## Shapiro-Wilk normality test
##
## data: birthwt$lwt
## W = 0.89331, p-value = 2.242e-10
②nonnormal = c(“lwt”)此处指定了偏态数据，那么在数据分析中会对偏态数据采用中位数，四分位数进行描述，未指定连续变量会采用均数±标准差，分类变量采用%进行描述。
vs<- c("low","age","lwt","race","smoke","ptd","ht","ui","ftv")
tableone <- CreateTableOne(vars = vs,data = birthwt)
## Warning in ModuleReturnVarsExist(vars, data): The data frame does not have: ptd
## Dropped
table1<-print(tableone, nonnormal = c("lwt"),showAllLevels = TRUE)
##
## level Overall
## n 189
## low (mean (SD)) 0.31 (0.46)
## age (mean (SD)) 23.24 (5.30)
## lwt (median [IQR]) 121.00 [110.00, 140.00]
## race (mean (SD)) 1.85 (0.92)
## smoke (mean (SD)) 0.39 (0.49)
## ht (mean (SD)) 0.06 (0.24)
## ui (mean (SD)) 0.15 (0.36)
## ftv (mean (SD)) 0.79 (1.06)
4.差异性分析
采用compareGroups包，也可以使用autoreg包
base_tab <- descrTable(low ~ age + lwt + race + smoke + ht + ui+ ftv,
data=birthwt,method = c (lwt=2))
print(base_tab)
##
## --------Summary descriptives table by 'low'---------
##
## ___________________________________________
## 0 1 p.overall
## N=130 N=59
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
## age 23.7 (5.58) 22.3 (4.51) 0.078
## lwt 124 [113;147] 120 [104;130] 0.013
## race 1.76 (0.91) 2.03 (0.91) 0.059
## smoke 0.34 (0.48) 0.51 (0.50) 0.031
## ht 0.04 (0.19) 0.12 (0.33) 0.083
## ui 0.11 (0.31) 0.24 (0.43) 0.040
## ftv 0.84 (1.07) 0.69 (1.04) 0.385
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
代码解读：descrTable(y~x1＋x2＋x3＋x4＋x……,data=数据集名，method = c (偏态数据=2))，指定因变量与自变量，设置数据集名。最后的method对指定偏态数据采用秩和检验，未指定连续变量采用t检验，分类数据采用卡方检验。
5.logisitic回归
采用autoReg包
使用glm() 函数构建回归模型，glm(y~x1＋x2＋x3＋x4＋x……,data=数据集名，family=“binomial”)，指定因变量与自变量，设置数据集名。
logfit<-glm(low ~ age + lwt + race + smoke + ht + ui + ftv,data=birthwt,family = "binomial")
summary(logfit)
##
## Call:
## glm(formula = low ~ age + lwt + race + smoke + ht + ui + ftv,
## family = "binomial", data = birthwt)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7426 -0.8398 -0.5698 1.0367 2.1293
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.117505 1.263702 -0.093 0.92592
## age -0.026944 0.035468 -0.760 0.44746
## lwt -0.013512 0.006547 -2.064 0.03901 *
## race 0.471209 0.213123 2.211 0.02704 *
## smoke 1.040777 0.391484 2.659 0.00785 **
## ht 1.851441 0.689782 2.684 0.00727 **
## ui 0.866535 0.451031 1.921 0.05470 .
## ftv 0.055545 0.169155 0.328 0.74263
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 206.72 on 181 degrees of freedom
## AIC: 222.72
##
## Number of Fisher Scoring iterations: 4
此处参数family规定了回归模型的类型:family=“binomial”指适用于二元离散因变量(binary)。
（1）显示单因素及多因素回归结果
logreg1<-autoReg(logfit,uni=TRUE)#显示单因素及多因素
logreg1
## —————————————————————————————————————————————————————————————————————————————————————————————————————————————
## Dependent: low 0 (N=130) 1 (N=59) OR (univariable) OR (multivariable)
## —————————————————————————————————————————————————————————————————————————————————————————————————————————————
## age Mean ± SD 23.7 ± 5.6 22.3 ± 4.5 0.95 (0.89-1.01, p=.104) 0.98 (0.91-1.04, p=.478)
## lwt Mean ± SD 133.3 ± 31.7 122.1 ± 26.6 0.99 (0.97-1.00, p=.023) 0.99 (0.97-1.00, p=.042)
## race Mean ± SD 1.8 ± 0.9 2.0 ± 0.9 1.38 (0.99-1.93, p=.060) 1.59 (1.05-2.41, p=.029)
## smoke Mean ± SD 0.3 ± 0.5 0.5 ± 0.5 2.02 (1.08-3.78, p=.028) 2.80 (1.30-6.03, p=.008)
## ht Mean ± SD 0.0 ± 0.2 0.1 ± 0.3 3.37 (1.02-11.09, p=.046) 6.20 (1.63-23.64, p=.008)
## ui Mean ± SD 0.1 ± 0.3 0.2 ± 0.4 2.58 (1.14-5.83, p=.023) 2.36 (0.98-5.72, p=.057)
## ftv Mean ± SD 0.8 ± 1.1 0.7 ± 1.0 0.87 (0.64-1.19, p=.389)
## —————————————————————————————————————————————————————————————————————————————————————————————————————————————
logtable1<-myft(logreg1)
（2）只显示多因素回归结果
logreg2<-autoReg(logfit,milti=TRUE)#只显示多因素
logreg2
## ——————————————————————————————————————————————————————————————————————————————————
## Dependent: low 0 (N=130) 1 (N=59) OR (multivariable)
## ——————————————————————————————————————————————————————————————————————————————————
## age Mean ± SD 23.7 ± 5.6 22.3 ± 4.5 0.97 (0.91-1.04, p=.448)
## lwt Mean ± SD 133.3 ± 31.7 122.1 ± 26.6 0.99 (0.97-1.00, p=.039)
## race Mean ± SD 1.8 ± 0.9 2.0 ± 0.9 1.60 (1.05-2.43, p=.027)
## smoke Mean ± SD 0.3 ± 0.5 0.5 ± 0.5 2.83 (1.31-6.10, p=.008)
## ht Mean ± SD 0.0 ± 0.2 0.1 ± 0.3 6.37 (1.65-24.62, p=.007)
## ui Mean ± SD 0.1 ± 0.3 0.2 ± 0.4 2.38 (0.98-5.76, p=.055)
## ftv Mean ± SD 0.8 ± 1.1 0.7 ± 1.0 1.06 (0.76-1.47, p=.743)
## ——————————————————————————————————————————————————————————————————————————————————
logtable2<-myft(logreg2)
（3）先单后多
logreg3<-autoReg(logfit,uni=TRUE,threshold=0.05)
logreg3
## —————————————————————————————————————————————————————————————————————————————————————————————————————————————
## Dependent: low 0 (N=130) 1 (N=59) OR (univariable) OR (multivariable)
## —————————————————————————————————————————————————————————————————————————————————————————————————————————————
## age Mean ± SD 23.7 ± 5.6 22.3 ± 4.5 0.95 (0.89-1.01, p=.104)
## lwt Mean ± SD 133.3 ± 31.7 122.1 ± 26.6 0.99 (0.97-1.00, p=.023) 0.98 (0.97-1.00, p=.013)
## race Mean ± SD 1.8 ± 0.9 2.0 ± 0.9 1.38 (0.99-1.93, p=.060)
## smoke Mean ± SD 0.3 ± 0.5 0.5 ± 0.5 2.02 (1.08-3.78, p=.028) 1.92 (1.00-3.71, p=.052)
## ht Mean ± SD 0.0 ± 0.2 0.1 ± 0.3 3.37 (1.02-11.09, p=.046) 6.84 (1.79-26.05, p=.005)
## ui Mean ± SD 0.1 ± 0.3 0.2 ± 0.4 2.58 (1.14-5.83, p=.023) 2.45 (1.03-5.84, p=.043)
## ftv Mean ± SD 0.8 ± 1.1 0.7 ± 1.0 0.87 (0.64-1.19, p=.389)
## —————————————————————————————————————————————————————————————————————————————————————————————————————————————
logtable3<-myft(logreg3)
（4）逐步回归法
logreg4<-autoReg(logfit,uni=TRUE,threshold=0.05, final=T) #final=T逐步回归 logreg4
logtable4<-myft(logreg4)
6.箱式图示例
Boxplot1<- ggboxplot(cancer,"race","size",
add = NULL,rug = TRUE,
color = "race",fill = NULL,width = 0.4,
palette = "npg", size = 0.4)
Boxplot<- Boxplot1+stat_compare_means(method = "kruskal.test",label.x.npc = "center")
Boxplot
