Logistic回归是最常用的多因素回归模型,在医学研究中,常用于研究疾病的危险因素,下面我们一起来看看,R语言是如何实现Logistic回归的。

1

第一步 导入数据

首先，在excel里全选数据集，右键复制。

然后，在Rstudio中，输入：

mydata "clipboard")

查看数据：

2

第二步 分类变量和等级变量转成因子型变量

mydata$sex$sex,levels=c(0,1),                  labels=c("Female","Male"))mydata$work$work,levels=c(0,1),                   labels=c("Non-work","Work"))mydata$disease$disease,levels=c(0,1),                       labels=c("Non-disease","Disease"))mydata$bmig$bmig,levels=c(1,2,3),                    labels=c("normalweight","overweight","obese"))mydata$bloodtype$bloodtype,levels=c(1,2,3,4),                      labels=c("O","A","B","AB"))

查看数据:

3

第三步进行单因素logistic回归

model1 <- glm(disease ~ bloodtype, data= mydata, family = binomial())

查看结果：

结果解读：

Estimate列表示回归系数，Std.Error表示标准误，z value表示统计量的值，Pr(>|z|)表示p值。我们可以看到bloodtypeA、bloodtypeB、bloodtypeC的结果，但是没看到bloodtypeO的结果。这是因为R语言logistic回归默认将分类变量的第一个factor设置为参照，通过前面的str(mydata)命令获得的数据集概况，可以看到bloodtype的levels顺序为：O、A、B、AB，所以在此回归模型中bloodtypeO当作参照。

获得OR值：

获得OR值95%可信区间：

科研论文通常需要我们提供OR值及其95%可信区间以及p值，虽然这些都得到了，但是貌似整理起来比较麻烦，别急，我们通过命令进行整理。具体而言：先通过broom包中tidy函数把model1的结果变成规范的数据框格式，提取出p值，再把OR值和OR的95%可信区间，组合到一起。

install.packages("broom")library(broom)y

zexp(coef(model1)),

结果解读：

exp.coef.model1表示OR值， X2.5和X97.5表示OR值的95%可信区间，y.p.value表示p值。Bloodtype作为无序多分类变量，需要设置为哑变量，一般哑变量的数目比分类变量的数目少一个，少掉的那个就作为参照(reference)。例如本例中，bloodtypeO就是参照，本文末尾有较为详细的说明。

4

第四步 多因素logistic回归

Model2 <- glm(disease ~sex+age+bmig+work+a+b+c+x+y, data = mydata, family = binomial())

整理成表格：

library(broom)mexp(coef(model2)),

结果输出到excel里：

write.csv(m,file="多因素logistic回归结果.csv")

在G盘的R文件夹中查看生成的多因素logistic回归结果，如图所示：

经过更改表头及简单的调整：

5

 进阶：哑变量及参照的设置

Logistic回归中一个重点内容，把无序多分类变量设置为哑变量，本文的例子中bloodtype是无序多分类变量，需要设置成哑变量，一个快捷的方法是把bmig设置为因子变量，做logistic时，默认为已经设置成哑变量。

还有一个重要的问题是如何设置哑变量的参照，本例中我们是以O型为参照，其他各型与它相比。如果我们想把A型设为参照，该如何操作呢？

我们先看下bloodtype的各level名称的顺序：

把A型作为参照，只需把A调到第一个位置：

mydata$bloodtype$bloodtype,

查看是否成功：

levels(mydata$bloodtype)

运行模型

model3 <- glm(disease ~sex+age+bmi+work++bloodtype+a+b+c+x+y, data = mydata, family = binomial())

关注我的，数据分析都不求人了不信你试试