科研星球

R语言统计与绘图:R-逻辑回归概论与结果解读

当结果变量为分类变量时,我们该如何创建回归模型呢?这个时候,结果变量不再是连续型变量,也无法与预测变量呈“线性关系“了。比如,如何预测检查者是否得病等“是否”的问题。这时,数据科学家们发明了逻辑回归(Logistic regression)。它的原理很简单,既然“是/否“(二分类)变量没法跟像连续型变量一样产生”线性关系“,那么我们就把二分类变量变成连续型变量呗~



1. 逻辑回归


逻辑回归(又称Logistic回归)用于根据一个或多个预测变量(x)预测样本的类别。这个“类别”作为结果变量,只可能是两个值:0或者1,是或否,患病或未患病。我们也将这种只能取两个值(二进制)的变量称为二分类变量。


Logistic回归属于广义线性模型(Generalized Linear Model, GLM)范畴,它是对线性回归理论的扩展而开发出来的统计学方法。它有时候也被称为:二进制逻辑回归,二项逻辑回归和logit模型,logit回归。


逻辑回归并不直接返回观察样本的所属类别,而是该样本被划分到某个类别的概率(p),概率在0到1之间。您需要自己定义一个是0还是1的阈值概率。默认情况下,此设置为p = 0.5,即结果被划分到1的概率大于0.5时,我们定义该事件为“是”,小于0.5时,该事件为“否”。但实际上有时候需要您自己根据实际情况来确定这个阈值概率。






Logistic函数
 


在一般情况下,标准的逻辑回归可以画成一条S形曲线,它的x轴为预测变量(x),y轴为二分类结果变量发生事件“是”的概率(P)。


图片


公式为:

p = exp(y) / [1 + exp(y)] 
也可以简单地写成p = 1/[1 + exp(-y)],其中:


·y = b0 + b1*x,

·exp() 是指数函数

·p是事件“是(1)”在预测变量取x值时发生的概率。在数学上,写为p(event=1|x) ,简写为p(x)



所以,以上公式也可以写为
p(x)= 1 / [1 + exp(-(b0 + b1 * x))]
以上公式也可以写成
p/(1-p) = exp(b0 + b1*x)
通过取双方的对数,公式成为预测变量的线性组合
log[p/(1-p)] = b0 + b1*x

当我们有多个预测变量时,逻辑函数即为:
log[p/(1-p)] = b0 + b1*x1 + b2*x2 + ... + bn*xn

b0和b1等是回归的beta系数。正数b1表示增加x会增加事件“是”的概率p 。相反,负数b1表示增加x减小p。


该公式中log[p/(1-p)]称为比值(odd)的对数,也称为log-oddlogit

比值p/(1-p)反映事件”是“发生的可能性。可以将其视为“是”与“否”之比。该比值是事件发生的概率除以事件不会发生的概率。例如,如果糖尿病阳性的概率是0.5,则“不会”的机率是1-0.5 = 0.5,总概率是1.0。


此外,我们也可以通过比值计算事件概率

p = Odds/(1 + Odds)


2. 逻辑回归R实现


2.1 加载所需的R包


·tidyverse 便于数据操作和可视化

·caret 简化工作流程




library(tidyverse)library(caret)


2.2 准备数据


在进行逻辑回归分析前,执行以下数据清理步骤,可提高模型的准确性


  • 删除潜在的异常值


  • 确保预测变量呈正态分布。如果没有,则可以使用对数,平方根等方式进行数据转换。但是这一步不是必须的,这跟线性回归的假设检验相同,最终目的是保证方差齐性。我们也会在之后对此进行详细解说。


  • 删除高度相关的预测变量,以最大程度地减少过度拟合。防止多重共线性的发生。高度相关的预测变量的存在可能会导致模型不稳定。多重共线性的相关内容。


我们将使用一组糖尿病的数据,它包含768人糖尿病相关的数据,包括怀孕情况,血糖,血压,皮肤厚度,胰岛素水平,体重指数,糖尿病谱系功能,年龄和糖尿病诊断结果(Outcome)。这个数据可以实现以下研究问题:


根据多个临床变量预测糖尿病阳性的概率。



导入数据

my_data<-read.csv('diabetes.csv')

检查数据

dim(my_data)
[1] 768  9



head(my_data)


 输出结果

Pregnancies Glucose BloodPressure SkinThickness Insulin  BMI DiabetesPedigreeFunction Age Outcome1          6     148            72            35       0 33.6                    0.627  50      12          1      85            66            29       0 26.6                    0.351  31      03          8     183            64             0       0 23.3                    0.672  32      14          1      89            66            23      94 28.1                    0.167  21      05          0     137            40            35     168 43.1                    2.288  33      16          5     116            74             0       0 25.6                    0.201  30      0



检查数据

