R语言机器学习篇——逻辑回归1（二分类）

参考书籍：陈强.机器学习及R应用.北京:高等教育出版社,2020

一 逻辑回归概念

在监督学习中存在大量关于“是与否”的二分类问题，以过滤垃圾邮件为例，假设响应变量只有两种可能取值（既0和1），此时该变量称为虚拟变量或是哑变量。

线性概率模型一般并不适合作预测。这是因为虽然明知y的取值非0即1，但根据线性概率模型所作的预测值却可能出现y＞1或y＜0的不现实情形，对于二分类问题,机器学习一般不用线性概率模型。

为使y的预测值总是介于[0,1]之间，在给定x时，考虑y的两点分布概率：

其中称为连接函数，它将特征向量x与响应变量y连接起来。

在给定x的情况下，y的条件期望为

E（y｜x）＝l *P（y＝l｜x）＋0*P（y＝0｜x）＝P（y＝l｜x）

如果连接函数为逻辑分布的累计分布函数，则

其中的定义为。此模型称为逻辑回归，简记为Logit，对逻辑函数求导，可得到逻辑分布的密度函数：

在R语言中，画出逻辑分布的密度函数及累积分布函数，如下所示：

x<-seq(-5,5,by=0.01)
y<-dlogis(x)    #密度函数
z<-plogis(x)    #分布函数
par(mfrow=c(1,2))
plot(x,y,type = "l",main="Logistic Density",ylab = "")
abline(v=0,lty=2)
plot(x,z,type = "l",main="Logistic CDF",ylab = "")
abline(v=0,lty=2)

输出结果如下所示：

二 Logit模型的解释

（一）参数意义

Logit模型属于非线性模型，因此参数并不表示边际效应，记事件“y＝1”发生的条件概率为p≡

P（y＝1｜x），则该事件不发生的条件概率为l-P＝P（y＝0｜x）。对于Logit模型来说，

,而，因此发生与不发生的几率为：

几率（odds）也称为相对风险，假如几率为2，此时说明发生概率是不发生概率的2倍，也就说明发生概率为2/3。

若对上述几率公式取对数，则可得出

也称为对数几率，因为右式为线性函数，该变换也称为逻辑变换

观察上式，回归系数表示当变量xk增加一个微小量时，引起对数几率的边际变化：

，因此可将参数理解为半弹性，既当对应特征变量x增加一单位所引起的几率变化的百分比。定义如下

这个式子可理解为特征变量变动了一个单位之后的新几率与原来几率的比值。

这里说明的虽为Logit模型的，但也预示了probit模型的略势：连接函数较为复杂，无法使用几率对其系数进行类似的解释。

（二）拟合优度

1.对于非线性模型,由于不存在平方和分解公式,故一般无法使用R方度量拟合优度。对于

使用MLE进行估计的非线性模型，可使用准R方或伪R方度量模型的拟合优度。

准

其中,为原模型的对数似然函数之最大值,而为以常数项为唯—变量的对数似然函数

之最大值。似然函数的最大可能值为1（即取值概率为1）,故对数似然函数的最大可能值为0,记为。显然’0≥≥,而0≤准≤1,

由于=0，可将上式写为

准

2.在统计学中，还可用偏离度的概念，也称为偏离残差度

定义为：（残差）偏离度≡

进—步，对于仅包含常数项的零模型，可定义其相应的零偏离度：

定义为：零偏离度≡

在零模型中加人其他变量之后，其偏离度的改进程度，则反映了这些变量对于y的解释能力:

零偏离度-偏离度＝-2（-）

事实上, -2（-）正是检验原假设“除常数项外所有回归系数均为0”的似然比检验

的统计量。根据残差偏离度与零偏离度，很容易计算准R方：

准

三 Logit模型的预测

1.得到模型估计的系数后，可预测条件概率

此时预测概率>0.5，则可预测出响应变量值为1，反之则为0

但当预测概率=0.5时，可预测为1或0。对于二分类问题，在特征空间中所有满足预测概率为0.5的样本点集合也称为决策边界。

决策边界

2.使用对数几率预测响应变量类别

如果对数几率>0，则可预测响应变量为1，反之则为0，

