如何在R中执行Logistic回归？

通常，分类变量y可以采用不同的值。在最简单的情况下，y是二进制的，意味着它可以假设值为1或0.机器学习中使用的经典示例是电子邮件分类：给定每个电子邮件的一组属性，如单词数量，链接和图片，算法应该决定电子邮件是否为垃圾邮件（1）或不（0）。在这篇文章中，我们将模型称为“二项逻辑回归”，因为要预测的变量是二元的，然而，逻辑回归也可以用来预测一个可以假设多于2个值的因变量。在第二种情况下，我们称之为“多项逻辑回归”模型。例如，一个典型的例子就是将电影在“娱乐”，“边缘”或“无聊”之间进行分类。

R中的逻辑回归实现

R可以很容易地拟合逻辑回归模型。被调用的函数是glm（），拟合过程与线性回归中使用的不同。在这篇文章中，我将拟合一个二元逻辑回归模型并解释每一步。

数据集

我们将研究泰坦尼克号数据集。这个数据集有不同版本可以在线免费获得，但我建议使用Kaggle提供的数据集，因为它几乎可以使用（为了下载它，你需要注册Kaggle）。数据集（训练）是关于一些乘客的数据集合（准确地说是889），并且竞赛的目标是预测生存（如果乘客幸存，则为1，否则为0）基于某些诸如服务等级，性别，年龄等特征。正如您所看到的，我们将使用分类变量和连续变量。

数据清理过程

在处理真实数据集时，我们需要考虑到某些数据可能丢失或损坏的事实，因此我们需要为我们的分析准备数据集。作为第一步，我们使用read.csv（）函数加载csv数据。确保参数na.strings等于c（“”），以便将每个缺失值编码为NA。这将帮助我们接下来的步骤。