summary(my_data)

  

输出结果

Pregnancies        Glucose      BloodPressure    SkinThickness      Insulin           BMI        DiabetesPedigreeFunction Min.   :0.000   Min.   : 0.0   Min.   : 0.00   Min.   : 0.00  Min.   :  0.0  Min.   : 0.00   Min.  :0.0780           1st Qu.: 1.000   1st Qu.: 99.0   1st Qu.: 62.00   1st Qu.: 0.00   1st Qu.: 0.0   1st Qu.:27.30   1st Qu.:0.2437           Median : 3.000   Median :117.0   Median : 72.00   Median :23.00   Median : 30.5   Median :32.00   Median :0.3725           Mean   :3.845   Mean   :120.9  Mean   : 69.11   Mean  :20.54   Mean   : 79.8  Mean   :31.99   Mean  :0.4719           3rd Qu.: 6.000   3rd Qu.:140.2   3rd Qu.: 80.00   3rd Qu.:32.00   3rd Qu.:127.2   3rd Qu.:36.60   3rd Qu.:0.6262           Max.  :17.000   Max.   :199.0  Max.   :122.00   Max.  :99.00   Max.   :846.0  Max.   :67.10   Max.  :2.4200                Age           Outcome      Min.  :21.00   Min.   :0.000  1st Qu.:24.00  1st Qu.:0.000   Median :29.00  Median :0.000   Mean  :33.24   Mean   :0.349  3rd Qu.:41.00  3rd Qu.:1.000   Max.  :81.00   Max.   :1.000


我们可以看到,该数据中血糖,血压,皮肤厚度,胰岛素水平,体重指数都存在0值,但是实际情况告诉我们,这些值都不可能为0值的,说明这些值都是缺失值,我们需要讲有缺失值的数据删除。在实际情况中,回归分析也要求数据保证完整性,不能有缺失值,如果有缺失值,系统可能会报错。


清理数据,去除缺失值。

new_data<-my_data[my_data$Glucose>0& my_data$Insulin >0 & my_data$BMI >0 & my_data$BloodPressure>0 & my_data$SkinThickness > 0,]dim(new_data)
[1] 392  9


我们将数据随机分为训练集(80%用于构建预测模型)和测试集(20%用于评估模型)。设置可重复性的set.seed(),确保该“随机”之后可重复。

set.seed(123)training.samples <- new_data$Outcome %>% createDataPartition(p = 0.8, list = FALSE)train.data  <- new_data[training.samples, ]test.data <- new_data[-training.samples, ]



2.3逻辑回归函数


glm()函数为广义线性模型的R函数,可用于计算逻辑回归。您需要指定选项 family = binomial,它告诉R,我们要拟合逻辑回归。注意!在进行逻辑回归时,这个选项一定要写,否则进行的将不是逻辑回归,而是线性回归。


R代码如下

# 拟合模型

model<-glm( Outcome ~., data=train.data, family=binomial)


# 查看结果

summary(model)


 # 输出结果
Call:glm(formula = Outcome ~ .,family = binomial, data = train.data)
Deviance Residuals:    Min      1Q   Median       3Q     Max  -2.7107  -0.6749 -0.3785   0.6385   2.5480 Coefficients:                           Estimate Std. Errorz value Pr(>|z|)    (Intercept)              -9.8698967  1.3253634 -7.447 9.55e-14 ***Pregnancies               0.0787163  0.0604556  1.302   0.1929    Glucose                   0.0377628  0.0063444  5.952 2.65e-09 ***BloodPressure            -0.0039203  0.0128683 -0.305   0.7606    SkinThickness             0.0083579  0.0189184  0.442   0.6586    Insulin                  -0.0006349  0.0014448 -0.439   0.6603    BMI                       0.0721097  0.0298554  2.415   0.0157 *  DiabetesPedigreeFunction  1.0185073 0.4694417   2.170   0.0300 * Age                       0.0376236  0.0200916  1.873   0.0611 .  ---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’1 (Dispersion parameter forbinomial family taken to be 1)    Null deviance: 398.80  on 313 degrees of freedomResidual deviance:277.28  on 305  degrees of freedomAIC: 295.28 Number of Fisher Scoringiterations: 5


# 建立预测

probabilities<-model%>%predict(test.data, type="response")predicted.classes<-ifelse(probabilities>0.5, 1, 0)


# 查看模型的拟合准确性(预测结果与实际结果相等的百分比)


mean(predicted.classes==test.data$Outcome)


[1] 0.8076923


2.4 结果解读


从以上结果我们可以看到,血糖,体重指数,糖尿病谱系功能都对结果产生独立预测功能(Pr(>|z|) < 0.05)。改模型对结果的整体预测准确性达到了80%以上(0.8076923)。


