翻译 R语言逻辑回归分析

firedata 2016-05-25 10050 3 3

## 如何进行逻辑回归分析 逻辑回归是当y=f(x)，而y为分类变量的时候的逻辑曲线拟合的方法。这种模型通常的用法就是通过给定的一个x的预测值来预测y。这些预测值可以说连续的、分类的，或者是混合的。 通常来说，分类变量y有多种不同的假设值。其中，最简单的一个例子就是y为一个二元变量，这意味着我们可以假设每个y值为1或0。在机器学习里，一个非常常用的案例就是邮件分类：有一个针对每封邮件的给定属性集，比如字数、链接和图片。这时，我们要用的算法就是根据它判断这封邮件是垃圾邮件（1）或者不是非垃圾邮件（0）。在这篇博客当中，基于这个预测变量是二元的，我把它称之为“二元逻辑回归”。当然，逻辑回归同样也可以预测一个超过2个值的独立变量。在这篇博客中，我们称它为“多元逻辑回归”。在这里，一个典型的例子就是分类一部电影，它是“有趣”、“一般般”还是“无聊”。 ## R运行逻辑回归 R可以让逻辑回归建模过程变得很简单。我们可以使用glm()函数进行逻辑回归建模，而且，它的拟合过程和我们之前所做的线性回归没有很大的区别。在这篇博客中，我们将介绍如何一步一步的进行逻辑回归建模。 ## 数据集 我们将以Titanic数据集作为实例进行分析。关于这个数据集，在网上流传了好几个免费的版本，而我则建议大家使用在Kaggle上的那个版本，因为它已经被用过很多次了（如果你要下载这个数据集，你需要在Kaggle注册一个账号）。 这个（训练）数据集是关于一些乘客的信息（精确的数值为889），而这个比赛的目标就是基于诸如服务等级、性别、年龄等一些因素来预测一个顾客是否幸存（1表示乘客幸存，0则表示死亡）。在这里，你已经看到了，这个实例即有分类变量，也有连续变量。 ## 数据清洗过程 当我们拿到一个真实的数据集的时候，我们需要确认这个数据集有哪些缺失值或者异常值。因此，我们需要做这方面的准备工作以便我们后面对此进行数据分析。第一步则是我们使用read.csv()来读取数据集，可以在这里下载数据集。 确认参数na.strings相当于c(“”)，因此，每个缺失值可以用NA表示。这可以帮助我们进行下一步的工作。 ```r training.data.raw <- read.csv('train.csv',header=T,na.strings=c("")) ``` 现在，我们要查阅一下缺失值，并观察一下每个变量存在多少异常值，使用sapply()函数来展示数据框的每一列的值。 ```r sapply(training.data.raw,function(x) sum(is.na(x))) PassengerId Survived Pclass Name Sex 0 0 0 0 0 Age SibSp Parch Ticket Fare 177 0 0 0 0 Cabin Embarked 687 2 sapply(training.data.raw, function(x) length(unique(x))) PassengerId Survived Pclass Name Sex 891 2 3 891 2 Age SibSp Parch Ticket Fare 89 7 7 681 248 Cabin Embarked 148 4 ``` 在这里对缺失值进行可视化会有帮助：Amelia包有一个特别的作图函数missmap()来针对我们的数据集进行作图，并且对缺失值做标记。 ```r library(Amelia) missmap(training.data.raw, main = "Missing values vs observed") ```  这个变量有非常多的缺失值，我们都不会使用他们。我们也会把Passengerld去掉，因为它只有索引和票。 使用subset()函数，我们可以从原始数据集中只选择与之相关的列。 ```r 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()函数，而默认情况下是分类变量被转成因子类型。一个因子就是R处理分类变量的方式。 我们可以使用下面的代码来查阅它： ```r is.factor(data$Sex) TRUE is.factor(data$Embarked) TRUE ``` 为了能更容易的理解R是怎么处理分类变量的，我们可以使用contrasts()函数。这个函数向我们展示了R是如何优化变量的，并把它转译成模型。 ```r contrasts(data$Sex) male female 0 male 1 contrasts(data$Embarked) Q S C 0 0 Q 1 0 S 0 1 ``` 比如，你可以看见在变量sex，female可被作为参考。在乘船上的缺失值，由于这里只有两行，我们可以把这两行删掉(我们也可以用众数取代缺失值，然后保留这些数据点）。 ```r data <- data[!is.na(data$Embarked),] rownames(data) <- NULL ``` 在继续拟合过程之前，让我们回顾一下数据清洗和重构的重要性。这样的预处理通常对于获得一个好而且预测能力强的模型是非常重要的。 ## 模型拟合 我们把数据分成两个部分，训练和测试集。训练数据集可用于模型的拟合，而测试数据集则对此进行测试。 ```r train <- data[1:800,] test <- data[801:889,] ``` 现在，让我们拟合一下这个模型。确认我们在glm()函数使用参数family=binomial。 ```r model <- glm(Survived ~.,family=binomial(link='logit'),data=train) ``` 现在，使用summary()函数来获得这个模型的结果。 ``` r summary(model) Call: glm(formula = Survived ~ ., family = binomial(link = "logit"), data = train) Deviance Residuals: Min 1Q Median 3Q Max -2.6064 -0.5954 -0.4254 0.6220 2.4165 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 5.137627 0.594998 8.635 < 2e-16 *** Pclass -1.087156 0.151168 -7.192 6.40e-13 *** Sexmale -2.756819 0.212026 -13.002 < 2e-16 *** Age -0.037267 0.008195 -4.547 5.43e-06 *** SibSp -0.292920 0.114642 -2.555 0.0106 * Parch -0.116576 0.128127 -0.910 0.3629 Fare 0.001528 0.002353 0.649 0.5160 EmbarkedQ -0.002656 0.400882 -0.007 0.9947 EmbarkedS -0.318786 0.252960 -1.260 0.2076 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 1065.39 on 799 degrees of freedom Residual deviance: 709.39 on 791 degrees of freedom AIC: 727.39 Number of Fisher Scoring iterations: 5 ``` ## 解读逻辑回归模型的结果 现在，我们可以分析这个拟合的模型，阅读一下它的结果。首先，我们可以看到SibSp，Fare和Embarked并不是重要的变量。对于那些比较重要的变量，sex的p值最低，这说明sex与乘客的生存几率的关系是最强的。这个负系数预测变量表明其它变量都一样，男乘客生存的机率更低。记住，逻辑模型的因变量是对数机率：ln(odds) = ln(p/(1-p)) = a*x1 + b*x2 + … + z*xn。male是一个优化变量，男性的生还机率下载2.75个对数机率，而年龄的增大则下降0.037个对数机率。 现在，我们可以使用anova()函数来分析这个数据表的偏差了： ```r anova(model, test="Chisq") Analysis of Deviance Table Model: binomial, link: logit Response: Survived Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev Pr(>Chi) NULL 799 1065.39 Pclass 1 83.607 798 981.79 < 2.2e-16 *** Sex 1 240.014 797 741.77 < 2.2e-16 *** Age 1 17.495 796 724.28 2.881e-05 *** SibSp 1 10.842 795 713.43 0.000992 *** Parch 1 0.863 794 712.57 0.352873 Fare 1 0.994 793 711.58 0.318717 Embarked 2 2.187 791 709.39 0.334990 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ``` 空偏差和残差的不同展示了我们的模型是空模型（只有截距）。这个沟阅读，效果越好。分析这个图表，我们可以看到每增加一个变量，同时也会降低其中的偏差。再一次，添加Pclass、Sex和Age能显著的降低残差。其它变量对于模型的提高没什么用，尽管SibSp的p值很低。一个p值较大的值预示着一个没有变量的模型多少解释了模型的变化。最终，你想要看到的就是偏差的降低和AIC。 然而，这里并不存在线性回归的残差R^2,McFaddenR^2指数可用于模型的拟合。 ```r library(pscl) pR2(model) llh llhNull G2 McFadden r2ML r2CU -354.6950111 -532.6961008 356.0021794 0.3341513 0.35917 ``` ## 评估模型的预测能力 在上述的步骤中，我们简单的评估一下模型的拟合效果。现在，我们看一下使用心得数据集，它的预测结果y是什么样的。设定参数type=’response’，R可以输出形如P(Y=1|X）的概率。我们的预测边界就是0.5。如果P(Y=1|X)>0.5，y=1或y=0。记住，其它的决策边界可能是更好的选择。 ```r 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)) "Accuracy 0.842696629213483" ``` 精度为0.84，说明这个模型拟合效果不错。然而，你要记住，结果多少都依赖于手动改变数据集的数据。因此，如果你想要得到更精确的精度，你最好执行一些交叉检验，例如k次交叉检验。 作为最后一步，我们会做ROC曲线并计算AUC（曲线下的面积），它常用于预测二元分类器的模型表现。 ROC曲线是一种曲线，它可以通过设定各种极值来让正例律（TPR)来抵消反正例律(FPR)，它就在ROC曲线之下。通常来说，一个预测能力强的模型应当能让ROC接近1（1是理想的）而不是0.5。 ```r 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图：  我希望这篇文章对于你来说有用。这篇文章的详细代码和案例在这里。 感谢您能阅读这篇文章，如果你有什么问题，都可以在后面留言。 >严禁修改，可以转载，请注明出自[数据人网](http://shujuren.org)和原文链接。