R-机器学习秘籍-二-

布客飞龙
115 阅读39分钟

R 机器学习秘籍（二）

原文：annas-archive.org/md5/167951e4c9873bc214756c25a63518c8

译者：飞龙

协议：CC BY-NC-SA 4.0

第四章：理解回归分析

在本章中，我们将介绍以下食谱：

  • 使用 lm 函数拟合线性回归模型

  • 总结线性模型拟合结果

  • 使用线性回归预测未知值

  • 生成拟合模型的诊断图

  • 使用 lm 函数拟合多项式回归模型

  • 使用 rlm 函数拟合稳健线性回归模型

  • 研究 SLID 数据上的线性回归案例

  • 应用高斯模型进行广义线性回归

  • 应用泊松模型进行广义线性回归

  • 应用二项式模型进行广义线性回归

  • 将广义加性模型拟合到数据中

  • 可视化广义加性模型

  • 诊断广义加性模型

简介

回归是一种监督学习方法，用于建模和分析一个因变量（响应变量）与一个或多个自变量（预测变量）之间的关系。可以使用回归来构建预测模型，首先用于找到具有最小拟合值平方误差的最佳拟合模型。然后可以将拟合模型进一步应用于数据进行连续值预测。

有许多种类的回归。如果只有一个预测变量，响应变量和独立变量之间的关系是线性的，我们可以应用线性模型。然而，如果有多个预测变量，应使用多重线性回归方法。当关系是非线性的，可以使用非线性模型来模拟预测变量和响应变量之间的关系。

在本章中，我们介绍如何使用lm函数将线性模型拟合到数据中。接下来，对于非正态高斯模型（例如泊松或二项式）的分布，我们使用带有适当对应数据分布的链接函数的glm函数。最后，我们介绍如何使用gam函数将广义加性模型拟合到数据中。

使用 lm 函数拟合线性回归模型

回归中最简单的模型是线性回归，当只有一个预测变量，响应变量和独立变量之间的关系是线性的时，最适合使用。在 R 中，我们可以使用lm函数将线性模型拟合到数据中。

准备工作

我们需要准备一个预测变量和一个响应变量，以及两个变量之间的线性关系的数据。

如何做...

执行以下步骤以使用lm进行线性回归：

  1. 您应该安装car包并加载其库：

    > install.packages("car")
> library(car)

  2. 从包中，您可以加载Quartet数据集：

    > data(Quartet)

  3. 您可以使用str函数来显示Quartet数据集的结构：

    > str(Quartet)
'data.frame':   11 obs. of  6 variables:
 $ x : int  10 8 13 9 11 14 6 4 12 7 ...
 $ y1: num  8.04 6.95 7.58 8.81 8.33 ...
 $ y2: num  9.14 8.14 8.74 8.77 9.26 8.1 6.13 3.1 9.13 7.26 ...
 $ y3: num  7.46 6.77 12.74 7.11 7.81 ...
 $ x4: int  8 8 8 8 8 8 8 19 8 8 ...
 $ y4: num  6.58 5.76 7.71 8.84 8.47 7.04 5.25 12.5 5.56 7.91 ...

  4. 使用plot函数绘制 x 和 y 变量的散点图，并通过lm和abline函数添加拟合线：

    > plot(Quartet$x, Quartet$y1)
> lmfit = lm(y1~x, Quartet) 
> abline(lmfit, col="red")

    如何做...

    使用 lm 拟合的简单回归图

  5. 要查看拟合模型，执行以下操作：

    > lmfit

Call:
lm(formula = y1 ~ x, data = Quartet)

Coefficients:
(Intercept)            x 
 3.0001       0.5001

它是如何工作的...

回归模型具有response ~ terms形式，其中response是响应向量，terms是一系列指定预测器的项。我们可以用一个简单的回归模型公式y=α+βx来表示，其中α是截距，而斜率β描述了x变化时y的变化。通过使用最小二乘法，我们可以估计如何工作...和如何工作...（其中如何工作...表示y的平均值，如何工作...表示x的平均值）。

要执行线性回归，我们首先准备具有预测变量和响应变量之间线性关系的资料。在这个例子中，我们从 car 包中加载了 Anscombe 的四重奏数据集。在数据集中，x和y1变量之间存在线性关系，我们为这些变量准备了一个散点图。为了生成回归线，我们使用 lm 函数生成两个变量的模型。进一步，我们使用abline在图上绘制回归线。根据前面的截图，回归线说明了x和y1变量的线性关系。我们可以看到，拟合模型的系数显示了截距等于 3.0001，系数等于 0.5001。因此，我们可以使用截距和系数来推断响应值。例如，我们可以推断当x等于 3 时的响应值等于 4.5103（3 * 0.5001 + 3.0001）。

还有更多...

除了lm函数外，您还可以使用lsfit函数执行简单线性回归。例如：

> plot(Quartet$x, Quartet$y1)
> lmfit2 = lsfit(Quartet$x,Quartet$y1)
> abline(lmfit2, col="red")

还有更多...

使用 lsfit 函数拟合的简单回归。

总结线性模型拟合

summary函数可以用来获取拟合模型的格式化系数、标准误差、自由度和其他总结信息。本食谱介绍了如何通过使用summary函数来获取模型的整体信息。

准备工作

您需要完成之前的食谱，通过从四重奏中计算x和y1变量的线性模型，并将拟合模型分配给lmfit变量。

如何做...

执行以下步骤以总结线性模型拟合：

  1. 计算拟合模型的详细总结：

    > summary(lmfit)

Call:
lm(formula = y1 ~ x)

Residuals:
 Min       1Q   Median       3Q      Max 
-1.92127 -0.45577 -0.04136  0.70941  1.83882 

Coefficients:
 Estimate Std. Error t value Pr(>|t|) 
(Intercept)   3.0001     1.1247   2.667  0.02573 * 
Quartet$x     0.5001     0.1179   4.241  0.00217 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.237 on 9 degrees of freedom
Multiple R-squared:  0.6665,    Adjusted R-squared:  0.6295 
F-statistic: 17.99 on 1 and 9 DF,  p-value: 0.00217

如何工作...

summary函数是一个通用函数，用于生成总结统计量。在这种情况下，它计算并返回拟合线性模型的总结统计量列表。这里，它将输出诸如残差、系数标准误差、R 平方、f 统计量和自由度等信息。在Call部分，显示了用于生成拟合模型的函数。在Residuals部分，它提供了分布的快速总结（最小值、1Q、中位数、3Q、最大值）。

在 系数 部分，每个系数都是一个高斯随机变量。在此部分中，估计值 代表变量的均值分布；标准误差 显示变量的标准误差；t 值是 估计值 除以 标准误差，而 p 值表示得到大于 t 值的概率。在此样本中，截距（0.002573）和 x（0.00217）的 p 值都有 95% 的置信水平。

残差标准误差输出残差的方差，而自由度表示训练样本中的观测值与模型中使用的数量之间的差异。多重判定系数是通过除以平方和得到的。一个人可以使用判定系数来衡量数据与回归线拟合的接近程度。通常，判定系数越高，模型与数据的拟合度越好。然而，这并不一定意味着回归模型是充分的。这意味着您可能得到一个低判定系数的好模型，或者您可能得到一个高判定系数的差模型。由于多重判定系数忽略了自由度，计算出的分数是有偏的。为了使计算公平，调整判定系数（0.6295）使用无偏估计，并将略低于多重判定系数（0.6665）。F 统计量是通过在模型上执行 f 检验得到的。一个等于 0.00217（< 0.05）的 p 值拒绝零假设（变量之间没有线性相关性）并表明观察到的 F 值大于临界 F 值。换句话说，结果表明变量之间存在显著的线性正相关。

相关阅读

  • 要获取拟合模型参数的更多信息，您可以使用 help 函数或 ? 来查看帮助页面：

    > ?summary.lm

  • 或者，您可以使用以下函数来显示模型的属性：

    >  coefficients(lmfit) # Extract model coefficients
>  confint(lmfit, level=0.95)  # Computes confidence intervals for model parameters.
>  fitted(lmfit) # Extract model fitted values
>  residuals(lmfit) # Extract model residuals 
>  anova(lmfit) # Compute analysis of variance tables for fitted model object
>  vcov(lmfit) # Calculate variance-covariance matrix for a fitted model object
>  influence(lmfit) # Diagnose quality of regression fits

使用线性回归预测未知值

使用拟合的回归模型，我们可以将模型应用于预测未知值。对于回归模型，我们可以使用预测区间和置信区间来表示预测的精度。在下面的配方中，我们介绍如何在两种测量下预测未知值。

准备工作

您需要完成上一个配方，通过从 quartet 数据集中计算 x 和 y1 变量的线性模型来完成。

如何操作...

按以下步骤使用线性回归预测值：

  1. 使用 x 和 y1 变量拟合线性模型：

    > lmfit = lm(y1~x, Quartet)

  2. 将要预测的值分配到 newdata：

    > newdata = data.frame(x = c(3,6,15))

  3. 使用置信水平设置为 0.95 的置信区间来计算预测结果：

    > predict(lmfit, newdata, interval="confidence", level=0.95)
 fit      lwr       upr
1  4.500364 2.691375  6.309352
2  6.000636 4.838027  7.163245
3 10.501455 8.692466 12.310443

  4. 使用此预测区间来计算预测结果：

    > predict(lmfit, newdata, interval="predict")
 fit      lwr       upr
1  4.500364 1.169022  7.831705
2  6.000636 2.971271  9.030002
3 10.501455 7.170113 13.832796

工作原理...

我们首先使用 x 和 y1 变量构建一个线性拟合模型。接下来，我们将要预测的值分配到一个数据框 newdata 中。重要的是要注意，生成的模型形式为 y1 ~ x。

接下来，我们通过在interval参数中指定confidence来使用置信区间计算预测结果。从第 1 行的输出中，我们得到x=3输入的拟合y1等于4.500364，以及x=3的y1均值的 95%置信区间（在level参数中设置为 0.95）介于2.691375和6.309352之间。此外，第 2 行和第 3 行给出了x=6和x=15输入的y1预测结果。

接下来，我们通过在interval参数中指定prediction来使用预测区间计算预测结果。从第 1 行的输出中，我们可以看到x=3输入的拟合y1等于4.500364，以及x=3的y1的 95%预测区间介于1.169022和7.831705之间。第 2 行和第 3 行输出了x=6和x=15输入的y1预测结果。

参考以下内容

生成拟合模型的诊断图

诊断是评估回归假设的方法，可用于确定拟合模型是否充分代表数据。在以下菜谱中，我们介绍如何通过使用诊断图来诊断回归模型。

准备工作

您需要通过从四元组中计算 x 和 y1 变量的线性模型来完成前面的菜谱，并将模型分配给lmfit变量。

如何操作...

执行以下步骤以生成拟合模型的诊断图：

  1. 绘制回归模型的诊断图：

    > par(mfrow=c(2,2))
> plot(lmfit)

    如何操作...

    回归模型的诊断图

工作原理...

绘图函数生成回归模型的四个诊断图：

  • 左上角的图显示了残差与拟合值之间的关系。在图中，残差代表一个点到回归线的垂直距离。如果所有点都正好落在回归线上，所有残差都将正好落在虚线灰线上。图中的红色线是关于残差的光滑曲线，如果所有点都正好落在回归线上，红色线的位置应该正好与虚线灰线相匹配。

  • 右上角显示了残差的正态分布。此图验证了残差正态分布的假设。因此，如果残差是正态分布的，它们应该正好位于灰色虚线上。

  • 左下角的尺度-位置图用于测量标准化残差的平方根与拟合值之间的关系。因此，如果所有点都位于回归线上，则y的值应该接近零。由于假设残差的方差不会显著改变分布，如果假设正确，红色线应该相对平坦。

  • 右下角的图显示了标准化残差与杠杆作用的关系。杠杆作用是衡量每个数据点如何影响回归的度量。它是回归质心距离和隔离程度的度量（通过是否有邻居来衡量）。此外，您还可以找到受高杠杆和大量残差影响的 Cook 距离的轮廓。您可以使用此方法来衡量如果删除一个点，回归将如何改变。红色线在标准化残差方面很平滑。对于完美的拟合回归，红色线应该接近虚线，Cook 距离没有超过 0.5 的点。

还有更多...

要查看更多的诊断图功能，可以使用help函数获取更多信息：

> ?plot.lm

为了发现是否存在具有大 Cook 距离的点，可以使用cooks.distance函数计算每个点的 Cook 距离，并通过可视化分析距离的分布：

> plot(cooks.distance(lmfit))

还有更多...

Cook 距离图

在这种情况下，如果索引为 3 的点显示的 Cook 距离大于其他点，可以调查这个点是否可能是异常值。

使用 lm 拟合多项式回归模型

一些预测变量和响应变量可能存在非线性关系，它们的关系可以建模为n次多项式。在这个菜谱中，我们将介绍如何使用lm和poly函数处理多项式回归。

准备工作

准备包含可以建模为n次多项式的预测变量和响应变量之间关系的数据库。在这个菜谱中，我们将继续使用car包中的Quartet数据库。

如何操作...

执行以下步骤以使用lm拟合多项式回归模型：

  1. 首先，我们绘制x和y2变量的散点图：

    > plot(Quartet$x, Quartet$y2)

    如何操作...

    变量 x 和 y2 的散点图

  2. 您可以通过指定参数中的2来应用poly函数：

    > lmfit = lm(Quartet$y2~poly(Quartet$x,2))
> lines(sort(Quartet$x), lmfit$fit[order(Quartet$x)], col = "red")

    如何操作...

    变量 x 和 y2 的回归图的一个二次拟合示例

工作原理

我们可以在公式工作原理中说明二次多项式回归模型，其中α是截距，而β表示回归系数。

在前面的截图（步骤 1）中，x和y2变量的散点图不呈线性关系，而是显示一个向下凹的曲线（或向上凸）的形状，转折点在 x=11。为了建模非线性关系，我们使用带有 2 作为参数的poly函数来拟合独立的x变量和依赖的y2变量。截图中的红线显示模型完美地拟合了数据。

更多...

您还可以使用一个独立变量等于一阶x变量和二阶x变量公式的二阶多项式模型：

> plot(Quartet$x, Quartet$y2)
> lmfit = lm(Quartet$y2~ I(Quartet$x)+I(Quartet$x²))

使用 rlm 拟合稳健线性回归模型

数据集中的异常值会将回归线从主流移动开去。除了移除它之外，我们还可以对包含异常值的数据集应用稳健线性回归。在本配方中，我们介绍了如何将rlm应用于对包含异常值的数据集进行稳健线性回归。

准备工作

准备一个可能将回归线从主流移动开去的异常值数据集。在这里，我们使用从先前配方中加载的Quartet数据集。

如何操作...

