第三章逻辑回归（下）

2018-09-30 本文已影响93人 Stone_Stan4d

利用lasso进行正则化

前一章线性回归，利用glmnet包，通过岭回归和lasso进行正则化。去除模型的某些特征也许是好思路，我们尝试对数据集运用lasso并评估结果。首先用glmnet()训练一系列正则化的模型，然后我们调用cv.glmnet()来估计合适的 $\lambda$ 值。之后利用 $\lambda$ 检验我们正则化的模型的系数：

library(glmnet)
heart_train_mat <- model.matrix(OUTPUT ~ ., heart_train)[, -1]
lambdas <- 10^seq(8, -4, length = 250)
heart_models_lasso <- glmnet(heart_train_mat, heart_train$OUTPUT, alpha = 1, 
                             lambda = lambdas, family = 'binomial')
lasso.cv <- cv.glmnet(heart_train_mat, heart_train$OUTPUT, lambda = lambdas, alpha = 1,
                      family = 'binomial')
lambda_lasso <- lasso.cv$lambda.min
lambda_lasso

predict(heart_models_lasso, type = 'coefficients',s = lambda_lasso)

我们不一样，同样的代码，四次运行出来的结果不一样。神奇！重复运行代码如下：

heart_models_lasso <- glmnet(heart_train_mat, heart_train$OUTPUT, alpha = 1, 
                             lambda = lambdas, family = 'binomial')
lasso.cv <- cv.glmnet(heart_train_mat, heart_train$OUTPUT, lambda = lambdas, alpha = 1,
                      family = 'binomial')
lambda_lasso <- lasso.cv$lambda.min
lambda_lasso

predict(heart_models_lasso, type = 'coefficients',s = lambda_lasso)

一次结果	二次结果
三次结果	四次结果

有一些特征的系数为0，它们被移除模型了。再次衡量训练集和测试集的分类精确度，结果再一次与原著不一样。。。模型少用了三个特征，精确度在训练集上略微提升，测试集没有变化。

lasso_train_predictions <- predict(heart_models_lasso, s = lambda_lasso, 
                                   newx = heart_train_mat, type = 'response')
lasso_train_class_predictions <- as.numeric(lasso_train_predictions > 0.5)
mean(lasso_train_class_predictions == heart_train$OUTPUT)
###
heart_test_mat <- model.matrix(OUTPUT ~ ., heart_test)[, -1]
lasso_test_predictions <- predict(heart_models_lasso, s = lambda_lasso, 
                                  newx = heart_test_mat, type = 'response')
lasso_test_class_predictions <- as.numeric(lasso_test_predictions > 0.5)
mean(lasso_test_class_predictions == heart_test$OUTPUT)

我的笔记本精分了！原著，测试集的精确率是0.925

3.6 分类指标

对训练集计算查准率、查全率和F1分数：

(confusion_matrix <- table(predict = train_class_predictions, actual = heart_train$OUTPUT))
(precision <- confusion_matrix[2, 2] / sum(confusion_matrix[2,]))
(recall <- confusion_matrix[2, 2] / sum(confusion_matrix[, 2]))
(f <- 2 * precision * recall /(precision + recall))

在医学里，查全率也叫做灵敏度（sensitivity），也就是被识别为1的数据占实际类别为1的数据的比例，也就是真阳性样本被识别出来的比例。查全率又叫做真阳性率。特异度是真阴性率，是阴性样本被正确识别为阴性的比例。

(specificity <- confusion_matrix[1,1] / sum(confusion_matrix[1,]))

image.png

阈值的重要性

library(ROCR)
train_predictions <- predict(heart_model, newdata = heart_train, type = 'response')
pred <- prediction(train_predictions, heart_train$OUTPUT)
perf <- performance(pred, measure = 'prec', x.measure = 'rec')
plot(perf, main = '心脏病模型查全-查准曲线')

在这些衡量指标中，我们先来看设置阈值为0.5的重要性。如果要选择一个不同的阈值，就必须清楚的了解，所有之前的衡量指标都会改变。具体说，有很多情况下，我们会想要调整阈值，让它偏向于把预测数据类别识别为类别1，当前的医疗场景就是一个主要的示例。例如，假定模型会被一位临床医生用来确定是否让一位病人接受一种详细也更昂贵的心脏病检查。我们很可能会认为把一位心脏有问题的病人错误标记为健康和要求健康的病人接受进一步的检查相比是更严重的错误，因为他们其实是不健康的。例如，为了达到这种偏误，我们可以把分类阈值降低到0.3到0.2。
理想情况下，我们想要的是一种评价改变阈值对性能指标影响效果的视觉方式，而查准率-查全率曲线就是这样的一种有用的图形，可以利用ROCR包来获得查准率-查全率曲线：

library(ROCR)
train_predictions <- predict(heart_model, newdata = heart_train, type = 'response')
pred <- prediction(train_predictions, heart_train$OUTPUT)
perf <- performance(pred, measure = 'prec', x.measure = 'rec')
plot(perf, main = '心脏病模型查全-查准曲线')

然后绘制perf对象来获得我们的查准率-查全率曲线。