诺对数几率恰好为0，则预测值落在决策边界上。

四 二分类模型的评估

对于监督学习来说，一般用预测效果评估模型性能，具体到分类问题上，常用指标如下：

1.准确率：也称”正确预测的百分比“（将样本数据实际值与预测值比较得出，公式为：

注：I（·）为示性函数，n为样本容量。

2.错误率：也叫错分率，是指错误预测的百分比

当出现”类别不平衡现象“时，可用混淆矩阵，计算其中部分指标

3.混淆矩阵指标：

真阳性：TP ， 假阳性：FP ， 假阴性：FN ，真阴性：TN

灵敏度：预测正确的比例，也称为真阳率

灵敏度=真阳率

特异度：也称为真阴率

特异度=真阴率

注意：假阳率：1-特异度

从横向角度考察混淆矩阵时，可分析召回率：

召回率=查全率

4.ROC与AUC：

在分类时，门槛值的确定时极为重要的指标，在门槛值更低时，会导致预测更多的正例，更少的反例，因此灵敏度与“1-特异度”均为门槛值的函数，因此把灵敏度作为横轴，另外一个作为纵轴，再让门槛值从0取到1，得到一条曲线，称为接收器工作特征曲线，也称ROC曲线，而ROC曲线下的面积称为AUC，如下所示：

AUC的取值一般为0.5-1，为1的时候一般是无法达到的理想程度，表示对所有正例和反例的预测都是正确的。

5.科恩的kappa

对于类别不平衡的样本，同样可使用Kappa指标，已去掉准确率指标中可能的虚高“水分”，将混淆矩阵的结果转化为概率如下所示：

假设实际值与真实值分别来自两个不同的”评分者“，希望评估两个评分者的评分结果的一致性，因此观测的一致性，也就是准确率为

如果两个评分者的打分完全相互独立之时，因此此时期望的一致性为：

因此定义科恩的Kappa为：

对于kappa的取值所对应的解释如下所示，需要注意的是，他也可以用于多分类问题。

五 二分类问题的R案例

数据集源于R自带的Titanic，因为该数据集的类为表格，因此第一步就是将其变为数据框，如下所示：

titanic<-as.data.frame(Titanic)   #转化为数据框类型
str(titanic)    #考察数据框的结构

结果显示数据框包括32个观测，5个变量（其中响应变量定为：Survived)，但Freq表示该观测在数据框中出现的频率，因此展开后将不止32个观测。

接下来是将数据集随机分为训练集与测试集，先要展开它，根据变量Freq让不同观测值以对应频次出现，展开后再去掉Freq变量，最后观察剩下变量的统计特征，再用函数prop.table()函数，查看响应变量与特征变量之间的频率关系。（其中margin=1为设置参数）

#描述统计部分
titanic<-titanic[rep(1:nrow(titanic),titanic$Freq),]  #根据Freq列，重复观测个数  
titanic$Freq<-NULL  #令Freq列为0
dim(titanic)   #查看维度
summary(titanic)  #统计特征
prop.table(table(titanic$Sex,titanic$Survived),margin=1) 
prop.table(table(titanic$Age,titanic$Survived),margin=1)     

选择30%的样本作为测试集，剩下70%作为训练集，train_index表示从样本中随机挑选70%（1541）作为训练集

#训练集与测试集划分
set.seed(1)   #生成随机数种子
train_index<-sample(2201,1541)
train<-titanic[train_index,]
test<-titanic[-train_index,]

对训练集用模型进行拟合，其中参数“famliy = binomial”表示连接函数为二项分布，默认就是Logit模型，summary函数显示逻辑回归的主要结果。

#用训练集进行逻辑回归
fit<-glm(Survived~.,data = train,family = binomial) #参数“famliy”表示连接函数为二项分布
summary(fit)

结果显示经过4次“费雪得分跌代”得到估计结果，可看到零偏离度（仅包含常数项的残差平方和）为1939.3，残差偏离度（完整的）为1564.2。

另外在进行逻辑回归时，会将所有因子变为虚拟变量，如因子Sex（male＝l，female＝2）变为虚拟变量Sex（male＝0，female＝1）;但classCrew放入会导致严重的多重共线性，就没放入。