按照以下步骤使用rlm拟合稳健线性回归模型：

  1. 生成x变量对y3的散点图：

    > plot(Quartet$x, Quartet$y3)

    如何操作...

    变量 x 和 y3 的散点图

  2. 接下来，您应该首先导入MASS库。然后，您可以使用rlm函数来拟合模型，并使用abline函数可视化拟合的线：

    > library(MASS)
> lmfit = rlm(Quartet$y3~Quartet$x)
> abline(lmfit, col="red")

    如何操作...

    对变量 x 和 y3 的稳健线性回归

工作原理

根据前面的截图（步骤 1），您可能会遇到包含远离主流的异常值的数据集。为了消除异常值的影响，我们展示了如何应用稳健线性回归（rlm）来拟合数据。在第二个截图（步骤 2）中，稳健回归线忽略了异常值并匹配了主流。

更多...

为了看到异常值如何将回归线从主流移动开去，您可以将此配方中使用的rlm函数替换为lm，并重新绘制图形：

> plot(Quartet$x, Quartet$y3)
> lmfit = lm(Quartet$y3~Quartet$x)
> abline(lmfit, col="red")

更多...

变量 x 和 y3 的线性回归

显然，异常值（x=13）将回归线从主流移动开去。

研究 SLID 数据上的线性回归案例

为了总结上一节的内容，我们使用线性回归探索更复杂的数据。在本配方中，我们展示了如何将线性回归应用于分析劳动与收入动态调查（SLID）数据集。

准备工作

检查是否已安装和加载car库，因为它需要访问 SLID 数据集。

如何操作...

按照以下步骤在 SLID 数据上执行线性回归：

  1. 您可以使用str函数来获取数据的概览：

    > str(SLID)
'data.frame':  7425 obs. of  5 variables:
 $ wages    : num  10.6 11 NA 17.8 NA ...
 $ education: num  15 13.2 16 14 8 16 12 14.5 15 10 ...
 $ age      : int  40 19 49 46 71 50 70 42 31 56 ...
 $ sex      : Factor w/ 2 levels "Female","Male": 2 2 2 2 2 1 1 1 2 1 ...
 $ language : Factor w/ 3 levels "English","French",..: 1 1 3 3 1 1 1 1 1 1 ..

  2. 首先，我们将变量工资与语言、年龄、教育和性别可视化：

    > par(mfrow=c(2,2))
> plot(SLID$wages ~ SLID$language)
> plot(SLID$wages ~ SLID$age)
> plot(SLID$wages ~ SLID$education)
> plot(SLID$wages ~ SLID$sex)

    如何操作...

    工资与多种组合的绘图

  3. 然后，我们可以使用lm来拟合模型：

    > lmfit = lm(wages ~ ., data = SLID)

  4. 您可以通过summary函数检查拟合模型的摘要：

    > summary(lmfit)

Call:
lm(formula = wages ~ ., data = SLID)

Residuals:
 Min      1Q  Median      3Q     Max 
-26.062  -4.347  -0.797   3.237  35.908 

Coefficients:
 Estimate Std. Error t value Pr(>|t|) 
(Intercept)    -7.888779   0.612263 -12.885   <2e-16 ***
education       0.916614   0.034762  26.368   <2e-16 ***
age             0.255137   0.008714  29.278   <2e-16 ***
sexMale         3.455411   0.209195  16.518   <2e-16 ***
languageFrench -0.015223   0.426732  -0.036    0.972 
languageOther   0.142605   0.325058   0.439    0.661 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.6 on 3981 degrees of freedom
 (3438 observations deleted due to missingness)
Multiple R-squared:  0.2973,	Adjusted R-squared:  0.2964 
F-statistic: 336.8 on 5 and 3981 DF,  p-value: < 2.2e-16

  5. 删除language属性，并使用lm函数重新拟合模型：

    > lmfit = lm(wages ~ age + sex + education, data = SLID)
> summary(lmfit)

Call:
lm(formula = wages ~ age + sex + education, data = SLID)

Residuals:
 Min      1Q  Median      3Q     Max 
-26.111  -4.328  -0.792   3.243  35.892 

Coefficients:
 Estimate Std. Error t value Pr(>|t|) 
(Intercept) -7.905243   0.607771  -13.01   <2e-16 ***
age          0.255101   0.008634   29.55   <2e-16 ***
sexMale      3.465251   0.208494   16.62   <2e-16 ***
education    0.918735   0.034514   26.62   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.602 on 4010 degrees of freedom
 (3411 observations deleted due to missingness)
Multiple R-squared:  0.2972,	Adjusted R-squared:  0.2967 
F-statistic: 565.3 on 3 and 4010 DF,  p-value: < 2.2e-16

  6. 然后，我们可以绘制lmfit的诊断图：

    > par(mfrow=c(2,2))
> plot(lmfit)

    如何操作...

    拟合模型的诊断图

  7. 接下来，我们对工资取对数并重新绘制诊断图：

    > lmfit = lm(log(wages) ~ age + sex + education, data = SLID)
> plot(lmfit)

    如何操作...

    调整后拟合模型的诊断图

  8. 接下来，您可以使用vif函数诊断回归模型的多重共线性：

    > vif(lmfit)
 age       sex education 
 1.011613  1.000834  1.012179 
> sqrt(vif(lmfit)) > 2
 age       sex education 
 FALSE     FALSE     FALSE

  9. 然后，您可以安装并加载lmtest包，并使用bptest函数诊断回归模型的异方差性：

    > install.packages("lmtest")
> library(lmtest)
> bptest(lmfit)

 studentized Breusch-Pagan test

data:  lmfit
BP = 29.0311, df = 3, p-value = 2.206e-06

  10. 最后，您可以安装并加载rms包。然后，您可以使用robcov来纠正标准误差：

    > install.packages("rms")
> library(rms)
> olsfit = ols(log(wages) ~ age + sex + education, data= SLID, x= TRUE, y= TRUE)
> robcov(olsfit)

Linear Regression Model

ols(formula = log(wages) ~ age + sex + education, data = SLID, 
 x = TRUE, y = TRUE)

Frequencies of Missing Values Due to Each Variable
log(wages)        age        sex  education 
 3278          0          0        249 

 Model Likelihood     Discrimination 
 Ratio Test           Indexes 
Obs     4014    LR chi2   1486.08    R2       0.309 
sigma 0.4187    d.f.            3    R2 adj   0.309 
d.f.    4010    Pr(> chi2) 0.0000    g        0.315 

Residuals

 Min       1Q   Median       3Q      Max 
-2.36252 -0.27716  0.01428  0.28625  1.56588 

 Coef   S.E.   t     Pr(>|t|)
Intercept 1.1169 0.0387 28.90 <0.0001 
age       0.0176 0.0006 30.15 <0.0001 
sex=Male  0.2244 0.0132 16.96 <0.0001 
education 0.0552 0.0022 24.82 <0.0001

它是如何工作的...

本菜谱演示了如何在 SLID 数据集上执行线性回归分析。首先，我们加载 SLID 数据并通过使用str函数显示其结构。从数据结构中，我们知道有四个独立变量将影响因变量工资。

接下来，我们通过可视化探索每个独立变量与因变量wages之间的关系；可视化结果如图所示（步骤 2）。在此截图的左上角，您可以找到三种不同语言与工资的箱线图；语言与工资之间的相关性不明显。截图右上角显示，年龄似乎与因变量wages呈正相关。在截图的左下角，显示教育似乎也与工资呈正相关。最后，截图右下角的箱线图显示，男性的工资略高于女性。

接下来，我们将除了wages以外的所有属性拟合到模型中作为预测变量。通过总结模型，可以看出教育、年龄和性别具有显著性（p-value < 0.05）。因此，我们删除了不显著的language属性（其 p 值大于 0.05），并在线性模型中将三个独立变量（教育、性别和年龄）与因变量（wages）拟合。这相应地将 f 统计量从 336.8 提高到 565.3。

接下来，我们生成拟合模型的诊断图。在诊断图中，所有四个图都表明回归模型遵循回归假设。然而，从残差与拟合值和尺度-位置图中可以看出，较小拟合值的残差偏向回归模型。由于工资跨越几个数量级，为了诱导对称性，我们对工资应用对数变换，并将数据重新拟合到回归模型中。现在，残差与拟合值图和尺度-位置图的红线与灰色虚线更接近。

接下来，我们想测试模型中是否存在多重共线性。当预测变量与其他变量高度相关时，就会发生多重共线性。如果模型中存在多重共线性，你可能会看到一些变量具有高的 R-squared 值，但显示为不显著的变量。为了检测多重共线性，我们可以使用vif函数计算线性模型和广义线性模型的方差膨胀因子和广义方差膨胀因子。如果存在多重共线性，我们应该找到方差膨胀因子平方根大于 2 的预测变量。然后，我们可以移除冗余预测变量或使用主成分分析将预测变量转换为较小的一组不相关成分。

最后，我们想测试模型中是否存在异方差性。在讨论异方差性的定义之前，我们首先要知道，在经典假设中，普通回归模型假设误差的方差在观测值之间是恒定的或同质的。相反，异方差性意味着方差在观测值之间是不等的。因此，异方差性可能会偏向我们估计的标准误差，从而误导假设的检验。为了检测和测试异方差性，我们可以在lmtest包中的bptest函数执行异方差性的Breusch-Pagan检验。在这种情况下，p 值显示为 2.206e-06（<0.5），这拒绝了同方差性的零假设（没有异方差性）。在这里，它意味着参数估计的标准误差是不正确的。然而，我们可以使用rms包中的robcov来使用稳健标准误差纠正标准误差（不消除异方差性）并提高真正显著参数的显著性。然而，由于它只接受rms系列中的拟合模型作为输入，我们必须先拟合普通最小二乘模型。

参考信息

  • 如需了解更多关于 SLID 数据集的信息，您可以使用help函数查看相关文档：

    >  ?SLID

应用高斯模型进行广义线性回归

广义线性模型（GLM）是线性回归的推广，可以包括连接函数以进行线性预测。作为默认设置，glm 的家族对象是高斯，这使得 glm 函数的表现与 lm 完全相同。在本例中，我们首先演示如何使用 glm 函数将模型拟合到数据中，然后展示具有高斯模型的 glm 与 lm 的表现完全相同。

准备工作

检查 car 库是否已安装并加载，因为我们需要从该包中获取 SLID 数据集。

如何操作...

执行以下步骤以拟合具有高斯模型的广义线性回归模型：

  1. 将独立变量 age、sex 和 education 以及因变量 wages 拟合到 glm：

    > lmfit1 = glm(wages ~ age + sex + education, data = SLID, family=gaussian)
> summary(lmfit1)

Call:
glm(formula = wages ~ age + sex + education, family = gaussian, 
 data = SLID)

Deviance Residuals: 
 Min       1Q   Median       3Q      Max 
-26.111   -4.328   -0.792    3.243   35.892 

Coefficients:
 Estimate Std. Error t value Pr(>|t|) 
(Intercept) -7.905243   0.607771  -13.01   <2e-16 ***
age          0.255101   0.008634   29.55   <2e-16 ***
sexMale      3.465251   0.208494   16.62   <2e-16 ***
education    0.918735   0.034514   26.62   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Gaussian family taken to be 43.58492)

 Null deviance: 248686  on 4013  degrees of freedom
Residual deviance: 174776  on 4010  degrees of freedom
 (3411 observations deleted due to missingness)
AIC: 26549

Number of Fisher Scoring iterations: 2

  2. 将独立变量 age、sex 和 education 以及因变量 wages 拟合到 lm：

    > lmfit2 = lm(wages ~ age + sex + education, data = SLID)
> summary(lmfit2)

Call:
lm(formula = wages ~ age + sex + education, data = SLID)

Residuals:
 Min      1Q  Median      3Q     Max 
-26.111  -4.328  -0.792   3.243  35.892 

Coefficients:
 Estimate Std. Error t value Pr(>|t|) 
(Intercept) -7.905243   0.607771  -13.01   <2e-16 ***
age          0.255101   0.008634   29.55   <2e-16 ***
sexMale      3.465251   0.208494   16.62   <2e-16 ***
education    0.918735   0.034514   26.62   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.602 on 4010 degrees of freedom
 (3411 observations deleted due to missingness)
Multiple R-squared:  0.2972,	Adjusted R-squared:  0.2967 
F-statistic: 565.3 on 3 and 4010 DF,  p-value: < 2.2e-16

  3. 使用 anova 比较两个拟合模型：

    > anova(lmfit1, lmfit2)
Analysis of Deviance Table

Model: gaussian, link: identity

Response: wages

Terms added sequentially (first to last)

 Df Deviance Resid. Df Resid. Dev
NULL                       4013     248686
age        1    31953      4012     216733
sex        1    11074      4011     205659
education  1    30883      4010     174776

工作原理...

glm 函数以与 lm 函数类似的方式拟合数据模型。唯一的区别是您可以在参数 family 中指定不同的连接函数（您可以在控制台中使用 ?family 来查找不同类型的连接函数）。在本例中，我们首先将独立变量 age、sex 和 education 以及因变量 wages 输入到 glm 函数中，并将构建的模型赋值给 lmfit1。您可以使用构建的模型进行进一步的预测。

接下来，为了确定具有高斯模型的 glm 是否与 lm 完全相同，我们将独立变量 age、sex 和 education 以及因变量 wages 拟合到 lm 模型中。通过将 summary 函数应用于两个不同的模型，它揭示了两个输出摘要的残差和系数完全相同。

最后，我们使用 anova 函数进一步比较两个拟合模型。anova 函数的结果显示，两个模型相似，具有相同的 残差自由度（Res.DF）和 残差平方和（RSS Df）。

参见

  • 对于广义线性模型与线性模型的比较，您可以参考 Venables, W. N., & Ripley, B. D. (2002). Modern applied statistics with S. Springer。

应用泊松模型进行广义线性回归

广义线性模型允许响应变量具有除正态分布（高斯）之外的错误分布模型。在本例中，我们展示了如何将泊松作为 glm 中的家族对象应用于计数数据。

准备工作

本任务的先决条件是准备计数数据，所有输入数据值均为整数。

如何操作...

执行以下步骤以拟合具有泊松模型的广义线性回归模型：

  1. 加载 warpbreaks 数据，并使用 head 查看前几行：

    > data(warpbreaks)
> head(warpbreaks)
 breaks wool tension
1     26    A       L
2     30    A       L
3     54    A       L
4     25    A       L
5     70    A       L
6     52    A       L

  2. 我们将泊松作为家族对象应用于独立变量 tension 和因变量 breaks：

    > rs1 = glm(breaks ~ tension, data=warpbreaks, family="poisson")
> summary(rs1)

Call:
glm(formula = breaks ~ tension, family = "poisson", data = warpbreaks)