解读R运行结果的内容,以及如何在论文写作中,解释这些变量。


1 数据前期准备

library(tidyverse)library(caret)


#导入数据:(如需获取数据Outcome.csv,请关注投必得医学公众号,后台回复“Outcome.csv”获取数据。)



my_data<-read.csv('Outcome.csv')


#检查数据





dim(my_data)head(my_data)summary(my_data)


#清理数据,去除缺失值。




new_data<-my_data[my_data$Glucose>0& my_data$Insulin >0 & my_data$BMI >0 & my_data$BloodPressure>0 & my_data$SkinThickness > 0,]dim(new_data)


#将数据随机分为训练集和测试集

set.seed(123)training.samples <- new_data$Outcome %>% createDataPartition(p = 0.8, list = FALSE)train.data  <- new_data[training.samples, ]test.data <- new_data[-training.samples, ]



2 一元逻辑回归


2.1 逻辑回归模型建立


一元逻辑回归为用一个预测变量来预测事件概率的逻辑回归,预测变量只有一个。

以下R代码建立了一个模型:

根据血糖浓度预测糖尿病阳性的概率:


model <- glm( Outcome ~ Glucose, data = train.data, family = binomial)summary(model)$coef


输出结果

             Estimate Std. Error z value Pr(>|z|)(Intercept) -6.0593351 0.704588302 -8.5998247.983874e-18Glucose     0.0419662 0.005307411  7.9070952.634654e-15


上面的输出显示了回归系数β的估计值及其显着性水平。截距(b0)为-6.06,变量血糖的系数为0.042。


逻辑方程可以写成

p = exp(-6.06 + 0.042*glucose)/ [1 + exp(-6.06 + 0.042*glucose)]

使用此公式,对于每个新获得的血糖值,您都可以预测个体出现糖尿病阳性的概率。


2.2 逻辑回归预测功能


使用预测函数predict()来获取。使用选项type=“response”直接获得概率


newdata <- data.frame(Glucose = c(120,  180)) #设置血糖值两个:120和180probabilities <- model %>% predict(newdata, type = "response")predicted.classes <- ifelse(probabilities > 0.5, 1, 0)predicted.classes


输出结果



#1 2 #0 1


上面2X2表格结果解读结果为:第一个血糖值(1,120)为糖尿病阴性(0),第二个血糖值(2,180)为糖尿病阳性(1)。



2.3 逻辑回归结果可视化


逻辑函数给出了S形概率曲线,我们可以通过如下代码画出:


640 (2).png

train.data %>%  mutate(prob = ifelse(Outcome == 1, 1, 0)) %>%  ggplot(aes(Glucose, prob)) +  geom_point(alpha = 0.2) +  geom_smooth(method = "glm", method.args = list(family = "binomial")) +  labs(    title = "Logistic Regression Model",    x = "Plasma Glucose Concentration",        y = "Probability of being diabete-pos"      )





3. 多元逻辑回归


3.1 多元逻辑回归模型建立


当有多个预测变量在模型中用来预测二分类的结果事件发生概率时,即可用多元逻辑回归。

例如,用血糖、体重指数和怀孕次数来预测糖尿病阳性概率,R代码如下所示:


model <- glm( Outcome ~ Glucose + BMI + Pregnancies,                 data = train.data, family = binomial)summary(model)$coef




     输出结果  


             Estimate  Std. Error   z value    Pr(>|z|)(Intercept) -8.76578780 1.072319167 -8.1746072.968317e-16Glucose     0.03811238 0.005376785  7.0883221.357479e-12BMI         0.07784041 0.021922802  3.5506603.842671e-04Pregnancies 0.15395137 0.044067743  3.4935164.767042e-04


如果我们要包括数据集中所有可用的预测变量,即数据集中除了结果变量意外所有的变量都被纳入预测变量的范畴,则可用使用~.以下命令完成此操作:


model <- glm( Outcome ~., data = train.data, family = binomial)summary(model)$coef




    输出结果


                             Estimate  Std. Error    z value     Pr(>|z|)   (Intercept)              -9.8698967344 1.325363416 -7.4469362 9.553282e-14   Pregnancies               0.0787162818 0.060455599  1.3020511 1.928989e-01   Glucose                   0.0377628160 0.006344416  5.9521339 2.646686e-09   BloodPressure            -0.0039203075 0.012868285 -0.3046488 7.606337e-01   SkinThickness             0.0083578960 0.018918377  0.4417871 6.586432e-01   Insulin                  -0.0006348993 0.001444800 -0.4394375 6.603446e-01   BMI                       0.0721097400 0.029855378  2.4153015 1.572219e-02   DiabetesPedigreeFunction  1.0185073117 0.469441703  2.1696140 3.003610e-02   Age                       0.0376235670 0.020091611  1.8726008 6.112353e-02