例如，这个图形为我们显示，要得到0.8以上的查全率，必须在查准率上做出相当大的牺牲。要对阈值进行调优，我们会需要看每一个用老计算这个绘图的阈值。一个游泳的联系是创建一个阈值的数据框，里面是在数据中查准率和查全率产生变化的阈值，以及它他们对应的查准率和查全率值。然后可以对这个数据框划分自己来逐个检查我们感兴趣的阈值。

thresholds <- data.frame(cutoffs = perf@alpha.values[[1]], recall = perf@x.values[[1]],
                         precision = perf@y.values[[1]])
subset(thresholds, (recall > 0.9)&(precision > 0.8))
[out] cutoffs    recall precision
112 0.3491857 0.9019608 0.8288288
113 0.3472740 0.9019608 0.8214286
114 0.3428354 0.9019608 0.8141593
115 0.3421438 0.9019608 0.8070175

正如我们所见，一个大约为0.35的阈值就可以满足我们的需求。

这里用@来访问perf对象的某些属性，该对象是S4类的特殊类型的对象。S4类在R语言里面提供面向对象的特性。面向对象编程与R语言可参考《Advanced R》

3.7 二元逻辑分类器的扩展

我们现在要处理多个类别的预测问题。在第一章里，我们研究了鸢尾花数据集（iris），当时的目标是根据描述鸢尾花的外观来区分三种不同的鸢尾花品种。在讲解其他多元分类的示例之前，有一个重要的提示，就是我们在本书中还会学习到其他几种分类方法，如神经网络和决策树，它们处理多类别的分类，比逻辑回归更自然也更常用。

多元逻辑回归