观察各参数的p值，可发现这些虚拟变量都在统计意义上高度显著。接下来计算准R方，计算拟合优度的大小，在观察参数值以及置信区间，

(fit$null.deviance-fit$deviance)/fit$null.deviance  #根据公式计算，下面为结果
#[1] 0.1934151
coef(fit)
confint(fit)   
exp(coef(fit))  #显示几率比
#结果 (Intercept)    Class2nd    Class3rd   ClassCrew   SexFemale   AgeAdult  
       2.0504460   0.3767243   0.1800478   0.4330681  10.5127092  0.3341689 

根据exp(coef(fit))可计算几率比，得出结果，此结果显示，如对于性别（分类型）来说，以响应变量取0（死亡）为参照，在给定其他特征变量时，若乘客性别由男性（SexFemale=0）变为女性（SexFemale=1）时，其存活的几率是其旧几率的10.51倍，如舱位(class)从一等变为三等时，存活几率将只有旧几率的0.18倍。

接下来使用margins包计算每个变量的平均边际效应，展示结果，并画图:

#平均边际效应AME
library(margins)
effects<-margins(fit)
summary(effects)
plot(effects,main="AME and Confidence Intervals")

至此模型的建模完成，使用该模型进行预测，并计算混淆矩阵

#模型预测
prob_train<-predict(fit,type="response")   #参数表示预测事件发生的条件概率,对数几率
pred_train<-prob_train>0.5
table<-table(predicted=pred_train,Actual=train$Survived)
table

#指标判断
Accuracy<-(table[1,1]+table[2,2])/sum(table)  #准确率
#[1] 0.7748215
Error_rate<-(table[2,1]+table[1,2])/sum(table)   #错分率
#[1] 0.2251785
Sensitivity<-table[2,2]/(table[1,2]+table[2,2]) #灵敏度
#[1] 0.4779116
Specificity<-table[1,1]/(table[1,1]+table[2,1])  #特异度
#[1] 0.9165868
Recall<-table[2,2]/(table[2,1]+table[2,2])   #召回率
#[1] 0.7323077

结果显示训练样本的预测准确率为0.77，而正例的预测准确率仅为0.47，而反例为0.91，故更易预测反例，较难预测正例。

接下来查看测试误差，用测试集进行预测

#测试集查看测试误差
prob_test<-predict(fit,type="response",newdata=test)  
pred_test<-prob_test>0.5
table<-table(predicted=pred_test,Actual=test$Survived)
table #混淆矩阵

查看测试集混淆矩阵的指标

Accuracy<-(table[1,1]+table[2,2])/sum(table)  
#[1] 0.7863636
Error_rate<-(table[2,1]+table[1,2])/sum(table)  
#[1] 0.2136364
Sensitivity<-table[2,2]/(table[1,2]+table[2,2]) 
#[1] 0.5211268
Specificity<-table[1,1]/(table[1,1]+table[2,1])  
#[1] 0.9127517
Recall<-table[2,2]/(table[2,1]+table[2,2])   
#[1] 0.74

可说明对于该案例，测试误差与训练误差比较接近。

最后画ROC曲线，并计算AUC面积，其中ROCR包的函数performance（）用来建立一个“性能对象”，其参数“y.measure＝＂tpr＂,表示y轴为“真阳率”，即灵敏度;而参数“x.measure＝

＂fpr＂则表示x轴为“假阳率”,即1-特异度。

#ROC曲线
library(ROCR)
pred_object<-prediction(prob_test,test$Survived)  #建立预测对象
perf<-performance(pred_object,measure = "tpr",x.measure = "fpr")
plot(perf,main="ROC curve(Test Set)",lwd=2,col="blue",
     xlab="1-Specificity",ylab="Sensitivity")
abline(0,1)
#计算AUC面积
performance(pred_object, "auc")@y.values    #用预测对象计算
#[1] 0.7733298
#计算Kappa值
library(vcd)
Kappa(table)

结果显示，在测试集中AUC面积为0.77，kappa值为0.47，这表明预测值与实际值之间具有中等一致性，但由于案例中的特征变量并不多且还是虚拟变量，这已经是不错的预测效果。