当结果变量为分类变量时,我们该如何创建回归模型呢?这个时候,结果变量不再是连续型变量,也无法与预测变量呈“线性关系“了。比如,如何预测检查者是否得病等“是否”的问题。这时,数据科学家们发明了逻辑回归(Logistic regression)。它的原理很简单,既然“是/否“(二分类)变量没法跟像连续型变量一样产生”线性关系“,那么我们就把二分类变量变成连续型变量呗~
逻辑回归(又称Logistic回归)用于根据一个或多个预测变量(x)预测样本的类别。这个“类别”作为结果变量,只可能是两个值:0或者1,是或否,患病或未患病。我们也将这种只能取两个值(二进制)的变量称为二分类变量。
Logistic回归属于广义线性模型(Generalized Linear Model, GLM)范畴,它是对线性回归理论的扩展而开发出来的统计学方法。它有时候也被称为:二进制逻辑回归,二项逻辑回归和logit模型,logit回归。
逻辑回归并不直接返回观察样本的所属类别,而是该样本被划分到某个类别的概率(p),概率在0到1之间。您需要自己定义一个是0还是1的阈值概率。默认情况下,此设置为p = 0.5,即结果被划分到1的概率大于0.5时,我们定义该事件为“是”,小于0.5时,该事件为“否”。但实际上有时候需要您自己根据实际情况来确定这个阈值概率。
在一般情况下,标准的逻辑回归可以画成一条S形曲线,它的x轴为预测变量(x),y轴为二分类结果变量发生事件“是”的概率(P)。
公式为:
·y = b0 + b1*x, ·exp() 是指数函数 ·p是事件“是(1)”在预测变量取x值时发生的概率。在数学上,写为p(event=1|x) ,简写为p(x)。
b0和b1等是回归的beta系数。正数b1表示增加x会增加事件“是”的概率p 。相反,负数b1表示增加x减小p。
该公式中log[p/(1-p)]称为比值(odd)的对数,也称为log-odd或logit。
该比值p/(1-p)反映事件”是“发生的可能性。可以将其视为“是”与“否”之比。该比值是事件发生的概率除以事件不会发生的概率。例如,如果糖尿病阳性的概率是0.5,则“不会”的机率是1-0.5 = 0.5,总概率是1.0。
此外,我们也可以通过比值计算事件概率
·tidyverse 便于数据操作和可视化
·caret 简化工作流程
library(tidyverse)library(caret)
在进行逻辑回归分析前,执行以下数据清理步骤,可提高模型的准确性
删除潜在的异常值
确保预测变量呈正态分布。如果没有,则可以使用对数,平方根等方式进行数据转换。但是这一步不是必须的,这跟线性回归的假设检验相同,最终目的是保证方差齐性。我们也会在之后对此进行详细解说。
删除高度相关的预测变量,以最大程度地减少过度拟合。防止多重共线性的发生。高度相关的预测变量的存在可能会导致模型不稳定。多重共线性的相关内容。
我们将使用一组糖尿病的数据,它包含768人糖尿病相关的数据,包括怀孕情况,血糖,血压,皮肤厚度,胰岛素水平,体重指数,糖尿病谱系功能,年龄和糖尿病诊断结果(Outcome)。这个数据可以实现以下研究问题:
根据多个临床变量预测糖尿病阳性的概率。
导入数据
my_data<-read.csv('diabetes.csv')
检查数据
dim(my_data)
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)
我们将数据随机分为训练集(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, ]
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 freedom
Residual deviance:277.28 on 305 degrees of freedom
AIC: 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)
从以上结果我们可以看到,血糖,体重指数,糖尿病谱系功能都对结果产生独立预测功能(Pr(>|z|) < 0.05)。改模型对结果的整体预测准确性达到了80%以上(0.8076923)。
解读R运行结果的内容,以及如何在论文写作中,解释这些变量。
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, ]
一元逻辑回归为用一个预测变量来预测事件概率的逻辑回归,预测变量只有一个。
以下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)]
使用此公式,对于每个新获得的血糖值,您都可以预测个体出现糖尿病阳性的概率。
使用预测函数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)。
逻辑函数给出了S形概率曲线,我们可以通过如下代码画出:
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" )
当有多个预测变量在模型中用来预测二分类的结果事件发生概率时,即可用多元逻辑回归。
例如,用血糖、体重指数和怀孕次数来预测糖尿病阳性概率,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
可以看出,在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)。
OR值表示疾病(结果变量)与暴露(预测变量)之间关联强度的指标,与相对危险度(RR)类似,指暴露者的疾病危险性为非暴露者的多少倍。
当我们使用病例对照研究时,不能计算发病率,所以也不能计算相对危险度(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)。
从逻辑回归结果中可以注意到,一些变量(如胰岛素水平和血压值)在统计学上不显着。将它们保留在模型中可能会导致过拟合。因此,在之后的模型中,应该把它们删除,以达到选择一个具有减少的变量集的最佳模型,而不影响模型的准确性。
我们将使用测试数据进行预测,以评估逻辑回归模型的性能。
步骤如下:
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
模型准确性的衡量标准是预测值与实际观测值(正确分类)进行比较,获取比值。分类误差定义为被错误分类的预测值与实际观测值的比例。
正确分类的比例:
mean(predicted.classes == test.data$Outcome)
分类预测的准确性约为81%,即错误分类错误率为19%。
参考内容:
1. Alboukadel Kassambara, Machine Learning Essentials: Practical Guide in R