假设我们的目标变量是由 $K$ 个类别组合而成，如在iris数据集中， $K = 3$ 。多元逻辑回归（Multinomial logistic regression）通过拟合 $K-1$ 个独立的二元回归分类器模型来处理这种。这是通过随意选择输出其中一个类别作为参考类别，并拟合 $K-1$ 个把每一个其他类别和这个类别比较的回归模型。例如，如果我们有两个特征 $X_1$ 和 $X_2$ ，以及可以称为0、1、2的三个类别，我们就构建如下的两个模型：
$\ln(\frac{P(Y=1|X)}{P(Y=0|X)}) = \beta_{10} + \beta_{11}X_1 + \beta_{12} X_2$
$\ln(\frac{P(Y=2|X)}{P(Y=0|X)}) = \beta_{20} + \beta_{21}X_1 + \beta_{22} X_2$
在这里以类别0做基线，构建两个二元逻辑回归模型。第一个模型用类别1和类别0作比较，第二个模型是类别2和类别0的比较。模型系数有两个下标，第一个表示模型，第二个表示某个特征的系数。当总类别为K，编号依次为类别0到k，通用表达式为：
$\begin{equation} P(Y=1|X)=\left \{ \begin{aligned} &\frac{1}{1 + \displaystyle\sum_{j=1}^{K-1} e^{f_j(x)}}, &k=0 \\ &\frac{ e^{f_j(x)}}{1 + \displaystyle\sum_{j=1}^{K-1} e^{f_j(x)}}, &k>0\\ \end{aligned} \right . \end{equation} \\f_j(x) = \beta_0 + \beta_{j1}X_1+\cdots+\beta_{jp}X_p$
读者应该检验一下，根据规则，所有输出类别概率的总和是1。这个由一个指数除以一组指数之和的数学公式就称为softmax函数。对于前面讨论的3个类别的问题，我们就可以直接把前面的公式用 $K=3$ 进行替换。到这里，我们应该提到这个方法的某些重要特点。

首先，训练集的模型数量比输出变量里的类别总数少1个，因此很容易就这种方法在输出具有较大候选类别的情况下，是不太好扩展的。如果要创建和训练相当多的模型，这就意味着我们往往需要一个大得多的数据集才能产生有合理精确度的结果。最后，由于是独立的用每一个输出类别去比较一个参考，所以我们要做出一个叫做无关选项的独立性(Independence of Irrelevant Alternatives， IIA)的假设。

简而言之，IIA假设表达的是，预测某个输出而不是其他输出的机率并不依赖于我们是否通过加入新的类别来增大可能的输出类别数量k。要说明这个假设，为了简单起见，我们要利用多元逻辑回归对iris数据集进行建模，使三种不同品种对应的输出类别是0.33:0.33:0.33，因此每个品种和其他品种的比例都是1:1。IIA假设表达的是，如果我们重新拟合一个包含了某个鸢尾花的新品种（例如，日本鸢尾花ensata）的模型，原来三种鸢尾花品种出现的机率会保持不变。如这4个品种之间的新的总体机率比为0.2:0.2:0.2:0.4就是合法的，原来3个品种之间保持1:1的关系不变。

预测玻璃模型

在本节，我们要通过一个示例数据集来讲解如何在R语言里训练多院逻辑回归模型。我们这次要分析的数据来自法医学邻域。在这个示例中，我们的目标是考察玻璃碎片的来源（比如车前灯）。这个玻璃识别数据集(glass identification data set)存放在UCI数据库（Glass + Identification）。我们首先要把这个数据集加载到一个数据框，利用来自网站的信息对各列进行重命名，然后丢弃样本识别码--第一列，因为它是任意分配的，无特别含义：

glass <- read.csv('YuCeFenXiRawData/glass.data', header = F)
names(glass) <- c('id', 'RI', 'Na', 'Mg', 'Al', 'Si', 'K', 'Ca', 'Ba', 'Fe', 'Type')
galss <- glass[, -1]

下一步，我们来看一个表格，它显示了数据框每一列的意义：

列名	类型	定义
Ri	数值	折射率
Na	数值	氧化钠重量的百分比
Mg	数值	氧化镁重量的百分比
Al	数值	氧化铝重量的百分比
Si	数值	二氧化硅重量的百分比
K	数值	氧化钾重量的百分比
Ca	数值	氧化钙重量的百分比
Ba	数值	氧化钡重量的百分比
Fe	数值	氧化铁重量的百分比
Type	类别	玻璃类别(1.浮法生产的建筑物窗户，2.非法浮法生产的建筑物，3.浮法生产的车辆窗户，4.非浮法生产的车辆窗户，5.玻璃容器，6.餐具，7.车前灯)

照例，我们会先为玻璃数据准备一个训练集和测试集：

set.seed(4365677)
glass_sampling_vector <- createDataPartition(glass$Type, p = 0.8, list = F)
glass_train <- glass[glass_sampling_vector, ]
glass_test <- glass[-glass_sampling_vector, ]

现在用nnet扩展包。这个包还有用在神经网络上的函数，所以我们在下一章中会再次用这个包。multinom()是用于多元逻辑回归的。它需要给定一个公式和一个数据框来调用，所以它有一个我们熟悉的数据接口。此外，我们可以指定maxit参数，它决定了要运行的相关优化程序的最大迭代次数。有时候，我们会发现某个模型的训练会返回一个“未达到收敛”的错误。在这种情况下，一种可能的方法是增大这个参数，让模型经过更多的迭代次数来训练。不过，如果这样做的话，我们必须注意到一个事实，那就是模型会需要更多的训练时间：

library(nnet)
glass_model <- multinom(Type ~ ., data = glass_train, maxit = 1000)
summary(glass_model)
[out]
Call:
multinom(formula = Type ~ ., data = glass_train, maxit = 1000)

Coefficients:
  (Intercept)        id         RI         Na         Mg       Al        Si          K         Ca         Ba         Fe
2   4.5071292  9.406355  20.273529  71.824943  -39.19729 210.1157 -25.19036  179.66305  -3.890667   80.13429  139.96585
3  11.9239996 16.059512  -9.038460   7.759758  105.14286 247.4771 -33.78759 -180.44821  19.827318   21.88280  -82.09221
5   0.5806134 15.618658   1.683779   6.406146 -105.39458 141.0011 -21.99433  234.02978 -14.317176  -56.88068  171.29687
6  -8.8594453 14.660756 -13.700699 170.947923  -69.08264 220.9859 -44.71005 -336.95778 -51.746540 -490.16299 -211.30828
7  -1.5278885 17.993106   9.842946  81.473920  -70.99193 129.1253 -34.64918   85.89812 -60.620948   16.69416 -198.39591

Std. Errors:
  (Intercept)        id        RI        Na         Mg         Al        Si            K        Ca           Ba           Fe
2   13.353250 1332.2895 20.619982 469.85029 328.406126 495.818262 1129.7073 1.873162e+02 175.76720 1.479062e+01 2.046871e+02
3   15.982381 1341.2168 21.960056 501.35812 306.374853 499.757725 1123.0411 1.851590e+02 209.00633 5.114700e-30 2.046871e+02
5    4.169490  691.3656  6.344835  48.19161   7.811845   6.492859  303.7922 1.952961e+00  47.62448 1.197345e-38 1.296635e-31
6    5.223616  929.8036  7.936605  72.03366  12.588914   6.216103  380.0703 6.469344e-25  51.03473 6.848498e-24 6.949558e-33
7   51.727755 1105.8908 76.543448 688.85030 191.882894  60.439775 3842.8012 2.785129e+01 294.30291 6.848498e-24 1.015335e-49

Residual Deviance: 6.812649e-05 
AIC: 110.0001

模型模型摘要表明，我们有5组系数。这是因为TYPE输出变量有6个水平，也就是说我们要在6种不同的玻璃来源中预测出一种。在数据里没有Type取值为4的示例。该模型还为我们显示了标准误，但是没有显著性检验结果。总体而言，对于系数的显著性检验比在二元逻辑回归里

第三章逻辑回归（下）

利用lasso进行正则化

3.6 分类指标

3.7 二元逻辑分类器的扩展

多元逻辑回归

猜你喜欢

热点阅读

第三章 逻辑回归 （下）

利用lasso进行正则化

3.6 分类指标

3.7 二元逻辑分类器的扩展

多元逻辑回归

猜你喜欢

热点阅读

第三章逻辑回归（下）