Deviance Residuals: 
 Min       1Q   Median       3Q      Max 
-4.2464  -1.6031  -0.5872   1.2813   4.9366 

Coefficients:
 Estimate Std. Error z value Pr(>|z|) 
(Intercept)  3.59426    0.03907  91.988  < 2e-16 ***
tensionM    -0.32132    0.06027  -5.332 9.73e-08 ***
tensionH    -0.51849    0.06396  -8.107 5.21e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Poisson family taken to be 1)

 Null deviance: 297.37  on 53  degrees of freedom
Residual deviance: 226.43  on 51  degrees of freedom
AIC: 507.09

Number of Fisher Scoring iterations: 4

工作原理...

在泊松分布的假设下，计数数据可以拟合到对数线性模型。在这个例子中，我们首先从warpbreaks数据集中加载了一个样本计数数据，其中包含了每台织机的断裂次数。接下来，我们使用glm函数，将断裂次数作为因变量，tension作为自变量，泊松作为家族对象。最后，我们使用摘要函数查看拟合的对数线性模型。

参考以下内容

  • 要了解更多关于泊松模型与计数数据的关系，您可以参考Cameron, A. C., & Trivedi, P. K. (2013). Regression analysis of count data (No. 53). Cambridge university press.

应用广义线性回归的二项模型

对于二元因变量，可以在glm函数中将二项模型作为家族对象应用。

准备工作

此任务的先决条件是准备一个二元因变量。在这里，我们使用vs变量（V 型发动机或直列发动机）作为因变量。

如何操作...

执行以下步骤以使用二项模型拟合广义线性回归模型：

  1. 首先，我们检查mtcars数据集中vs的前六个元素：

    > head(mtcars$vs)
[1] 0 0 1 1 0 1

  2. 我们使用glm函数，将binomial作为家族对象应用：

    > lm1 = glm(vs ~ hp+mpg+gear,data=mtcars, family=binomial)
> summary(lm1)

Call:
glm(formula = vs ~ hp + mpg + gear, family = binomial, data = mtcars)

Deviance Residuals: 
 Min        1Q    Median        3Q       Max 
-1.68166  -0.23743  -0.00945   0.30884   1.55688 

Coefficients:
 Estimate Std. Error z value Pr(>|z|) 
(Intercept) 11.95183    8.00322   1.493   0.1353 
hp          -0.07322    0.03440  -2.129   0.0333 *
mpg          0.16051    0.27538   0.583   0.5600 
gear        -1.66526    1.76407  -0.944   0.3452 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

 Null deviance: 43.860  on 31  degrees of freedom
Residual deviance: 15.651  on 28  degrees of freedom
AIC: 23.651

Number of Fisher Scoring iterations: 7

工作原理...

在二元数据中，每个响应值的观测值被编码为0或1。将二元数据拟合到回归模型需要二项分布函数。在这个例子中，我们首先从mtcars数据集中加载二元因变量vs。vs适合二项模型，因为它包含二元数据。接下来，我们使用glm函数通过指定binomial作为家族对象将模型拟合到二元数据。最后，通过查看摘要，我们可以获得拟合模型的描述。

参考以下内容

  • 如果您仅在参数中指定家族对象，您将使用默认的链接来拟合模型。但是，要使用替代链接函数，您可以添加一个链接参数。例如：

    > lm1 = glm(vs ~ hp+mpg+gear,data=mtcars, family=binomial(link="probit"))

  • 如果您想知道可以使用多少个替代链接，请通过帮助功能参考家族文档：

     > ?family

将广义加性模型拟合到数据中

广义加性模型（GAM），用于拟合广义加性模型，可以被视为 GLM 的半参数扩展。虽然 GLM 假设因变量和自变量之间存在线性关系，但 GAM 根据数据的局部行为拟合模型。因此，GAM 具有处理因变量和自变量之间高度非线性关系的能力。在下面的步骤中，我们介绍如何使用广义加性模型拟合回归。

准备工作

我们需要准备一个包含变量的数据框，其中一个变量是响应变量，其他变量可能是预测变量。

如何操作...

执行以下步骤以将广义加性模型拟合到数据中：

  1. 首先，加载包含gam函数的mgcv包：

    > install.packages("mgcv")
> library(mgcv)

  2. 然后，安装MASS包并加载Boston数据集：

    > install.packages("MASS")
> library(MASS)
> attach(Boston)
> str(Boston)

  3. 使用gam拟合回归：

    > fit = gam(dis ~ s(nox))

  4. 获取拟合模型的摘要信息：

    > summary(fit)
Family: gaussian 
Link function: identity 

Formula:
dis ~ s(nox)

Parametric coefficients:
 Estimate Std. Error t value Pr(>|t|) 
(Intercept)  3.79504    0.04507   84.21   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
 edf Ref.df   F p-value 
s(nox) 8.434  8.893 189  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.768   Deviance explained = 77.2%
GCV = 1.0472  Scale est. = 1.0277    n = 506

工作原理

GAM 旨在通过估计与因变量y通过链接函数连接的预测者的非参数函数来最大化从各种分布中预测因变量y。GAM 的概念工作原理，其中指定了指数族E用于y，以及g链接函数；f表示预测者的链接函数。

gam函数包含在mgcv包中，因此，首先安装此包并将其加载到 R 会话中。接下来，从MASS包中加载波士顿数据集（波士顿郊区的住房价值）。从数据集中，我们使用dis（到五个波士顿就业中心的加权平均距离）作为因变量，nox（氮氧化物浓度）作为自变量，然后将它们输入到gam函数中以生成拟合模型。

与glm类似，gam允许用户总结gam拟合。从总结中，可以找到参数参数、平滑项的显著性以及其他有用信息。

相关内容

  • 除了gam之外，mgcv包还提供另一个用于大数据集的广义加性模型bam。bam包与gam非常相似，但使用更少的内存，效率相对更高。请使用help函数获取有关此模型的更多信息：

     > ? bam

  • 有关 R 中广义加性模型的更多信息，请参阅 Wood, S. (2006). 广义加性模型：R 语言导论. CRC 出版社。

可视化广义加性模型

在这个配方中，我们演示了如何将gam拟合回归线添加到散点图中。此外，我们使用plot函数可视化gam拟合。

准备工作

通过将gam拟合模型分配给fit变量来完成前面的配方。

如何做...

执行以下步骤以可视化广义加性模型：

  1. 使用nox和dis变量生成散点图：

    > plot(nox, dis)

    如何做...

    变量 nox 对 dis 的散点图

  2. 将回归添加到散点图：

    > x = seq(0, 1, length = 500)
> y = predict(fit, data.frame(nox = x))
> lines(x, y, col = "red", lwd = 2)

    如何做...

    散点图上的gam拟合回归

  3. 或者，您可以使用plot函数绘制拟合模型：

    > plot(fit)

    如何做...

    拟合gam的图

工作原理...

要可视化拟合回归，我们首先使用dis和nox变量生成散点图。然后，我们生成x轴的序列，并通过在拟合模型fit上使用predict函数来响应y。最后，我们使用lines函数将回归线添加到散点图中。

除了在散点图上使用线条添加拟合回归线外，gam还有一个plot函数来可视化包含置信区域的拟合回归线。为了阴影置信区域，我们在函数中指定shade = TRUE。

更多...

vis.gam函数用于生成gam模型预测的透视或等高线图视图。这有助于观察响应变量如何与两个预测变量相互作用。以下是在Boston数据集上的等高线图示例：

> fit2=gam(medv~crim+zn+crim:zn, data=Boston)
> vis.gam(fit2)

更多...

由 vis.gam 生成的样本等高线图

诊断广义加性模型

GAM 还提供了关于拟合过程和广义加性模型结果的诊断信息。在这个菜谱中，我们展示了如何通过gam.check函数绘制诊断图。

准备工作

确保之前的菜谱已完成，并将gam拟合模型分配给fit变量。

如何做...

执行以下步骤以诊断广义加性模型：

  1. 使用gam.check在拟合模型上生成诊断图：

    > gam.check(fit)

Method: GCV   Optimizer: magic
Smoothing parameter selection converged after 7 iterations.
The RMS GCV score gradient at convergence was 8.79622e-06 .
The Hessian was positive definite.
The estimated model rank was 10 (maximum possible: 10)
Model rank =  10 / 10 

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

 k'   edf k-index p-value
s(nox) 9.000 8.434   0.397       0

    如何做...

    拟合的 gam 诊断图

它是如何工作的...

gam.check 函数首先生成平滑参数估计收敛信息。在这个例子中，平滑参数，GCV/UBRE（广义交叉验证/无偏风险估计器）得分在七次迭代后收敛。GCV/UBRE 函数在最小值处的平均绝对梯度为 8.79622e-06，估计的秩为10。维度检查是为了测试平滑函数的基函数维度是否足够。从这个例子中可以看出，低 p 值表明 k 设置得太低。可以通过指定参数 k，通过拟合gam到数据中来调整平滑的维度选择。

除了提供有关平滑参数估计收敛的信息外，该函数还返回四个诊断图。截图中的图的上左部分显示了一个分位数比较图。此图有助于识别异常值和重尾。图的右上部分显示了残差与线性预测器的对比，这在寻找非恒定误差方差时很有用。图的左下部分显示了残差的直方图，有助于检测非正态性。图的右下部分显示了响应值与拟合值的关系。

更多...

您可以通过help函数获取有关gam.check的更多信息。特别是，这包括平滑参数估计收敛的详细说明和四个返回的图：

> ?gam.check

此外，可以通过以下命令获取choose.k的更多信息：

> ?choose.k

第五章. 分类（I）- 树、懒惰和概率

在本章中，我们将介绍以下步骤：

  • 准备训练和测试数据集

  • 使用递归分区树构建分类模型

  • 可视化递归分区树

  • 测量递归分区树的预测性能

  • 剪枝递归分区树

  • 使用条件推断树构建分类模型

  • 可视化条件推断树

  • 测量条件推断树的预测性能

  • 使用 k 近邻分类器对数据进行分类

  • 使用逻辑回归对数据进行分类

  • 使用朴素贝叶斯分类器对数据进行分类

简介

分类用于根据从训练数据集构建的分类模型识别新观察值的类别（测试数据集），其中类别已经已知。与回归类似，分类被归类为监督学习方法，因为它使用训练数据集的已知答案（标签）来预测测试数据集的答案（标签）。回归和分类之间的主要区别在于，回归用于预测连续值。

与此相反，分类用于识别给定观察值的类别。例如，一个人可能使用回归来根据历史价格预测给定股票的未来价格。然而，应该使用分类方法来预测股票价格是上涨还是下跌。

在本章中，我们将说明如何使用 R 进行分类。我们首先从客户流失数据集中构建训练数据集和测试数据集，然后应用不同的分类方法对客户流失数据集进行分类。在下面的步骤中，我们将介绍使用传统分类树和条件推断树、基于懒惰的算法以及使用基于概率的方法（使用训练数据集构建分类模型），然后使用该模型预测测试数据集的类别（类别标签）。我们还将使用混淆矩阵来衡量性能。

准备训练和测试数据集

构建分类模型需要训练数据集来训练分类模型，并且需要测试数据来验证预测性能。在下面的步骤中，我们将演示如何将电信客户流失数据集分割成训练数据集和测试数据集。

准备工作

在这个步骤中，我们将使用电信客户流失数据集作为输入数据源，并将数据分割成训练数据集和测试数据集。

如何做...

执行以下步骤将客户流失数据集分割成训练数据集和测试数据集：

  1. 您可以从C50包中检索客户流失数据集：

    > install.packages("C50")
> library(C50)
> data(churn)

  2. 使用str读取数据集的结构：

    > str(churnTrain)

  3. 我们可以移除state、area_code和account_length属性，这些属性不适合作为分类特征：

    > churnTrain = churnTrain[,! names(churnTrain) %in% c("state", "area_code", "account_length") ]

  4. 然后，将 70%的数据分割为训练集，30%的数据分割为测试集：

    > set.seed(2)
> ind = sample(2, nrow(churnTrain), replace = TRUE, prob=c(0.7, 0.3))
> trainset = churnTrain[ind == 1,]
> testset = churnTrain[ind == 2,]

  5. 最后，使用dim来探索训练集和测试集的维度：

    > dim(trainset)
[1] 2315   17
> dim(testset)
[1] 1018   17

工作原理...

在本食谱中，我们使用电信流失数据集作为示例数据源。该数据集包含 20 个变量和 3,333 个观测值。我们希望构建一个分类模型来预测客户是否会流失，这对电信公司来说非常重要，因为获取新客户的成本显著高于保留现有客户。

在构建分类模型之前，我们首先需要预处理数据。因此，我们使用变量名churn将流失数据从C50包加载到 R 会话中。由于我们确定属性如state、area_code和account_length对于构建分类模型不是有用的特征，我们移除了这些属性。

在预处理数据后，我们将数据分别分割为训练集和测试集。然后，我们使用一个样本函数随机生成一个包含 70%训练集和 30%测试集的序列，其大小等于观测值的数量。然后，我们使用生成的序列将流失数据集分割为训练集trainset和测试集testset。最后，通过使用dim函数，我们发现 3,333 个观测值中有 2,315 个被归类到训练集trainset，而其他 1,018 个被归类到测试集testset。

更多内容...

您可以将训练集和测试集的分割过程合并到split.data函数中。因此，您可以通过调用此函数并指定参数中的比例和种子来轻松地将数据分割成两个数据集：

> split.data = function(data, p = 0.7, s = 666){
+   set.seed(s)
+   index = sample(1:dim(data)[1])
+   train = data[index[1:floor(dim(data)[1] * p)], ]
+   test = data[index[((ceiling(dim(data)[1] * p)) + 1):dim(data)[1]], ]
+   return(list(train = train, test = test))
+ }

使用递归分割树构建分类模型

分类树使用分割条件根据一个或多个输入变量预测类标签。分类过程从树的根节点开始；在每个节点，过程将检查输入值是否应该根据分割条件递归地继续到右子分支或左子分支，并在遇到决策树的任何叶（终端）节点时停止。在本食谱中，我们将介绍如何将递归分割树应用于客户流失数据集。

准备工作

您需要完成之前的食谱，将流失数据集分割为训练集（trainset）和测试集（testset），并且每个数据集应恰好包含 17 个变量。

如何操作...

执行以下步骤将流失数据集分割为训练集和测试集：

  1. 加载rpart包：

    > library(rpart)

  2. 使用rpart函数构建分类树模型：

    > churn.rp = rpart(churn ~ ., data=trainset)

  3. 输入churn.rp以检索分类树的节点细节：

    > churn.rp

  4. 接下来，使用printcp函数来检查复杂度参数：

    > printcp(churn.rp)

Classification tree:
rpart(formula = churn ~ ., data = trainset)