从上面的输出中,系数表显示β系数估计值及其显着性水平:



·Estimate:截距(b0,以上显示为-9.87),即与每个预测变量相关的β系数估计(如Pregnancies的系数为0.078)

·Std.Error:系数估计值的标准误。这代表了系数的准确性。标准误差越大,估计值的准确性越低。

·z value:z统计量,即系数估计值(第2列Estimate)除以估算值的标准误(第3列Std.Error)

·Pr(>|z|):对应于z统计量的p值。p值越小,估计值对结果变量来说越重要。



请注意,这些函数coef()和summary()可用于提取系数,如下所示:

coef(model)



    输出结果

 (Intercept)              Pregnancies                  Glucose            BloodPressure            SkinThickness             -9.8698967344             0.0787162818             0.0377628160            -0.0039203075             0.0083578960                   Insulin                      BMI DiabetesPedigreeFunction                      Age             -0.0006348993             0.0721097400             1.0185073117             0.037623567



summary(model )$coef



3.2 结果解释


可以看出,在8个预测变量中,只有3个与结果显着相关。这些包括:血糖, 体重指数和血糖谱系功能。


Glucose的系数估计为b = 0.038,为正。这意味着血糖的增加与糖尿病阳性概率的增加有关。BloodPressure的系数b = -0.004,这是负数。这意味着血压升高与糖尿病阳性的可能性降低有关,但是由于p = 0.76,p > 0.05,不具有统计显着性,所以不能认为这种负相关有实际统计学意义。


与线性回归的直观性不同,逻辑回归是广义的线性回归模型。预测变量本身与结果变量不成线性相关,但是经过变换以后,

log[p/(1-p)] = b0 + b1*x1 + b2*x2 + ... + bn*xn

预测变量的系数(b1)与变换后的结果变量log[p/(1-p)]线性相关。


为了解释逻辑回归中的系数b,我们首先要理解的一个重要概念比值比(odds ratio,OR)。


640.png


OR值表示疾病(结果变量)与暴露(预测变量)之间关联强度的指标,与相对危险度(RR)类似,指暴露者的疾病危险性为非暴露者的多少倍。


640 (1).png


当我们使用病例对照研究时,不能计算发病率,所以也不能计算相对危险度(relative risk,RR),只能用OR作为反映关联强度的指标。OR>1说明疾病的危险度因暴露而增加,与疾病之间为“正”关联;OR<1说明疾病的危险度因暴露而减少,暴露与疾病之间为“负”关联;or=1说明疾病与暴露无关联。


在逻辑回归中,我们得到的系数,即为OR值的对数。对于连续型变量,比如血糖值,即表示改变一个单位的情况下,比值比的改变。


例如,血糖的回归系数为0.038,表明在其他预测变量都相同的情况下,血糖浓度每改变1个单位,患糖尿病的比值比的对数就会增加0.038,即比值比改变为改变前的exp(0.038) = 1.39倍。如果血糖改变2个单位,则比值比改变为改变前的exp(0.038X2)=1.079倍。在实际情况中,总体率较低时(如病例发生率小于10%),病例对照研究的OR值近似相对危险度(RR)。


从逻辑回归结果中可以注意到,一些变量(如胰岛素水平和血压值)在统计学上不显着。将它们保留在模型中可能会导致过拟合。因此,在之后的模型中,应该把它们删除,以达到选择一个具有减少的变量集的最佳模型,而不影响模型的准确性。


4. 预测


我们将使用测试数据进行预测,以评估逻辑回归模型的性能。

步骤如下:



1.根据预测变量实际观察值,通过模型对事件发生概率进行计算

2.定义概率大于0.5为事件发生“1”,小于0.5为事件不发生“0”。



predict()函数可用于计算预测的糖尿病阳性的概率。

预测患糖尿病的可能性:


probabilities <- model %>% predict(test.data, type = "response")head(probabilities)



    输出结果     


   5          9         26         28         32         36 0.88210454 0.86271962 0.37788421 0.03377201 0.593642740.14298930


这些概率指的是哪些类?我们可以通过以下代码对其进行分类,概率大于0.5,Outcome定义为1,否则定义为0。

predicted.classes <- ifelse(probabilities > 0.5, 1,0)head(predicted.classes)


 输出结果    

5 9 26 28 32 36  1 1  0  0 1  0



5. 评估模型准确性


模型准确性的衡量标准是预测值与实际观测值(正确分类)进行比较,获取比值。分类误差定义为被错误分类的预测值与实际观测值的比例。

正确分类的比例:


mean(predicted.classes == test.data$Outcome)
[1]0.8076923

分类预测的准确性约为81%,即错误分类错误率为19%。





参考内容:

1. Alboukadel Kassambara, Machine Learning Essentials: Practical Guide in R


没有账号?