training.data.raw <- read.csv(
‘train.csv’,
header = TRUE,
na.strings = c(“”)

现在我们需要检查缺失的值，并使用sapply（）函数查看每个变量有多少个唯一值，该函数将作为参数传递的函数应用于数据帧的每一列。

sapply(
training.data.raw,
function(x) sum(is.na(x))
)

sapply(
training.data.raw,
function(x) length(unique(x))
)

对缺失值进行可视化处理可能会有所帮助：Amelia包具有特殊的绘图函数missmap（），它将绘制数据集并突出显示缺失值：

变量carbin有太多的缺失值，我们不会使用它。我们也会放弃PassengerId，因为它只是一个索引和Ticket。使用subset（）函数，我们将原始数据集进行子集化，只选择相关列。

data <- subset(training.data.raw,select=c(2,3,5,6,7,8,10,12))

处理缺失值

现在我们需要解释其他缺失的值。通过在拟合函数内设置参数来拟合广义线性模型时，R可以很容易地处理它们。但是，我个人更倾向于在可能的情况下“手动”处理缺失值。有不同的方法可以做到这一点，一种典型的方法是用现有的平均值，中位数或模式代替缺失值。我将使用平均值。

data$Age[is.na(data$Age)] <- mean(data$Age,na.rm=T)

就分类变量而言，默认情况下使用read.table（）或read.csv（）会将分类变量编码为因子。一个factor是R如何处理分类变量。
我们可以使用以下代码行来检查编码

is.factor(data$Sex)
is.factor(data$Embarked)

为了更好地理解R如何处理分类变量，我们可以使用contrasts（）函数。这个函数将向我们展示变量如何被R虚拟化，以及如何在模型中解释它们。

contrasts(data$Sex)
contrasts(data$Embarked)

例如，你可以看到变量sex时，女性将被用作参考。至于Embark中的缺失值，由于只有两个值，我们将丢弃这两行（我们也可以用模式替换缺失的值并保留数据点）。

data <- data[!is.na(data$Embarked),]
rownames(data) <- NULL

在进行拟合过程之前，让我提醒您清洁和格式化数据的重要性。这个预处理步骤对于获得模型的良好拟合和更好的预测能力通常是至关重要的。

模型拟合

我们将数据分成两部分：训练和测试集。训练集将用于拟合我们将在测试集上进行测试的模型。

train <- data[1:800,]
test <- data[801:889,]

现在，我们来拟合模型。一定要在glm（）函数中指定参数family = binomial。

model <- glm(Survived ~.,family=binomial(link=’logit’),data=train)

通过使用函数summary（），我们获得了我们模型的结果：

summary(model)

部分结果如下图：

解释我们的逻辑回归模型的结果

现在我们可以分析拟合并解释模型告诉我们什么。
首先，我们可以看到SibSp，Fare和Embarked没有统计学意义。至于统计上显着的变量，性别具有最低的p值，这表明乘客的性别与存活的可能性有很强的关联。该预测因子的负系数表明所有其他变量相同，男性乘客不太可能存活下来。请记住，在Logit模型中，响应变量是对数可能性：ln（odds）= ln（p /（1-p））= a x1 + b x2 + … + z * xn。由于男性是虚拟变量，因此男性将对数赔率降低2.75，而单位年龄增加则将对数赔率降低0.037。

现在我们可以在模型上运行anova（）函数来分析偏差表

anova(model, test=”Chisq”)

零偏差和残差偏差之间的差异显示了我们的模型如何对付零模型（仅包含截距的模型）。这个差距越大越好。通过分析表格，我们可以看到每次添加一个变量时偏差下降的情况。同样，增加Pclass，Sex and Age可以显着减少残余偏差。尽管SibSp具有较低的p值，但其他变量似乎可以改进模型。这里的大p值表示没有变量的模型差不多也能解释相同的变化量。最终你想看到的是偏差和AIC的显着下降。

尽管不存在与线性回归R2的精确等价，但McFadden R2指数可用于评估模型拟合。

library(pscl)
pR2(model)

评估模型的预测能力

在上面的步骤中，我们简要地评估了模型的拟合，现在我们想看看模型在预测新的一组数据时是如何进行的。通过设置参数type =’response’，R将以P（y = 1 | X）的形式输出概率。我们的决策边界将是0.5。如果P（y = 1 | X）> 0.5，则y = 1，否则y = 0。请注意，对于某些应用程序，不同的阈值可能是更好的选择。

fitted.results <- predict(model,newdata=subset(test,select=c(2,3,4,5,6,7,8)),type=’response’)
fitted.results <- ifelse(fitted.results > 0.5,1,0)
misClasificError <- mean(fitted.results != test$Survived)
print(paste(‘Accuracy’,1-misClasificError))

“准确性0.842696629213483”
测试集上的0.84精度是相当不错的结果。但是，请记住，这个结果有点依赖于我之前做的数据的手动分割，因此如果您希望得到更精确的分数，最好运行某种交叉验证，如k折交叉验证。

作为最后一步，我们将绘制ROC曲线并计算二元分类器典型性能测量的AUC（曲线下面积）。
ROC是通过在各种阈值设置下绘制真阳性率（TPR）与假阳性率（FPR）而产生的曲线，而AUC是ROC曲线下的面积。作为一个经验法则，具有良好预测能力的模型应该具有接近于1（1是理想的）的AUC。

library(ROCR)
p <- predict(model, newdata=subset(test,select=c(2,3,4,5,6,7,8)), type=”response”)
pr <- prediction(p, test$Survived)
prf <- performance(pr, measure = “tpr”, x.measure = “fpr”)
plot(prf)
auc <- performance(pr, measure = “auc”)
auc <- auc@y.values[[1]]
auc

0.8647186
ROC如下图所示：

完整代码

# 任务描述
# R执行逻辑回归算法
training.data.raw <- read.csv(
‘train.csv’,
header = TRUE,
na.strings = c(“”)
sapply(
training.data.raw,
function(x) sum(is.na(x))
sapply(
training.data.raw,
function(x) length(unique(x))
library(Amelia)
missmap(training.data.raw, main = “Missing values vs observed”)
data <- subset(training.data.raw,select=c(2,3,5,6,7,8,10,12))
data$Age[is.na(data$Age)] <- mean(data$Age,na.rm=T)
is.factor(data$Sex)
is.factor(data$Embarked)
contrasts(data$Sex)
contrasts(data$Embarked)
data <- data[!is.na(data$Embarked),]
rownames(data) <- NULL
train <- data[1:800,]
test <- data[801:889,]
model <- glm(Survived ~.,family=binomial(link=’logit’),data=train)
summary(model)
anova(model, test=”Chisq”)
library(pscl)
pR2(model)
fitted.results <- predict(model,newdata=subset(test,select=c(2,3,4,5,6,7,8)),type=’response’)
fitted.results <- ifelse(fitted.results > 0.5,1,0)
misClasificError <- mean(fitted.results != test$Survived)
print(paste(‘Accuracy’,1-misClasificError))
library(ROCR)
p <- predict(model, newdata=subset(test,select=c(2,3,4,5,6,7,8)), type=”response”)
pr <- prediction(p, test$Survived)
prf <- performance(pr, measure = “tpr”, x.measure = “fpr”)
plot(prf)
auc <- performance(pr, measure = “auc”)
auc <- auc@y.values[[1]]
auc