Variables actually used in tree construction:
[1] international_plan            number_customer_service_calls
[3] total_day_minutes             total_eve_minutes 
[5] total_intl_calls              total_intl_minutes 
[7] voice_mail_plan 

Root node error: 342/2315 = 0.14773

n= 2315 

 CP nsplit rel error  xerror     xstd
1 0.076023      0   1.00000 1.00000 0.049920
2 0.074561      2   0.84795 0.99708 0.049860
3 0.055556      4   0.69883 0.76023 0.044421
4 0.026316      7   0.49415 0.52632 0.037673
5 0.023392      8   0.46784 0.52047 0.037481
6 0.020468     10   0.42105 0.50877 0.037092
7 0.017544     11   0.40058 0.47076 0.035788
8 0.010000     12   0.38304 0.47661 0.035993

  5. 接下来，使用plotcp函数来绘制成本复杂度参数：

    > plotcp(churn.rp)

    如何操作...

    图 1：成本复杂度参数图

  6. 最后，使用summary函数来检查构建的模型：

    > summary(churn.rp)

它是如何工作的...

在这个菜谱中，我们使用rpart包中的递归分割树来构建基于树的分类模型。递归分割树包括两个过程：递归和分割。在决策诱导的过程中，我们必须考虑一个统计评估问题（或者简单地说是一个是/否问题），根据评估结果将数据分割成不同的分区。然后，当我们确定了子节点后，我们可以重复执行分割，直到满足停止标准。

例如，根节点中的数据（如下所示）可以根据f1是否小于X的问题分为两组。如果是，数据被分割到左侧。否则，它被分割到右侧。然后，我们可以继续使用f2是否小于Y的问题来分割左侧数据：

如何操作...

图 2：递归分割树

在第一步中，我们使用library函数加载rpart包。接下来，我们使用churn变量作为分类类别（类别标签）和剩余变量作为输入特征来构建分类模型。

模型构建完成后，你可以输入构建的模型变量名churn.rp来显示树节点细节。在打印的节点细节中，n表示样本大小，loss表示误分类成本，yval表示分类成员（在这种情况下为no或yes），而yprob表示两个类的概率（左值表示达到标签no的概率，右值表示达到标签yes的概率）。

然后，我们使用printcp函数来打印构建的树模型的复杂度参数。从printcp的输出中，应该找到一个复杂度参数 CP 的值，它作为惩罚来控制树的大小。简而言之，CP 值越大，分割的数量（nsplit）就越少。输出值（rel误差）表示当前树的平均偏差除以空树的平均偏差。xerror值表示由 10 折分类估计的相对误差。xstd表示相对误差的标准误差。

为了使CP（成本复杂度参数）表更易于阅读，我们使用plotcp生成 CP 表的信息图形。根据截图（步骤 5），底部 x 轴表示cp值，y 轴表示相对误差，上 x 轴显示树的大小。虚线表示标准差的上线。从截图可以确定，当树的大小为 12 时，最小交叉验证误差发生。

我们还可以使用summary函数显示函数调用、拟合树模型的复杂度参数表、变量重要性，这有助于识别对树分类最重要的变量（总和为 100），以及每个节点的详细信息。

使用决策树的优点在于它非常灵活且易于解释。它适用于分类和回归问题，以及更多；它是非参数的。因此，无需担心数据是否线性可分。至于使用决策树的缺点，它倾向于有偏差且过拟合。然而，你可以通过使用条件推断树来克服偏差问题，并通过随机森林方法或树剪枝来解决过拟合问题。

参见

  • 关于rpart、printcp和summary函数的更多信息，请使用help函数：

    > ?rpart
> ?printcp
> ?summary.rpart

  • C50是另一个提供决策树和基于规则的模型的包。如果你对这个包感兴趣，可以参考cran.r-project.org/web/packages/C50/C50.pdf文档。

可视化递归分割树

从最后一个配方中，我们学习了如何以文本格式打印分类树。为了使树更易于阅读，我们可以使用plot函数来获取已建分类树的图形显示。

准备工作

需要完成上一个配方，通过生成分类模型，并将模型分配给churn.rp变量。

如何做...

执行以下步骤以可视化分类树：

  1. 使用plot函数和text函数绘制分类树：

    > plot(churn.rp, margin= 0.1)
> text(churn.rp, all=TRUE, use.n = TRUE)

    如何做...

    图 3：分类树的图形显示

  2. 你还可以指定uniform、branch和margin参数来调整布局：

    > plot(churn.rp, uniform=TRUE, branch=0.6, margin=0.1)
> text(churn.rp, all=TRUE, use.n = TRUE)

    如何做...

    图 4：调整分类树的布局

它是如何工作的...

在这里，我们演示如何使用plot函数图形化显示分类树。plot函数可以简单地可视化分类树，然后你可以使用text函数向图中添加文本。

在图 3中，我们将margin参数设置为 0.1，以在边界周围添加额外的空白，防止显示的文本被边界截断。它显示分支的长度表示偏差下降的相对幅度。然后我们使用文本函数为节点和分支添加标签。默认情况下，文本函数将在每个分割处添加一个分割条件，并在每个终端节点添加一个类别标签。为了在树图中添加更多信息，我们将参数设置为所有等于TRUE以添加所有节点的标签。此外，我们通过指定use.n = TRUE添加一个参数，以添加额外信息，这表明实际观测数分为两个不同的类别（是和否）。

在图 4中，我们将分支选项设置为 0.6，为每个绘制的分支添加一个肩部。此外，为了显示等长的分支而不是偏差下降的相对幅度，我们将选项uniform设置为TRUE。因此，图 4显示了一个具有短肩部和等长分支的分类树。

参见

  • 您可以使用?plot.rpart来了解更多关于分类树绘制的相关信息。此文档还包括如何指定参数uniform、branch、compress、nspace、margin和minbranch以调整分类树布局的信息。

测量递归分割树的预测性能

由于我们在之前的配方中构建了一个分类树，我们可以用它来预测新观测的类别（类别标签）。在做出预测之前，我们首先验证分类树的预测能力，这可以通过在测试数据集上生成分类表来完成。在本配方中，我们将介绍如何使用predict函数和table函数生成预测标签与真实标签表，并解释如何生成混淆矩阵以衡量性能。

准备工作

您需要完成之前的配方，生成分类模型churn.rp。此外，您还需要准备训练数据集trainset和测试数据集testset，这些数据集是在本章第一节的第一个配方中生成的。

如何操作...

执行以下步骤以验证分类树的预测性能：

  1. 您可以使用predict函数为测试数据集生成预测标签：

    > predictions = predict(churn.rp, testset, type="class")

  2. 使用table函数为测试数据集生成分类表：

    > table(testset$churn, predictions)
 predictions
 yes  no
 yes 100  41
 no   18 859

  3. 可以进一步使用caret包中提供的confusionMatrix函数生成混淆矩阵：

    > library(caret)
> confusionMatrix(table(predictions, testset$churn))
Confusion Matrix and Statistics

predictions yes  no
 yes 100  18
 no   41 859

 Accuracy : 0.942 
 95% CI : (0.9259, 0.9556)
 No Information Rate : 0.8615 
 P-Value [Acc > NIR] : < 2.2e-16 

 Kappa : 0.7393 
 Mcnemar's Test P-Value : 0.004181 

 Sensitivity : 0.70922 
 Specificity : 0.97948 
 Pos Pred Value : 0.84746 
 Neg Pred Value : 0.95444 
 Prevalence : 0.13851 
 Detection Rate : 0.09823 
 Detection Prevalence : 0.11591 
 Balanced Accuracy : 0.84435 

 'Positive' Class : yes

工作原理...

在这个菜谱中，我们使用 predict 函数和构建的分类模型 churn.rp 来预测测试数据集 testset 的可能类别标签。预测的类别（类别标签）编码为是或否。然后，我们使用 table 函数在测试数据集上生成一个分类表。从表中，我们发现 859 个被正确预测为否，而 18 个被错误分类为是。100 个是的预测被正确预测，但有 41 个观测值被错误分类为否。进一步，我们使用来自 caret 包的 confusionMatrix 函数来产生分类模型的测量值。

相关内容

  • 你可以使用 ?confusionMatrix 了解更多关于使用混淆矩阵进行性能测量的信息

  • 对于那些对混淆矩阵定义输出感兴趣的人，请参阅维基百科条目，混淆矩阵 (en.wikipedia.org/wiki/Confusion_matrix)

修剪递归分割树

在之前的菜谱中，我们已经为 churn 数据集构建了一个复杂的决策树。然而，有时我们必须移除在分类实例中不起作用的部分，以避免过拟合，并提高预测精度。因此，在这个菜谱中，我们介绍了成本复杂度修剪方法来修剪分类树。

准备工作

你需要完成之前的菜谱，生成一个分类模型，并将模型分配给 churn.rp 变量。

如何做...

执行以下步骤来修剪分类树：

  1. 找到分类树模型的最低交叉验证误差：

    > min(churn.rp$cptable[,"xerror"])
[1] 0.4707602

  2. 定位具有最低交叉验证误差的记录：

    > which.min(churn.rp$cptable[,"xerror"])
7

  3. 获取具有最低交叉验证误差的记录的成本复杂度参数：

    > churn.cp = churn.rp$cptable[7,"CP"]
> churn.cp
[1] 0.01754386

  4. 通过将 cp 参数设置为具有最低交叉验证误差的记录的 CP 值来修剪树：

    > prune.tree = prune(churn.rp, cp= churn.cp)

  5. 使用 plot 和 text 函数可视化分类树：

    > plot(prune.tree, margin= 0.1)
> text(prune.tree, all=TRUE , use.n=TRUE)

    如何做...

    图 5：修剪后的分类树

  6. 接下来，你可以根据修剪后的分类树模型生成一个分类表：

    > predictions = predict(prune.tree, testset, type="class")
> table(testset$churn, predictions)
 predictions
 yes  no
 yes  95  46
 no   14 863

  7. 最后，你可以根据分类表生成一个混淆矩阵：

    > confusionMatrix(table(predictions, testset$churn))
Confusion Matrix and Statistics

predictions yes  no
 yes  95  14
 no   46 863

 Accuracy : 0.9411 
 95% CI : (0.9248, 0.9547)
 No Information Rate : 0.8615 
 P-Value [Acc > NIR] : 2.786e-16 

 Kappa : 0.727 
 Mcnemar's Test P-Value : 6.279e-05 

 Sensitivity : 0.67376 
 Specificity : 0.98404 
 Pos Pred Value : 0.87156 
 Neg Pred Value : 0.94939 
 Prevalence : 0.13851 
 Detection Rate : 0.09332 
 Detection Prevalence : 0.10707 
 Balanced Accuracy : 0.82890 

 'Positive' Class : yes

它是如何工作的...

在本食谱中，我们讨论了剪枝分类树以避免过拟合并产生更健壮的分类模型。我们首先在cptable中找到具有最小交叉验证错误的记录，然后提取该记录的 CP 并将其值分配给churn.cp。接下来，我们使用prune函数以churn.cp作为参数剪枝分类树。然后，通过使用plot函数，我们图形化地显示了剪枝后的分类树。从图 5中可以看出，树的分割小于原始分类树（图 3）。最后，我们生成了一个分类表，并使用混淆矩阵来验证剪枝树的性能。结果显示，准确率（0.9411）略低于原始模型（0.942），同时也表明剪枝树可能不如原始分类树表现好，因为我们已经剪掉了一些分割条件（尽管如此，人们应该检查敏感性和特异性的变化）。然而，剪枝后的树模型更健壮，因为它去除了可能导致过拟合的一些分割条件。

相关内容

  • 对于那些想了解更多关于成本复杂度剪枝的人来说，请参阅维基百科关于剪枝（决策树）的文章：[en.wikipedia.org/wiki/Pruning_(decision_trees](en.wikipedia.org/wiki/Prunin…

使用条件推理树构建分类模型

除了传统的决策树（rpart）之外，条件推理树（ctree）是另一种流行的基于树的分类方法。与传统决策树类似，条件推理树也通过对因变量进行单变量分割来递归地分割数据。然而，使条件推理树与传统决策树不同的地方在于，条件推理树将显著性检验程序适应于选择变量，而不是通过最大化信息度量（rpart使用基尼系数）来选择变量。在本食谱中，我们将介绍如何适应条件推理树以构建分类模型。

准备工作

您需要完成第一个步骤，通过生成训练数据集trainset和测试数据集testset。

如何操作...

执行以下步骤以构建条件推理树：

  1. 首先，我们使用party包中的ctree构建分类模型：

    > library(party)
> ctree.model = ctree(churn ~ . , data = trainset)

  2. 然后，我们检查构建的树模型：

    > ctree.model

它是如何工作的...

在这个菜谱中，我们使用条件推理树构建了一个分类树。ctree的使用与rpart类似。因此，在面临分类问题时，你可以轻松地使用传统的决策树或条件推理树来测试分类能力。接下来，我们通过检查构建的模型来获取分类树的节点细节。在模型中，我们发现ctree提供的信息类似于分割条件、标准（1 – p 值）、统计（测试统计量）和权重（与节点对应的案例权重）。然而，它提供的信息不如rpart通过使用summary函数提供的信息多。

参见

  • 你可以使用help函数查阅二叉树类的定义，并了解更多关于二叉树属性的信息：

     > help("BinaryTree-class")

可视化条件推理树

与rpart类似，party包也为用户提供了一种可视化条件推理树的方法。在下面的菜谱中，我们将介绍如何使用plot函数来可视化条件推理树。

准备工作

你需要完成第一个菜谱，生成条件推理树模型ctree.model。此外，你还需要在 R 会话中加载trainset和testset。

如何做...

执行以下步骤以可视化条件推理树：

  1. 使用plot函数绘制在上一菜谱中构建的ctree.model：

    > plot(ctree.model)

    如何做...

    图 6：客户流失数据的条件推理树

  2. 要获得一个简单的条件推理树，可以通过减少输入特征来简化构建的模型，并重新绘制分类树：

    > daycharge.model = ctree(churn ~ total_day_charge, data = trainset)
> plot(daycharge.model)

    如何做...

    图 7：使用 total_day_charge 变量作为唯一分割条件的条件推理树

工作原理...

要可视化条件推理树的节点细节，我们可以对构建的分类模型应用plot函数。输出图显示，每个中间节点都显示了因变量名称和 p 值。分割条件显示在左右分支上。终端节点显示了分类观察数的数量，n，以及类别标签为 0 或 1 的概率。

以图 7为例，我们首先使用total_day_charge作为唯一特征，churn作为类别标签构建一个分类模型。构建的分类树显示，当total_day_charge超过 48.18 时，节点 9 中较浅的灰色区域大于较深的灰色区域，这表明日收费超过 48.18 的客户更有可能流失（标签=是）。

参见

  • 条件推理树的可视化来自plot.BinaryTree函数。如果你对调整分类树的布局感兴趣，可以使用help函数阅读以下文档：

    > ?plot.BinaryTree

测量条件推理树的预测性能

在构建条件推理树作为分类模型之后，我们可以使用treeresponse和predict函数来预测测试数据集testset的类别，并进一步通过分类表和混淆矩阵验证预测能力。

准备工作

您需要完成上一个菜谱，生成条件推理树模型ctree.model。此外，您还需要在 R 会话中加载trainset和testset。

如何操作...

执行以下步骤以测量条件推理树的预测性能：

  1. 您可以使用predict函数来预测测试数据集testset的类别：

    > ctree.predict = predict(ctree.model ,testset)
> table(ctree.predict, testset$churn)

ctree.predict yes  no
 yes  99  15
 no   42 862

  2. 此外，您可以使用来自 caret 包的confusionMatrix函数来生成预测结果的性能度量：

    > confusionMatrix(table(ctree.predict, testset$churn))
Confusion Matrix and Statistics

ctree.predict yes  no
 yes  99  15
 no   42 862

 Accuracy : 0.944 
 95% CI : (0.9281, 0.9573)
 No Information Rate : 0.8615 
 P-Value [Acc > NIR] : < 2.2e-16 

 Kappa : 0.7449 
 Mcnemar's Test P-Value : 0.0005736 

 Sensitivity : 0.70213 
 Specificity : 0.98290 
 Pos Pred Value : 0.86842 
 Neg Pred Value : 0.95354 
 Prevalence : 0.13851 
 Detection Rate : 0.09725 
 Detection Prevalence : 0.11198 
 Balanced Accuracy : 0.84251 

 'Positive' Class : yes

  3. 您还可以使用treeresponse函数，它将告诉您类概率列表：

    > tr = treeresponse(ctree.model, newdata = testset[1:5,])
> tr
[[1]]
[1] 0.03497409 0.96502591

[[2]]
[1] 0.02586207 0.97413793

[[3]]
[1] 0.02586207 0.97413793

[[4]]
[1] 0.02586207 0.97413793

[[5]]
[1] 0.03497409 0.96502591

工作原理...

在本菜谱中，我们首先演示了如何使用prediction函数预测测试数据集testset的类别（类标签），然后使用table函数生成分类表。接下来，您可以使用内置在 caret 包中的confusionMatrix函数来确定性能度量。

除了predict函数外，treeresponse还可以估计类概率，这通常会将具有更高概率的标签分类。在本例中，我们演示了如何使用测试数据集testset的前五条记录来获取估计的类概率。treeresponse函数返回一个包含五个概率的列表。您可以使用该列表来确定实例的标签。

相关内容

  • 对于predict函数，您可以指定类型为response、prob或node。如果您在调用predict函数时指定类型为prob（例如，predict(… type="prob")），您将得到与treeresponse返回的完全相同的结果。

使用 k 近邻分类器对数据进行分类

K 近邻（knn）是一种非参数的懒惰学习方法。从非参数的角度来看，它不对数据分布做出任何假设。在懒惰学习方面，它不需要一个显式的学习阶段来进行泛化。以下菜谱将介绍如何在流失数据集上应用 k 近邻算法。

准备工作

您需要完成上一个菜谱，生成训练和测试数据集。

如何操作...

执行以下步骤，使用 k 近邻算法对流失数据进行分类：

  1. 首先，必须安装class包并在 R 会话中加载它：

    > install.packages("class")
> library(class)

  2. 将训练数据集和测试数据集中voice_mail_plan和international_plan属性的yes和no替换为 1 和 0：

    > levels(trainset$international_plan) = list("0"="no", "1"="yes")
> levels(trainset$voice_mail_plan) = list("0"="no", "1"="yes")
> levels(testset$international_plan) = list("0"="no", "1"="yes")
> levels(testset$voice_mail_plan) = list("0"="no", "1"="yes")

  3. 在训练数据集和测试数据集上使用 knn 分类方法：

    > churn.knn  = knn(trainset[,! names(trainset) %in% c("churn")], testset[,! names(testset) %in% c("churn")], trainset$churn, k=3)

  4. 然后，你可以使用 summary 函数检索预测标签的数量：

    > summary(churn.knn)
yes  no
 77 941

  5. 接下来，你可以使用 table 函数生成分类矩阵：

    > table(testset$churn, churn.knn)
 churn.knn
 yes  no
 yes  44  97
 no   33 844

  6. 最后，你可以使用 confusionMatrix 函数生成混淆矩阵：

    > confusionMatrix(table(testset$churn, churn.knn))
Confusion Matrix and Statistics

 churn.knn
 yes  no
 yes  44  97
 no   33 844

 Accuracy : 0.8723 
 95% CI : (0.8502, 0.8922)
 No Information Rate : 0.9244 
 P-Value [Acc > NIR] : 1 

 Kappa : 0.339 
 Mcnemar's Test P-Value : 3.286e-08 

 Sensitivity : 0.57143 
 Specificity : 0.89692 
 Pos Pred Value : 0.31206 
 Neg Pred Value : 0.96237 
 Prevalence : 0.07564 
 Detection Rate : 0.04322 
 Detection Prevalence : 0.13851 
 Balanced Accuracy : 0.73417 

 'Positive' Class : yes

如何工作...

knn 通过训练所有样本并根据相似性（距离）度量对新实例进行分类。例如，相似性度量可以表示如下：

  • 欧几里得距离: 如何工作...

  • 曼哈顿距离: 如何工作...

在 knn 中，一个新的实例被分类到一个标签（类别），这个标签在 k 个最近邻中是共同的。如果 k = 1，那么新的实例将被分配到其最近邻所属的类别。算法的唯一输入是 k。如果我们给出小的 k 输入，可能会导致过拟合。另一方面，如果我们给出大的 k 输入，可能会导致欠拟合。为了选择合适的 k 值，可以依靠交叉验证。

knn 的优点是：

  • 学习过程没有成本

  • 它是非参数的，这意味着你不需要对数据分布做出假设。

  • 当你能找到给定实例的相似性度量时，你可以对任何数据进行分类

knn 的主要缺点是：

  • 解释分类结果很困难。

  • 对于大数据集来说，这是一个昂贵的计算。

  • 性能依赖于维数的数量。因此，对于高维问题，你应该首先降低维度以提高处理性能。

knn 的使用与之前菜谱中提到的基于树的算法应用没有显著差异。然而，虽然基于树的算法可能会展示决策树模型，但 knn 生成的输出仅揭示分类类别因素。然而，在构建分类模型之前，应该将属性从字符串类型替换为整数，因为 k 近邻算法需要计算观测之间的距离。然后，我们通过指定 k=3 来构建分类模型，这意味着选择最近的三个邻居。在分类模型构建完成后，我们可以使用预测因素和测试数据集标签作为输入生成分类表。最后，我们可以从分类表中生成混淆矩阵。混淆矩阵的输出显示准确率为 (0.8723)，这表明之前菜谱中提到的基于树的算法在此情况下优于 k 近邻分类方法的准确率。尽管如此，我们无法仅根据准确率来确定哪个模型更好，还应该检查输出的特异性和敏感性。

参见

使用逻辑回归进行数据分类

逻辑回归是一种概率统计分类模型，可以用来根据一个或多个特征预测类别标签。分类是通过使用logit函数来估计结果概率来完成的。可以通过指定家族为二项式并使用glm函数来使用逻辑回归。在本菜谱中，我们将介绍如何使用逻辑回归进行数据分类。

准备工作

您需要完成第一个菜谱，通过生成训练集和测试集。

如何操作...

执行以下步骤以使用逻辑回归对流失数据进行分类：

  1. 在指定家族为二项式的情况下，我们通过使用glm函数对数据集trainset应用，以churn作为类别标签，其余变量作为输入特征：

    > fit = glm(churn ~ ., data = trainset, family=binomial)

  2. 使用summary函数来获取构建的逻辑回归模型的摘要信息：

    > summary(fit)

Call:
glm(formula = churn ~ ., family = binomial, data = trainset)

Deviance Residuals: 
 Min       1Q   Median       3Q      Max 
-3.1519   0.1983   0.3460   0.5186   2.1284 

Coefficients:
 Estimate Std. Error z value Pr(>|z|)
(Intercept)                    8.3462866  0.8364914   9.978  < 2e-16
international_planyes         -2.0534243  0.1726694 -11.892  < 2e-16
voice_mail_planyes             1.3445887  0.6618905   2.031 0.042211
number_vmail_messages         -0.0155101  0.0209220  -0.741 0.458496
total_day_minutes              0.2398946  3.9168466   0.061 0.951163
total_day_calls               -0.0014003  0.0032769  -0.427 0.669141
total_day_charge              -1.4855284 23.0402950  -0.064 0.948592
total_eve_minutes              0.3600678  1.9349825   0.186 0.852379
total_eve_calls               -0.0028484  0.0033061  -0.862 0.388928
total_eve_charge              -4.3204432 22.7644698  -0.190 0.849475
total_night_minutes            0.4431210  1.0478105   0.423 0.672367
total_night_calls              0.0003978  0.0033188   0.120 0.904588
total_night_charge            -9.9162795 23.2836376  -0.426 0.670188
total_intl_minutes             0.4587114  6.3524560   0.072 0.942435
total_intl_calls               0.1065264  0.0304318   3.500 0.000464
total_intl_charge             -2.0803428 23.5262100  -0.088 0.929538
number_customer_service_calls -0.5109077  0.0476289 -10.727  < 2e-16

(Intercept)                   ***
international_planyes         ***
voice_mail_planyes            * 
number_vmail_messages 
total_day_minutes 
total_day_calls 
total_day_charge 
total_eve_minutes 
total_eve_calls 
total_eve_charge 
total_night_minutes 
total_night_calls 
total_night_charge 
total_intl_minutes 
total_intl_calls              ***
total_intl_charge 
number_customer_service_calls ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

 Null deviance: 1938.8  on 2314  degrees of freedom
Residual deviance: 1515.3  on 2298  degrees of freedom
AIC: 1549.3

Number of Fisher Scoring iterations: 6

  3. 然后，我们发现构建的模型包含不显著的变量，这可能导致误分类。因此，我们只使用显著变量来训练分类模型：

    > fit = glm(churn ~ international_plan + voice_mail_plan+total_intl_calls+number_customer_service_calls, data = trainset, family=binomial)
> summary(fit)

Call:
glm(formula = churn ~ international_plan + voice_mail_plan + 
 total_intl_calls + number_customer_service_calls, family = binomial, 
 data = trainset)

Deviance Residuals: 
 Min       1Q   Median       3Q      Max 
-2.7308   0.3103   0.4196   0.5381   1.6716 

Coefficients:
 Estimate Std. Error z value
(Intercept)                    2.32304    0.16770  13.852
international_planyes         -2.00346    0.16096 -12.447
voice_mail_planyes             0.79228    0.16380   4.837
total_intl_calls               0.08414    0.02862   2.939
number_customer_service_calls -0.44227    0.04451  -9.937
 Pr(>|z|) 
(Intercept)                    < 2e-16 ***
international_planyes          < 2e-16 ***
voice_mail_planyes            1.32e-06 ***
total_intl_calls               0.00329 ** 
number_customer_service_calls  < 2e-16 ***
---
Signif. codes: 
0  es:    des:  **rvice_calls  < '.  es:    de

(Dispersion parameter for binomial family taken to be 1)

 Null deviance: 1938.8  on 2314  degrees of freedom
Residual deviance: 1669.4  on 2310  degrees of freedom
AIC: 1679.4

Number of Fisher Scoring iterations: 5

  4. 然后，您可以使用拟合模型fit来预测testset的结果。您也可以通过判断概率是否高于 0.5 来确定类别：

    > pred = predict(fit,testset, type="response")
> Class = pred >.5

  5. 接下来，使用summary函数将显示二元结果计数，并揭示概率是否高于 0.5：

    > summary(Class)
 Mode   FALSE    TRUE    NA's 
logical      29     989       0

  6. 您可以根据测试数据集的标签和预测结果生成计数统计信息：

    > tb = table(testset$churn,Class)
> tb
     Class
 FALSE TRUE
 yes    18  123
 no     11  866

  7. 您可以将上一步的统计信息转换为分类表，然后生成混淆矩阵：

    > churn.mod = ifelse(testset$churn == "yes", 1, 0)
> pred_class = churn.mod
> pred_class[pred<=.5] = 1- pred_class[pred<=.5]
> ctb = table(churn.mod, pred_class)
> ctb
 pred_class
churn.mod   0   1
 0 866  11
 1  18 123
> confusionMatrix(ctb)
Confusion Matrix and Statistics

 pred_class
churn.mod   0   1
 0 866  11
 1  18 123

 Accuracy : 0.9715 
 95% CI : (0.9593, 0.9808)
 No Information Rate : 0.8684 
 P-Value [Acc > NIR] : <2e-16 

 Kappa : 0.8781 
 Mcnemar's Test P-Value : 0.2652 

 Sensitivity : 0.9796 
 Specificity : 0.9179 
 Pos Pred Value : 0.9875 
 Neg Pred Value : 0.8723 
 Prevalence : 0.8684 
 Detection Rate : 0.8507 
 Detection Prevalence : 0.8615 
 Balanced Accuracy : 0.9488 

 'Positive' Class : 0

它是如何工作的...

逻辑回归与线性回归非常相似；主要区别在于线性回归中的因变量是连续的，而逻辑回归中的因变量是二元的（或名义的）。逻辑回归的主要目标是使用 logit 来得出名义变量与测量变量相关的概率。我们可以用以下方程表示 logit：ln(P/(1-P))，其中 P 是某个事件发生的概率。

逻辑回归的优势在于它易于解释，它指导模型逻辑概率，并为结果提供置信区间。与难以更新模型的决策树不同，您可以在逻辑回归中快速更新分类模型以包含新数据。该算法的主要缺点是它受到多重共线性问题的影响，因此解释变量必须是线性独立的。glm 提供了一个广义线性回归模型，允许在选项中指定模型。如果将家族指定为二项逻辑，则可以将家族设置为二项以对分类因变量进行分类。

分类过程从使用训练数据集生成逻辑回归模型开始，指定 Churn 作为类标签，其他变量作为训练特征，并将家族设置为二项。然后我们使用 summary 函数生成模型的摘要信息。从摘要信息中，我们可能会发现一些不显著的变量（p 值 > 0.05），这可能导致误分类。因此，我们应该只考虑显著变量来构建模型。

接下来，我们使用 fit 函数来预测测试数据集 testset 的分类因变量。fit 函数输出一个类标签的概率，结果等于或低于 0.5 表示预测的标签与测试数据集的标签不匹配，而概率高于 0.5 则表示预测的标签与测试数据集的标签匹配。此外，我们可以使用 summary 函数来获取预测标签是否与测试数据集标签匹配的统计信息。最后，为了生成混淆矩阵，我们首先生成一个分类表，然后使用 confusionMatrix 生成性能度量。

参见

  • 有关如何使用 glm 函数的更多信息，请参阅第四章，理解回归分析，其中涵盖了如何解释 glm 函数的输出。

使用 Naïve Bayes 分类器进行数据分类

Naïve Bayes 分类器也是一种基于概率的分类器，它基于应用贝叶斯定理并假设强独立性。在本食谱中，我们将介绍如何使用 Naïve Bayes 分类器对数据进行分类。

准备工作

您需要完成第一个食谱，生成训练和测试数据集。

如何操作...

使用 Naïve Bayes 分类器对流失数据进行分类的以下步骤：

  1. 加载 e1071 库并使用 naiveBayes 函数构建分类器：

    > library(e1071) 
> classifier=naiveBayes(trainset[, !names(trainset) %in% c("churn")], trainset$churn)

  2. 输入 classifier 以检查函数调用、先验概率和条件概率：

    > classifier

Naive Bayes Classifier for Discrete Predictors

Call:
naiveBayes.default(x = trainset[, !names(trainset) %in% c("churn")], 
 y = trainset$churn)

A-priori probabilities:
trainset$churn
 yes        no 
0.1477322 0.8522678 

Conditional probabilities:
 international_plan
trainset$churn         no        yes
 yes 0.70467836 0.29532164
 no  0.93512418 0.06487582

  3. 接下来，您可以生成测试数据集的分类表：

    > bayes.table = table(predict(classifier, testset[, !names(testset) %in% c("churn")]), testset$churn)
> bayes.table

 yes  no
 yes  68  45
 no   73 832

  4. 最后，你可以从分类表中生成混淆矩阵：

    > confusionMatrix(bayes.table)
Confusion Matrix and Statistics

 yes  no
 yes  68  45
 no   73 832

 Accuracy : 0.8841 
 95% CI : (0.8628, 0.9031)
 No Information Rate : 0.8615 
 P-Value [Acc > NIR] : 0.01880 

 Kappa : 0.4701 
 Mcnemar's Test P-Value : 0.01294 

 Sensitivity : 0.4823 
 Specificity : 0.9487 
 Pos Pred Value : 0.6018 
 Neg Pred Value : 0.9193 
 Prevalence : 0.1385 
 Detection Rate : 0.0668 
 Detection Prevalence : 0.1110 
 Balanced Accuracy : 0.7155 

 'Positive' Class : yes

如何工作...

朴素贝叶斯假设特征在条件上是独立的，即预测变量（x）对类别（c）的影响独立于其他预测变量对类别（c）的影响。它计算后验概率 P(c|x)，如下公式所示：

如何工作...

其中 P(x|c) 被称为似然，p(x) 被称为边缘似然，p(c) 被称为先验概率。如果有许多预测变量，我们可以将后验概率如下公式化：

如何工作...

朴素贝叶斯的优势在于它相对简单且易于使用。当训练集相对较小，可能包含一些噪声和缺失数据时，它很适用。此外，你可以轻松地获得预测的概率。朴素贝叶斯的缺点在于它假设所有特征都是独立的且同等重要，这在现实世界中是非常不可能的。

在这个配方中，我们使用来自 e1071 包的朴素贝叶斯分类器来构建分类模型。首先，我们将所有变量（不包括 churn 类标签）指定为第一个输入参数，并将 churn 类标签指定为 naiveBayes 函数调用中的第二个参数。然后，我们将分类模型分配给变量分类器。接下来，我们打印变量分类器以获取信息，例如函数调用、先验概率和条件概率。我们还可以使用 predict 函数获取预测结果，以及使用 table 函数检索测试数据集的分类表。最后，我们使用混淆矩阵来计算分类模型的性能度量。

最后，我们列出本章中提到的所有算法的比较表：

算法优点缺点
递归分割树

  • 非常灵活且易于解释

  • 适用于分类和回归问题

  • 非参数

|

  • 容易产生偏差和过拟合

|

条件推断树

  • 非常灵活且易于解释

  • 适用于分类和回归问题

  • 非参数

  • 比递归分割树更不容易产生偏差

|

  • 容易过拟合

|

K 最近邻分类器

  • 学习过程成本为零

  • 非参数方法

  • 当你能找到任何给定实例的相似度度量时，你可以对任何数据进行分类

|

  • 分类结果难以解释

  • 对于大数据集，计算成本高昂

  • 性能依赖于维度数量

|

逻辑回归

  • 易于解释

  • 提供模型逻辑概率

  • 提供置信区间

  • 你可以快速更新分类模型以包含新数据

|

  • 患有多重共线性

  • 无法处理连续变量的缺失值

  • 对连续变量的极端值敏感

|

朴素贝叶斯

  • 相对简单直观易用

  • 当训练集相对较小时适用

  • 可以处理一些噪声和缺失数据

  • 可以轻松获得预测的概率

|

  • 假设所有特征都是独立且同等重要的，这在现实世界中是非常不可能的

  • 当训练集数量增加时，容易出现偏差

|

参见

第六章. 分类（II）- 神经网络和 SVM

在本章中，我们将涵盖以下食谱：

  • 使用支持向量机对数据进行分类

  • 选择支持向量机的成本

  • 可视化支持向量机拟合

  • 基于由支持向量机训练的模型预测标签

  • 调整支持向量机

  • 使用 neuralnet 训练神经网络

  • 可视化由 neuralnet 训练的神经网络

  • 基于由 neuralnet 训练的模型预测标签

  • 使用 nnet 训练神经网络

  • 基于由 nnet 训练的模型预测标签

简介

大多数研究表明，支持向量机（SVM）和神经网络（NN）是强大的分类工具，可以应用于多个不同领域。与上一章中提到的基于树或基于概率的方法不同，支持向量机和神经网络从输入到输出的转换过程不太清晰，可能难以解释。因此，支持向量机和神经网络都被称为黑盒方法。

神经网络的发展灵感来源于人类大脑活动。因此，这类网络是一种模仿人类心智模式的计算模型。相比之下，支持向量机首先将输入数据映射到由核函数定义的高维特征空间，并找到通过最大间隔分离训练数据的最佳超平面。简而言之，我们可以将支持向量机视为高维空间中的线性算法。

这两种方法在解决分类问题时都有优点和缺点。例如，支持向量机解决方案是全局最优解，而神经网络可能会遭受多个局部最优解。因此，选择哪种方法取决于数据源的特征。在本章中，我们将说明以下内容：

  • 如何训练支持向量机

  • 观察成本选择如何影响支持向量机分类器

  • 可视化支持向量机拟合

  • 基于由 SVM 训练的模型预测测试数据集的标签

  • 调整 SVM

在神经网络部分，我们将涵盖以下内容：

  • 如何训练神经网络

  • 如何可视化神经网络模型

  • 基于由neuralnet训练的模型预测测试数据集的标签

  • 最后，我们将展示如何使用nnet训练神经网络，以及如何使用它来预测测试数据集的标签

使用支持向量机对数据进行分类

最知名且最受欢迎的支持向量机工具是libsvm和SVMLite。对于 R 语言用户，你可以在e1071包中找到libsvm的实现，在klaR包中找到SVMLite。因此，你可以使用这两个包中实现的功能来训练支持向量机。在本例中，我们将重点关注使用来自e1071包的svm函数（libsvm实现版本）来训练基于电信客户流失数据训练数据集的支持向量机。

准备工作

在本例中，我们将继续使用电信客户流失数据集作为输入数据源来训练支持向量机。对于那些尚未准备数据集的用户，请参阅第五章，分类（I）-树、懒惰和概率性，以获取详细信息。

如何操作...

执行以下步骤来训练 SVM：

  1. 加载e1071包：

    > library(e1071)

  2. 使用svm函数并通过trainset作为输入数据集来训练支持向量机，并使用churn作为分类类别：

    > model  = svm(churn~., data = trainset, kernel="radial", cost=1, gamma = 1/ncol(trainset))

  3. 最后，你可以使用summary函数获取关于构建的模型的整体信息：

    > summary(model)

Call:
svm(formula = churn ~ ., data = trainset, kernel = "radial", cost = 1, gamma = 1/ncol(trainset))

Parameters:
 SVM-Type:  C-classification 
 SVM-Kernel:  radial 
 cost:  1 
 gamma:  0.05882353 

Number of Support Vectors:  691

 ( 394 297 )

Number of Classes:  2 

Levels: 
 yes no

它是如何工作的...

支持向量机构建一个超平面（或一组超平面），在多维空间中最大化两个类别之间的边缘宽度。在这些情况下，定义超平面的点是支持向量，如图所示：

它是如何工作的...

图 1：支持向量机

支持向量机从构建一个最大化边缘宽度的超平面开始。然后，它将定义扩展到非线性可分问题。最后，它将数据映射到一个高维空间，在那里数据可以更容易地用线性边界分离。

使用 SVM 的优势在于它通过一个面向工程问题的核函数构建了一个高度准确模型。此外，它利用正则化项来避免过拟合。它也不受局部最优和多重共线性影响。SVM 的主要局限性在于其训练和测试过程中的速度和大小。因此，它不适合或效率不足以构建大型数据集的分类模型。此外，由于 SVM 难以解释，核函数的确定是如何进行的？正则化是我们需要解决的问题之一。

在本菜谱中，我们继续使用电信 churn 数据集作为我们的示例数据源。我们开始使用 e1071 软件包中提供的 libsvm 训练支持向量机。在训练函数 svm 中，可以指定 kernel 函数、成本和 gamma 函数。对于 kernel 参数，默认值是径向的，可以将核指定为线性、多项式、径向基和 sigmoid。至于 gamma 参数，默认值等于（1/数据维度），它控制分离超平面的形状。增加 gamma 参数通常会增加支持向量的数量。

关于成本，默认值设置为 1，这表示正则化项是常数，值越大，边界越小。我们将在下一菜谱中进一步讨论成本如何影响 SVM 分类器。一旦构建了支持向量机，可以使用 summary 函数获取信息，例如调用次数、参数、类别数量和标签类型。

参见

另一个流行的支持向量机工具是 SVMLight。与提供 libsvm 完整实现的 e1071 软件包不同，klaR 软件包仅提供对 SVMLight 的接口。要使用 SVMLight，可以执行以下步骤：

  1. 安装 klaR 软件包：

    > install.packages("klaR")
> library(klaR)

  2. 从 svmlight.joachims.org/ 下载您平台上的 SVMLight 源代码和二进制文件。例如，如果您的虚拟操作系统是 Windows 64 位，您应该从 download.joachims.org/svm_light/current/svm_light_windows64.zip 下载文件。

  3. 然后，您应该解压文件并将可工作的二进制文件放入工作目录；您可以使用 getwd 函数检查您的工作目录：

    > getwd()

  4. 使用 svmlight 函数训练支持向量机：

    > model.light  = svmlight(churn~., data = trainset, kernel="radial", cost=1, gamma = 1/ncol(trainset))

选择支持向量机的成本

支持向量机通过最大边界创建一个最优的超平面来分离训练数据。然而，有时我们希望在分离类别时允许一些误分类。SVM 模型有一个成本函数，它控制训练错误和边界。例如，小的成本创建大的边界（软边界）并允许更多的误分类。另一方面，大的成本创建窄的边界（硬边界）并允许较少的误分类。在本菜谱中，我们将说明大成本和小成本如何影响 SVM 分类器。

准备工作

在本菜谱中，我们将使用 iris 数据集作为我们的示例数据源。

如何操作...

执行以下步骤以生成两个具有不同成本的不同的分类示例：

  1. 使用列名为 Sepal.Length、Sepal.Width 和 Species 的 iris 数据集子集，其中物种为 setosa 和 virginica：

    > iris.subset = subset(iris, select=c("Sepal.Length", "Sepal.Width", "Species"), Species %in% c("setosa","virginica"))

  2. 然后，你可以生成一个散点图，其中Sepal.Length作为 x 轴，Sepal.Width作为 y 轴：

    > plot(x=iris.subset$Sepal.Length,y=iris.subset$Sepal.Width, col=iris.subset$Species, pch=19)

    如何操作...

    图 2：鸢尾花数据集子集的Sepal.Length和Sepal.Width的散点图

  3. 接下来，你可以使用iris.subset训练成本为 1 的 SVM：

    > svm.model = svm(Species ~ ., data=iris.subset, kernel='linear', cost=1, scale=FALSE)

  4. 然后，我们可以用蓝色圆圈圈出支持向量：

    > points(iris.subset[svm.model$index,c(1,2)],col="blue",cex=2)

    如何操作...

    图 3：用蓝色圆圈圈出支持向量

  5. 最后，我们可以在图上添加分离线：

    > w = t(svm.model$coefs) %*% svm.model$SV
> b = -svm.model$rho
> abline(a=-b/w[1,2], b=-w[1,1]/w[1,2], col="red", lty=5)

    如何操作...

    图 4：在散点图上添加分离线

  6. 此外，我们创建另一个成本为10,000的 SVM 分类器：

    > plot(x=iris.subset$Sepal.Length,y=iris.subset$Sepal.Width, col=iris.subset$Species, pch=19)
> svm.model = svm(Species ~ ., data=iris.subset, type='C-classification', kernel='linear', cost=10000, scale=FALSE)
> points(iris.subset[svm.model$index,c(1,2)],col="blue",cex=2)
> w = t(svm.model$coefs) %*% svm.model$SV
> b = -svm.model$rho
> abline(a=-b/w[1,2], b=-w[1,1]/w[1,2], col="red", lty=5)

    如何操作...

    图 5：大成本分类示例

工作原理...

在这个菜谱中，我们展示了不同的成本如何影响 SVM 分类器。首先，我们创建了一个包含物种setosa和virginica的Sepal.Length、Sepal.Width和Species列的鸢尾花子集。然后，为了创建软边界并允许一些误分类，我们使用成本较小的 SVM（其中cost = 1）来训练支持向量机。接下来，我们用蓝色圆圈圈出支持向量并添加分离线。根据图 5，一个绿色点（virginica）由于成本选择较小而被误分类（被分类为setosa）到分离线的另一侧。

此外，我们还想确定大成本如何影响 SVM 分类器。因此，我们选择了一个大成本（其中cost = 10,000）。从图 5 中，我们可以看到创建的边界很窄（硬边界）且没有误分类案例。因此，这两个示例表明，不同成本的选择可能会影响创建的边界，并影响误分类的可能性。

参见

  • 允许误分类的软边界概念是由 Corinna Cortes 和 Vladimir N. Vapnik 在 1995 年以下论文中提出的：Cortes, C.，and Vapnik, V. (1995). 支持向量机. 机器学习，20(3)，273-297。

可视化 SVM 拟合

为了可视化构建的模型，可以使用绘图函数首先生成数据输入和 SVM 拟合的散点图。在这个图中，支持向量和类别通过颜色符号突出显示。此外，还可以绘制类区域的填充轮廓图，以便从图中轻松识别误分类样本。

准备工作

在这个菜谱中，我们将使用两个数据集：iris数据集和电信churn数据集。对于电信churn数据集，需要完成之前的菜谱，通过训练 SVM 来训练支持向量机，并保存 SVM 拟合模型。

如何操作...

执行以下步骤以可视化 SVM 拟合对象：

  1. 使用基于鸢尾花数据集的 SVM 训练支持向量机，并使用plot函数可视化拟合的模型：

    > data(iris)
> model.iris  = svm(Species~., iris)
> plot(model.iris, iris, Petal.Width ~ Petal.Length, slice = list(Sepal.Width = 3, Sepal.Length = 4))

    如何操作...

    图 6：基于 iris 数据集训练的 SVM 分类图

  2. 使用plot函数，以total_day_minutes和total_intl_charge的维度可视化 SVM 拟合对象model：

    > plot(model, trainset, total_day_minutes ~ total_intl_charge)

    如何操作...

    图 7：基于流失数据集训练的 SVM 分类图

工作原理...

在本配方中，我们展示了如何使用plot函数来可视化 SVM 拟合。在第一个图中，我们使用iris数据集训练了一个支持向量机。然后，我们使用plot函数来可视化拟合的 SVM。

在参数列表中，我们将拟合的模型指定为第一个参数，将数据集（这应该是用于构建模型的数据）作为第二个参数。第三个参数表示用于生成分类图的维度。默认情况下，plot函数只能根据两个维度（用于 x 轴和 y 轴）生成散点图。因此，我们选择变量Petal.Length和Petal.Width作为两个维度来生成散点图。

从图 6中，我们发现Petal.Length被分配到 x 轴，Petal.Width被分配到 y 轴，带有X和O符号的数据点散布在图上。在散点图中，X符号表示支持向量，O符号表示数据点。这两个符号可以通过配置svSymbol和dataSymbol选项来更改。支持向量和真实类别根据它们的标签（绿色代表 virginica，红色代表 versicolor，黑色代表 setosa）突出显示并着色。最后一个参数slice在存在多于两个变量时设置。因此，在这个例子中，我们使用额外的变量Sepal.width和Sepal.length，通过分配常数3和4。

接下来，我们采用相同的方法来绘制基于客户流失数据的 SVM 拟合。在这个例子中，我们使用total_day_minutes和total_intl_charge作为绘制散点图的两个维度。根据图 7，支持向量和红色、黑色数据点在图的中央区域紧密散布，没有简单的方法可以将它们分开。

参见

  • 还有其他参数，如fill、grid、symbolPalette等，可以配置以改变图的布局。您可以使用help函数查看以下文档以获取更多信息：

    > ?svm.plot

基于支持向量机训练的模型预测标签

在前面的配方中，我们基于训练数据集训练了一个 SVM。训练过程通过最大间隔找到将训练数据分开的最优超平面。然后，我们可以利用 SVM 拟合来预测新观察值的标签（类别）。在本配方中，我们将演示如何使用predict函数根据 SVM 训练的模型预测值。

准备工作

你需要完成之前的菜谱，通过生成拟合 SVM，并将拟合模型保存在 model 中。

如何操作...

执行以下步骤以预测测试数据集的标签：

  1. 根据拟合的 SVM 和测试数据集的属性预测测试数据集的标签：

    > svm.pred = predict(model, testset[, !names(testset) %in% c("churn")])

  2. 然后，你可以使用table函数来生成一个包含测试数据集预测结果和标签的分类表：

    > svm.table=table(svm.pred, testset$churn)
> svm.table

svm.pred yes  no
 yes  70  12
 no   71 865

  3. 接下来，你可以使用classAgreement来计算与分类一致性相比的系数：

    > classAgreement(svm.table)
$diag
[1] 0.9184676

$kappa
[1] 0.5855903

$rand
[1] 0.850083

$crand
[1] 0.5260472

  4. 现在，你可以使用confusionMatrix来根据分类表衡量预测性能：

    > library(caret)
> confusionMatrix(svm.table)
Confusion Matrix and Statistics

svm.pred yes  no
 yes  70  12
 no   71 865

 Accuracy : 0.9185 
 95% CI : (0.8999, 0.9345)
 No Information Rate : 0.8615 
 P-Value [Acc > NIR] : 1.251e-08 

 Kappa : 0.5856 
 Mcnemar's Test P-Value : 1.936e-10 

 Sensitivity : 0.49645 
 Specificity : 0.98632 
 Pos Pred Value : 0.85366 
 Neg Pred Value : 0.92415 
 Prevalence : 0.13851 
 Detection Rate : 0.06876 
 Detection Prevalence : 0.08055 
 Balanced Accuracy : 0.74139 

 'Positive' Class : yes

它是如何工作的...

在这个菜谱中，我们首先使用predict函数获取测试数据集的预测标签。然后，我们使用table函数根据测试数据集的预测标签生成分类表。到目前为止，评估过程与上一章中提到的评估过程非常相似。

我们随后引入了一个新的函数，classAgreement，它计算了二维列联表中行和列之间的一致性系数。这些系数包括 diag、kappa、rand 和 crand。diag系数表示分类表主对角线上数据点的百分比，kappa指的是diag，它通过一个变化（随机一致性的概率）进行了校正，rand代表 Rand 指数，它衡量两个数据簇之间的相似性，而crand表示调整了元素随机分组机会的 Rand 指数。

最后，我们使用了caret包中的confusionMatrix来衡量分类模型的表现。准确率 0.9185 表明训练好的支持向量机可以正确分类大多数观测值。然而，准确率本身并不能很好地衡量一个分类模型。还应参考敏感性和特异性。

更多内容...

除了使用 SVM 来预测新观测值的类别外，你还可以使用 SVM 来预测连续值。换句话说，可以使用 SVM 进行回归分析。

在以下示例中，我们将展示如何基于指定为eps-regression类型的拟合 SVM 执行简单的回归预测：

执行以下步骤以使用 SVM 训练回归模型：

  1. 基于一个Quartet数据集训练支持向量机：

    > library(car)
> data(Quartet)
> model.regression = svm(Quartet$y1~Quartet$x,type="eps-regression")

  2. 使用predict函数获取预测结果：

    > predict.y = predict(model.regression, Quartet$x) 
> predict.y
 1        2        3        4        5        6        7        8 
8.196894 7.152946 8.807471 7.713099 8.533578 8.774046 6.186349 5.763689 
 9       10       11 
8.726925 6.621373 5.882946

  3. 在同一张图上绘制预测点为正方形，训练数据点为圆圈：

    > plot(Quartet$x, Quartet$y1, pch=19)
> points(Quartet$x, predict.y, pch=15, col="red")

    更多内容...

    图 8：散点图包含预测数据点和训练数据点

调整支持向量机

除了使用不同的特征集和支持向量机中的kernel函数外，你可以调整配置在参数中的 gamma 和成本来调整其性能。测试不同 gamma 和成本组合值性能的一个可能方法是编写一个for循环来生成所有 gamma 和成本的组合，作为训练不同支持向量机的输入。幸运的是，SVM 提供了一个调整函数tune.svm，这使得调整变得容易得多。在这个配方中，我们将演示如何通过使用tune.svm来调整支持向量机。

准备工作

你需要完成之前的配方，准备一个训练数据集trainset。

如何操作...

执行以下步骤来调整支持向量机：

  1. 首先，使用tune.svm调整支持向量机：

    > tuned = tune.svm(churn~., data = trainset, gamma = 10^(-6:-1), cost = 10^(1:2))

  2. 接下来，你可以使用summary函数获取调整结果：

    > summary(tuned)

Parameter tuning of 'svm':

- sampling method: 10-fold cross validation 

- best parameters:
 gamma cost
 0.01  100

- best performance: 0.08077885 

- Detailed performance results:
 gamma cost      error dispersion
1  1e-06   10 0.14774780 0.02399512
2  1e-05   10 0.14774780 0.02399512
3  1e-04   10 0.14774780 0.02399512
4  1e-03   10 0.14774780 0.02399512
5  1e-02   10 0.09245223 0.02046032
6  1e-01   10 0.09202306 0.01938475
7  1e-06  100 0.14774780 0.02399512
8  1e-05  100 0.14774780 0.02399512
9  1e-04  100 0.14774780 0.02399512
10 1e-03  100 0.11794484 0.02368343
11 1e-02  100 0.08077885 0.01858195
12 1e-01  100 0.12356135 0.01661508

  3. 在调整结果中检索最佳性能参数后，你可以使用最佳性能参数重新训练支持向量机：

    > model.tuned = svm(churn~., data = trainset, gamma = tuned$best.parameters$gamma, cost = tuned$best.parameters$cost)
> summary(model.tuned)

Call:
svm(formula = churn ~ ., data = trainset, gamma = 10^-2, cost = 100)

Parameters:
 SVM-Type:  C-classification 
 SVM-Kernel:  radial 
 cost:  100 
 gamma:  0.01 

Number of Support Vectors:  547

 ( 304 243 )

Number of Classes:  2 

Levels: 
 yes no

  4. 然后，你可以使用predict函数根据拟合的支持向量机预测标签：

    > svm.tuned.pred = predict(model.tuned, testset[, !names(testset) %in% c("churn")])

  5. 接下来，根据测试数据集的预测标签和原始标签生成一个分类表：

    > svm.tuned.table=table(svm.tuned.pred, testset$churn)
> svm.tuned.table

svm.tuned.pred yes  no
 yes  95  24
 no   46 853

  6. 此外，生成一个类一致性来衡量性能：

    > classAgreement(svm.tuned.table)
$diag
[1] 0.9312377

$kappa
[1] 0.691678

$rand
[1] 0.871806

$crand
[1] 0.6303615

  7. 最后，你可以使用混淆矩阵来衡量重新训练的模型性能：

    > confusionMatrix(svm.tuned.table)
Confusion Matrix and Statistics

svm.tuned.pred yes  no
 yes  95  24
 no   46 853

 Accuracy : 0.9312 
 95% CI : (0.9139, 0.946)
 No Information Rate : 0.8615 
 P-Value [Acc > NIR] : 1.56e-12 

 Kappa : 0.6917 
 Mcnemar's Test P-Value : 0.01207 

 Sensitivity : 0.67376 
 Specificity : 0.97263 
 Pos Pred Value : 0.79832 
 Neg Pred Value : 0.94883 
 Prevalence : 0.13851 
 Detection Rate : 0.09332 
 Detection Prevalence : 0.11690 
 Balanced Accuracy : 0.82320 

 'Positive' Class : yes

它是如何工作的...

为了调整支持向量机，你可以使用试错法来找到最佳的 gamma 和成本参数。换句话说，必须生成各种 gamma 和成本的组合，以训练不同的支持向量机。

在这个例子中，我们生成从10^-6到10^-1的不同 gamma 值，以及值为 10 或 100 的成本。因此，你可以使用调整函数svm.tune生成 12 组参数。该函数然后进行 10 次交叉验证，并输出每个组合的错误分散度。结果，错误分散度最低的组合被认为是最佳参数集。从摘要表中，我们发现 gamma 值为 0.01 且成本值为 100 是 SVM 拟合的最佳参数。

在获得最佳参数后，我们可以使用 gamma 等于 0.01 且成本等于 100 的新支持向量机进行训练。此外，我们还可以根据预测标签和测试数据集的标签获得一个分类表。我们还可以从分类表中获得一个混淆矩阵。从混淆矩阵的输出中，你可以确定新训练的模型与原始模型相比的准确性。

参考信息

  • 关于如何使用svm.tune调整 SVM 的更多信息，你可以使用help函数访问此文档：

    > ?svm.tune

使用 neuralnet 训练神经网络

神经网络是由相互连接的节点组构建的，涉及输入、连接权重、处理元素和输出。神经网络可以应用于许多领域，如分类、聚类和预测。要在 R 中训练神经网络，您可以使用 neuralnet，它是在回归分析背景下构建的，用于训练多层感知器，并包含许多灵活的函数来训练前向神经网络。在这个菜谱中，我们将介绍如何使用 neuralnet 来训练神经网络。

准备工作

在这个菜谱中，我们将使用iris数据集作为我们的示例数据集。我们首先将iris数据集分为训练集和测试集。

如何操作...

执行以下步骤以使用 neuralnet 训练神经网络：

  1. 首先加载iris数据集并将数据分为训练集和测试集：

    > data(iris)
> ind = sample(2, nrow(iris), replace = TRUE, prob=c(0.7, 0.3))
> trainset = iris[ind == 1,]
> testset = iris[ind == 2,]

  2. 然后，安装并加载neuralnet包：

    > install.packages("neuralnet")
> library(neuralnet)

  3. 根据名称匹配值在Species列中添加 versicolor、setosa 和 virginica 列：

    > trainset$setosa = trainset$Species == "setosa"
> trainset$virginica = trainset$Species == "virginica"
> trainset$versicolor = trainset$Species == "versicolor"

  4. 接下来，使用具有每层三个隐藏神经元的neuralnet函数训练神经网络。请注意，每次训练的结果可能会有所不同，因此您可能不会得到相同的结果。然而，您可以在开始时使用 set.seed，这样您可以在每次训练过程中得到相同的结果。

    > network = neuralnet(versicolor + virginica + setosa~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, trainset, hidden=3)
> network
Call: neuralnet(formula = versicolor + virginica + setosa ~ Sepal.Length +     Sepal.Width + Petal.Length + Petal.Width, data = trainset,     hidden = 3)

1 repetition was calculated.

 Error Reached Threshold Steps
1 0.8156100175    0.009994274769 11063

  5. 现在，您可以通过访问构建的神经网络模型的result.matrix属性来查看summary信息：

    > network$result.matrix
 1
error                        0.815610017474
reached.threshold            0.009994274769
steps                    11063.000000000000
Intercept.to.1layhid1        1.686593311644
Sepal.Length.to.1layhid1     0.947415215237
Sepal.Width.to.1layhid1     -7.220058260187
Petal.Length.to.1layhid1     1.790333443486
Petal.Width.to.1layhid1      9.943109233330
Intercept.to.1layhid2        1.411026063895
Sepal.Length.to.1layhid2     0.240309549505
Sepal.Width.to.1layhid2      0.480654059973
Petal.Length.to.1layhid2     2.221435192437
Petal.Width.to.1layhid2      0.154879347818
Intercept.to.1layhid3       24.399329878242
Sepal.Length.to.1layhid3     3.313958088512
Sepal.Width.to.1layhid3      5.845670010464
Petal.Length.to.1layhid3    -6.337082722485
Petal.Width.to.1layhid3    -17.990352566695
Intercept.to.versicolor     -1.959842102421
1layhid.1.to.versicolor      1.010292389835
1layhid.2.to.versicolor      0.936519720978
1layhid.3.to.versicolor      1.023305801833
Intercept.to.virginica      -0.908909982893
1layhid.1.to.virginica      -0.009904635231
1layhid.2.to.virginica       1.931747950462
1layhid.3.to.virginica      -1.021438938226
Intercept.to.setosa          1.500533827729
1layhid.1.to.setosa         -1.001683936613
1layhid.2.to.setosa         -0.498758815934
1layhid.3.to.setosa         -0.001881935696

  6. 最后，您可以通过在网络上访问它来查看广义权重：

    > head(network$generalized.weights[[1]])

它是如何工作的...

神经网络是由人工神经元（或节点）组成的网络。网络中有三种类型的神经元：输入神经元、隐藏神经元和输出神经元。在网络中，神经元是相互连接的；神经元之间的连接强度称为权重。如果权重大于零，则处于兴奋状态。否则，处于抑制状态。输入神经元接收输入信息；输入值越高，激活度越大。然后，激活值根据图中的权重和传递函数在网络中传递。隐藏神经元（或输出神经元）随后将激活值相加，并使用传递函数修改总和值。激活值随后通过隐藏神经元流动，并在达到输出节点时停止。因此，可以使用输出神经元的输出值来分类数据。

它是如何工作的...

图 9：人工神经网络

神经网络的优点是：首先，它可以检测因变量和自变量之间的非线性关系。其次，可以使用并行架构高效地训练大型数据集。第三，它是一个非参数模型，因此可以消除参数估计中的误差。神经网络的主要缺点是它往往收敛到局部最小值而不是全局最小值。此外，当训练过程过长时，可能会出现过拟合。

在这个菜谱中，我们展示了如何训练一个神经网络。首先，我们将iris数据集分为训练集和测试集，然后安装neuralnet包并将库加载到 R 会话中。接下来，我们根据Species列中匹配的名称分别添加versicolor、setosa和virginica列。然后，我们使用neuralnet函数来训练网络模型。在函数中指定标签（名称等于 versicolor、virginica 和 setosa 的列）和训练属性之外，我们还配置了每层的隐藏神经元（顶点）数量为三个。

然后，我们检查训练过程和保存在网络中的训练网络的基本信息。从输出信息中，可以看出训练过程需要 11,063 步才能使所有误差函数的绝对偏导数低于 0.01（指定在阈值中）。误差指的是计算赤池信息量准则（AIC）的可能性。要查看详细信息，您可以访问构建的神经网络的result.matrix以查看估计的权重。输出显示估计的权重范围从-18 到 24.40；第一隐藏层的截距为 1.69、1.41 和 24.40，连接到第一隐藏神经元的两个权重估计为 0.95（Sepal.Length）、-7.22（Sepal.Width）、1.79（Petal.Length）和 9.94（Petal.Width）。我们最后可以确定训练的神经网络信息包括广义权重，这些权重表示每个协变量的影响。在这个菜谱中，模型生成了 12 个广义权重，这是四个协变量（Sepal.Length、Sepal.Width、Petal.Length、Petal.Width）与三个响应（setosa、virginica、versicolor）的组合。

参见

  • 对于神经网络的更详细介绍，可以参考以下论文：Günther, F. 和 Fritsch, S. (2010). neuralnet: 神经网络的训练. R 期刊, 2(1), 30-38。

可视化由 neuralnet 训练的神经网络

neuralnet包提供了plot函数来可视化构建的神经网络，以及gwplot函数来可视化广义权重。在接下来的菜谱中，我们将介绍如何使用这两个函数。

准备工作

您需要完成之前的菜谱，通过训练神经网络并保存所有基本信息到网络中。

如何操作...

执行以下步骤以可视化神经网络和广义权重：

  1. 您可以使用plot函数可视化训练好的神经网络：

    > plot(network)

    如何操作...

    图 10：训练好的神经网络的图示

  2. 此外，您还可以使用gwplot可视化广义权重：

    > par(mfrow=c(2,2))
> gwplot(network,selected.covariate="Petal.Width")
> gwplot(network,selected.covariate="Sepal.Width")
> gwplot(network,selected.covariate="Petal.Length")
> gwplot(network,selected.covariate="Petal.Width")

    如何操作...

    图 11：广义权重的图示

工作原理...

在本食谱中，我们展示了如何可视化训练好的神经网络以及每个训练属性的广义权重。根据图 10，该图显示了训练好的神经网络的网络拓扑结构。此外，该图还包括估计的权重、截距以及训练过程的基本信息。在图的下部，可以找到总体误差和收敛所需的步数。

图 11展示了与network$generalized.weights相关的广义权重图。图 11中的四个图显示了四个协变量：Petal.Width、Sepal.Width、Petal.Length和Petal.Width，相对于 versicolor 响应。如果图上的所有广义权重都接近于零，这意味着协变量几乎没有影响。然而，如果整体方差大于一，这意味着协变量有非线性影响。

参见

  • 关于gwplot的更多信息，可以使用help函数访问以下文档：

    > ?gwplot

基于 neuralnet 模型训练预测标签

与其他分类方法类似，我们可以根据训练好的神经网络预测新观察值的标签。此外，我们可以通过使用混淆矩阵来验证这些网络的表现。在接下来的食谱中，我们将介绍如何使用神经网络中的compute函数获取测试数据集标签的概率矩阵，并使用表格和混淆矩阵来衡量预测性能。

准备工作

您需要通过生成训练数据集trainset和测试数据集testset来完成前面的食谱。训练好的神经网络需要保存在网络中。

如何操作...

执行以下步骤以衡量训练好的神经网络的预测性能：

  1. 首先，基于训练好的神经网络和测试数据集testset生成一个预测概率矩阵：

    > net.predict = compute(network, testset[-5])$net.result

  2. 然后，通过找到概率最大的列来获取其他可能的标签：

    > net.prediction = c("versicolor", "virginica", "setosa")[apply(net.predict, 1, which.max)]

  3. 根据预测标签和测试数据集的标签生成一个分类表：

    > predict.table = table(testset$Species, net.prediction)
> predict.table
 prediction
 setosa versicolor virginica
 setosa         20          0         0
 versicolor      0         19         1
 virginica       0          2        16

  4. 接下来，从分类表中生成classAgreement：

    > classAgreement(predict.table)
$diag
[1] 0.9444444444

$kappa
[1] 0.9154488518

$rand
[1] 0.9224318658

$crand
[1] 0.8248251737

  5. 最后，使用confusionMatrix来衡量预测性能：

    > confusionMatrix(predict.table)
Confusion Matrix and Statistics

 prediction
 setosa versicolor virginica
 setosa         20          0         0
 versicolor      0         19         1
 virginica       0          2        16

Overall Statistics

 Accuracy : 0.9482759 
 95% CI : (0.8561954, 0.9892035) 
 No Information Rate : 0.362069 
 P-Value [Acc > NIR] : < 0.00000000000000022204

 Kappa : 0.922252 
 Mcnemar's Test P-Value : NA 

Statistics by Class:

 Class: setosa Class: versicolor Class: virginica
Sensitivity              1.0000000         0.9047619        0.9411765
Specificity              1.0000000         0.9729730        0.9512195
Pos Pred Value           1.0000000         0.9500000        0.8888889
Neg Pred Value           1.0000000         0.9473684        0.9750000
Prevalence               0.3448276         0.3620690        0.2931034
Detection Rate           0.3448276         0.3275862        0.2758621
Detection Prevalence     0.3448276         0.3448276        0.3103448
Balanced Accuracy        1.0000000         0.9388674        0.9461980

工作原理...

在本食谱中，我们展示了如何根据由 neuralnet 训练的模型预测标签。最初，我们使用compute函数根据训练好的神经网络和测试数据集创建输出概率矩阵。然后，为了将概率矩阵转换为类别标签，我们使用which.max函数通过选择行内概率最大的列来确定类别标签。接下来，我们使用一个表格根据测试数据集的标签和预测标签生成分类矩阵。由于我们已经创建了分类表，我们可以使用混淆矩阵来衡量构建的神经网络的预测性能。

参见

  • 在本食谱中，我们使用net.result函数，这是神经网络的总体结果，用于预测测试数据集的标签。除了通过访问net.result来检查总体结果之外，compute函数还生成了每一层的神经元输出。你可以检查神经元的输出，以更好地理解compute是如何工作的：

    > compute(network, testset[-5])

使用 nnet 训练神经网络

nnet包是另一个可以处理人工神经网络的包。此包提供了使用传统反向传播训练前馈神经网络的函数。正如你可以在neuralnet包中找到的大多数神经网络函数一样，在本食谱中，我们提供了一个简短的概述，说明如何使用nnet训练神经网络。

准备工作

在本食谱中，我们不使用从上一步生成的trainset和trainset；请重新加载iris数据集。

如何做...

执行以下步骤以使用nnet训练神经网络：

  1. 首先，安装并加载nnet包：

    > install.packages("nnet")
> library(nnet)

  2. 接下来，将数据集分为训练集和测试集：

    > data(iris)
> set.seed(2)
> ind = sample(2, nrow(iris), replace = TRUE, prob=c(0.7, 0.3))
> trainset = iris[ind == 1,]
> testset = iris[ind == 2,]

  3. 然后，使用nnet训练神经网络：

    > iris.nn = nnet(Species ~ ., data = trainset, size = 2, rang = 0.1, decay = 5e-4, maxit = 200)
# weights:  19
initial  value 165.086674 
iter  10 value 70.447976
iter  20 value 69.667465
iter  30 value 69.505739
iter  40 value 21.588943
iter  50 value 8.691760
iter  60 value 8.521214
iter  70 value 8.138961
iter  80 value 7.291365
iter  90 value 7.039209
iter 100 value 6.570987
iter 110 value 6.355346
iter 120 value 6.345511
iter 130 value 6.340208
iter 140 value 6.337271
iter 150 value 6.334285
iter 160 value 6.333792
iter 170 value 6.333578
iter 180 value 6.333498
final  value 6.333471 
converged

  4. 使用summary获取已训练神经网络的详细信息：

    > summary(iris.nn)
a 4-2-3 network with 19 weights
options were - softmax modelling  decay=0.0005
 b->h1 i1->h1 i2->h1 i3->h1 i4->h1 
 -0.38  -0.63  -1.96   3.13   1.53 
 b->h2 i1->h2 i2->h2 i3->h2 i4->h2 
 8.95   0.52   1.42  -1.98  -3.85 
 b->o1 h1->o1 h2->o1 
 3.08 -10.78   4.99 
 b->o2 h1->o2 h2->o2 
 -7.41   6.37   7.18 
 b->o3 h1->o3 h2->o3 
 4.33   4.42 -12.16

它是如何工作的...

在本食谱中，我们展示了使用nnet包训练神经网络模型的步骤。我们首先使用nnet来训练神经网络。通过这个函数，我们可以设置分类公式、数据源、size参数中的隐藏单元数量、rang参数中的初始随机权重、decay参数中的权重衰减参数以及maxit参数中的最大迭代次数。由于我们将maxit设置为 200，训练过程会反复运行，直到拟合准则值加上衰减项收敛。最后，我们使用summary函数获取构建的神经网络的详细信息，这表明模型使用 4-2-3 网络构建，共有 19 个权重。此外，模型在打印消息的底部显示了一个从节点到另一个节点的权重转换列表。

参见

对于对nnet的背景理论和其制作方式感兴趣的人，请参阅以下文章：

  • Ripley, B. D. (1996) Pattern Recognition and Neural Networks. Cambridge

  • Venables, W. N., and Ripley, B. D. (2002). Modern applied statistics with S. Fourth edition. Springer

基于 nnet 训练的模型进行标签预测

在前一个菜谱中，我们已经使用nnet训练了一个神经网络，现在我们可以根据训练好的神经网络预测测试数据集的标签。此外，我们可以使用来自caret包的混淆矩阵来评估模型。

准备工作

您需要完成上一个菜谱，通过从iris数据集中生成训练数据集trainset和测试数据集testset。训练好的神经网络也需要保存为iris.nn。

如何做...

执行以下步骤以基于训练好的神经网络进行标签预测：

  1. 基于模型iris.nn生成测试数据集的预测：

    > iris.predict = predict(iris.nn, testset, type="class")

  2. 根据预测标签和测试数据集的标签生成分类表：

    > nn.table = table(testset$Species, iris.predict)
 iris.predict
 setosa versicolor virginica
 setosa         17          0         0
 versicolor      0         14         0
 virginica       0          1        14

  3. 最后，根据分类表生成混淆矩阵：

    > confusionMatrix(nn.table)
Confusion Matrix and Statistics

 iris.predict
 setosa versicolor virginica
 setosa         17          0         0
 versicolor      0         14         0
 virginica       0          1        14

Overall Statistics

 Accuracy : 0.9782609 
 95% CI : (0.8847282, 0.9994498) 
 No Information Rate : 0.3695652 
 P-Value [Acc > NIR] : < 0.00000000000000022204

 Kappa : 0.9673063 
 Mcnemar's Test P-Value : NA 

Statistics by Class:

 Class: setosa Class: versicolor
Sensitivity              1.0000000         0.9333333
Specificity              1.0000000         1.0000000
Pos Pred Value           1.0000000         1.0000000
Neg Pred Value           1.0000000         0.9687500
Prevalence               0.3695652         0.3260870
Detection Rate           0.3695652         0.3043478
Detection Prevalence     0.3695652         0.3043478
Balanced Accuracy        1.0000000         0.9666667
 Class: virginica
Sensitivity                 1.0000000
Specificity                 0.9687500
Pos Pred Value              0.9333333
Neg Pred Value              1.0000000
Prevalence                  0.3043478
Detection Rate              0.3043478
Detection Prevalence        0.3260870
Balanced Accuracy           0.9843750

它是如何工作的...

与其他分类方法类似，也可以基于由nnet训练的神经网络进行标签预测。首先，我们使用predict函数根据测试数据集testset生成预测标签。在predict函数中，我们指定类的type参数，以便输出将是类标签而不是概率矩阵。接下来，我们使用table函数根据预测标签和测试数据集中的标签生成分类表。最后，由于我们已经创建了分类表，我们可以使用caret包中的混淆矩阵来衡量训练好的神经网络的预测性能。

参见

  • 对于predict函数，如果未指定class的type参数，默认情况下，它将生成一个概率矩阵作为预测结果，这与neuralnet包内compute函数生成的net.result非常相似：

    > head(predict(iris.nn, testset))
avatar
文章
阅读
粉丝