六、ML：简介
机器学习(ML)是一个相当无定形的，至少在作者看来，是计算机辅助统计的工具包。虽然我们的最终目标将是支持向量机、分类和回归树，以及使用一些最新R包的人工神经网络，但在它们的核心，机器学习只是各种味道的模式识别。
这一介绍性章节的两个焦点，理解样本结构和介绍并行处理，试图为机器学习奠定基础。每个工作流和项目都是不同的，因此我们的意图是深入第六章中的数据将会发生什么，一旦我们的脚湿了，就在下一章中后退一步，适当地清理并设置和减少我们的数据，最后在第三幕中进行一些真正的机器学习。
亲爱的读者，如果你熟悉训练/验证/测试数据集、引导和并行/多核处理，请随意跳到下一章。另一方面，如果你是一名受过更传统训练的统计学家，来自一位作者深情地回忆的舒适的查找表世界，那么准备好以一种与以前截然不同的方式思考数据。
在任何情况下，在这一章中，我们假装数据干净、小而整洁地来到我们面前——选择这个词是因为我们使用了tibbles和名为tidyverse的新 uberpackage 的其他元素，它组合了几个对数据管理有用的包[111]。也许在我们机器学习的探索中，它会更有用。现在，看一下设置和库调用，我们将介绍本章中使用的其他包。
library(checkpoint)
checkpoint("2018-09-28", R.version = "3.5.1",
project = book_directory,
checkpointLocation = checkpoint_directory,
scanForPackages = FALSE,
scan.rnw.with.knitr = TRUE, use.knitr = TRUE)
library(knitr)
library(tidyverse)
library(rsample)
library(data.table)
library(boot)
library(parallel)
library(foreach)
library(doParallel)
options(width = 70, digits = 3)
像往常一样，我们使用几个包，并在这里简要讨论新增加的内容。rsample包[52]允许简单的重采样，而boot包为我们处理通用引导程序[26]。parallel包提供了在多个内核上并行运行计算的功能，从而缩短了得出结果的时间[76]。foreach包是众所周知的 for 循环的变体，它有助于并行化[63]。最后(目前)，doParallel包允许foreach实际使用并行化。
6.1 培训和验证数据
我们从iris样本和一个关于数据的词开始。iris 数据集有 150 个完美的观察值。在机器学习中，一个风险是机器对我们的样本数据学习得太好，然后在真实世界的测试中不太准确。这被称为过度训练或过度适应。换句话说，我们用来估计模型和现实生活之间可能的误差的标准可能太乐观，太小了。一种解决方案是将数据分成两部分。一个是训练片，展示给机器看。另一个是测试部分，作为储备，用于以后评估该模型对新数据或“野生”数据的处理效果。没有什么是不需要权衡或成本的，使用这种技术，挑战是一样的，提供的数据越多，机器应该做得越好。常见的比率包括 80/20、75/25 或 70/30，大部分数据在训练样本中，而较小的一组数据则用于测试。事实上，为了阻止甚至是人为的偏见，在进行探索性数据分析之前进行分割可能是推荐的。
这里值得花点时间来讨论一下这个理论框架。从技术上讲，如果我们遵循这个逻辑结论，有三个数据是有意义的。一组训练集是数据的主体，用于探索性数据分析和模型训练。第二部分，一个较小的集合，将用于验证模型，看看它做得有多好。由于我们可能会从几个模型中进行选择，此验证集(不用于训练模型)可用于在线性拟合和二次拟合之间或者在人工神经网络和随机森林之间做出决定。由于这种验证数据本质上仍将用于模型选择，它在技术上仍是拟合过程的一部分，因此不能真正给出新数据的真实性能。最后，第三部分将被保留到最后，只用于评估模型在面对真实世界的新数据时的误差。
这当然需要足够的数据来牺牲。对于我们在机器学习中的例子，我们并没有忽视这种关注。相反，我们是在本地假设下展示任何单一计算或技术背后的机制，即所研究的技术已经通过文献或一些其他过程被确定为最佳前进方式。未能保留最终的第三个测试集可能是机器学习技术在面对新数据时表现不佳的一个原因(当然，这不是唯一的原因)。
在我们的示例中，使用 80/20 分割可能是有意义的。我们使用基数为R的函数set.seed(5)来确定rsample重采样函数使用的随机数发生器。这允许感兴趣的读者精确复制值。使用的函数是initial_split()，它接受要分割的数据的参数以及要设置到训练集中的数据的比例。鉴于我们的具体情况，在iris的 150 行中，我们保留 29 行作为我们的测试集，保留 80%用于模型构建和训练。
set.seed(5)
case_data <- initial_split(data = iris, prop = 0.8)
case_data
## <121/29/150>
这种重采样实际上是随机抽取数据——这很重要，因为如果数据以某种方式排序，它可能会过度影响模型，尤其是因为我们拆分了数据。我们没有选择在这个重采样分割中使用任何strata。如果这对数据很重要，那么这些层将根据给定的比例进行分割。从case_data对象中，我们通过rsample包中的相关函数调用来提取我们的训练和测试数据。此外，我们使用来自tidyverse的glimpse()来查看我们的训练数据集。
data_train <- training(case_data)
data_test <- testing(case_data)
glimpse(data_train)
## Observations: 121
## Variables: 5
## $ Sepal.Length <dbl> 5.1, 4.9, 4.7, 4.6, 5.0, 5.4, 4.6, 5.0, 4.4,...
## $ Sepal.Width <dbl> 3.5, 3.0, 3.2, 3.1, 3.6, 3.9, 3.4, 3.4, 2.9,...
## $ Petal.Length <dbl> 1.4, 1.4, 1.3, 1.5, 1.4, 1.7, 1.4, 1.5, 1.4,...
## $ Petal.Width <dbl> 0.2, 0.2, 0.2, 0.2, 0.2, 0.4, 0.3, 0.2, 0.2,...
## $ Species <fct> setosa, setosa, setosa, setosa, setosa, seto...
注意我们有四个类型为double的列，我们确实看到有一个用于species的因子列。经过进一步检查，这些竟是三种鸢尾花。
unique(data_train$Species)
## [1] setosa versicolor virginica
## Levels: setosa versicolor virginica
为了将重点放在 ML 和支持合理快速 ML 的计算机环境中使用样本的方式上，我们将一个简单的线性模型拟合到我们的数值数据中，以展示考虑样本数据结构的三种不同方式。目的是提供一个可操作的例子，作为理解这三种方法的一种透镜。
本文之前已经拟合过线性模型，因此我们只注意到该模型是线性拟合，其中数据集中可用的四个数字变量中的三个用于预测第四个变量Petal.Length。
length.lm = lm(Petal.Length ~ Sepal.Length +
Sepal.Width + Petal.Width,
data = data_train)
length.lm
##
## Call:
## lm(formula = Petal.Length ~ Sepal.Length + Sepal.Width + Petal.Width,
## data = data_train)
##
## Coefficients:
## (Intercept) Sepal.Length Sepal.Width Petal.Width
## -0.274 0.723 -0.630 1.466
summary(length.lm)
##
## Call:
## lm(formula = Petal.Length ~ Sepal.Length + Sepal.Width + Petal.Width,
## data = data_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.0349 -0.1699 -0.0061 0.1976 0.5751
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2735 0.3091 -0.88 0.38
## Sepal.Length 0.7230 0.0609 11.87 < 2e-16 ***
## Sepal.Width -0.6298 0.0715 -8.81 1.3e-14 ***
## Petal.Width 1.4661 0.0700 20.93 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1
##
## Residual standard error: 0.299 on 117 degrees of freedom
## Multiple R-squared: 0.971, Adjusted R-squared: 0.97
## F-statistic: 1.32e+03 on 3 and 117 DF, p-value: <2e-16
剩余标准误差(RSE)为 0.299，RSE 可视为回归线的 y 高度与定型集中各点的实际 y 高度之间的一种平均值。接近零的值比远离零的值更可取(就告诉我们我们模型的一些功效而言)。形式上，方程[47]，这是一个真正伟大的文本背后的所有这些更深层次的理论，给出如下:
)
(6.1)
翻译成R，我们在这里展示代码:
sqrt(
sum(
(fitted(length.lm)-data_train$Petal.Length)²
)/(nrow(data_train)-4)
)
## [1] 0.299
在这两种情况下，请注意除法部分的数据集大小减去了变量的数量(在我们的例子中是 4)。虽然免费为我们提供这个值有点儿R,但均方差(MSE)可能是一种更常见的测量拟合优度或质量的方法。这两个等式本质上是相同的，因为在本质上，它们都涉及由模型生成的 y 高度减去在给定成对输入值的数据中找到的 y 高度。由于R已经为我们存储了模型的这些残差，代码看起来有点不同，但是一个快速的实验会显示前面代码的第三行和后面的mean()中的参数是相同的。
mse_train<- mean(length.lm$residuals²)
mse_train
## [1] 0.0862
这两种度量都是以原始数据为单位的，并且这两种度量都有一个关键缺陷。它们是在训练数据上测量的，并且模型被专门训练为在这些相同的数据上表现良好。这是样本内精度测量的一个示例，它不太可能是我们的模型在野生或真实世界数据上表现如何的一个很好的估计器。不是从训练数据计算 MSE(或者实际上任何类型的拟合优度度量),同时拥有训练和测试数据的目的正是为了让我们能够根据训练中没有使用的数据来估计模型性能。我们很自然地转向我们的data_test。
如果我们在测试数据上测量相同类型的 y 值差异，我们会看到两种情况下的数字都更大。这给了我们一些证据，我们有一个模型犯了过度拟合。这也是一个很大的激励因素，说明考虑模型时不仅要考虑样本内数据，还要考虑样本外数据是多么重要。
sqrt(
sum(
(predict(length.lm, data_test)-data_test$Petal.Length)²
)/(nrow(data_test)-2)
)
## [1] 0.41
mse_test <- mean((predict(length.lm, data_test) -
data_test$Petal.Length)²)
mse_test
## [1] 0.156
这种训练与测试过程是关键。任何好的模型的目标都是具有小的残差。事实上，作为模型系数选择的一部分，训练数据的残差通常被最小化，因此，我们不是从训练样本内部的数据来看 MSE，而是从样本外部的数据来看测试数据 MSE。这可以帮助我们理解我们的模型可能如何执行。请记住，正如在开始时提到的，我们目前并不确定这个模型是好是坏。我们真的在使用我们的测试数据来理解我们的模型在部署后如何支持真实世界的使用。如果我们要根据测试数据修改我们的模型，同样，我们实际上需要保留三组数据。
记住，MSE 是原始数据单位的平方。因此，虽然测试数据的 MSE 高于训练数据——实际上几乎是两倍——但它仍然足够好，可以继续进行。有时，任何种类的可操作数据都比没有好，测试 MSE 较高的事实可能只是警告我们的用户在使用模型数据来通知决策时要小心。它还取决于所使用的数据类型。例如，大学入学趋势可能比健康结果更能容忍错误。在这些花数据的情况下，这些数据是以平方厘米为单位的，我们的测试数据表明我们可能会偏离大约半厘米，而训练数据声称可能只偏离四分之一厘米。
我们有一个模型，我们知道它在现实世界中的表现。然而，我们确实损失了 29 个数据点的价值。如果我们能够得到样本外的 MSE 估计值，并且仍然拥有全部 150 个数据点，那就太好了。损失的一部分不仅仅是我们的模型可能在多 20%的数据下做得更好。那就是特定的数据分割，虽然是随机的，但可能仍然有点极端，我们的意思是我们选择在小数据集上分割数据的位置，即使是随机的，也可能影响我们的最终结果。
仅在两个维度上考虑我们的测试与训练数据。在图的顶行，可以看到我们的输入仅仅是Sepal.Length与我们的输出Petal.Length在训练与测试中有一些不同的范围。Q-Q 图的第二行证实了这些样品在图 6-1 中看起来确实略有不同。
par(mfrow=c(2,2))
plot(data_train$Sepal.Length, data_train$Petal.Length)
plot(data_test$Sepal.Length, data_test$Petal.Length)
qqnorm(data_train$Petal.Length,
xlab = "Theoretical Quantiles Train")
qqnorm(data_test$Petal.Length,
xlab = "Theoretical Quantiles Test")
即使是六位数的总结也显示测试数据与我们的样本数据略有不同:
summary(data_train$Sepal.Length)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.3 5.1 5.7 5.8 6.3 7.9
summary(data_test$Sepal.Length)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.40 5.20 6.10 6.02 6.50 7.70
因此，虽然我们可能有比只看训练数据 MSE 更好的想法，我们的真实世界 MSE 可能是什么，但我们受set.seed(5)数据拉动的支配。我们这样做的两个风险是基于样本分割的潜在高方差，在其他条件相同的情况下，如果在 121 和 150 个观察值之间进行选择，我们预计基于 150 个观察值训练的模型会更好一些。
我们在这一部分的目标是激发对单独的训练和测试数据的需求。在这个过于简单的案例中，我们的测试数据可能会让我们对我们训练过的模型在现实世界中的表现有所了解。但是一般来说，我们会因为丢失数据而牺牲很多。我们如何做得更好，尤其是对于较小的样本？
图 6-1
一些数据的 2D 和 Q-Q 图
6.2 重采样和交叉验证
解决这个问题的方法是所谓的交叉验证。我们不是只做一次训练/测试分割，而是做多次。
有几种方法可以完成这个过程。与其将我们的数据分成两部分，不如让我们设想 5 部分，每部分包含我们的 iris 数据集的 30 个观察值。我们将前四个箱合并在一起以形成我们的训练集，同时保留箱 5 作为测试。将模型拟合到训练集，并计算测试集的 MSE。重复这个过程五次，接下来我们只保留 bin 4 作为我们的测试集，合并 1-3 和 5 作为我们的训练集，依此类推。我们依次从每个箱中获得五个不同的 MSE 值。如果我们平均这些，我们现在有所谓的交叉验证 MSE(有时称为 CV)。
虽然可能有一些软件包为我们做了这种计算，但偶尔让我们亲自动手看看到底发生了什么也是有帮助的。首先，我们将我们的虹膜样本随机放入一个变量中，我们称之为crossData，没有替换(这确保我们的箱不以任何方式排序)。我们通过tidyverse包中的sample_n()函数来实现。通过使用%>%管道操作符，我们不必将数据集作为函数的第一个参数。一旦我们的数据被打乱，我们就添加一个列，把它分成五个箱。选择 5 是因为它很小，能很好地分成 150，我们喜欢它。另外，调用一个store变量来保存 MSEs。
crossData <- iris %>%
sample_n(nrow(iris), replace = FALSE)
crossData <- add_column(crossData,
Bin = cut(1:150, breaks = 5, labels = c(1:5)))
store <- tibble(Fold=1:5, MSE=NA_integer_)
现在我们使用一个 for 循环来遍历我们的五个 bin。我们这样做五次，因为我们有五个箱，每次我们都留出一个箱作为测试数据，如上所述。在交叉验证的说法中，折叠的数量是我们的 K 数，这被称为 K 折叠(或在我们的情况下为 5 折叠)。这个 for 循环的每一次传递都代表我们数据的一个“折叠”。我们系统地通过我们的数据来训练我们的线性模型，计算保留的测试箱上的 MSE，将该 MSE 存储在我们的store变量中，并移动到下一个折叠。
for(i in 1:5){
data_train<-crossData %>% filter(Bin != i)
data_test<-crossData %>% filter(Bin == i)
lengthFold.lm = lm(Petal.Length ~ Sepal.Length +
Sepal.Width + Petal.Width,
data = data_train)
store[i,]$MSE <- mean((predict(lengthFold.lm, data_test) -
data_test$Petal.Length)²)
}
现在我们得到了所有五次折叠的平均标准误差。请注意，这个 MSE 来自我们的折叠测试部分，因此这可能会给我们一个很好的想法，即我们可能会期望一个根据所有数据训练的模型在现实生活中对新输入的数据进行处理。
mse_k <- mean(store$MSE)
mse_k
## [1] 0.109
在这种情况下，MSE 约为 0.109，这小于我们在前面章节中只有一组训练/测试数据时的估计值。
当进行 K 折叠和交叉验证时，我们使用的模型不是单个折叠模型，而是一个基于所有数据训练的模型，我们看到样本外估计值与样本内估计值的 MSEs 相当接近。因此，虽然我们的交叉验证 MSE 比我们的完整数据集 MSE 大一点，但它们彼此之间并不是不可能远离的，至少对于花瓣长度这样的东西来说不是太远。
lengthFold.lm <- lm(Petal.Length ~ Sepal.Length +
Sepal.Width + Petal.Width,
data = iris)
lengthFold.lm
##
## Call:
## lm(formula = Petal.Length ~ Sepal.Length + Sepal.Width + Petal.Width,
## data = iris)
##
## Coefficients:
## (Intercept) Sepal.Length Sepal.Width Petal.Width
## -0.263 0.729 -0.646 1.447
mse_ALL <- mean(lengthFold.lm$residuals²)
mse_ALL
## [1] 0.099
请注意，在存储的每个折叠的 MSE 值中(其中模型适合于其他箱，因此这是 for 循环每次迭代的真正样本外测试)，我们有一个从 0.075 到 0.18 的很大范围，这只是显示了单个测试/训练数据集的风险！这又回到了我们之前的警告，当将太多的数据集划分给测试数据集时，较小的数据集会有一点风险。
store
## # A tibble: 5 x 2
## Fold MSE
## <int> <dbl>
## 1 1 0.180
## 2 2 0.109
## 3 3 0.0822
## 4 4 0.0746
## 5 5 0.101
现在，我们没有理由不能只做 4 次折叠，甚至 100 次折叠(进行一些代码编辑),因此一般称为 K-Fold，K 可以是从 1(没有数据分割)、2(我们的第一个例子是一个特殊的加权情况)到 n(这也是一个特殊情况，称为留一个交叉验证或 LOOCV)中选择的任何数字。最后一件要注意的事，好吧，注意到一段时间后就忘了，是关于我们的 for 循环。虽然我们是按顺序运行的，即 1、2、3、4、5，但我们没有理由必须这样运行。换句话说，迭代 2 不依赖于迭代 1。如果我们要拟合一个复杂的模型或者拟合更多的数据，我们的 for 循环每次运行都需要一些时间。尽管如此，只要我们得到 fold-MSEs 表，每次运行都是独立的。再次注意到这一点，现在马上忘记它。
k 倍似乎是一个体面的技术，它是。我们能够在一个完整的数据集上训练一个模型，但我们仍然对我们的模型在样本之外的数据上的表现有所了解。然而，虽然这对于我们正在使用的线性模型来说可能是很好的，但取决于我们真正部署的实际机器学习模型，K-Fold 模型和在完整数据集上使用的模型之间可能没有明显的联系。或者至少，我们可能不想处理这样一个事实，即我们在整个数据集上部署的模型很可能是基于与任何单个 K 倍模型相当不同的系数权重构建的。我们该如何解决这个问题？
还要注意，虽然 K-Fold 模型可以工作，但它确实比我们的简单训练/测试模型需要更多的计算能力。我们从只拟合一个模型到拟合六次模型。就工作量而言，这可能是一个相当大的增长。因此，当我们考虑允许我们的模型在整个数据集上训练并从完全训练的模型中获得我们的误差估计的方法时，要认识到每一种改进都取决于机器技术的提高。
6.3 引导
我们现在把注意力转向另一项技术，即自举。在我们刚刚使用的 K-Fold 或交叉验证技术中，虽然最终我们确实从样本外数据中获得了我们的“真实”MSE 估计，但每个模型都是使用不到一个完整的数据集样本构建的，这意味着每个模型可能都有点弱。因此，我们的误差估计虽然比 pass 好，但可能仍然比它需要的更强。
Bootstrap 用不同的方式解决这个问题。假设我们的样本，最初的iris样本是从鸢尾花群体中精心抽取的样本，我们可以假设，如果我们从iris重新取样，替换，并做很多很多次，我们应该会看到一个看起来非常接近群体的整体。
首先，我们编写一个函数，它是我们自举中感兴趣的统计量。我们之前已经知道 MSE 是什么了。注意我们每次都符合模型。我们的函数将获取一个数据集，计算线性模型，然后计算样本内 MSE 并返回该值。想必我们会在别处收集，事实上我们会！
mse <- function(data, i) {
lengthBoot.lm <- lm(Petal.Length ~ Sepal.Length +
Sepal.Width + Petal.Width,
data=data[i,])
return(mean(lengthBoot.lm$residuals²))
}
现在，使用来自library(boot),的函数boot()，我们运行引导程序。正在发生的是，这个函数通过替换，从我们的虹膜数据样本中进行重采样。事实上，它这样做了一万次。每次，重采样都被发送到我们的mse()函数，该模型的 MSE 返回值存储在bootResults中。
bootResults <- boot(data=iris, statistic=mse, R=10000)
bootResults
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = iris, statistic = mse, R = 10000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.099 -0.00263 0.0129
我们可以用数字和图形来查看结果。图 6-2 显示了所有这些样本内 MSE 的总体分布。
图 6-2
自举结果。
plot(bootResults)
从这个分布中，我们看到 MSEs 的 95%置信区间(注意 0.95 也是默认值)在 0.08 和 0.13 之间。
boot.ci(bootResults, conf = 0.95, type="bca")
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = bootResults, conf = 0.95, type = "bca")
##
## Intervals :
## Level BCa
## 95% ( 0.0794, 0.1343 )
## Calculations and Intervals on Original Scale
就像 K-Fold 方法一样，对于实际模型，我们将使用整个数据集lengthFull.lm的线性拟合来获得可变系数。这只是给我们一种感觉，我们可能期望现实世界的 MSE 生活在哪里。
我们可能认为这个置信区间有点太低，因为它是基于样本内均值标准误差的。人们可以想象对 bootstrap 的改进，其中使用重采样完成后剩余的数据点来计算样本外 bootstrap MSE。当然，这种估计可能有点过高，因为这些数据点是从训练中删除的。虽然每个启动的线性模型都是针对 150 个值进行训练的，但不可能提取所有 150 个原始值-事实上，这就是为什么会得到直方图中显示的正态分布。这种分布正是因为不是所有的值都被拉出来，所以 MSE 是变化的。使用剩余的数据值当然会让我们到达与 K 倍相同的地方，因为我们预计这种 bootstrap 技术的 MSE 可能会比真实世界的数据在完整模型上经历的稍大一些。尽管如此，通过多次进行 bootstrap，我们也许会对我们的模型如何工作有所了解。聪明的研究人员经常开发新技术，目前.632 bootstrap值得探索。然而，我们现在停留在我们的基本引导。这一介绍性章节是关于在执行机器学习时要考虑的三个简单结构。
K-Fold 和 bootstrap 的另一个有趣之处是它需要一些时间来运行。很可能该数据集的 K 倍太小，很难发现。然而，10，000 个引导可能确实变得引人注目。现在，时间不会太长，可能只有 6 秒钟，但我们确实抽取了 10，000 个大小为 150 的样本，总计 1，500，000 个点！
如果我们的原始数据有 1000 个点呢？还是一万？这很可能会降低我们使用 bootstrap 的能力。至少，在我们对如何使用电脑变得更聪明之前是这样。
6.4 并行处理和随机数
至少就计算环境而言，它有两个主要特征。它存在于内存中，而不是硬盘中。因此，任何关于拥有数百万或数十亿数据点与数百或数千数据点的讨论都需要考虑随机存取存储器的大小。尽管如此，由于现代图形贪婪，大多数计算机都有相当多的这种能力。另一个特性是，默认情况下，R 在单个处理器上工作。大多数现代计算机至少有两个。因为我们已经加载了parallel库，所以继续在您的系统上运行代码片段并比较结果:
detectCores()
## [1] 12
detectCores(logical = TRUE)
## [1] 12
detectCores(logical = FALSE)
## [1] 6
我们写这篇文章的计算机有四个物理核心，每个都有两个逻辑“核心”现在，这不是一本关于计算机硬件的书，所以让我们忽略物理和逻辑内核的细微差别。事实上，根据您的操作系统和硬件系统，您可能无法区分这三者。在任何情况下，如果返回的数字大于 1，那么如果将更多的处理器放入工作池，bootstrapping 和 K-Fold 可以运行得更快。
还记得我们在 K-Fold 中说过，我们的 for 循环每次都是完全相同的操作，迭代 c 不依赖于要计算的迭代 a。自举也是一样。我们正在进行 10，000 次随机拉动，并将其输入到我们的函数mse(),中，如果我们正在做一个比简单线性拟合更复杂的模型，或者如果我们的数据变得更大，我们可能会发现自己需要等待相当长的时间。如果同样的过程是独立完成的，那么本节的方法可以很好地工作。如果这个过程不是独立的，那就超出了本文的范围(当然，尽管有方法来适应这个过程)。
现在，多核与集群和并行之间存在细微差别。如果您的数据生活在真正的大数据世界中，这些细微差别可能对您有意义。但是，我们冒险记忆也可能成为一个问题，这些也超出了本文的范围。然而，如果我们对 10，000 次重采样的 6 秒估计是正确的，并且我们得到了两个内核而不是一个，我们应该能够将执行时间减少一半。
这在实践中并不完全正确，因为启动和停止多核环境会产生一些开销。然而，回想一下，我们的分析只是一个简单的线性模型，然后是一个平均值。这些间接成本大多是固定的，因此如果模型拟合更奇特，并且引导花费一天或半天，那么将该过程削减大约一半开始变得有意义。而且，由于我们的系统有几个核心，我们实际上可以做得比一半更好！一个警告:无论你的系统中有多少个内核，至少留一个给系统使用——否则一切都会停滞不前，直到R完成。最有可能的是，你会听到风扇运转的声音——毕竟你要做一些计算工作了！
在下面的代码中，我们将只在 2 个内核上运行，因此我们将 10，000 次启动分为 5，000 次。为了简单起见，我们构建一个函数runP()来处理并行流程运行。接下来，我们使用makeCluster()函数创建R的两个并行副本，在其上运行我们的函数。我们的R副本需要与我们的机器拥有相同的环境，所以我们使用clusterEvalQ()函数来设置我们的cl集群，使其拥有我们需要的库结构。虽然我们使用checkpoint()来确保版本是正确的，但是R的每个副本只需要boot库。在这里吝啬可以节省开销时间，所以在传递信息时要明智。从那里，我们还使用clusterExport()从我们的全球环境中导出一些已经构建的函数。我们做的最后一点集群开销是设置我们的随机数流，以便我们的结果通过clusterSetRNGStream()与你的相匹配。
从这里开始，其余的大部分应该对引导很熟悉，所以我们花一点时间来剖析实际的并行调用——记住前面的步骤是开销。从 base R调用apply()函数对某种数组应用函数。函数parLapply()是相同的，除了第一个参数链接到我们刚刚建立的集群，第二个参数只是说我们将运行我们的函数两次，第三个当然是我们之前构建的 bootstrapping 5,000 runs 函数。do.call()函数也是一个基本函数，它只是在我们的并行函数上调用组合函数c(),以便将两个结果合并到一个引导对象中。这样做是为了让熟悉的boot.ci()像以前一样工作。当然，我们确实花了一些时间通过stopCluster()关闭我们的集群。
## notice 10000/2 = 5000
runP <- function(...) boot(data=iris, statistic=mse, R=5000)
## makes a cluster with 2 cores as 10000/5000 = 2
cl<-makeCluster(2)
## passes along parts of the global environment
## to each node / part of the cluster
## again, base is a file path variable to our book's path
## set book_directory <- "C:/YourPathHere/"
clusterExport(cl, c("runP", "mse", "book_directory", "checkpoint_directory" ))
## creates the library and some environment on
## each of the parts of the cluster
clusterEvalQ(cl, {
library(checkpoint)
checkpoint("2018-09-28", R.version = "3.5.1",
project = book_directory,
checkpointLocation = checkpoint_directory,
scanForPackages = FALSE,
scan.rnw.with.knitr = TRUE, use.knitr = TRUE)
library(boot)
})
## [[1]]
## [1] "boot" "checkpoint" "RevoUtils" "stats"
## [5] "graphics" "grDevices" "utils" "datasets"
## [9] "RevoUtilsMath" "methods" "base"
##
## [[2]]
## [1] "boot" "checkpoint" "RevoUtils" "stats"
## [5] "graphics" "grDevices" "utils" "datasets"
## [9] "RevoUtilsMath" "methods" "base"
## similar to set.seed() except for clusters
clusterSetRNGStream(cl, 5)
## uses the parLapply() function
which works on windows too
pBootResults <- do.call(c, parLapply(cl, seq_len(2), runP))
#stop the cluster
stopCluster(cl)
# view results
pBootResults
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = iris, statistic = mse, R = 5000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.099 -0.00295 0.0128
## get 95% confidence interval of the MSEs
## (note 0.95 is the default)
boot.ci(pBootResults, conf = 0.95, type="bca")
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = pBootResults, conf = 0.95, type = "bca")
##
## Intervals :
## Level BCa
## 95% ( 0.0795, 0.1352 )
## Calculations and Intervals on Original Scale
就是这样。根据我们的计算，这还不到使用一个处理器所用时间的一半，但这没关系。请注意，这其中的大部分可以根据您所处的环境进行修改。由于我们有了更多的内核，我们可以增加激活的R实例的数量，并相应地减少runP()函数中的引导运行次数。或者，通过将seq_len()编辑为四个而不是两个，我们可以在两个内核上运行不止一个周期的runP()。在这种情况下，我们让引导函数运行四次，在每个 R 实例上运行两次，我们得到 20，000 次引导运行。因此，您可以看到，它基于可用的硬件和数据细微差别，非常可定制和灵活。图 6-2 中的结果与图 6-3 中的结果相似:
plot(pBootResults)
不过，值得花点时间承认这两个数字略有不同——尽管幸运的是，您的机器在每种情况下都与我们的相符。原因在于随机数及其生成方式。要使研究具有可重复性，有助于set.seed()或clusterSetRNGStream()。这不仅有助于我们在各种机器上的输出与您的相匹配，还有助于我们今天的输出与明天的输出相匹配。随机数是从一个起点或种子生成的，这就是这两种技术的作用。然而，在集群的情况下，虽然我们希望在全局级别上有相同的种子，但对于每个处理器，我们确实希望有稍微不同的起点，以便我们真正获得 10，000 个唯一的重采样，而不是 5，000 个重复的副本。这是通过使用种子生成器的集群变体来实现的，它在幕后相当巧妙地处理这个过程，以便恰好出现可重复的随机性。尽管如此，这是一个不同于之前使用的随机性，所以数字会略有不同。
图 6-3
引导结果
在本节中，我们看到了如何为任何流程构建一个集群，并合并集群的结果。我们使用了我们的boot()函数(不可否认，现在它已经内置了一些并行功能)。虽然引导过程可以通过 for 循环手工编码，但我们选择使用一个包。这在R中经常是可能的，在CRAN上有许多包和库。这种不断更新和改进的生态系统是checkpoint包的原因之一。适用于一个版本库的代码可能在下个月或明年就不再适用了。对于任何处理不断变化的数据的人来说，这是机器学习的一个功能，其中真实世界的数据进来并需要分析，这是非常有用的，既可以将代码锁定在特定的时间点，以便它在未来继续工作，也可以偶尔重新访问代码并改进写作，以利用更新或更新的库。
为每一个
我们的最后一节是关于foreach()函数和同名的库。现在，虽然我们在前面介绍的 K-Fold for()函数的上下文中介绍了这一点，但这当然适用于任何循环，只要每个迭代都是独立的。然而，一般来说，在R中经常有对数据矢量化循环函数的方法，并且这些方法会比这个过程更快。此外，对于这种特定情况，由于创建集群的开销，这实际上并没有节省时间，因为我们在简单的线性模型上只处理了 150 位数据 5 次。我们的 K-Fold 在计算上比我们的 bootstrap 简单得多！
接下来的大部分代码都很熟悉。我们从建立集群的代码开始。这里唯一的新函数是registerDoParallel()，它将我们的foreach()函数连接到正确的框架。和以前一样，我们选择在两个内核上完成这项工作。因为我们不使用boot，我们不叫那个库。然而，我们确实调用了tidyverse库，因为我们在 tibbles 上使用了管道。
cl <- makeCluster(2)
registerDoParallel(cl)
clusterExport(cl, c("book_directory", "checkpoint_directory"))
clusterEvalQ(cl, {
library(checkpoint)
checkpoint("2018-09-28", R.version = "3.5.1",
project = book_directory,
checkpointLocation = checkpoint_directory,
scanForPackages = FALSE,
scan.rnw.with.knitr = TRUE, use.knitr = TRUE)
library("tidyverse")
})
## [[1]]
## [1] "forcats" "stringr" "dplyr" "purrr"
## [5] "readr" "tidyr" "tibble" "ggplot2"
## [9] "tidyverse" "checkpoint" "RevoUtils" "stats"
## [13] "graphics" "grDevices" "utils" "datasets"
## [17] "RevoUtilsMath" "methods" "base"
##
## [[2]]
## [1] "forcats" "stringr" "dplyr" "purrr"
## [5] "readr" "tidyr" "tibble" "ggplot2"
## [9] "tidyverse" "checkpoint" "RevoUtils" "stats"
## [13] "graphics" "grDevices" "utils" "datasets"
## [17] "RevoUtilsMath" "methods" "base"
感兴趣的主要代码在这里。同样，这几乎正是我们第一次做 K 折叠时发生的事情，这就是重点。唯一的变化是foreach()函数和参数。我们再次通过折叠 1 到 5 运行我们的函数。我们确实使用了一个新的参数，.combine,并把它设置到我们的组合函数中，这个函数我们在R和之前的do.call()函数中很熟悉。然后，我们使用%dopar%通过我们在前面的设置代码中设置的注册后端来调用并行处理。这些包处理剩下的，将这五次运行发送到两个内核并合并结果。
k <- foreach(i=1:5, .combine = c) %dopar% {
data_train <- crossData %>% filter(Bin != i)
data_test <- crossData %>% filter(Bin == i)
lengthFold.lm <- lm(Petal.Length ~ Sepal.Length +
Sepal.Width + Petal.Width,
data = data_train)
mean((predict(lengthFold.lm, data_test) -
data_test$Petal.Length)²)
}
请注意，由于要计算五次折叠，而只有两个处理器，所以有点不匹配。根据我们的循环，可以想象这些迭代中的一个或多个可能或多或少是计算密集型的。在这种情况下，当然不可能。尽管如此，还是有可能以这样一种方式建立并行处理，即将花费很长时间的交互首先发送到一个内核，而其余的发送到另一个内核。默认值是它们先入先出，按顺序运行。
stopCluster(cl)
mse_Pk<-mean(k)
mse_Pk
## [1] 0.109
我们看到，之前 for 循环中的 MSE 0.109 与现在的 MSE 0.109 相匹配。当对相同的数据运行时，这两个过程是相同的。事实上，当第一次(甚至是最后一次)通过foreach()让代码运行得更快以确保没有犯愚蠢的错误时，这可能是明智的。如果我们的读者发现自己处于这样一种情景中，作者只能增加他们冰冷的安慰，即我们过去和现在都比期望的更频繁。
6.5 总结
这就完成了机器学习的设置。有了三种好的方法来训练数据，然后获得对现实世界性能的一些估计，我们现在有了思维导图来了解我们的数据将如何被格式化和使用。有了几个好的多核计算包和代码在手，我们就有了让我们的机器学习不需要太长时间的技术工具。
事实上，并行计算技术——尤其是foreach()——非常容易实现，可以节省大量时间。现在，有很多方法可以绕过 for 循环。矢量化计算就是其中之一。举一个个人的例子，一个作者每天处理相当大的数据(大约 200，000 行和 70 列)。从数据库中提取数据后，基本的数据转换至少需要半天时间。代码是几年前写的，场景最初对时间并不敏感。一般的用例是在周末运行。然而，操作发生了变化，代码开始需要更频繁地运行。在写这段文字的时候，你的作者决定清理代码。在一台小型的四核办公计算机上，通过对一些过程进行矢量化，对其他过程进行并行化，时间从半天缩短到不到 2 小时！
这也许是看待这一章的一个很好的方式。它确实让我们可以用机器学习做一些非常棒的事情。然而，这些技术可以混合搭配，让生活和数据变得更简单、更友好。就像一些工作的旧代码一样，只是现在有点太慢了，不方便，有可能通过一些非常小的改变获得一些非常重要的突破！
像往常一样，我们通过调用表 6-1 中的主要函数来结束。
表 6-1
本章中描述的关键功能列表及其功能摘要
|功能
|
它的作用
|
| — | — |
| boot() | 从引导库并运行引导程序 |
| boot.ci() | 计算引导输出的置信区间 |
| clusterEvalQ() | 从并行库中设置要复制到每个集群实例的环境 |
| clusterExport() | 从全局环境中导出值，供每个集群环境使用 |
| clusterSetRNGStream() | 创建一个随机数种子，并适当地渗透到每个集群 |
| detectCores() | 确定 R 中有多少内核可用 |
| do.call() | 调用函数并通过列表传递参数的基本函数 |
| filter() | tidyverse 函数返回某些行 |
| fitted() | 从线性模型中提取拟合的 y 值 |
| foreach() | for的平行版本(大部分) |
| glimpse() | tidyverse 函数，水平显示列名，后跟前几个值 |
| initial_split() | 将数据分成训练和测试用例的 rsample 函数 |
| lm() | 基于提供的关系创建线性模型 |
| makeCluster() | 形成一个集群 |
| parLapply() | 除适用于并行计算之外的应用功能 |
| predict() | 显示线性模型的预测 y 值 |
| registerDoParallel() | 向并行集群注册foreach()函数的 doParallel 包操作符 |
| seq_len() | 创建一定长度的序列 |
| set.seed() | 为可重现的结果创建随机数种子 |
| stopCluster() | 停止由makeCluster创建的集群 |
| testing() | 提取测试数据的 rsample 函数 |
| training() | 提取列车数据的 rsample 函数 |
| unique() | 仅返回唯一值 |
七、ML：无监督
这一章关注无监督机器学习，它通常处理未标记的数据。目标是基于共同的特征将这些数据分类到相似的组中。通常，尽管不总是，无监督的机器学习也被用作一种降维。例如，如果您得到的数据集包含数百或数千个要素，但只有几千个案例，您可能希望首先利用无监督学习将大量要素提取到数量较少的维度中，这些维度仍然可以从较大的集合中捕获大部分信息。无监督机器学习也是探索性数据分析阶段的最后一步。可以利用无监督机器学习中的部分排序或聚类来了解您的数据有多少“唯一”的组或维度。想象一个由来自几个不同地理区域的各种指标组成的数据集。人们可能期望无监督的分组技术来指示关于地理区域的一些事情。或者，人们可能会发现相距遥远的位置有几个高度共同的特征。
无论有无标签，数据总是杂乱无章地出现在我们面前。因此，预处理阶段通常需要付出大量的努力。通常，数据不符合机器学习算法预期的格式。通常，需要的是一系列列，其中前几列在某种意义上是关键列(可能是地理区域和收集某些测量值的年份)，而后面的列表示对特定特征或变量的观察。可能需要进行操作才能将数据转换成正确的格式，其中每一行都是特定时间点的唯一观察值。我们将在本章的示例数据中看到一些数据操作。
像往常一样，我们使用几个包，并在这里简要讨论新增加的内容。readxl包[112]允许从 Excel 数据文件中快速读取我们的示例数据，而ape包为我们处理树状图选项[72]。MASS包提供了处理非线性降维的功能[99]。matrixStats包提供了矩阵和向量的功能。最后，viridis软件包提供了一个更好的调色板——特别是对于色盲来说。
library(checkpoint)
checkpoint("2018-09-28", R.version = "3.5.1",
project = book_directory,
checkpointLocation = checkpoint_directory,
scanForPackages = FALSE,
scan.rnw.with.knitr = TRUE, use.knitr = TRUE)
library(ggplot2)
library(cowplot)
library(viridis)
library(scales)
library(readxl)
library(data.table)
library(ape)
library(MASS)
library(matrixStats)
options(width = 70, digits = 2)
本章讨论的一种方法是主成分分析(PCA)。虽然R有一些 PCA 的内置函数，但是在包pcaMethods中有更广泛的选项，它是Bioconductor的一部分，是 CRAN 包库的替代，有许多用于生物信息学的R包。我们可以使用下面的代码安装pcaMethods。注意到
source("https://bioconductor.org/biocLite.R")
biocLite("pcaMethods")
一旦安装了pcaMethods，我们就可以像加载任何其他R包一样加载它。
library(pcaMethods)
7.1 数据背景和探索性分析
所使用的样本数据来自世界银行，在 CC-BY 4.0 下获得许可，并已被修改为仅包括来自某些地区的某些列的数据[2]。数据在Gender_StatsData_worldbank.org_ccby40.xlsx里。特别是，我们将数据简化为近年来的一些关键信息，这样就没有遗漏变量。所使用的文件可以从出版商的网站或 GitHub 存储库中下载。
应该注意的是，read_excel()默认为stringsAsFactors = FALSE，这通常是人们希望读入数据的方式。如果需要因子，很有可能在以后通过适当的函数调用来控制它。
## Note: download Excel file from publisher website first
dRaw <- read_excel("Gender_StatsData_worldbank.org_ccby40.xlsx")
dRaw <- as.data.table(dRaw) # convert data to data.table format.
丢失的变量可以用不同的方法处理。处理缺失数据的最简单方法是删除所有缺失的信息。代价是将我们的数据缩减为全部可用信息的一个非常小的子集。更好地处理缺失数据的低成本方法将在我们关于缺失数据的章节中讨论。现在，只需注意在本章的示例数据中没有丢失数据。
理解数据的结构是任何分析的关键的第一步。例如，至少，知道存在什么类型的数据以及数据是如何组织的是很重要的。对我们的数据使用str()函数告诉我们，该结构具有按年份组织的字符和数字数据，从 1997 年到 2014 年。事实上，大多数数据都是数字。
summary()函数按列生成一个简短的摘要，并显示我们的数据变化很大。在这一点上，我们开始意识到我们的数据对于机器学习格式来说不是最佳的。换句话说，每一列都不是单一的、唯一的度量。此外，时间信息不是在 year 变量中捕获的，而是在变量名本身中捕获的。
深入研究我们数据的标签，我们发现通过unique()函数，我们的数据来自几个大的地理区域，似乎围绕着某些特定的指标。
str(dRaw)
## Classes 'data.table' and 'data.frame': 99 obs. of 21 variables:
## $ CountryName : chr "Sub-Saharan Africa" "Sub-Saharan Africa" "Sub-Saharan Africa" "Sub-Saharan Africa" ...
## $ Indicator Name: chr "Adolescent fertility rate (births per 1,000 women ages 15-19)" "Age dependency ratio (% of working-age population)" "Children out of school, primary, female" "Children out of school, primary, male" ...
## $ IndicatorCode : chr "SP.ADO.TFRT" "SP.POP.DPND" "SE.PRM.UNER.FE" "SE.PRM.UNER.MA" ...
## $ 1997 : num 1.32e+02 9.17e+01 2.44e+07 2.03e+07 1.55e+01 ...
## $ 1998 : num 1.31e+02 9.13e+01 2.44e+07 2.05e+07 1.53e+01 ...
## $ 1999 : num 1.30e+02 9.09e+01 2.43e+07 2.08e+07 1.52e+01 ...
## $ 2000 : num 1.28e+02 9.04e+01 2.37e+07 2.00e+07 1.49e+01 ...
## $ 2001 : num 1.27e+02 9.04e+01 2.31e+07 1.95e+07 1.47e+01 ...
## $ 2002 : num 1.26e+02 9.02e+01 2.28e+07 1.91e+07 1.44e+01 ...
## $ 2003 : num 124 90 21938840 18230741 14 ...
## $ 2004 : num 1.23e+02 8.97e+01 2.14e+07 1.79e+07 1.36e+01 ...
## $ 2005 : num 1.21e+02 8.94e+01 2.06e+07 1.71e+07 1.32e+01 ...
## $ 2006 : num 1.20e+02 8.94e+01 1.99e+07 1.67e+07 1.28e+01 ...
## $ 2007 : num 1.18e+02 8.94e+01 1.94e+07 1.52e+07 1.23e+01 ...
## $ 2008 : num 1.16e+02 8.92e+01 1.90e+07 1.52e+07 1.19e+01 ...
## $ 2009 : num 1.15e+02 8.89e+01 1.92e+07 1.56e+07 1.15e+01 ...
## $ 2010 : num 1.13e+02 8.85e+01 1.98e+07 1.61e+07 1.10e+01 ...
## $ 2011 : num 1.11e+02 8.83e+01 1.92e+07 1.54e+07 1.07e+01 ...
## $ 2012 : num 1.09e+02 8.80e+01 1.91e+07 1.55e+07 1.03e+01 ...
## $ 2013 : num 1.07e+02 8.75e+01 1.91e+07 1.53e+07 1.00e+01 ...
## $ 2014 : num 1.05e+02 8.69e+01 1.92e+07 1.56e+07 9.73 ...
## - attr(*, ".internal.selfref")=<externalptr>
summary(dRaw)
## CountryName Indicator Name IndicatorCode
## Length:99 Length:99 Length:99
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
## 1997 1998 1999
## Min. : 6 Min. : 6 Min. : 6
## 1st Qu.: 16 1st Qu.: 15 1st Qu.: 15
## Median : 71 Median : 71 Median : 71
## Mean : 949741 Mean : 872971 Mean : 817616
## 3rd Qu.: 4366 3rd Qu.: 4338 3rd Qu.: 3993
## Max. :24371987 Max. :24437801 Max. :24292225
## 2000 2001 2002
## Min. : 6 Min. : 6 Min. : 5
## 1st Qu.: 15 1st Qu.: 16 1st Qu.: 16
## Median : 70 Median : 70 Median : 70
## Mean : 781078 Mean : 736806 Mean : 674889
## 3rd Qu.: 4224 3rd Qu.: 4275 3rd Qu.: 4672
## Max. :23672959 Max. :23125633 Max. :22795557
## 2003 2004 2005
## Min. : 5 Min. : 5 Min. : 5
## 1st Qu.: 16 1st Qu.: 15 1st Qu.: 16
## Median : 70 Median : 71 Median : 71
## Mean : 651075 Mean : 637985 Mean : 659420
## 3rd Qu.: 5568 3rd Qu.: 6772 3rd Qu.: 8042
## Max. :21938840 Max. :21350198 Max. :20582825
## 2006 2007 2008
## Min. : 5 Min. : 5 Min. : 5
## 1st Qu.: 15 1st Qu.: 16 1st Qu.: 16
## Median : 71 Median : 71 Median : 72
## Mean : 653180 Mean : 597847 Mean : 573176
## 3rd Qu.: 9166 3rd Qu.: 11168 3rd Qu.: 13452
## Max. :19904220 Max. :19402096 Max. :19015196
## 2009 2010 2011
## Min. : 5 Min. : 5 Min. : 5
## 1st Qu.: 16 1st Qu.: 16 1st Qu.: 16
## Median : 72 Median : 72 Median : 72
## Mean : 569320 Mean : 569669 Mean : 561551
## 3rd Qu.: 12484 3rd Qu.: 12654 3rd Qu.: 13404
## Max. :19209252 Max. :19774011 Max. :19191406
## 2012 2013 2014
## Min. : 5 Min. : 5 Min. : 5
## 1st Qu.: 16 1st Qu.: 16 1st Qu.: 16
## Median : 72 Median : 73 Median : 73
## Mean : 567238 Mean : 592806 Mean : 610288
## 3rd Qu.: 13047 3rd Qu.: 13574 3rd Qu.: 13852
## Max. :19068296 Max. :19092876 Max. :19207489
unique(dRaw$CountryName)
## [1] "Sub-Saharan Africa" "North America"
## [3] "Middle East & North Africa" "Latin America & Caribbean"
## [5] "European Union" "Europe & Central Asia"
## [7] "East Asia & Pacific" "Central Europe and the Baltics"
## [9] "Arab World"
unique(dRaw$IndicatorCode)
## [1] "SP.ADO.TFRT" "SP.POP.DPND" "SE.PRM.UNER.FE"
## [4] "SE.PRM.UNER.MA" "SP.DYN.CDRT.IN" "SE.SCH.LIFE.FE"
## [7] "SE.SCH.LIFE.MA" "NY.GDP.PCAP.CD" "NY.GNP.PCAP.CD"
## [10] "SP.DYN.LE00.FE.IN" "SP.DYN.LE00.MA.IN"
在我们可以对数据使用机器学习算法之前，我们需要将数据重新组织为一种格式，其中每一列都只有与一个指标相关的数据。
我们将使用 IndicatorCodes 作为新的列名，并删除对人友好的描述，尽管一旦我们重新组织了数据，让这些描述便于解释列标签是很好的。因为指示器名称列有一个空格，所以我们用刻度线(在键盘的波浪号键上可以找到)来描述列名的开始和结束，并通过赋值 null 来完全删除该列。
dRaw[,'Indicator Name':= NULL]
为了转换数据，我们对原始数据使用melt()函数将年份列折叠成一个名为 year 的变量。这样，关于时间的信息被捕获到一个变量中，而不是不同变量的名称中。这使得所有数值都在一个名为“值”的列中。情况似乎更糟，因为现在要确定任何单个值的含义，必须同时检查指标代码列和年份列。
但是，使用dcast()，我们可以将数据“转换”成正确的结构，其中指标代码是值变量的列，而我们的CountryName和Year列表示随着时间的推移对国家地区的观察。
## collapse columns into a super long dataset
## with Year as a new variable
d <- melt(dRaw, measure.vars = 3:20, variable.name = "Year")
head(d)
## CountryName IndicatorCode Year value
## 1: Sub-Saharan Africa SP.ADO.TFRT 1997 1.3e+02
## 2: Sub-Saharan Africa SP.POP.DPND 1997 9.2e+01
## 3: Sub-Saharan Africa SE.PRM.UNER.FE 1997 2.4e+07
## 4: Sub-Saharan Africa SE.PRM.UNER.MA 1997 2.0e+07
## 5: Sub-Saharan Africa SP.DYN.CDRT.IN 1997 1.6e+01
## 6: Sub-Saharan Africa SE.SCH.LIFE.FE 1997 5.7e+00
str(d)
## Classes 'data.table' and 'data.frame': 1782 obs. of 4 variables:
## $ CountryName : chr "Sub-Saharan Africa" "Sub-Saharan Africa" "Sub-Saharan Africa" "Sub-Saharan Africa" ...
## $ IndicatorCode: chr "SP.ADO.TFRT" "SP.POP.DPND" "SE.PRM.UNER.FE" "SE.PRM.UNER.MA" ...
## $ Year : Factor w/ 18 levels "1997","1998",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ value : num 1.32e+02 9.17e+01 2.44e+07 2.03e+07 1.55e+01 ...
## - attr(*, ".internal.selfref")=<externalptr>
## finally cast the data wide again
## this time with separate variables by indicator code
## keeping a country and time (Year) variable
d <- dcast(d, CountryName + Year ~ IndicatorCode)
head(d)
## CountryName Year NY.GDP.PCAP.CD NY.GNP.PCAP.CD SE.PRM.UNER.FE
## 1: Arab World 1997 2299 2310 6078141
## 2: Arab World 1998 2170 2311 5961001
## 3: Arab World 1999 2314 2288 5684714
## 4: Arab World 2000 2589 2410 5425963
## 5: Arab World 2001 2495 2496 5087547
## 6: Arab World 2002 2463 2476 4813368
## SE.PRM.UNER.MA SE.SCH.LIFE.FE SE.SCH.LIFE.MA SP.ADO.TFRT
## 1: 4181176 8.1 9.7 57
## 2: 4222039 8.3 9.8 56
## 3: 4131775 8.5 10.0 55
## 4: 3955257 8.7 10.0 54
## 5: 3726838 8.8 10.1 53
## 6: 3534138 9.1 10.2 52
## SP.DYN.CDRT.IN SP.DYN.LE00.FE.IN SP.DYN.LE00.MA.IN SP.POP.DPND
## 1: 6.8 69 65 79
## 2: 6.7 69 65 78
## 3: 6.6 69 66 76
## 4: 6.5 70 66 75
## 5: 6.4 70 66 73
## 6: 6.3 70 66 72
str(d)
## Classes 'data.table' and 'data.frame': 162 obs. of 13 variables:
## $ CountryName : chr "Arab World" "Arab World" "Arab World" "Arab World" ...
## $ Year : Factor w/ 18 levels "1997","1998",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ NY.GDP.PCAP.CD : num 2299 2170 2314 2589 2495 ...
## $ NY.GNP.PCAP.CD : num 2310 2311 2288 2410 2496 ...
## $ SE.PRM.UNER.FE : num 6078141 5961001 5684714 5425963 5087547 ...
## $ SE.PRM.UNER.MA : num 4181176 4222039 4131775 3955257 3726838 ...
## $ SE.SCH.LIFE.FE : num 8.08 8.27 8.5 8.65 8.84 ...
## $ SE.SCH.LIFE.MA : num 9.73 9.82 9.97 10.02 10.12 ...
## $ SP.ADO.TFRT : num 56.6 55.7 54.9 54.2 53.3 ...
## $ SP.DYN.CDRT.IN : num 6.8 6.68 6.57 6.48 6.4 ...
## $ SP.DYN.LE00.FE.IN: num 68.7 69 69.3 69.6 69.8 ...
## $ SP.DYN.LE00.MA.IN: num 65 65.3 65.7 65.9 66.2 ...
## $ SP.POP.DPND : num 79.1 77.7 76.2 74.7 73.2 ...
## - attr(*, ".internal.selfref")=<externalptr>
## - attr(*, "sorted")= chr "CountryName" "Year"
现在，数据采用了适合机器学习的格式，其中每一行都是按年份对一个全球区域的观察，所有信息都是在一个变量中捕获的，而不是在变量(列)名称中。这种格式的数据有时被称为“整洁”数据，这是 Hadley Wickham 在 2014 年的一篇文章[113]中深入描述的一个概念。
由于我们的列名本身和我们的区域名可能会被图形化，因此缩短这些名称的长度对于提高视觉清晰度是有价值的。此外，一些算法可能有保留字符(如完整的点或其他标点符号)。在这种情况下，我们的数据的列名被分配给x，然后使用gsub函数删除所有标点符号。接下来，使用函数abbreviate将每个列名减少到四个字符。这些名称通过names函数分配给我们的数据集。最后，国名本身也被缩短了。我们希望我们的读者原谅我们在图形易读性和理解性之间的权衡。注意，在这种情况下，我们使用了left.kept选项来增强理解。
## rename columns with shortened, unique names
x<-colnames(d)
x<-gsub("[[:punct:]]", "", x)
(y <- abbreviate(x, minlength = 4, method = "both.sides"))
## CountryName Year NYGDPPCAPCD NYGNPPCAPCD SEPRMUNERFE
## "CntN" "Year" "NYGD" "NYGN" "SEPR"
## SEPRMUNERMA SESCHLIFEFE SESCHLIFEMA SPADOTFRT SPDYNCDRTIN
## "ERMA" "SESC" "FEMA" "SPAD" "SPDY"
## SPDYNLE00FEIN SPDYNLE00MAIN SPPOPDPND
## "FEIN" "MAIN" "SPPO"
names(d) <- y
## shorten regional names to abbreviations.
d$CntN<-abbreviate(d$CntN, minlength = 5,
method = "left.kept")
我们简要描述了表 7-1 中每列数据所代表的含义。我们首先显示原始的列名，接着是一个|，然后是我们的缩写。
表 7-1
性别数据中的列列表
|变量(特征)
|
描述
|
| — | — |
| CountryName | CntN | 地理区域或国家集团的简称 |
| Year | Year | 每个数据的来源年份 |
| SP.ADO.TFRT | SPAD | 青少年生育率(每 1，000 名 15-19 岁女性的出生率) |
| SP.POP.DPND | SPPO | 受抚养年龄比率(占工作年龄人口的百分比) |
| SE.PRM.UNER.FE | SEPR | 失学儿童，小学，女性 |
| SE.PRM.UNER.MA | ERMA | 失学儿童，小学，男性 |
| SP.DYN.CDRT.IN | SPDY | 粗死亡率(每千人) |
| SE.SCH.LIFE.FE | SESC | 预期受教育年限，女性 |
| SE.SCH.LIFE.MA | FEMA | 预期受教育年限，男性 |
| NY.GDP.PCAP.CD | NYGD | 人均国内生产总值(现值美元) |
| NY.GNP.PCAP.CD | NYGN | 人均国民总收入，阿特拉斯法(现值美元) |
| SP.DYN.LE00.FE.IN | FEIN | 女性出生时预期寿命(岁) |
| SP.DYN.LE00.MA.IN | MAIN | 出生时预期寿命，男性(岁) |
既然数据已经有了合适的结构，summary()函数将为我们提供一些关于各种度量单位的信息。有趣的是，在某些情况下，每列中的数据范围似乎变化很大。另一项需要注意的是，Year 列不是数字，而是一个因子。虽然把它作为一个因素可能有意义，但也可能没有意义。现在，我们用as.character()函数将Year转换成字符串。
summary(d)
## CntN Year NYGD NYGN
## Length:162 1997 : 9 Min. : 496 Min. : 487
## Class :character 1998 : 9 1st Qu.: 3761 1st Qu.: 3839
## Mode :character 1999 : 9 Median : 7458 Median : 7060
## 2000 : 9 Mean :13616 Mean :13453
## 2001 : 9 3rd Qu.:19708 3rd Qu.:19747
## 2002 : 9 Max. :54295 Max. :55010
## (Other):108
## SEPR ERMA SESC FEMA
## Min. : 100024 Min. : 109075 Min. : 5.7 Min. : 7.0
## 1st Qu.: 482710 1st Qu.: 563119 1st Qu.:10.3 1st Qu.:11.2
## Median : 1338898 Median : 1195360 Median :13.3 Median :13.1
## Mean : 3992637 Mean : 3360191 Mean :12.8 Mean :12.8
## 3rd Qu.: 3936040 3rd Qu.: 3339679 3rd Qu.:15.7 3rd Qu.:14.9
## Max. :24437801 Max. :20766960 Max. :17.3 Max. :16.5
##
## SPAD SPDY FEIN MAIN SPPO
## Min. : 11 Min. : 5.0 Min. :52 Min. :48 Min. :41
## 1st Qu.: 21 1st Qu.: 6.0 1st Qu.:72 1st Qu.:68 1st Qu.:49
## Median : 38 Median : 8.1 Median :77 Median :70 Median :51
## Mean : 45 Mean : 8.5 Mean :74 Mean :69 Mean :57
## 3rd Qu.: 53 3rd Qu.:10.6 3rd Qu.:80 3rd Qu.:73 3rd Qu.:63
## Max. :132 Max. :15.5 Max. :84 Max. :78 Max. :92
##
str(d)
## Classes 'data.table' and 'data.frame': 162 obs. of 13 variables:
## $ CntN: chr "ArbWr" "ArbWr" "ArbWr" "ArbWr" ...
## $ Year: Factor w/ 18 levels "1997","1998",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ NYGD: num 2299 2170 2314 2589 2495 ...
## $ NYGN: num 2310 2311 2288 2410 2496 ...
## $ SEPR: num 6078141 5961001 5684714 5425963 5087547 ...
## $ ERMA: num 4181176 4222039 4131775 3955257 3726838 ...
## $ SESC: num 8.08 8.27 8.5 8.65 8.84 ...
## $ FEMA: num 9.73 9.82 9.97 10.02 10.12 ...
## $ SPAD: num 56.6 55.7 54.9 54.2 53.3 ...
## $ SPDY: num 6.8 6.68 6.57 6.48 6.4 ...
## $ FEIN: num 68.7 69 69.3 69.6 69.8 ...
## $ MAIN: num 65 65.3 65.7 65.9 66.2 ...
## $ SPPO: num 79.1 77.7 76.2 74.7 73.2 ...
## - attr(*, ".internal.selfref")=<externalptr>
d[, Year := as.character(Year)]
接下来，我们开始探索我们各种特征之间的关系。我们使用ggplot2包中的ggplot()函数制作一些图表，将人均国民总收入作为输入(x 轴)，青少年生育率作为响应变量(y 轴)。一个是人均美元，另一个是每 1000 名 15-19 岁女性的出生率。plot_grid()函数帮助我们制作了一个图表面板来显示这两个变量和另一个图表，其中我们按数据年份给数据点着色。
## ggplot2 plot object indicating x and y variables
p1 <- ggplot(d, aes(NYGN, SPAD))
## make a grid of two plots
plot_grid(
## first plot data points only
p1 + geom_point(),
## data poins colored by year
p1 + geom_point(aes(colour = Year)) +
scale_colour_viridis(discrete = TRUE),
ncol = 1
)
看图 7-1 ，似乎有明显的群体。下一部分转向机器学习方法，试图准确地确定有多少不同的群体。人们认为，并非我们所有的地区都是截然不同的(例如，欧盟与欧洲和中亚之间可能有重叠)。因此，简单地按地理区域划分每个群体可能是不必要的。相反，我们将检查一些算法，这些算法试图根据经验确定需要多少个聚类或分组来捕获观察到的数据。
图 7-1
每千名妇女的人均国民生产总值和青少年生育率图
720 万美元
比较常见的分组算法之一是kmeans。我们可以使用stats包中名副其实的kmeans()函数在R中执行kmeans集群，该包是少数几个“附带”base R的包之一。
kmeans ()有两个主要参数:包含连续数字数据的数据集x和告诉函数我们希望从数据中提取多少个聚类的中心数k。简而言之，该算法所做的是创建一组良好的 k 中心点，使得组内平方和(欧几里德距离的代表)最小化，给定我们指定的要提取的聚类数量的约束。
该算法的启动方式有一定的可变性，因此如果您希望您的分析具有可重复性，在R中设置随机启动种子是很重要的。我们通过使用值为2468的set.seed()函数来做到这一点，如果您想重现我们的结果，就应该使用这个函数。
接下来，初始化变量wgss以包含组内平方和。对于任何给定数量的聚类，这是算法寻求最小化的度量。
案例k = 1会很无聊，所以我们从两个中心开始，让我们的for()循环在 2 到 9 个中心上迭代算法，用于我们的NYGN对SPAD的数据。对于每次迭代，我们将组内平方和的总和存储在变量中，并绘制我们熟悉的数据图，这一次通过分配的分类成员关系来给点着色。
最后，我们在最后一个图块中添加了一个 scree 图，以显示聚类数与组内平方和的减少量相比的情况。最终的数字如图 7-2 所示。
set.seed(2468)
wgss <- vector("numeric", 8)
plots <- vector("list", 9)
p1 <- ggplot(d, aes(NYGN, SPAD))
for(i in 2:9) {
km <- kmeans(d[, .(NYGN, SPAD)],
centers = i)
wgss[i - 1] <- km$tot.withinss
plots[[i - 1]] <- p1 +
geom_point(aes_(colour = factor(km$cluster))) +
scale_color_viridis(discrete = TRUE) +
theme(legend.position = "none") +
ggtitle(paste("kmeans centers = ", i))
}
plots[[9]] <- ggplot() +
geom_point(aes(x = 2:9, y = wgss)) +
xlab("Number of Clusters") +
ylab("Within SS") +
ggtitle("Scree Plot")
do.call(plot_grid, c(plots, ncol = 3))
在图 7-2 中，我们看到了各种可能的集群数量。请记住，在这种情况下，我们碰巧知道这些数据来自九个地理区域。然而，观察右下方的 scree 图，我们看到在 3 或 4 次之后，组内平方和的减少没有显著改善。作为一种无监督的学习技术，kmeans 向我们展示了这些数据的变化远小于仅基于地理的预期。
图 7-2
不同数量 k 组的人均国民生产总值和每千名妇女的青少年生育率图
注意图 7-2 中的主要因素似乎是 GNP。因为那是最远的距离，所以说它是决定因素是有道理的。kmeans 的目标是最小化每个组的中心和组成该组的点之间的距离。更仔细地看我们的数据，我们看到我们的两个变量之间的规模相当大的差异。也就是说，一个范围超过一百个左右，另一个超过几千个。
summary(d[,.(NYGN, SPAD)])
## NYGN SPAD
## Min. : 487 Min. : 11
## 1st Qu.: 3839 1st Qu.: 21
## Median : 7060 Median : 38
## Mean :13453 Mean : 45
## 3rd Qu.:19747 3rd Qu.: 53
## Max. :55010 Max. :132
通常，数据的规模可能是相对任意的，并不表明在分析中哪个变量应该给予更大的权重。解决这类问题的一种方法是通过居中和缩放。中心数据从每个特征中减去数据的平均值，因此每个特征的平均值等于零。缩放将我们的数据除以每个特征的标准偏差，因此也有一个标准范围。base R中的scale()函数为我们处理这两种操作，并处理矩阵或数据框和数据表样式的数据。
x <- scale(d[,.(NYGN, SPAD)])
summary(x)
## NYGN SPAD
## Min. :-0.92 Min. :-1.04
## 1st Qu.:-0.68 1st Qu.:-0.74
## Median :-0.46 Median :-0.21
## Mean : 0.00 Mean : 0.00
## 3rd Qu.: 0.45 3rd Qu.: 0.26
## Max. : 2.96 Max. : 2.72
现在我们可以重现上次绘制的图表，但使用的是我们缩放后的数据。这表明现在 y 轴高度是分组中的一个区分符(图 7-3 )。
图 7-3
不同数量 k 组的人均国民生产总值和每千名妇女的青少年生育率图
set.seed(2468)
wgss <- vector("numeric", 8)
plots <- vector("list", 9)
p1 <- ggplot(d, aes(NYGN, SPAD))
for(i in 2:9) {
km <- kmeans(x, centers = i)
wgss[i - 1] <- km$tot.withinss
plots[[i - 1]] <- p1 +
geom_point(aes_(colour = factor(km$cluster))) +
scale_color_viridis(discrete = TRUE) +
theme(legend.position = "none") +
ggtitle(paste("kmeans centers = ", i))
}
plots[[9]] <- ggplot() +
geom_point(aes(x = 2:9, y = wgss)) +
xlab("Number of Clusters") +
ylab("Within SS") +
ggtitle("Scree Plot")
do.call(plot_grid, c(plots, ncol = 3))
通过缩放我们的数据，我们使 kmeans 更有可能根据 y 轴的生育率高度来识别各组，尽管在我们图表的左侧仍然很清楚，除了 GNP 之外，肯定还有一些因素可以解释各种生育率。同样，尽管规模很大，但我们从 scree 图中可以看出，“逻辑”组可能比我们选择的九个区域要少。事实上，如果我们稍微思考一下现实世界的地理情况，我们可能会决定六个组可能是最大值，并且查看各种 kmeans 图和 scree 图可能会向我们显示，在 6 之后，组内总平方和值的改善相对较小，甚至四个组也可能足够接近数据。
然而，假设我们决定将中心数设为 6，让我们看看将默认的最大迭代次数从 10 改变会得到什么。虽然在分组中有一些细微的差异，但我们真正注意到的是我们的视觉检查告诉我们的，无论如何，最右边的那个组是相当稀疏定义的(因此它切换组，因为它缺乏强信号)。最大迭代次数控制算法在排序过程中运行的次数。
kmeans取一组被称为质心的 k 随机起始点，然后根据哪个质心与该观测值的欧氏距离最小，为每个值(观测值)分配聚类成员。从那里，每个聚类中的所有点被用于计算新的质心，该质心现在位于该聚类的中心。所以现在有了第二代 k 形心，这意味着每一个点都要再次对照每个形心进行检查，并根据最小欧几里得距离分配聚类成员。如果组成员发生变化，那么每个聚类中的点将用于计算该组的中心，这将成为新的质心，我们现在有了第三代质心。重复该过程，直到组成员不再有变化或者达到最大迭代次数。在这种情况下，似乎由于弱信号，我们有一些点可以向任何一个方向移动，因此连续的迭代会导致一些成员关系的转换。结果如图 7-4 所示。
图 7-4
不同迭代次数下每 1 000 名妇女的人均国民生产总值和青少年生育率图
set.seed(2468)
plots <- vector("list", 9)
p1 <- ggplot(d, aes(NYGN, SPAD))
for(i in 6:14) {
km <- kmeans(x, centers = 6, iter.max = i)
plots[[i - 5]] <- p1 +
geom_point(aes_(colour = factor(km$cluster))) +
scale_color_viridis(discrete = TRUE) +
theme(legend.position = "none") +
ggtitle(paste("kmeans iters = ", i))
}
do.call(plot_grid, c(plots, ncol = 3))
最后一个在kmeans()中定期修改的形式参数是nstart。kmeans 通过测量每个点到组中心的距离来计算组，并试图在每次迭代中使距离测量值变小。然而，对于第一代，在开始时，它会随机设置 k 个中心。因此，nstart决定了有多少不同的随机化竞争成为第一代 k 质心。该算法选择最佳选择作为开始第一代和运行下一次迭代的位置。尽管如此，在本例中，我们没有发现任何重大变化，如图 7-5 所示。
图 7-5
不同 nstart 值的人均国民生产总值和每 1 000 名妇女的青少年生育率图
set.seed(2468)
plots <- vector("list", 9)
p1 <- ggplot(d, aes(NYGN, SPAD))
for(i in 1:9) {
km <- kmeans(x, centers = 6, iter.max = 10, nstart = i)
plots[[i]] <- p1 +
geom_point(aes_(colour = factor(km$cluster))) +
scale_color_viridis(discrete = TRUE) +
theme(legend.position = "none") +
ggtitle(paste("kmeans nstarts = ", i))
}
do.call(plot_grid, c(plots, ncol = 3))
虽然这些示例是在二维空间中完成的，但是我们可以在更高维的数据上执行 kmeans 聚类。虽然不可能通过二维图显示所有 11 个变量，但我们仍然可以观察 scree 图并确定最佳组数。我们执行与之前相同的计算，只是这次只绘制了碎石图。
x <- scale(d[,-c(1,2)])
wgss<-0
set.seed(2468)
for( i in 1:11){
km <- kmeans(x, centers = i)
wgss[i]<-km$tot.withinss
}
ggplot() +
geom_point(aes(x = 1:11, y = wgss)) +
xlab("Number of Clusters") +
ylab("Within SS") +
ggtitle("Scree Plot - All Variables")
根据图 7-6 中的碎石图，我们确定了完整数据集的最佳组数。在这种情况下，我们可能选择 4，现在最后一次运行算法，看看分组会发生什么。结果存储在kmAll变量中，我们使用cbind()以及我们的国家名称和年份将这些结果绑定到我们的数据。
图 7-6
所有人的屏幕图
kmAll <- kmeans(x, centers = 4, nstart = 25)
x <- cbind(d[, c(1,2)], x,
Cluster = kmAll$cluster)
tail(x)
## CntN Year NYGD NYGN SEPR ERMA SESC FEMA SPAD SPDY FEIN MAIN
## 1: Sb-SA 2009 -0.89 -0.87 2.4 2.3 -1.5 -1.5 2.2 1.11 -2.2 -2.2
## 2: Sb-SA 2010 -0.87 -0.87 2.5 2.4 -1.5 -1.5 2.1 0.95 -2.1 -2.1
## 3: Sb-SA 2011 -0.86 -0.86 2.4 2.3 -1.4 -1.4 2.1 0.81 -2.0 -2.0
## 4: Sb-SA 2012 -0.85 -0.84 2.4 2.3 -1.3 -1.4 2.0 0.68 -1.9 -1.9
## 5: Sb-SA 2013 -0.85 -0.84 2.4 2.3 -1.3 -1.3 1.9 0.56 -1.8 -1.8
## 6: Sb-SA 2014 -0.85 -0.83 2.4 2.3 -1.3 -1.4 1.9 0.45 -1.7 -1.7
## SPPO Cluster
## 1: 2.2 1
## 2: 2.2 1
## 3: 2.2 1
## 4: 2.2 1
## 5: 2.1 1
## 6: 2.1 1
如果我们交叉列表显示病例属于一个国家和一个特定群集的频率，我们可以看到一些区域相当一致，而其他区域似乎偶尔改变群集。
xtabs(~ CntN + Cluster, data = x)
## Cluster
## CntN 1 2 3 4
## ArbWr 0 18 0 0
## CEatB 0 0 0 18
## Er&CA 0 0 1 17
## ErpnU 0 0 11 7
## EsA&P 0 14 0 4
## LtA&C 0 18 0 0
## ME&NA 0 18 0 0
## NrthA 0 0 18 0
## Sb-SA 18 0 0 0
这是否意味着我们的算法失败了？首先，由于这是无监督的学习，答案很可能是，即使特定的地理区域没有被数据完全再现，经验结构仍然可以告诉我们一些有用的东西。此外，在这种情况下，我们的数据确实有时间成分。当我们按时间检查结果时，我们看到随着时间的推移，这种变化是一致的，一个区域切换集群成员。
unique(x[
order(CntN, Year, Cluster),
.(CntN, Year, Cluster)][
CntN=="EsA&P"])
## CntN Year Cluster
## 1: EsA&P 1997 2
## 2: EsA&P 1998 2
## 3: EsA&P 1999 2
## 4: EsA&P 2000 2
## 5: EsA&P 2001 2
## 6: EsA&P 2002 2
## 7: EsA&P 2003 2
## 8: EsA&P 2004 2
## 9: EsA&P 2005 2
## 10: EsA&P 2006 2
## 11: EsA&P 2007 2
## 12: EsA&P 2008 2
## 13: EsA&P 2009 2
## 14: EsA&P 2010 2
## 15: EsA&P 2011 4
## 16: EsA&P 2012 4
## 17: EsA&P 2013 4
## 18: EsA&P 2014 4
unique(x[
order(CntN, Year, Cluster),
.(CntN, Year, Cluster)][
CntN == "ErpnU"])
## CntN Year Cluster
## 1: ErpnU 1997 4
## 2: ErpnU 1998 4
## 3: ErpnU 1999 4
## 4: ErpnU 2000 4
## 5: ErpnU 2001 4
## 6: ErpnU 2002 4
## 7: ErpnU 2003 4
## 8: ErpnU 2004 3
## 9: ErpnU 2005 3
## 10: ErpnU 2006 3
## 11: ErpnU 2007 3
## 12: ErpnU 2008 3
## 13: ErpnU 2009 3
## 14: ErpnU 2010 3
## 15: ErpnU 2011 3
## 16: ErpnU 2012 3
## 17: ErpnU 2013 3
## 18: ErpnU 2014 3
正如我们所看到的，kmeans()是一种将相似案例分组在一起的方法，尽管我们的数据出于教学目的而被标记，但并不要求我们的数据有标签。使用 scree 图和for循环，可以确定合理的最佳聚类，以便将数据组织到相似的组中。在某种程度上，这是探索性数据分析的最后一步，因为结果可能会对您组织的客户或您研究的参与者有所启发。但是，请记住，未缩放的数据可能会对某些维度赋予不同的权重，并且算法的目标是最小化欧几里德距离，这可能不会产生人眼可能看到的相同的明显分组。
7.3 分级集群
虽然kmeans开始时选择k个任意中心，并将每个点分配到最近的中心，然后重复该过程一定次数，但层次聚类是不同的。取而代之的是，每个点被分配到它自己唯一的簇中——你有多少个点就有多少个组！然后，确定每个点之间的距离，并找到最近的邻居。邻居然后加入一个更大的群体。这个过程一直重复，直到最后两个超群合并成一个大群。
因为这依赖于距离，所以第一步是计算每个点与其他每个点的距离。距离函数再次默认为欧几里德距离。因此，在最基本的情况下，这也需要连续的数字数据。然而，可以编写定制的距离函数，允许其他类型的“距离”只要输出是距离矩阵，其中较大的数字表示较大的距离，该算法将进行处理(尽管该处理的功效当然不能保证)。
第一阶段是对我们的数据使用dist()函数。为了建立我们的直觉，我们只使用二维数据。
hdist <- dist(d[,.(NYGN, SPAD)])
str(hdist)
## 'dist' num [1:13041] 1.13 22.04 100.03 186.5 166.08 ...
## - attr(*, "Size")= int 162
## - attr(*, "Diag")= logi FALSE
## - attr(*, "Upper")= logi FALSE
## - attr(*, "method")= chr "euclidean"
## - attr(*, "call")= language dist(x = d[, .(NYGN, SPAD)])
生成的对象编码了 2D 数据集中每个点之间的相对距离。该信息被传递给hclust()函数，该函数创建了我们的层次集群，可以通过调用plot()来绘制该集群。这是一个树状图，每个线段的高度显示了我们数据的任意两行/观察值之间的距离。该图如图 7-7 所示。
图 7-7
带行号的聚类树状图
hclust <- hclust(hdist)
plot(hclust)
请注意，这里使用了行名，在这种情况下，行名是数字，没有太大的帮助。我们可以设置一个关键列，以便更好地理解这个树状图。虽然这不是最清晰的图，但人们确实注意到，在图 7-8 中，相似的区域似乎经常彼此紧密相连。
图 7-8
通过 ape 包改变树状图
x <- d[, .(CntN, Year, NYGN, SPAD)]
x[, Key := paste(CntN, Year)]
x[, CntN := NULL]
x[, Year := NULL]
hdist <- dist(x[,.(NYGN, SPAD)])
hclust <- hclust(hdist)
plot(hclust, labels = x$Key)
现在，虽然kmeans可能更简单，但它显示了一个最终结果，即一个人选择多少组就有多少组。因此，scree 图对于确定多少个集群可能有意义是有价值的。在这里，每一行都从自己的组开始，然后每一个邻居都被系统地连接成对。
高度向我们显示了任何两个相连组之间的距离。此外，请注意，对于给定的高度，我们可以了解我们的数据适合多少“组”。
plot(hclust, labels = x$Key)
abline(h = 30000, col = "blue")
在h = 30,000的高度，如图 7-9 所示，我们只有两组。一组包括北美(NrtA)和 2000 年代的大部分欧盟国家(普尔)，而另一组包括其余国家。看一下摘要以及对我们的父数据的几个数据表调用表明，主要高度(记住，这是基于欧几里德距离的)可能是 GNP (NYGN)。请注意，这很好地说明了为什么一旦我们完成了探索性的工作，扩展我们的数据仍然很重要。如果我们不努力确保数据列大致相等，GNP (NYGN)之类的东西可能会淹没其他显著特征。
图 7-9
带有国家名称、年份和高度线的聚类树状图
summary(x)
## NYGN SPAD Key
## Min. : 487 Min. : 11 Length:162
## 1st Qu.: 3839 1st Qu.: 21 Class :character
## Median : 7060 Median : 38 Mode :character
## Mean :13453 Mean : 45
## 3rd Qu.:19747 3rd Qu.: 53
## Max. :55010 Max. :132
d[, mean(NYGN), by = CntN][order(V1)]
## CntN V1
## 1: Sb-SA 953
## 2: ArbWr 4260
## 3: ME&NA 5045
## 4: EsA&P 5801
## 5: LtA&C 6004
## 6: CEatB 8531
## 7: Er&CA 19021
## 8: ErpnU 28278
## 9: NrthA 43188
d[, mean(SPAD), by = CntN][order(V1)]
## CntN V1
## 1: ErpnU 15
## 2: EsA&P 20
## 3: CEatB 23
## 4: Er&CA 23
## 5: NrthA 37
## 6: ME&NA 40
## 7: ArbWr 52
## 8: LtA&C 73
## 9: Sb-SA 120
调整我们的abline()的高度允许我们改变我们查看数据的方式，现在它在四个集群中。现在，特别是在书本形式中，很难清晰地打印出如此大的图表，尽管我们最好的尝试是在图 7-10 中。试着在你自己的机器上运行代码，记住我们还处于探索阶段。这有助于我们理解国家之间的相似性，以及这两个维度的结构。
图 7-10
带有国家名称、年份和另一条高度线的聚类树状图
plot(hclust, labels = x$Key)
abline(h = 20000, col = "blue")
考虑到维度，认识到这个模型非常适合可视化整个数据集——在所有列上。它只需要对我们的代码做一点小小的调整，而不是将变量限制为两个。当然，高度测量将发生巨大变化，因为我们现在有一个更大的空间来填充我们的欧几里德距离。此外，如图 7-11 所示，最后两个分组将撒哈拉以南非洲(S-SA)与我们的其他地区联系在一起。
图 7-11
带有国家名称、年份和所有维度数据的聚类树状图
x <- copy(d)
x[, Key := paste(CntN, Year)]
x[, CntN := NULL]
x[, Year := NULL]
hdist <- dist(x[, -12])
hclust <- hclust(hdist)
plot(hclust, labels = x$Key)
现在，hclust函数将距离矩阵作为第一个形式参数，将方法类型作为第二个形式参数。第二个形式参数默认为“complete ”,它使用聚类点之间的最大距离来确定组距离，然后总是根据该最大距离对最近的组进行分组。其他方法将产生不同的结果，因为它们选择使用其他方法(例如最小化类内方差的迭代增加)。请注意，这里我们没有刷新我们的hdist矩阵——我们只是简单地改变了使用的方法。结果如图 7-12 所示。
图 7-12
使用 ward 聚类树图。D2 方法
hclust <- hclust(hdist, method = "ward.D2")
plot(hclust, labels = x$Key)
值得注意的是，欧几里德距离可以被替换为一些其他的距离度量；事实上，有几个内置选项。除了内置函数之外，在不太深奥的情况下，使用一些对推文进行情感分析的函数来确定某个特定趋势标签有多少组意见可能是有意义的。当然，一旦情感分析被转换成一些数字，在这些数字上使用欧几里德距离也许是有意义的。在任何情况下，关键是，不要求数据从数字级别开始。
如前所述，尤其是当我们在分析中增加更多维度时，如果不调整数据，范围最大的列将对我们的分组产生巨大的影响。缩放后，使用所有变量并重新绘图，结果如图 7-13 所示。
图 7-13
带标度的聚类树状图
x <- scale(d[,-c(1,2)])
row.names(x) <- paste(d$CntN, d$Year)
hdist <- dist(x)
hclust <- hclust(hdist)
plot(hclust, labels = paste(d$CntN, d$Year))
abline(h = 6, col = "blue")
现在，虽然我们可以选择根据一定的高度来砍伐我们的树，但这只是直观地向我们展示了我们的群体是如何排列的。我们还可以使用cutree()函数在某个高度对我们的数据进行聚类。在 h = 6 的情况下，这将我们的树分成三组，我们可以看到这个函数给了我们。
cut_hclust <- cutree(hclust, h = 6)
unique(cut_hclust)
## [1] 1 2 3
或者，不是在特定的高度上切割，而是在一定数量的簇上切割。通过创建数据的副本，我们可以在名为 cluster 的新列中记录集群分配。
dcopy <- as.data.table(copy(d))
dcopy[, cluster:= NA_integer_]
dcopy$cluster <- cutree(hclust, k = 3)
tail(dcopy)
## CntN Year NYGD NYGN SEPR ERMA SESC FEMA SPAD SPDY FEIN MAIN
## 1: Sb-SA 2009 1198 1186 1.9e+07 1.6e+07 8.2 9.4 115 11.5 58 55
## 2: Sb-SA 2010 1555 1287 2.0e+07 1.6e+07 8.3 9.4 113 11.0 58 55
## 3: Sb-SA 2011 1706 1412 1.9e+07 1.5e+07 8.4 9.6 111 10.7 59 56
## 4: Sb-SA 2012 1740 1631 1.9e+07 1.6e+07 8.6 9.7 109 10.3 60 57
## 5: Sb-SA 2013 1787 1686 1.9e+07 1.5e+07 8.9 9.9 107 10.0 61 57
## 6: Sb-SA 2014 1822 1751 1.9e+07 1.6e+07 8.7 9.7 105 9.7 61 58
## SPPO cluster
## 1: 89 3
## 2: 89 3
## 3: 88 3
## 4: 88 3
## 5: 87 3
## 6: 87 3
为了结束我们关于层次化集群的部分，我们提到了ape包，它有几个可视化选项，可以显示默认情况下由plot()提供的树状图。诚然，用较少的类别来可视化模型更容易。请记住，在无监督学习的通常应用中，一个自然的目标是确定合理存在多少个类别。所以我们的教学例子将不仅仅是有点笨拙，正是因为我们保持了一个相当广泛的，固定的类别集。我们在图 7-14 中展示了结果。
图 7-14
通过 ape 包改变树状图
plot(as.phylo(hclust), type = "cladogram")
plot(as.phylo(hclust), type = "fan")
plot(as.phylo(hclust), type = "radial")
分级聚类的部分优势在于，它们可以是一种可视化观察结果之间的相似性和不相似性的方式。作为这个过程的一部分，将树分成清晰的组可能是一个有用的步骤。现在，基于我们早期的 kmeans 分析，我们相信可能有四个重要的组，所以我们尝试可视化。这里，cutree()的第二个形式参数被设置为用k = 4期望的集群数。此外，as.phylo()函数用于将我们的hclust对象转换成供ape包使用的对象类型(一个phylo对象)。我们展示了一种新的图表类型— unrooted—和一个label.offset，它在我们的标签和图表之间提供了一些距离(从而在我们的人类友好但适合图表不友好标签的字符之间创建了一些非常重要的空间)。最后，我们在tip.color参数中引入cutree()信息，并用cex缩小文本的放大倍数。最终结果如图 7-15 所示。
图 7-15
通过 ape 包改变树状图
hclust4 <- cutree(hclust, k = 4)
plot(as.phylo(hclust), type = "unrooted", label.offset = 1,
tip.color = hclust4, cex = 0.8)
请注意，该图本身是基于原始的层次聚类模型的—只有尖端的颜色用于将我们的数据分组为配对(但是这些配对在该图中是有意义的)。
图 7-16
通过 ape 包改变树状图
7.4 主成分分析
到目前为止，在无监督学习中，我们已经看到了两种技术来了解我们的数据中有多少组或簇。一个是kmeans，它被预先告知应该找到多少个组(可能基于 scree plot 分析)。另一个是hclust，将每个观察结果放在单个组中，然后最终将这些点连接起来，直到只有一个组。分析员要看哪些观察值最接近，哪些子群最接近。主成分分析(PCA)也试图以某种方式确定组。PCA 将数据分解成独特的(即不相关的、独立的、正交的)分量。例如，在我们的数据集中，我们可以想象 GDP 和 GNP 在某种程度上是相关的。其实它们是高度相关的！这些本质上都是一样的。甚至他们的范围和手段都相当接近，如图 7-18 所示。
图 7-17
四个集群上的无根类型
cor(d$NYGD, d$NYGN)
## [1] 1
summary(d[,.(NYGD, NYGN)])
## NYGD NYGN
## Min. : 496 Min. : 487
## 1st Qu.: 3761 1st Qu.: 3839
## Median : 7458 Median : 7060
## Mean :13616 Mean :13453
## 3rd Qu.:19708 3rd Qu.:19747
## Max. :54295 Max. :55010
ggplot(d, aes(NYGD, NYGN)) +
geom_point()
从视觉上看，图 7-18 本质上是直线 y = x 。主成分分析可以认为是一种分组操作。分组的目的是看看我们有多少真正的独特维度。在这种情况下，虽然看起来我们在数据集中有两个不同的列，但事实是我们实际上只有一个唯一信息列(维度)似乎是合理的。
图 7-18
一个高度相关的图——这里真的有二维吗？
主成分分析(PCA)允许我们做的是确定我们不需要两列，并以有原则的方式将它们组合成一个维度。这将简化我们的特征空间，正如我们所看到的，不会导致太多的信息损失。另一方面，如果我们看看其他一些变量，我们会发现简单地删除一列以支持另一列是没有意义的。即使它们高度相关，仍然清楚的是，有时，特别是在某一点上，SE.SCH.Life.FE (SESC)有相当多的变化，而 GDP(纽约市)变化很小，如图 7-19 所示。显然，这里有不止一个维度的信息。
图 7-19
一个高度相关的图——这里真的有二维吗？
ggplot(d, aes(NYGD, SESC)) +
geom_point()
cor(d$NYGD, d$SESC)
## [1] 0.79
与 PCA 相关的技术是因子分析(FA)。这两种方法都具有探索性，因为数据的真实维度是未知的，尽管存在验证性 FA，其目标是测试假设的维度。此外，PCA 和 FA 都试图找到一个更低维的空间来提供数据的合理近似。但是，PCA 和 FA 也有很多不同之处。它们来自不同的理论背景，对主成分(PCA)和因子(FA)的解释也不相同，对数据的基本假设也不同。一个显著的区别是 FA 通常关注变量间的共享方差，而 PCA 包含共享和唯一方差。这种差异是由于 PCA 和 FA 的基础和目标不同。FA 来源于心理测量学的传统，通常与专门设计的测试一起使用，例如，分析考试中的不同问题或旨在测量智商或评估某些心理结构的多个问题。在一次考试中，所有问题的目标是提供一个学生对课程整体理解的指数，这是无法直接观察到的。与任何给定测试问题特有的方差相比，问题间共享的方差被认为是整体理解的更好指标，任何给定测试问题特有的方差可能代表措辞不当的问题、某个特定概念的糟糕教学或对特定概念缺乏理解。在 PCA 中，目标通常是用尽可能少的维度再现更高维的空间。因此，如果考试中的一个问题与任何其他问题没有太多重叠，PCA 不会认为这是整体表现的潜在不良指标，而是认为这是另一个需要的独特维度。一般来说，PCA 在机器学习中使用得更多，因为它较少建立在任何特定的理论或信念上，即一组项目应该有一些共享的重叠，并且因为使用 PCA 通常更容易包括足够多的成分，使得原始数据几乎可以完美地恢复，但是是从更小的维度集恢复。
PCA 所做的是查看我们的数据，并将左右移动数据的部分与上下移动数据的部分分开，换句话说，分成主要成分。如果你有线性代数的背景，标准 PCA 是协方差(或者如果标准化，相关性)矩阵的特征值分解。无论如何，我们都要回到我们的主题，看看我们感兴趣的两列。再次参考这两部的原著剧情。
pca()函数既有对scale数据的参数，如果设置为“uv ”,则将方差设置为 1(单位方差),如果设置为TRUE,则将数据设置为零center。传统类型的 PCA 基于特征值和奇异值分解(SVD)，因此对于传统 PCA，我们告诉pca()函数使用method = "svd"。
首先，我们只收集国民生产总值的原始数据。PCAP 对阿多。TFRT 进入我们的工作数据集，x。计算主成分分析非常简单，我们可以立即缩放和集中我们的数据，并使用奇异值分解来估计传统的主成分分析。summary()函数向我们显示了每个主成分可以单独解释总方差的多少，另一行显示了解释的累积方差。在只有两个变量的情况下，100%的方差， R 2 = 1，可以用两个主成分来解释。
x <- d[,.( NYGN, SPAD)]
res <- pca(x, method="svd", center=TRUE, scale = "uv")
summary(res)
## svd calculated PCA
## Importance of component(s):
## PC1 PC2
## R2 0.7213 0.2787
## Cumulative R2 0.7213 1.0000
为了将我们的结果可视化，我们可以创建一个图表，显示原始数据是如何转换为双标图的，如图 7-20 所示。你能认出正交旋转吗？当然，不仅仅是轮换在起作用。PC1 变量是可以用一条直线解释的最大方差。差异的剩余部分预计到 PC2。biplot()函数显示新空间中原始数据值的向量。
图 7-20
原始数据和主成分分析的比较
biplot(res, main = "Biplot of PCA")
接下来，我们将检查当我们使用所有变量进行 PCA 时会发生什么。我们可以绘制一个 scree 图，看看当我们添加额外的主成分时，精度如何提高，作为一种尝试确定需要多少成分才能充分代表数据中的基本维度的方法。我们为数据中的列数设置了一个额外的参数,nPcs,,这是我们可能需要的主成分的最大可能数量，尽管如果我们数据的真实维度更少，那么很少的成分可能足以代表数据。
我们使用图 7-21 中的plot()函数绘制了从较少维度恢复的原始信息的比例图，这是一种反向 scree 图。再次查看 scree 图，我们会发现这里实际上只有 4 个维度的独特数据，即使只有一个维度也能够捕捉 11 个变量中超过一半的方差。
图 7-21
数据中所有特征的传统 PCA 的 Scree 图
x <- d[, -c(1,2)]
res <- pca(x, method="svd", center=TRUE, scale = "uv",
nPcs = ncol(x))
summary(res)
## svd calculated PCA
## Importance of component(s):
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## R2 0.7278 0.1614 0.06685 0.02554 0.01226 0.00447 0.0011
## Cumulative R2 0.7278 0.8892 0.95607 0.98161 0.99387 0.99834 0.9994
## PC8 PC9 PC10 PC11
## R2 0.00021 0.0002 1e-04 5e-05
## Cumulative R2 0.99965 0.9999 1e+00 1e+00
## reverse scree plot
ggplot() +
geom_bar(aes(1:11, cumsum(res@R2)),
stat = "identity") +
scale_x_continuous("Principal Component", 1:11) +
scale_y_continuous(expression(R²), labels = percent) +
ggtitle("Scree Plot") +
coord_cartesian(xlim = c(.5, 11.5), ylim = c(.5, 1),
expand = FALSE)
我们也可以从这个 PCA 创建双绘图，但是一次只能绘制两个组件。图 7-22 显示了前两个组成部分的结果，它们共同解释了大部分的差异。
图 7-22
前两个主成分的双标图
biplot(res, choices = c(1, 2))
顺便提一下，现在我们在res变量中有了我们的数据，我们可以使用该数据(以及尽可能多的列，比如前四列)来代替我们的原始数据。特别是，如果我们希望执行监督机器学习技术，正如下一章所讨论的，使用 PCA 数据可能会有优势。例如，来自我们的 PCA 的结果表明，我们可以用四个分量捕获所有 11 个原始特征中除 2%之外的所有变化。
使用主成分代替原始特征的优点是主成分是正交的。缺点是，由于每一个主成分分析都可能包含几个原始变量，因此很难准确解释每一个特征意味着什么或代表什么。但是，如果模型的预测能力得到提高或者是唯一的目标，这可能是我们愿意做出的权衡。在某些情况下，它还可以加速模型计算，因为后期分析只需管理较少的要素，从而减少分析所需的内存和计算。
我们使用scores()函数提取主成分得分，查看前几行，然后使用相关矩阵显示它们确实是线性独立的。
head(scores(res))
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
## [1,] -2.5 -1.2 0.74 -0.31 0.99 -0.075 -0.22 0.0046 0.0015 -0.00660
## [2,] -2.4 -1.2 0.72 -0.31 0.89 -0.057 -0.21 0.0075 -0.0035 0.00095
## [3,] -2.3 -1.3 0.70 -0.29 0.82 -0.047 -0.19 0.0121 -0.0046 -0.00521
## [4,] -2.2 -1.3 0.69 -0.29 0.75 -0.074 -0.18 0.0202 -0.0086 -0.01269
## [5,] -2.0 -1.3 0.66 -0.27 0.70 -0.095 -0.15 0.0290 -0.0149 -0.00243
## [6,] -1.9 -1.4 0.62 -0.25 0.64 -0.104 -0.13 0.0355 -0.0237 -0.00071
## PC11
## [1,] -9.9e-03
## [2,] 4.7e-06
## [3,] 1.3e-02
## [4,] 1.1e-02
## [5,] 1.3e-02
## [6,] 1.3e-02
round(cor(scores(res)),2)
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11
## PC1 1 0 0 0 0 0 0 0 0 0 0
## PC2 0 1 0 0 0 0 0 0 0 0 0
## PC3 0 0 1 0 0 0 0 0 0 0 0
## PC4 0 0 0 1 0 0 0 0 0 0 0
## PC5 0 0 0 0 1 0 0 0 0 0 0
## PC6 0 0 0 0 0 1 0 0 0 0 0
## PC7 0 0 0 0 0 0 1 0 0 0 0
## PC8 0 0 0 0 0 0 0 1 0 0 0
## PC9 0 0 0 0 0 0 0 0 1 0 0
## PC10 0 0 0 0 0 0 0 0 0 1 0
## PC11 0 0 0 0 0 0 0 0 0 0 1
除了特征值和奇异值分解的传统 PCA 方法外，pca()函数还有多种其他方法可用于各种其他情况。一种这样的方法是robustPca，它改变算法，试图对极端值或异常值更加鲁棒。为了查看其效果，我们将使用相同的数据，使用prep()函数对其进行预缩放，然后创建两个版本，一个没有异常值，一个有异常值。
接下来，我们运行四个主成分分析模型，传统的奇异值分解主成分分析和鲁棒主成分分析。最后，我们使用常规和稳健的 PCA 方法，从无异常值和有异常值的数据中，绘制负载，即用于将原始数据投影到主成分上的值。结果如图 7-23 所示。在该图中，您可以清楚地看到传统 PCA 中的负载如何显著变化，但在鲁棒 PCA 中，异常值的添加对负载几乎没有影响。一般来说，如果您担心您的数据可能有极值或异常值，那么首先清理或删除它们是值得的，或者至少比较传统和稳健 PCA 方法的结果，以确保结果不会对少数极值分数的存在过于敏感。
图 7-23
使用传统的奇异值分解主成分分析和稳健主成分分析绘制有和没有异常值的主成分分析模型的负荷图
x <- d[, -c(1,2)]
x <- prep(x, center = TRUE, scale = "uv")
xout <- copy(x)
xout[1:5, "NYGD"] <- (-10)
res1 <- pca(x, method = "svd",
center = FALSE, nPcs = 4)
res2 <- pca(xout, method = "svd",
center = FALSE, nPcs = 4)
res1rob <- pca(x, method = "robustPca",
center = FALSE, nPcs = 4)
res2rob <- pca(xout, method = "robustPca",
center = FALSE, nPcs = 4)
plot_grid(
ggplot() +
geom_point(aes(
x = as.numeric(loadings(res1)),
y = as.numeric(loadings(res2)))) +
xlab("Loadings, SVD, No Outliers") +
ylab("Loadings, SVD, Outliers"),
ggplot() +
geom_point(aes(
x = as.numeric(loadings(res1rob)),
y = as.numeric(loadings(res2rob)))) +
xlab("Loadings, Robust PCA, No Outliers") +
ylab("Loadings, Robust PCA, Outliers"),
ncol = 1)
7.5 非线性聚类分析
主成分分析的一个常见目标是减少维数，从而简化模型并提高预测精度。PCA 的一个基本假设是在这些较少的维度中存在可感测的线性关系，换句话说，对于 n 个标称维度，存在一些基本的正交向量要被投影，而不会有大的信号损失。如果有理由怀疑情况并非如此，那么非线性方法可能是合适的。
虽然前四行代码在分层集群中很常见，但是MASS包中的sammon()函数获取这些距离，并尝试将更高维度(在我们的示例中为 11 维)映射到k维(在本例中，我们设置 k = 2 以绘制图表)。我们提到这一点的主要目的只是为了观察，在明显非线性关系的情况下，PCA 可能不如放松线性假设的其他方法有效。
x <- scale(d[, -c(1,2)])
row.names(x) <- paste(d$CntN, d$Year)
head(x)
## NYGD NYGN SEPR ERMA SESC FEMA SPAD SPDY FEIN MAIN
## ArbWr 1997 -0.81 -0.79 0.33 0.156 -1.5 -1.4 0.37 -0.66 -0.74 -0.62
## ArbWr 1998 -0.82 -0.79 0.31 0.164 -1.5 -1.3 0.34 -0.71 -0.70 -0.57
## ArbWr 1999 -0.81 -0.79 0.27 0.147 -1.4 -1.2 0.32 -0.75 -0.67 -0.53
## ArbWr 2000 -0.79 -0.79 0.22 0.113 -1.3 -1.2 0.30 -0.78 -0.63 -0.48
## ArbWr 2001 -0.80 -0.78 0.17 0.070 -1.3 -1.2 0.27 -0.81 -0.60 -0.44
## ArbWr 2002 -0.80 -0.78 0.13 0.033 -1.2 -1.1 0.24 -0.84 -0.57 -0.41
## SPPO
## ArbWr 1997 1.6
## ArbWr 1998 1.4
## ArbWr 1999 1.3
## ArbWr 2000 1.2
## ArbWr 2001 1.1
## ArbWr 2002 1.0
sdist <- dist(x)
xSammon <- sammon(sdist, k = 2)
## Initial stress : 0.04343
## stress after 7 iters: 0.03619
head(xSammon$points)
## [,1] [,2]
## ArbWr 1997 -2.6 -1.2
## ArbWr 1998 -2.5 -1.2
## ArbWr 1999 -2.3 -1.3
## ArbWr 2000 -2.2 -1.3
## ArbWr 2001 -2.1 -1.4
## ArbWr 2002 -2.0 -1.4
Sammon 技术试图最小化一个名为stress的度量，该度量试图测量高维对象被“挤压”到低维空间的效率。
如果将前面的代码调整为 k = 3，我们可以预期应力会降低。图 7-24 显示了使用 Sammon 将 11 维减少到只有 2 维的结果。
图 7-24
带文本标签的 Sammon 点图
plot(xSammon$points, type = "n")
text(xSammon$points, labels = row.names(x) )
7.6 摘要
在这一章中，我们发展了无监督机器学习的概念。虽然采用这种技术的数据通常是未标记的，但在这种情况下，我们用标记数据进行了演示，以便对我们可能预期的聚类或分组类型有所了解。常用的技术kmeans和hclust与标准欧几里德距离一起使用(从而将我们限制在数字数据)。然而，这些强大的技术使我们能够理解在一个复杂的多维数据集中可能有多少真正的组。最后，介绍了降维的概念，使用主成分分析作为最常用的方法，同时以简单的观察结束，即主成分分析的线性假设可以通过更复杂的技术放宽。表 7-2 中概述了本章使用的一些功能，并简要介绍了它们的作用。
表 7-2
本章中描述的关键功能列表及其功能摘要
|功能
|
它的作用
|
| — | — |
| read_excel() | 将文件路径作为字符串，并读入 Excel 文件 |
| as.data.table() | 将 data.frame 转换为 data.table 对象 |
| str() | 显示了 R 中对象的底层结构 |
| summary() | 尝试对数据或模型对象进行统计汇总 |
| melt() | 将宽数据转换为长数据 |
| dcast() | 将熔化的长数据转换为宽格式 |
| sort() | 按升序对数据进行排序 |
| unique() | 删除重复项 |
| plot() | 数据和模型结果的一般绘图 |
| set.seed() | 允许为伪随机算法复制代码 |
| kmeans() | 运行 kmeans 算法 |
| scale() | 在平均值= 0 时将数据按列居中，并调整到单位标准偏差 |
| for() | 调用迭代器上的 for 循环操作 |
| cbind() | 通过列将数据绑定在一起 |
| dist() | 创建一个距离矩阵，显示每个元素与其他元素之间的欧几里德距离 |
| hclust() | 创建分层集群对象 |
| abline() | 在绘图对象上绘制一条线 |
| row.names() | 允许访问 data.frame 行名属性 |
| cutree() | 按高度或簇切割层次树对象 |
| tail() | 显示最后六行数据 |
| copy() | 复制数据对象，而不仅仅是通过引用指定一个新名称 |
| as.phylo() | 创建一个 phylo 对象—用于 ape 包绘图 |
| cor() | 显示项目之间的相关性 |
| pca() | 执行 PCA 简化计算 |
| scores() | 提取主成分得分 |
| loadings() | 提取主成分载荷，用于将原始数据投影到主成分空间 |
| biplot() | 在 PC1 和 PC2 上绘制 PCA 图并显示原始数据向量 |
| sammon() | 非线性降维算法；第二种形式是新的维度 |
| text() | 地块对象上的文字标签 |
八、ML：监督机器学习和分类
在前一章中，数据被输入一个算法，然后该算法试图对常见类型的数据进行分组。虽然我们偶尔会保留标签，以便了解算法的执行情况，但在现实生活中，无监督的机器学习是一种探索性分析。有时称为预处理，它通常是一个初始阶段，包括数据的缩放和居中。
通常，目标不是简单的分组，而是使用当前数据预测未来的能力。机器学习的常见目标是训练具有良好预测能力的模型。什么是可接受的预测准确度水平将因应用而异。其中一位作者在一所社区大学工作，在那里，即使对学生的最终字母分数等事情的预测准确性水平很低，也可能允许更相关、更有针对性的干预发生(例如，关于辅导中心时间或分配给导师或教练的短信)。在那种情况下，假阳性的风险可能很小。相反，另一位作者研究健康结果数据，其中假阳性对患者整体健康带来更高水平的可能风险。
对于这一章，新的软件包包括caret [32]，其中包含了用于准备数据和模型开发的有用功能。kernlab包用于支持向量机【49】。结合caret，我们将使用DALEX【11】，这将有助于解释模型(并需要spdep【12】)。以类似的方式，rattle [119]允许额外的可视化。最后，我们还将包RSNNS [10]用于多层感知器模型——一种神经网络。
关于caret的一个注意事项是它不会自动安装特定分析方法所需的软件包。因此，当您逐步完成后面的分析时，预计需要安装一些包。为了减轻这种情况，我们在这里包括了包ranger【124】、e1071【62】、gbm【120】和plyr【110】。
library(checkpoint)
checkpoint("2018-09-28", R.version = "3.5.1",
project = book_directory,
checkpointLocation = checkpoint_directory,
scanForPackages = FALSE,
scan.rnw.with.knitr = TRUE, use.knitr = TRUE)
library(ggplot2)
library(cowplot)
library(data.table)
library(readxl)
library(viridis)
library(RSNNS)
library(kernlab)
library(rpart)
library(rattle)
library(DALEX)
library(caret)
library(spdep)
library(ranger)
library(e1071)
library(gbm)
library(plyr)
set.seed(1234)
options(width = 70, digits = 2)
8.1 数据准备
对于监督学习的大多数应用，数据需要某种程度的准备。我们使用与上一章相同的初始步骤开始这个过程，读入我们的样本数据集，并将其转换为我们现在熟悉的布局。回想一下,'是标签上方波浪号键上的勾号。
## Note: download Excel file from publisher website first
dRaw <- read_excel("Gender_StatsData_worldbank.org_ccby40.xlsx")
dRaw <- as.data.table(dRaw) # convert data to data.table format.
dRaw[,'Indicator Name':= NULL]
## collapse columns into a super long dataset
## with Year as a new variable
data <- melt(dRaw, measure.vars = 3:20, variable.name = "Year", variable.factor = FALSE)
## cast the data wide again
## this time with separate variables by indicator code
## keeping a country and time (Year) variable
data <- dcast(data, CountryName + Year ~ IndicatorCode)
rm(dRaw) #remove unneeded variable
#rename columns with shortened, unique names
x<-colnames(data)
x<-gsub("[[:punct:]]", "", x)
(y <- abbreviate(x, minlength = 4, method = "both.sides"))
## CountryName Year NYGDPPCAPCD NYGNPPCAPCD SEPRMUNERFE
## "CntN" "Year" "NYGD" "NYGN" "SEPR"
## SEPRMUNERMA SESCHLIFEFE SESCHLIFEMA SPADOTFRT SPDYNCDRTIN
## "ERMA" "SESC" "FEMA" "SPAD" "SPDY"
## SPDYNLE00FEIN SPDYNLE00MAIN SPPOPDPND
## "FEIN" "MAIN" "SPPO"
names(data) <- y
#shorten regional names to abbreviations.
data$CntN<-abbreviate(data$CntN, minlength = 5, method = "left.kept")
我们再次简要描述表 8-1 中每列数据所代表的内容。
警告:虽然许多现代软件包将执行我们在本节中试图概括和可视化的所有或部分步骤，但它有助于理解某些预处理技术背后的价值。某些度量标准的使用取决于所选的模型和输入模型的基础数据。例如，树方法可以很好地处理字符或因子数据，而传统的线性回归需要数字输入。因此，即使在本章中，期望讨论的每个模型以完全相同的方式使用我们的固定数据集也是没有意义的。因此，当我们通过这些方法工作时，我们将保留data作为上一章工作数据集的持有者。我们将对其副本执行额外的处理，并在合理的情况下返回到data。
表 8-1
性别数据中的列列表
|变量(特征)
|
描述
|
| — | — |
| CountryName | CntN | 地理区域或国家集团的简称 |
| Year | Year | 每个数据的来源年份 |
| SP.ADO.TFRT | SPAD | 青少年生育率(每 1，000 名 15-19 岁女性的出生率) |
| SP.POP.DPND | SPPO | 受抚养年龄比率(占工作年龄人口的百分比) |
| SE.PRM.UNER.FE | SEPR | 失学儿童，小学，女性 |
| SE.PRM.UNER.MA | ERMA | 失学儿童，小学，男性 |
| SP.DYN.CDRT.IN | SPDY | 粗死亡率(每千人) |
| SE.SCH.LIFE.FE | SESC | 预期受教育年限，女性 |
| SE.SCH.LIFE.MA | FEMA | 预期受教育年限，男性 |
| NY.GDP.PCAP.CD | NYGD | 人均国内生产总值(现值美元) |
| NY.GNP.PCAP.CD | NYGN | 人均国民总收入，阿特拉斯法(现值美元) |
| SP.DYN.LE00.FE.IN | FEIN | 女性出生时预期寿命(岁) |
| SP.DYN.LE00.MA.IN | MAIN | 出生时预期寿命，男性(岁) |
一个热编码
我们看到在CntN(地区名称)列中有九种独特的分类数据类型。要将名义数据转换成适用于回归或计算的数据(换句话说就是数字)，基本的技术是从名义类别中创建一组新的列，如果一行不适合，则包含值0，如果一行适合，则包含值1。如果我们有九个值，我们可以创建九列，每列以ArbWr(阿拉伯世界)开始，以S-SA(撒哈拉以南非洲)结束。在这种情况下，每行将有一个1，其余的留在0。这有时被称为one hot encoding。或者，由于我们只有九个选项，我们可以只对前八个选项进行编码。如果前八个都是0，这仍然会向我们的模型发出第九个选项的信号。这种方法有时被称为dummy coding。实际上，这两个术语似乎可以互换使用。不同的型号可能偏好不同的布局。像往常一样，花一些时间从可接受的输入方面理解特定模型的需求。
d <- copy(data)
sort(unique(d$CntN))
## [1] "ArbWr" "CEatB" "Er&CA" "ErpnU" "EsA&P" "LtA&C" "ME&NA" "NrthA"
## [9] "Sb-SA"
我们使用来自caret的dummyVars函数为我们执行一个热编码。这个函数接受一个常用的R公式的输入，这里我们要求它对整个数据集进行操作，而不是只对其中的一部分进行操作。该函数将忽略数字数据，但它将转换字符和因子数据。从str函数可以看出，我们的年份信息是一个字符类型。如果不加处理，年份也会被虚拟编码。在这种情况下，我们继续转换到数字，希望有一个体面的适合没有。不可否认，目前这仅仅是基于保持低的总列数的愿望。
str(d)
## Classes 'data.table' and 'data.frame': 162 obs. of 13 variables:
## $ CntN: chr "ArbWr" "ArbWr" "ArbWr" "ArbWr" ...
## $ Year: chr "1997" "1998" "1999" "2000" ...
## $ NYGD: num 2299 2170 2314 2589 2495 ...
## $ NYGN: num 2310 2311 2288 2410 2496 ...
## $ SEPR: num 6078141 5961001 5684714 5425963 5087547 ...
## $ ERMA: num 4181176 4222039 4131775 3955257 3726838 ...
## $ SESC: num 8.08 8.27 8.5 8.65 8.84 ...
## $ FEMA: num 9.73 9.82 9.97 10.02 10.12 ...
## $ SPAD: num 56.6 55.7 54.9 54.2 53.3 ...
## $ SPDY: num 6.8 6.68 6.57 6.48 6.4 ...
## $ FEIN: num 68.7 69 69.3 69.6 69.8 ...
## $ MAIN: num 65 65.3 65.7 65.9 66.2 ...
## $ SPPO: num 79.1 77.7 76.2 74.7 73.2 ...
## - attr(*, ".internal.selfref")=<externalptr>
d[,Year:=as.numeric(Year)]
ddum <- dummyVars("~.", data = d)
d <- data.table(predict(ddum, newdata = d))
rm(ddum) #remove ddum as unneeded
str(d)
## Classes 'data.table' and 'data.frame': 162 obs. of 21 variables:
## $ CntNArbWr: num 1 1 1 1 1 1 1 1 1 1 ...
## $ CntNCEatB: num 0 0 0 0 0 0 0 0 0 0 ...
## $ CntNEr&CA: num 0 0 0 0 0 0 0 0 0 0 ...
## $ CntNErpnU: num 0 0 0 0 0 0 0 0 0 0 ...
## $ CntNEsA&P: num 0 0 0 0 0 0 0 0 0 0 ...
## $ CntNLtA&C: num 0 0 0 0 0 0 0 0 0 0 ...
## $ CntNME&NA: num 0 0 0 0 0 0 0 0 0 0 ...
## $ CntNNrthA: num 0 0 0 0 0 0 0 0 0 0 ...
## $ CntNSb-SA: num 0 0 0 0 0 0 0 0 0 0 ...
## $ Year : num 1997 1998 1999 2000 2001 ...
## $ NYGD : num 2299 2170 2314 2589 2495 ...
## $ NYGN : num 2310 2311 2288 2410 2496 ...
## $ SEPR : num 6078141 5961001 5684714 5425963 5087547 ...
## $ ERMA : num 4181176 4222039 4131775 3955257 3726838 ...
## $ SESC : num 8.08 8.27 8.5 8.65 8.84 ...
## $ FEMA : num 9.73 9.82 9.97 10.02 10.12 ...
## $ SPAD : num 56.6 55.7 54.9 54.2 53.3 ...
## $ SPDY : num 6.8 6.68 6.57 6.48 6.4 ...
## $ FEIN : num 68.7 69 69.3 69.6 69.8 ...
## $ MAIN : num 65 65.3 65.7 65.9 66.2 ...
## $ SPPO : num 79.1 77.7 76.2 74.7 73.2 ...
## - attr(*, ".internal.selfref")=<externalptr>
缩放和居中
准备工作的第二阶段是根据上一章提到的原因对我们的数据进行缩放和集中。简而言之，这些原因是为了防止可能具有较长射程的某些组件意外超重。我们不缩放前九列，它们是我们的虚拟编码列。我们将使用cbind重新添加这些列。
dScaled<-scale(d[,-c(1:9)])
dScaled<-as.data.table(dScaled)
d <- cbind(d[,c(1:9)], dScaled)
rm(dScaled) #remove d2 as unneeded
str(d)
## Classes 'data.table' and 'data.frame': 162 obs. of 21 variables:
## $ CntNArbWr: num 1 1 1 1 1 1 1 1 1 1 ...
## $ CntNCEatB: num 0 0 0 0 0 0 0 0 0 0 ...
## $ CntNEr&CA: num 0 0 0 0 0 0 0 0 0 0 ...
## $ CntNErpnU: num 0 0 0 0 0 0 0 0 0 0 ...
## $ CntNEsA&P: num 0 0 0 0 0 0 0 0 0 0 ...
## $ CntNLtA&C: num 0 0 0 0 0 0 0 0 0 0 ...
## $ CntNME&NA: num 0 0 0 0 0 0 0 0 0 0 ...
## $ CntNNrthA: num 0 0 0 0 0 0 0 0 0 0 ...
## $ CntNSb-SA: num 0 0 0 0 0 0 0 0 0 0 ...
## $ Year : num -1.633 -1.441 -1.249 -1.057 -0.865 ...
## $ NYGD : num -0.813 -0.822 -0.812 -0.792 -0.799 ...
## $ NYGN : num -0.793 -0.793 -0.795 -0.786 -0.78 ...
## $ SEPR : num 0.327 0.309 0.266 0.225 0.172 ...
## $ ERMA : num 0.1565 0.1643 0.1471 0.1134 0.0699 ...
## $ SESC : num -1.52 -1.46 -1.39 -1.34 -1.28 ...
## $ FEMA : num -1.35 -1.31 -1.25 -1.22 -1.18 ...
## $ SPAD : num 0.37 0.344 0.319 0.295 0.269 ...
## $ SPDY : num -0.66 -0.708 -0.749 -0.784 -0.815 ...
## $ FEIN : num -0.744 -0.704 -0.668 -0.635 -0.603 ...
## $ MAIN : num -0.622 -0.573 -0.527 -0.484 -0.444 ...
## $ SPPO : num 1.55 1.45 1.34 1.24 1.13 ...
## - attr(*, ".internal.selfref")=<externalptr>
转换
这个缩放过程不会删除异常值。我们应该检查数据，看看是否存在显著的非正常行为，因为一些机器学习技术在正常数据上表现得更好。在这些情况下，执行其他变换(如对数缩放)可能会有所帮助。正如我们在图 8-1 中看到的，数据中存在一些异常值。
图 8-1
寻找重要的畸形数据
boxplot(d[,-c(1:9)], las = 2)
在图 8-2 中更仔细地观察小学毕业的男性儿童，可以发现非正常行为的证据。直方图不是正态分布的传统钟形，Q-Q 图不像直线y = x。
par(mfrow = c(1,2))
hist(d$ERMA, 100)
qqnorm(d$ERMA)
par(mfrow = c(1,1))
事实上，如果我们使用夏皮罗-维尔克正态性检验进行假设检验，我们看到 p 值是显著的，因此我们拒绝正态数据的零假设。
图 8-2
寻找重要的畸形数据
shapiro.test(d$ERMA)
##
## Shapiro-Wilk normality test
##
## data: d$ERMA
## W = 0.6, p-value <2e-16
基于图 8-2 中的直方图和 Q-Q，转换数据可能是有意义的，也许使用对数转换。首先，我们观察到缩放和居中数据的范围自然包括零，这对于对数变换来说是一个挑战。原始的data显示这不是底层信息的问题。然而，夏皮罗-维尔克正态性检验表明，即使是对数变换仍然不利于符合正态分布。
range(d$ERMA)
## [1] -0.62 3.32
range(data$ERMA)
## [1] 1.1e+05 2.1e+07
shapiro.test( log(data$ERMA) )
##
## Shapiro-Wilk normality test
##
## data: log(data$ERMA)
## W = 1, p-value = 3e-04
现在，从视觉上看，我们确实看到对数转换确实在一定程度上分隔了我们的数据，如图 8-3 所示。
图 8-3
寻找重要的畸形数据
par(mfrow = c(1,2))
hist(data$ERMA, 100)
hist( log (data$ERMA) , 100)
par(mfrow = c(1,1))
这个间距足以保证一个转换吗？一种确定方法是比较原始数据和转换数据与所需响应变量之间的相关性。我们的反应/输出变量将是SPAD，青春期生育率。我们创建一个数据集的副本，其中包括日志格式和原始格式的ERMA,并观察相关性。由于第一列中产生的相关性显示在未转换的列中相关性更强，所以我们最终选择不执行转换。
d2 <- copy(data[,.(SPAD, ERMA)])
d2[, Log.ERMA := log(ERMA)]
cor(d2)
## SPAD ERMA Log.ERMA
## SPAD 1.00 0.83 0.72
## ERMA 0.83 1.00 0.81
## Log.ERMA 0.72 0.81 1.00
rm(d2) #no longer needed
从技术上讲，应该检查每个考虑中的变量的正态性(同样，这只是因为我们的一些模型在使用正态预测时表现更好)。同样，一般来说，在应用特定模型时，理解最佳模型性能的预期输入是很重要的。首先，使用lapply，我们观察到所有的 p 值都是显著的。接下来，我们使用sapply对数据进行对数转换，并使用colnames添加后缀。最后，我们将原始数据相关性与SPAD和对数转换相关性进行比较和对比。在某些情况下，相关性稍好一些。在其他国家，情况并非如此。
lapply(data[,-c(1:2)], shapiro.test)
## $NYGD
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.8, p-value = 3e-13
##
##
## $NYGN
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.8, p-value = 1e-13
##
##
## $SEPR
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.6, p-value <2e-16
##
##
## $ERMA
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.6, p-value <2e-16
##
##
## $SESC
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.9, p-value = 7e-06
##
##
## $FEMA
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 1, p-value = 5e-05
##
##
## $SPAD
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.8, p-value = 5e-13
##
##
## $SPDY
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.9, p-value = 6e-08
##
##
## $FEIN
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.8, p-value = 1e-13
##
##
## $MAIN
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.8, p-value = 9e-14
##
##
## $SPPO
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.8, p-value = 6e-13
dlog <- copy(data)
dlog <- sapply(dlog[,-c(1:2)], log)
dlog<-as.data.table(dlog)
colnames(dlog) <- paste(colnames(dlog), "LOG", sep = ".")
dlog<-cbind(data, dlog)
View(cor(dlog[,-c(1:2)]))
rm(dlog) #remove as we will not use.
现在，我们在所有情况下都不进行对数转换。根据分析结果，根据所需的模型精度，甚至根据所选的具体模型，我们保留记录哪些变量在转换后具有更强相关性的权利。通常，各种级别的预处理成为模型调整过程的一部分。
培训与验证数据
我们到达预处理阶段，讨论训练和测试数据。回想一下引言章节，训练数据集的选择相当重要。对于任何单个模型，选择一个训练集和一个验证集(通常是 80/20 分割)是有意义的。这允许在训练集上创建模型，然后在验证集上获得一些“真实”世界的准确性(因此模型没有看到数据)。然而，在本章中，我们打算讨论几种可能的模型。如果我们选择这些模型中的一个，我们又一次(因为我们使用验证数据来选择“最佳”模型)冒了某种程度的过度拟合的风险。因此，在我们的模型用于现实生活之前，我们需要从训练和验证中保留最后一点测试数据，以用作我们最终的“真实”世界估计。这变得很复杂，因为 60/20/20 的分割现在只留给我们一半多一点的数据用于训练。
当然，交叉验证是避免留给训练的数据太少的一种方法。我们可以执行一个 80/20 的训练/测试分割，对我们的训练数据使用交叉验证，并节省我们自己对正式验证集的需求。这在计算上非常昂贵，尤其是在多个模型上(或者甚至是一个模型的多次迭代)。回想一下，最后，一旦选择了一个特定的模型，在现实世界中使用之前，正确的最后一步确实是对所有数据重新训练该模型——不再进行拆分。根据模型本身的计算强度，这可能要付出相当大的代价。
由于我们的数据集只有 162 个观察值，为了简单起见，我们将只进行标准的 80/20 分割。我们的目标是在我们的模型中使用交叉验证，这样保留的测试集可以在最后使用。
使用set.seed(1234)来允许再现性，我们引入一个新的函数来执行从caret包createDataPartition的分割。第一种形式允许函数了解我们的数据是如何分层的，第二种形式告诉我们分割的比率(在这个例子中是 80/20)。一旦index被填充，我们使用它来分离我们的训练和验证或测试数据。拆分完数据后，我们就为最后阶段的预处理做好了准备。
set.seed(1234)
index <- createDataPartition(data$CntN, p = 0.8, list = FALSE)
trainData <- data[index, ]
validationData <- data[-index, ]
我们应该注意，从技术上讲，对于训练数据、测试数据和验证数据，缩放和居中应该分开进行。事实上，音阶和中心公式必须在火车上开发；否则，这些集合之间可能会有信息泄漏。然而，由于我们一直在重新设置我们的数据，并且我们打算在即将到来的部分再次缩放，我们希望我们的读者原谅我们的订购。
主成分分析
我们在上一章看到了主成分分析。降维的优点包括减少处理时间、减少过拟合和模型简化。如果您已经在无监督学习一章中安装了pcaMethods包，则不需要再次安装。如果您需要安装，请删除注释掉代码的井号，然后运行。无论哪种情况，第三行代码肯定都需要运行。
#source("https://bioconductor.org/biocLite.R")
#biocLite("pcaMethods")
library(pcaMethods)
## Loading required package: Biobase
## Loading required package: BiocGenerics
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:Matrix':
##
## colMeans, colSums, rowMeans, rowSums, which
## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind,
## colMeans, colnames, colSums, dirname, do.call, duplicated,
## eval, evalq, Filter, Find, get, grep, grepl, intersect,
## is.unsorted, lapply, lengths, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position,
## rank, rbind, Reduce, rowMeans, rownames, rowSums, sapply,
## setdiff, sort, table, tapply, union, unique, unsplit,
## which, which.max, which.min
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages
## 'citation("pkgname")'.
##
## Attaching package: 'pcaMethods'
## The following object is masked from 'package:stats':
##
## loadings
已经分割了我们的数据，我们为 PCA 准备数据。在这里，我们的选择一部分是由我们的数据决定的，一部分是由我们的最终模型决定的。例如，各种类型的森林和基于树的方法可能不存在字符或因子数据问题。另一方面，各种类型的回归将需要前面提到的虚拟编码。PCA 不是特别适合伪编码的数据。然而，如果大量数据碰巧是分类的，那么首先对所有分类预测器进行虚拟编码，然后尝试各种类型的降维，这可能是有意义的，其中 PCA(又名svd)只是一种方法。
另一方面，我们的数据只有一位完全分类的数据—CntN。如果我们的数据集有几个伪编码列，然后我们将其他 10 个变量减少到 3 或 4 个，我们仍然可以获得可观的降维。
一般来说，PCA 是应用于数据的最后预处理步骤。当然，在 PCA 之前，应该对数值数据进行缩放和居中。至于分类数据是否也应该经过主成分分析？为简洁起见，我们将展示如何手动排除分类数据。如果您的数据需要不同的或更细微的技术，我们希望我们展示了足够的过程，以便您可以根据自己的目的修改代码。
首先，记住我们的数据的已知结构，我们用str函数验证我们已经选择了我们的训练集的正确行。
默认情况下，stats库带有基类R，prcomp是该库中的一个函数。它对我们的训练数据执行维度计算。请记住，任何符合我们 PCA 训练数据的模型都需要提供新数据或验证数据，这些数据执行了相同的 PCA 过程。为此，在通过summary(pc)查看 PCA 的结果后，我们使用predict功能。第一种形式是我们的 PCA 分析，在本例中存储在pc中。第二个是新数据，在本例中是我们的validationData。当然，对于实际生产中使用的模型，新数据将是我们希望放入模型的新数据。
然而，从上一章我们知道，除了“基本的”PCA 之外，还可以做更多的工作。因此，我们还展示了使用pcaMethods库中的pca函数所需的步骤。该功能的细节在无监督机器学习一章中有解释。现在我们注意到这个过程的结果与predict完美地工作在一起，我们展示了两种方法的第一维的打印输出，以表明它们确实是相同的。
当然，pca的优势在于通过method = " "形式的便利，方法转换成listPcaMethods()列出的任何方法都很容易。
#confirm structure
str(trainData[,c(3:8,10:13)])
## Classes 'data.table' and 'data.frame': 135 obs. of 10 variables:
## $ NYGD: num 2299 2170 2314 2589 2495 ...
## $ NYGN: num 2310 2311 2288 2410 2496 ...
## $ SEPR: num 6078141 5961001 5684714 5425963 5087547 ...
## $ ERMA: num 4181176 4222039 4131775 3955257 3726838 ...
## $ SESC: num 8.08 8.27 8.5 8.65 8.84 ...
## $ FEMA: num 9.73 9.82 9.97 10.02 10.12 ...
## $ SPDY: num 6.8 6.68 6.57 6.48 6.4 ...
## $ FEIN: num 68.7 69 69.3 69.6 69.8 ...
## $ MAIN: num 65 65.3 65.7 65.9 66.2 ...
## $ SPPO: num 79.1 77.7 76.2 74.7 73.2 ...
## - attr(*, ".internal.selfref")=<externalptr>
#base R / traditional method
pc <- prcomp(trainData[,c(3:8,10:13)], center = TRUE, scale. = TRUE)
summary(pc)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 2.703 1.321 0.8175 0.4029 0.27986 0.16986
## Proportion of Variance 0.731 0.175 0.0668 0.0162 0.00783 0.00289
## Cumulative Proportion 0.731 0.905 0.9719 0.9881 0.99597 0.99886
## PC7 PC8 PC9 PC10
## Standard deviation 0.08417 0.04726 0.03535 0.02894
## Proportion of Variance 0.00071 0.00022 0.00012 0.00008
## Cumulative Proportion 0.99957 0.99979 0.99992 1.00000
pcValidationData1 <- predict(pc, newdata = validationData[,c(3:8,10:13)])
#scalable method using PcaMethods
pc<-pca(trainData[,c(1:8,10:13)], method = "svd",nPcs = 4, scale = "uv", center = TRUE)
pc
## svd calculated PCA
## Importance of component(s):
## PC1 PC2 PC3 PC4
## R2 0.7306 0.1745 0.06683 0.01623
## Cumulative R2 0.7306 0.9051 0.97191 0.98814
## 10 Variables
## 135 Samples
## 0 NAs ( 0 %)
## 4 Calculated component(s)
## Data was mean centered before running PCA
## Data was scaled before running PCA
## Scores structure:
## [1] 135 4
## Loadings structure:
## [1] 10 4
summary(pc)
## svd calculated PCA
## Importance of component(s):
## PC1 PC2 PC3 PC4
## R2 0.7306 0.1745 0.06683 0.01623
## Cumulative R2 0.7306 0.9051 0.97191 0.98814
pcValidationData2 <- predict(pc, newdata = validationData[,c(3:8,10:13)])
#demonstration of how to access transformed validation data
pcValidationData1[,1]
## [1] -1.11 -0.87 -0.77 0.96 1.08 1.67 -1.52 -1.17 -0.82 0.37 1.06
## [12] 1.67 1.94 2.67 3.36 0.76 0.86 1.00 -1.14 -0.49 0.40 2.66
## [23] 3.72 3.80 -7.52 -5.80 -4.95
pcValidationData2$scores[,1]
## [1] -1.11 -0.87 -0.77 0.96 1.08 1.67 -1.52 -1.17 -0.82 0.37 1.06
## [12] 1.67 1.94 2.67 3.36 0.76 0.86 1.00 -1.14 -0.49 0.40 2.66
## [23] 3.72 3.80 -7.52 -5.80 -4.95
转换完我们的数据后，现在使用cbind将一个热编码的CntN和我们的响应变量SPAD放入就很简单了。
在我们结束预处理数据这一节时，记住几个关键原则很重要。
如果使用单一模型，那么训练和验证分离(通常是 80/20)是合适的，除非使用交叉验证或引导。请记住，后两者的计算量可能很大。如果从一个以上的模型中进行选择，则需要一个训练集来训练各种模型，需要一个验证集来估计模型在未知数据上的性能，最后需要一个测试集来估计最终选择的模型在真实数据上的性能(否则，我们实际上在验证步骤中有第二阶段过度拟合)。交叉验证可以消除对验证集的需求，但在使用多个模型时，计算成本会更高。另一方面，使用 70/20/10 进行训练、验证和测试的成本可能很高。无论选择哪种方法，都应该首先进行，以防止不同组之间的信息泄漏。
如果您的数据有丢失的值，首先阅读专门讨论丢失数据的章节可能会有所收获。插补通常是数据布局正确后的第一步。
分类数据的热编码或虚拟编码是一种选择。一些模型只需要数字数据，其他模型可以用因子格式处理分类数据。正如每种类型的模型可能并不适合所有类型的数据一样，某些类型的数据也不适合某些模型。因此，要同时考虑数据和模型！对于某些模型来说，保持数据的分类和因子格式，可以为模型提供有效使用数据所需的信息，同时不会对数据产生任意的线性影响。相比之下，其他模型可能无法处理这些因素，因此除了虚拟代码之外可能别无选择。
一般来说，数字数据应该缩放和居中。事实上，删除方差为零或接近零的列也可能是明智的。具有接近零的方差的列可以被保留，并且稍后在预处理中经受 PCA。
除了基本的居中和缩放之外，它还有助于转换数据以获得适当的数据分布。虽然我们只讨论了对数变换，但是对于特定类型和特征的数据以及特定的模型来说，有许多可能的变换是有价值的。同样，如果一个特定模型的第一遍不够精确，更高级的数学可能是合适的。
8.2 监督学习模式
现在我们了解了一些预处理数据的方法，我们的下一个任务是满足一些模型。记住，监督学习仅仅意味着我们头脑中有一个特定的反应变量。对于各种类型的响应或预测变量，不同的模型可能更适合于更好地分析特定类型的数据。例如，如果我们有更多的分类变量，那么我们可能更喜欢基于分类的监督学习。另一方面，如果我们的数据主要是数字，我们可能会在回归方法中发现更多的价值。一些模型生活在两个世界中，或者可以适应这样做。而且，正如我们在预处理中已经看到的，数据本身可以变异，例如，通过虚拟编码从分类到数字。当然，这是一个很差的方法，只能以一种方式工作——直方图条可能是一个很好的例子，可以将连续的数字数据重新编码为分类数据。
在研究这些模型时，需要记住的是，具体应用需要什么样的精度。找到一个提供比机会更好的结果的模型是合理可行的。仅仅比机会好就足够了吗？有时间(人和机器)测试更多的模型吗？大多数模型都有调整参数、允许更多交叉验证的额外输入或惩罚额外的模型复杂性。花多少精力调优才算现实？
在下面的章节中，我们采取了一种平衡的方法。我们准备数据，使其与特定模型最佳交互，我们假设硬件水平合理，并选择使用合理的计算时间，我们讨论各种调整参数选项，但不必穷尽这些选项。对于任何特定的数据应用，在数据收集、预处理和最终模型调整上花费比在模型选择步骤上更多的时间是有意义的。有了这个想法，然后安全地储存在我们的脑海中，我们继续前进。
支持向量机
支持向量机(SVM)在概念上有点类似于 k 近邻，因为其思想是找到我们的数据分组。这是通过寻找通过识别关键数据点(支持向量)确定的边界曲线来实现的。支持向量是那些具有一些最大间隔(在一些度量下)的向量，并且在这些向量之间绘制边界曲线。由于边界划分可以是线性的、多项式的或更奇异的，因此有几种方法来执行基本的 SVM 算法。界定曲线类型的确定被称为 SVM 的核心。自然地，更复杂的边界核往往在计算上更困难。与其他一些包一样，caret包本身支持多种内核风格的几种变体，包括线性、指数、多项式和径向。我们将把这里的讨论限于线性和多项式。
让我们从一个简单的数据集开始——不是我们完整的正式数据集——来直观地理解支持向量机是如何工作的。我们使用青少年生育率和人均国民生产总值作为预测指标，我们的分类回答将是国名。我们重新利用我们熟知的图表 8-4 并注意到，对于任何一个国家分组，通常都有可能存在明显的线性分离。诚然，我们并不期待左下角出现奇迹。
图 8-4
人均国民生产总值与青少年生育率
svmDataTrain <- trainData[,.(SPAD, NYGN)]
svmDataValidate <- validationData[,.(SPAD, NYGN)]
p1 <- ggplot(data = svmDataTrain,
aes(x = NYGN, y = SPAD))
## data poins colored by country
p1 + geom_point(aes(colour = trainData$CntN)) +
scale_colour_viridis(discrete = TRUE)
caret包由一个名为train的函数驱动。花点时间看看下面的代码布局。
第一种形式用于预测，第二种形式用于响应。在撰写本文时，第三种形式有超过 200 种可能的方法。在这种情况下，我们使用的是线性支持向量机，因此svmLinear。因为我们打算在原始数据上绘制我们的结果，所以现在我们选择保持preProcess为空。将来，它将采用center, scale, and pca等数值以及各种插补选项。正式参数metric = "Accuracy"表示我们希望我们的汇总输出让我们知道CountryName的训练区域与实际国家相比有多准确。同样，可以将 svm 看作是使用线性边界将数据分成不同的区域。我们想知道这条分界线多长时间划好一次。最后，trainControl功能用于控制模型的训练方式。在这种情况下，我们将使用交叉验证cv和five折叠。实际上有大量的数据包含在svm变量中。显示的打印输出只是一个小摘要。它确实提醒我们没有预处理，并且按照要求，告诉我们交叉验证的估计准确性。请记住，这种准确性是基于每个折叠的交叉验证过程。
set.seed(12345)
svm <- train(x = svmDataTrain,
y = trainData$CntN,
method = "svmLinear",
preProcess = NULL,
metric = "Accuracy",
trControl = trainControl(method = "cv",
number = 5,
seeds = c(123, 234, 345, 456, 567, 678)
)
)
svm
## Support Vector Machines with Linear Kernel
##
## 135 samples
## 2 predictor
## 9 classes: 'ArbWr', 'CEatB', 'Er&CA', 'ErpnU', 'EsA&P', 'LtA&C', 'ME&NA', 'NrthA', 'Sb-SA'
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 108, 108, 108, 108, 108
## Resampling results:
##
## Accuracy Kappa
## 0.82 0.8
##
## Tuning parameter 'C' was held constant at a value of 1
接下来，我们稍微深入研究一下模型的准确性。因为我们实际上使用了交叉验证，这被认为是一个很好的估计。首先，我们使用有两个输入的predict函数。一个是包含我们模型的新svm。另一个是我们用来训练那个模型的训练数据。与一次只使用108元素的交叉验证折叠不同，我们使用所有的135值。接下来，我们发现这些预测的国家名称何时与在我们的trainData变量中找到的正确国家名称相同。这创建了布尔值TRUE和FALSE，它们分别等于1和0,。因此，这些事件中的mean将会给我们true发生频率的比率。请注意，尽管对0.84的不可见数据的估计是准确的，但在我们的训练数据上，我们实际上看到了0.85。这是我们所说的过度拟合的一部分。
#predict the country name on our training data using our new model
predictOnTrain <- predict(svm, newdata = svmDataTrain)
mean( predictOnTrain == trainData$CntN)
## [1] 0.82
现在，使用我们提供的验证数据是不正确的。从技术上来说，我们必须把它保存到最后，而不是用它来通知我们的模型选择，然后在我们看到所有的模型之后才使用它来执行最后的验证步骤。然而，鉴于在第一个示例草图中展示了这个过程的样子，展示最后一步有一定的教学意义。我们执行与训练数据相同的计算(实际上，这个预测过程就是我们对真正新的真实世界数据所做的)。虽然交叉验证比实际训练数据可能显示的更保守，但它并不像它本来可能的那样保守。
同样，我们模型中完全看不见的数据上的0.78是模型过度拟合的一个例子。训练数据在某种程度上使我们的模型产生了偏差，在完全看不见的验证数据的情况下，它表现得不那么好。
predictOnTest <- predict(svm, newdata = svmDataValidate)
mean(predictOnTest == validationData$CntN)
## [1] 0.81
使用与上一章相同的过程，我们绘制了看不见的验证数据，并将正确答案的结果与图 8-5 中的预测进行比较。我们的模型很难理解线性分离的确切位置，也正是不准确性显现的地方。
图 8-5
训练与测试预测
p1 <- ggplot(data = validationData,
aes(x = NYGN, y = SPAD))
plot_grid(
## data poins colored by country
p1 + geom_point(aes(colour = validationData$CntN, size = validationData$CntN)) +
scale_colour_viridis(discrete = TRUE),
## data poins colored by predicted country
p1 + geom_point(aes(colour = predictOnTest, size = predictOnTest)) +
scale_colour_viridis(discrete = TRUE),
ncol = 1
)
## Warning: Using size for a discrete variable is not advised.
## Warning: Using size for a discrete variable is not advised.
预计前面的代码会生成两个警告—也就是说，对于区域名称这样的离散变量，不建议使用这样的大小。这是一个合理的警告，在现实生活中，我们永远不会使用这样的技术(一个重要的原因是人类倾向于将大小本身解释为重要的，在这种情况下没有这样的意义)。我们在这里使用它只是为了帮助区分类别。
既然我们已经看到并理解了caret模型构建的整体布局，以及如何在最后使用验证数据来理解选择哪个模型，我们将验证数据安全地保存起来，直到最后才会再次使用，以了解我们最终选择的模型可能会如何执行。我们还花了一些时间来清理我们的环境。
rm(p1)
rm(svm)
rm(svmDataTrain)
rm(svmDataValidate)
rm(pcValidationData1)
rm(pcValidationData2)
rm(predictOnTest)
rm(predictOnTrain)
rm(pc)
rm(d)
我们现在准备在我们的整个trainData布景上表演我们的 SVM 模型。我们留下的是我们原始的data的一个副本，我们用来将数据分成训练集和验证集的index，以及实际的训练集和验证集本身。
我们再次构建我们的模型训练数据减去我们的分类国家名称列。现在，我们将字符年份数据重新编码为数字，并为可再现性设置种子。除了包含更多的可变列之外，唯一的不同是设置了preProcess = c("scale", "center", "pca")。这默认为 PCA thresh = 0.95，如果需要，可以按照?preProcess中的说明进行调整。在这种情况下，我们使用默认值。我们看到估算模型的准确性有所提高。
# set up training & validation data
svmDataTrain <- trainData[,-1]
svmDataTrain[,Year:=as.numeric(Year)]
svmDataValidation <- validationData[,-1]
svmDataValidation[,Year:=as.numeric(Year)]
#run linear SVM on the full data set
set.seed(12345)
svmLinear <- train(x = svmDataTrain,
y = trainData$CntN,
method = "svmLinear",
preProcess = c("scale", "center", "pca"),
metric = "Accuracy",
trControl = trainControl(method = "cv",
number = 5,
seeds = c(123, 234, 345, 456, 567, 678)
)
)
svmLinear
## Support Vector Machines with Linear Kernel
##
## 135 samples
## 12 predictor
## 9 classes: 'ArbWr', 'CEatB', 'Er&CA', 'ErpnU', 'EsA&P', 'LtA&C', 'ME&NA', 'NrthA', 'Sb-SA'
##
## Pre-processing: scaled (12), centered (12), principal
## component signal extraction (12)
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 108, 108, 108, 108, 108
## Resampling results:
##
## Accuracy Kappa
## 0.99 0.98
##
## Tuning parameter 'C' was held constant at a value of 1
在这一点上，描述线性 SVM 是如何工作的是很重要的。该算法是为二元分类而设计的。在这种情况下，对于多个结果，它以一个对其余的方式在数据中循环。所以这个线性模型实际上是画了一条直线，一次将一个群体与其他群体隔离开来。目标是在数据之间画一条线，使这条线尽可能远离分界线两侧的点。这有助于一个模型预测未来数据应该走向何方。通过调整参数C，我们控制需要多宽的余量(如果数据过于接近，可能会使生产线选择不同的路径)。更高的C值导致更高的训练精度和更小的裕度——存在过度拟合的风险。
当我们从线性变为多项式时，它允许直线变成多项式曲线。此外，默认情况下，该算法循环遍历几个C值。注意(在运行下面的代码之后),一般来说，对于每组固定的度数和比例，C的值越大，产生的精确度越高。代价是，随着C的增加，我们曲线的边缘缩小，这意味着我们的模型过度拟合数据的可能性增加。
接下来，我们只对我们的模型做一个改变，从线性模型到多项式模型。在caret中，这很容易通过切换到method = "svmPoly"来完成。虽然代码的变化很小，但计算工作量的变化更小。注意svmPoly是一个完全不同于线性版本的模型。虽然算法和数学超出了本文的范围，但是即使在线性程度级别，过程也是不同的。虽然这目前看起来不如我们的线性模型好，但如果是真的，那么0.88的准确性仍然可能相当好。
#run polynomial SVM on the full data set
set.seed(12345)
svmPoly <- train(x = svmDataTrain,
y = trainData$CntN,
method = "svmPoly",
preProcess = c("scale", "center", "pca"),
metric = "Accuracy",
trControl = trainControl(method = "cv",
number = 5
)
)
svmPoly
## Support Vector Machines with Polynomial Kernel
##
## 135 samples
## 12 predictor
## 9 classes: 'ArbWr', 'CEatB', 'Er&CA', 'ErpnU', 'EsA&P', 'LtA&C', 'ME&NA', 'NrthA', 'Sb-SA'
##
## Pre-processing: scaled (12), centered (12), principal
## component signal extraction (12)
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 108, 108, 108, 108, 108
## Resampling results across tuning parameters:
##
## degree scale C Accuracy Kappa
## 1 0.001 0.25 0.76 0.73
## 1 0.001 0.50 0.76 0.73
## 1 0.001 1.00 0.76 0.73
## 1 0.010 0.25 0.76 0.73
## 1 0.010 0.50 0.76 0.73
## 1 0.010 1.00 0.76 0.73
## 1 0.100 0.25 0.76 0.73
## 1 0.100 0.50 0.78 0.75
## 1 0.100 1.00 0.81 0.79
## 2 0.001 0.25 0.76 0.73
## 2 0.001 0.50 0.76 0.73
## 2 0.001 1.00 0.76 0.73
## 2 0.010 0.25 0.76 0.72
## 2 0.010 0.50 0.76 0.72
## 2 0.010 1.00 0.76 0.72
## 2 0.100 0.25 0.79 0.76
## 2 0.100 0.50 0.81 0.79
## 2 0.100 1.00 0.87 0.85
## 3 0.001 0.25 0.76 0.73
## 3 0.001 0.50 0.76 0.73
## 3 0.001 1.00 0.76 0.73
## 3 0.010 0.25 0.76 0.72
## 3 0.010 0.50 0.76 0.72
## 3 0.010 1.00 0.76 0.72
## 3 0.100 0.25 0.80 0.78
## 3 0.100 0.50 0.86 0.84
## 3 0.100 1.00 0.90 0.88
##
## Accuracy was used to select the optimal model using the
## largest value.
## The final values used for the model were degree = 3, scale = 0.1
## and C = 1.
请注意，这两个模型完全符合我们的完整训练集。因此，尽管多项式预测的估计精度较低，但这两个模型似乎都非常符合我们的数据。
predictOnTrainL <- predict(svmLinear, newdata = svmDataTrain)
mean( predictOnTrainL == trainData$CntN)
## [1] 1
predictOnTrainP <- predict(svmPoly, newdata = svmDataTrain)
mean( predictOnTrainP == trainData$CntN)
## [1] 0.98
基于精度水平，如果在现实生活中面对这两个选项，我们会选择线性选项。在这个阶段，我们已经做出了选择，然后我们将只在我们选择的模型上运行最后一次验证数据。我们在下面的代码中做到了这一点，并取得了很好的结果。因此，我们变得相当确定，基于我们选择的数据，给定 2015 年对我们的其余数据点的新观察，我们将对我们正确排序到正确区域的能力相当有信心。当然，我们最好在未来跟踪我们模型的准确性。
predictOnTestL <- predict(svmLinear, newdata = svmDataValidation)
mean(predictOnTestL == validationData$CntN)
## [1] 1
总的来说，我们的数据集似乎相当完整。换句话说，这些不同的地理区域在最近几年被收集的一组预测因子很好地描述了。考虑到图 8-4 中非常清晰的分组，这并不特别令人震惊。事实上，意识到 PCA 不会对这两个变量有特别大的改变。有很好的垂直和水平分隔。
这就结束了关于支持向量机的部分。在接下来的章节中，我们将继续使用caret的高度一致的结构来利用其他模型。随着我们对整体模型结构越来越熟悉，我们试图增加额外的技能。请记住，在现实生活中，后面介绍的技巧很可能会用到支持向量机。
分类和回归树
分类和回归树(CART)是为数字、连续预测值和分类响应变量而设计的。这非常符合我们的数据集。我们刷新 CART 的训练和验证数据，使用set.seed来获得再现性，并使用我们现在熟悉的train函数。我们添加了一个新的特性，tuneLength = 10，它控制复杂度的迭代次数。增加该值会增加模型的总计算时间；因此，对于较大的数据集，模型精度的提高需要与训练模型的时间进行权衡。
cartDataTrain <- copy(trainData[,-1])
cartDataTrain[,Year:=as.numeric(Year)]
cartDataValidation <- copy(validationData[,-1])
cartDataValidation[,Year:=as.numeric(Year)]
set.seed(12345)
cartModel <- train(x = cartDataTrain,
y = trainData$CntN,
method = "rpart",
preProcess = c("scale", "center", "pca"),
metric = "Accuracy",
tuneLength = 10,
trControl = trainControl(method = "cv",
number = 5
)
)
cartModel
## CART
##
## 135 samples
## 12 predictor
## 9 classes: 'ArbWr', 'CEatB', 'Er&CA', 'ErpnU', 'EsA&P', 'LtA&C', 'ME&NA', 'NrthA', 'Sb-SA'
##
## Pre-processing: scaled (12), centered (12), principal
## component signal extraction (12)
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 108, 108, 108, 108, 108
## Resampling results across tuning parameters:
##
## cp Accuracy Kappa
## 0.000 0.84 0.83
## 0.014 0.84 0.83
## 0.028 0.84 0.83
## 0.042 0.85 0.83
## 0.056 0.85 0.83
## 0.069 0.85 0.83
## 0.083 0.85 0.83
## 0.097 0.76 0.73
## 0.111 0.73 0.70
## 0.125 0.11 0.00
##
## Accuracy was used to select the optimal model using the
## largest value.
## The final value used for the model was cp = 0.083.
在精确度降低之前，模型最终根据精确度(如我们所要求的)选择最大的复杂性值(这是交叉验证的一个例子，帮助我们避免过度拟合模型)。
可使用通用绘图功能绘制finalModel，如图 8-6 所示。这棵树没有太多的层次。请注意，虽然我们的估计准确性可能相当高(因此我们可能希望我们的模型具有可靠的预测值)，但该模型对于理解为什么应该做出特定预测并不那么有用。因此，虽然通常树模型可能被认为是相当容易解释的，但使用主成分分析会降低我们理解为什么我们的模型可能是正确的能力。
图 8-6
分类树形图
plot(cartModel$finalModel)
text(cartModel$finalModel, cex = 0.5)
虽然通用的plot函数很有帮助，但是fancyRpartPlot通常看起来更整洁，尽管在这种情况下，很难找到合适的文本缩放级别，如图 8-7 所示。在任何情况下，重要的是要注意，虽然模型的准确性可以通过主成分分析来提高，但模型的可解释性变得更加复杂。开发一种方法来更好地解释将是有帮助的。
图 8-7
花式分类树形图
fancyRpartPlot(cartModel$finalModel, cex = 0.4, main = "")
接下来，我们重复我们用于 SVM 的准确性度量。我们的最终模型在训练数据上是高度精确的，并且确实在完全看不见的数据上保持了这种精确度(嗯，几乎保持了)。同样，如果我们希望对模型在真实世界数据中的表现有一个真实的估计，当我们继续使用交叉验证时，正确的方法是等到本章的最后，通过交叉验证选择最准确的模型，然后最后只对选择的模型运行predictOnTestT。然而，作为一种教学方法，在完全看不见的数据上看到结果是有价值的，所以我们请求我们的读者继续原谅我们。
predictOnTrainT <- predict(cartModel, newdata = cartDataTrain)
mean( predictOnTrainT == trainData$CntN)
## [1] 0.77
predictOnTestT <- predict(cartModel, newdata = cartDataValidation)
mean(predictOnTestT == validationData$CntN)
## [1] 0.67
caret包有一个我们还没有介绍的功能，叫做confusionMatrix。虽然这种情况下的输出很长，但是我们可以看到预测与验证参考的结果。一点滚动和良好的记忆显示，有一个阿拉伯世界的预测是不正确的，实际的数据值是东亚和太平洋。混淆矩阵会有所帮助，因为它显示了超出简单准确性的细节。它可以允许检测显示错误发生位置的模式。如果该模型旨在处理真实世界的数据，则可能需要收集额外的信息来支持这些数据之间的分类。此外，它还提供了一个置信区间。
confusionMatrix(predictOnTestT, as.factor(validationData$CntN))
## Confusion Matrix and Statistics
##
## Reference
## Prediction ArbWr CEatB Er&CA ErpnU EsA&P LtA&C ME&NA NrthA Sb-SA
## ArbWr 3 0 0 0 2 3 3 0 0
## CEatB 0 3 0 0 0 0 0 0 0
## Er&CA 0 0 2 0 0 0 0 0 0
## ErpnU 0 0 1 3 0 0 0 0 0
## EsA&P 0 0 0 0 1 0 0 0 0
## LtA&C 0 0 0 0 0 0 0 0 0
## ME&NA 0 0 0 0 0 0 0 0 0
## NrthA 0 0 0 0 0 0 0 3 0
## Sb-SA 0 0 0 0 0 0 0 0 3
##
## Overall Statistics
##
## Accuracy : 0.667
## 95% CI : (0.46, 0.835)
## No Information Rate : 0.111
## P-Value [Acc > NIR] : 1.15e-11
##
## Kappa : 0.625
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: ArbWr Class: CEatB Class: Er&CA
## Sensitivity 1.000 1.000 0.6667
## Specificity 0.667 1.000 1.0000
## Pos Pred Value 0.273 1.000 1.0000
## Neg Pred Value 1.000 1.000 0.9600
## Prevalence 0.111 0.111 0.1111
## Detection Rate 0.111 0.111 0.0741
## Detection Prevalence 0.407 0.111 0.0741
## Balanced Accuracy 0.833 1.000 0.8333
## Class: ErpnU Class: EsA&P Class: LtA&C
## Sensitivity 1.000 0.333 0.000
## Specificity 0.958 1.000 1.000
## Pos Pred Value 0.750 1.000 NaN
## Neg Pred Value 1.000 0.923 0.889
## Prevalence 0.111 0.111 0.111
## Detection Rate 0.111 0.037 0.000
## Detection Prevalence 0.148 0.037 0.000
## Balanced Accuracy 0.979 0.667 0.500
## Class: ME&NA Class: NrthA Class: Sb-SA
## Sensitivity 0.000 1.000 1.000
## Specificity 1.000 1.000 1.000
## Pos Pred Value NaN 1.000 1.000
## Neg Pred Value 0.889 1.000 1.000
## Prevalence 0.111 0.111 0.111
## Detection Rate 0.000 0.111 0.111
## Detection Prevalence 0.000 0.111 0.111
## Balanced Accuracy 0.500 1.000 1.000
如您所见，在这个数据集中，分类和回归树工作得相当好。虽然 PCA 方面可能会使对哪些输入值驱动哪些决策的完美理解变得复杂，但树通常可以简洁地用图表表示。此外，如果愿意，可以尝试不使用五氯苯甲醚。使用原始数据(尽管可能经过标准化和缩放)可能会产生一个足够准确且更容易解释的模型。另一方面，单棵树在某些情况下可能不够准确。和模型总是有取舍的。
随机森林
随机森林比树的概念前进了一步。回想上一节，意识到我们的树真的很擅长识别中欧。虽然它在阿拉伯世界犯了一些错误，但它在中欧(CEatB)表现良好。如果我们不是训练一棵树，而是训练很多棵树，会怎么样？当然，简单地循环 CART 算法会产生一个克隆森林，这里不需要这样做。相反，随机森林概念上的第一步是随机选取某些预测值列，并从该列子集中随机选取某些观察值。这创建了训练数据集的随机子集，并且这些子集的每一个都用于训练树。一旦树木长成森林，模型就可以进行预测了。回归预测或数值数据是每棵树预测响应值的平均值。对我们的分类或分类数据的预测将是机器学习式的民主，其中多数投票获胜(如果需要，还附带一个概率)。
我们从设置现在常用的数据集副本开始。caret的一个非常有用的特性仍然是train函数的标准化，它允许我们非常容易地在这些不同的模型之间工作。对于随机森林，我们使用method = "ranger"。模型变量num.trees被设置为控制我们森林的大小，我们从森林中只有 20 棵树开始。默认情况下，该模型在一些调整参数上运行一个小的网格搜索，稍后我们将对此进行更深入的讨论。在这种情况下，基于交叉验证选择的模型估计具有很高的准确性，只有 20 棵树。请注意，finalModel中包含了几个特定的变量，我们确认模型中只有 20 棵树。在这种情况下，初始阶段将创建 20 个随机子集，每个子集用于训练 20 棵树中的一棵树。
rfDataTrain <- copy(trainData[,-1])
rfDataTrain[,Year:=as.numeric(Year)]
rfDataValidation <- copy(validationData[,-1])
rfDataValidation[,Year:=as.numeric(Year)]
set.seed(12345)
rfModel <- train(x = rfDataTrain,
y = trainData$CntN,
method = "ranger",
preProcess = c("scale", "center", "pca"),
metric = "Accuracy",
num.trees = 20,
trControl = trainControl(method = "cv",
number = 5
)
)
rfModel
## Random Forest
##
## 135 samples
## 12 predictor
## 9 classes: 'ArbWr', 'CEatB', 'Er&CA', 'ErpnU', 'EsA&P', 'LtA&C', 'ME&NA', 'NrthA', 'Sb-SA'
##
## Pre-processing: scaled (12), centered (12), principal
## component signal extraction (12)
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 108, 108, 108, 108, 108
## Resampling results across tuning parameters:
##
## mtry splitrule Accuracy Kappa
## 2 gini 0.93 0.92
## 2 extratrees 0.96 0.96
## 3 gini 0.93 0.92
## 3 extratrees 0.96 0.96
## 4 gini 0.93 0.92
## 4 extratrees 0.99 0.98
##
## Tuning parameter 'min.node.size' was held constant at a value of 1
## Accuracy was used to select the optimal model using the
## largest value.
## The final values used for the model were mtry = 4, splitrule
## = extratrees and min.node.size = 1.
rfModel$finalModel$num.trees
## [1] 20
自然地，我们运行我们现在通常的检查，看看最终的模型在我们的测试集和验证集上表现如何。基于validationData的比较，训练数据的保留部分的估计精度可能有点高。
predictOnTrainR <- predict(rfModel, newdata = rfDataTrain)
mean( predictOnTrainR == trainData$CntN)
## [1] 1
predictOnTestR <- predict(rfModel, newdata = rfDataValidation)
mean(predictOnTestR == validationData$CntN)
## [1] 1
随机森林的计算复杂性在一定程度上是被训练的树的数量的特征——这确实也开始对未来的预测产生计算成本(虽然通常不是一个巨大的负担，但预测并不仅仅通过单棵树进行过滤)。我们接下来调整我们的代码来拥有50树。
set.seed(12345)
rfModel <- train(x = rfDataTrain,
y = trainData$CntN,
method = "ranger",
preProcess = c("scale", "center", "pca"),
metric = "Accuracy",
num.trees = 50,
trControl = trainControl(method = "cv",
number = 5
)
)
rfModel
## Random Forest
##
## 135 samples
## 12 predictor
## 9 classes: 'ArbWr', 'CEatB', 'Er&CA', 'ErpnU', 'EsA&P', 'LtA&C', 'ME&NA', 'NrthA', 'Sb-SA'
##
## Pre-processing: scaled (12), centered (12), principal
## component signal extraction (12)
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 108, 108, 108, 108, 108
## Resampling results across tuning parameters:
##
## mtry splitrule Accuracy Kappa
## 2 gini 0.94 0.93
## 2 extratrees 0.99 0.99
## 3 gini 0.93 0.92
## 3 extratrees 0.98 0.97
## 4 gini 0.93 0.92
## 4 extratrees 0.98 0.97
##
## Tuning parameter 'min.node.size' was held constant at a value of 1
## Accuracy was used to select the optimal model using the
## largest value
.
## The final values used for the model were mtry = 2, splitrule
## = extratrees and min.node.size = 1.
rfModel$finalModel$num.trees
## [1] 50
这种增加的净效果是为我们的保留验证数据提供了一个完美的匹配。至于它是否经得起未来额外的真实世界测试，则是另一回事了。重要的是要认识到每个数据集都有各种各样的特征，这些特征使得数据集在不同的模型下表现不同。在这种情况下，就目前而言，模型给我们提供了高度准确的结果。如前所述，其中一位作者处理学生的成绩数据，在这种情况下，准确的分数预测会更加困难。
predictOnTrainR <- predict(rfModel, newdata = rfDataTrain)
mean( predictOnTrainR == trainData$CntN)
## [1] 1
predictOnTestR <- predict(rfModel, newdata = rfDataValidation)
mean(predictOnTestR == validationData$CntN)
## [1] 1
提高给定模型准确性的一种方法是通过调整。对于随机森林，调整参数包括num.trees、mtry、splitrule和min.node.size。虽然已经注意到num.trees可以控制树的数量，但是其他的调整参数值得讨论。基于我们之前对 PCA 的探索，有理由假设我们的数据仍然有四个主要成分。在这种情况下，当在任何特定节点时，随机森林算法将随机选择这些预测值中的一些，然后决定这些选择的预测值中的哪一个将给出最大排序增益(回忆树分支到两个子节点中的一个)。mtry变量设置每次随机选择四个预测值中的多少个。因此，一到四是我们的范围，并成为调整我们的模型的一种方式。splitrule允许选择算法的各种细微差别。在我们的例子中，当我们使用模型对数据进行分类时，有意义的选项是gini和extratrees。由于extratrees是迄今为止最佳选择的具体方法，我们将坚持使用它来简化调整过程。最后，min.node.size决定在树停止生长之前允许多少行数据。对于分类数据，默认值为 1。然而，我们的数据有 9 个区域的 135 个观察值，这给了我们每个区域 15 个实例。因此，虽然大于 15 的最小尺寸不太理想，但尝试多个尺寸可能是有意义的。为了用这些变量调整模型，我们使用 base R中的expand.grid函数从我们选择的范围的所有可能组合中创建一个数据框。这被传递到正式的tuneGride和模型运行。
set.seed(12345)
rfModel <- train(x = rfDataTrain,
y = trainData$CntN,
method = "ranger",
preProcess = c("scale", "center", "pca"),
metric = "Accuracy",
num.trees = 20,
trControl = trainControl(method = "cv",
number = 5
),
tuneGrid = expand.grid(mtry = c(1, 2, 3, 4),
splitrule = "extratrees",
min.node.size = c(1, 5, 10, 15))
)
rfModel
## Random Forest
##
## 135 samples
## 12 predictor
## 9 classes: 'ArbWr', 'CEatB', 'Er&CA', 'ErpnU', 'EsA&P', 'LtA&C', 'ME&NA', 'NrthA', 'Sb-SA'
##
## Pre-processing: scaled (12), centered (12), principal
## component signal extraction (12)
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 108, 108, 108, 108, 108
## Resampling results across tuning parameters:
##
## mtry min.node.size Accuracy Kappa
## 1 1 0.95 0.94
## 1 5 0.94 0.93
## 1 10 0.88 0.87
## 1 15 0.77 0.74
## 2 1 0.98 0.97
## 2 5 0.96 0.95
## 2 10 0.95 0.94
## 2 15 0.90 0.88
## 3 1 0.96 0.96
## 3 5 0.95 0.94
## 3 10 0.93 0.92
## 3 15 0.88 0.87
## 4 1 0.96 0.95
## 4 5 0.96 0.95
## 4 10 0.94 0.93
## 4 15 0.90 0.89
##
## Tuning parameter 'splitrule' was held constant at a value
## of extratrees
## Accuracy was used to select the optimal model using the
## largest value.
## The final values used for the model were mtry = 2, splitrule =
## extratrees and min.node.size = 1.
rfModel$finalModel$num.trees
## [1] 20
rfModel$finalModel$mtry
## [1] 2
rfModel$finalModel$splitrule
## [1] "extratrees"
rfModel$finalModel$min.node.size
## [1] 1
值得注意的是，在这些更广泛的选项上运行模型将模型在我们的系统上运行的时间增加了 1.9 倍。我们使用函数system.time()作为包含 20 棵树的原始运行的包装器，并使用tuneGride来确定时间成本。作为一般的工作流程，模型调优是最后一个步骤，一旦选择了特定的模型，就要进行调优，正是因为这个原因。
我们最后一次运行我们模型的验证，看到在 20 棵树的限制下，调整并没有增加我们验证数据的准确性。
predictOnTrainR <- predict(rfModel, newdata = rfDataTrain)
mean( predictOnTrainR == trainData$CntN)
## [1] 1
predictOnTestR <- predict(rfModel, newdata = rfDataValidation)
mean(predictOnTestR == validationData$CntN)
## [1] 1
随机梯度推进
随机梯度推进是一种迭代创建树木“森林”的方法。这种方法和随机森林的区别在于迭代。在训练第一个模型之后，计算预测输出和已知训练数据输出之间的误差。然后，误差被用作预测值的附加响应变量，以努力降低误差水平。这个过程重复进行，减少了错误率。该技术非常有效，既不需要删除/填补缺失数据，也不需要缩放/居中/pca。收益往往伴随着权衡，在没有交叉验证的情况下运行时，模型几乎肯定会过度拟合。迭代性质和减少误差的目标在计算时间和存储需求方面也可能相对昂贵。
这一次，我们将使用回归而不是分类模型。为此，我们对国家的训练和验证数据进行虚拟编码，并再次选择将变量Year设置为一个数值。
sgbDataTrain <- copy(trainData)
sgbDataTrain[,Year:=as.numeric(Year)]
sgbDataValidation <- copy(validationData)
sgbDataValidation[,Year:=as.numeric(Year)]
ddum <- dummyVars("˜.", data = sgbDataTrain)
sgbDataTrain <- data.table(predict(ddum, newdata = sgbDataTrain))
sgbDataValidation <- data.table(predict(ddum, newdata = sgbDataValidation))
rm(ddum)
这一次，我们的火车功能确实有些不同。因为我们正在执行回归，caret确实允许第一个形式是回归函数。在这种情况下，我们选择让SPAD作为因变量，所有其他变量都将是预测变量。由于第一个形式不再包含对数据集的提及，我们显式断言了在train函数的第二个形式参数中使用的数据。熟悉的缩放和居中已经完成，尽管我们现在放弃主成分分析。此外，由于这些是回归中的数字数据，我们设置metric = "RMSE"通过均方根误差检测最佳模型。我们继续使用交叉验证。这个模型将相当多的文本打印到屏幕上，显示许多迭代。通过设置verbose = FALSE，我们否定了这一点，尽管我们鼓励感兴趣的读者将设置改为TRUE。
tuneGride选项被设置为默认设置(因此不是严格要求的，尽管为了简洁起见我们在这里包括了它们)。interaction.depth的范围从 1 到 3，所以这些树的“高度”范围从单个节点到孙节点。shrinkage的值控制迭代的移动，较小的值允许更多的微调，但可能会花费模型更长的时间来找到一个“足够好”的位置。现在熟悉的是n.trees控制森林的最大尺寸，n.minobsinnode设置为同样熟悉的 10。一般来说，只有 150 棵树可能不够，价值超过 10，000 的树在野外很常见。
set.seed(12345)
sgbModel <- train(SPAD ~.,
data = sgbDataTrain,
method = "gbm",
preProcess = c("scale", "center"),
metric = "RMSE",
trControl = trainControl(method = "cv",
number = 5
),
tuneGrid = expand.grid(interaction.depth = 1:3,
shrinkage = 0.1,
n.trees = c(50, 100, 150),
n.minobsinnode = 10),
verbose = FALSE
)
sgbModel
## Stochastic Gradient Boosting
##
## 135 samples
## 20 predictor
##
## Pre-processing: scaled (20), centered (20)
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 108, 108, 108, 107, 109
## Resampling results across tuning parameters:
##
## interaction.depth n.trees RMSE Rsquared MAE
## 1 50 14.6 0.81 11.5
## 1 100 12.3 0.86 9.8
## 1 150 11.1 0.88 8.8
## 2 50 10.4 0.91 8.3
## 2 100 7.8 0.94 5.7
## 2 150 7.0 0.95 5.0
## 3 50 8.9 0.93 6.7
## 3 100 6.8 0.95 4.9
## 3 150 6.2 0.96 4.4
##
## Tuning parameter 'shrinkage' was held constant at a value of
## 0.1
## Tuning parameter 'n.minobsinnode' was held constant at a
## value of 10
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were n.trees =
## 150, interaction.depth = 3, shrinkage = 0.1 and n.minobsinnode = 10.
所选模型有 150 棵树，每个节点至少有 10 个观察值(与早期的分类选择形成对比，在早期的分类选择中，默认的最小观察值是 1)。这一模型的 RMSE 是每 1000 名 15-19 岁的妇女中有 5.9 个孩子。
在第一次尝试中避免 PCA 有助于理解是什么以某种程度的准确性和效率为代价驱动了模型。使用summary函数，我们可以看到最有影响力的预测值。影响最大的是SPPO，显示的是劳动年龄人口的百分比。紧随其后的是统计小学失学女童的SEPR，以及每千人死亡率的SPDY。在这种情况下，summary也会生成一个图形 8-8 ，这显然不是最佳的。然而，数据打印输出是有用的，并且可视化是有帮助的。
图 8-8
相对影响可视化
summary(sgbModel)
## var rel.inf
## SPPO SPPO 31.337
## SEPR SEPR 17.903
## SPDY SPDY 11.753
## MAIN MAIN 8.223
## NYGN NYGN 6.945
## ERMA ERMA 6.127
## NYGD NYGD 5.158
## SESC SESC 4.519
## `CntNLtA&C` `CntNLtA&C` 3.912
## FEMA FEMA 2.843
## Year Year 0.504
## CntNNrthA CntNNrthA 0.337
## CntNErpnU CntNErpnU 0.250
## FEIN FEIN 0.104
## `CntNME&NA` `CntNME&NA` 0.049
## CntNArbWr CntNArbWr 0.018
## `CntNSb-SA` `CntNSb-SA` 0.016
## CntNCEatB CntNCEatB 0.000
## `CntNEr&CA` `CntNEr&CA` 0.000
## `CntNEsA&P` `CntNEsA&P` 0.000
接下来，我们将注意力转向理解我们的模型的准确性。residuals函数用于计算训练模型对训练数据的 MSE(回想一下机器学习介绍章节中关于 MSE 的讨论)。我们将此与验证数据的预测值和实际值进行比较，发现模型很可能过度拟合。
mean(stats::residuals(sgbModel)²)
## [1] 3.6
mean((predict(sgbModel, sgbDataValidation) -
sgbDataValidation$SPAD)²)
## [1] 15
这种方法的一个挑战是，任何特定模型本身的精确机制都非常复杂，很难理解是什么驱动了预测。软件包DALEX有一个名为explain的函数，它试图帮助理解模型的更多内容，以及预测器如何与响应进行交互。虽然这个软件包非常适用于比较测试数据上的各种模型，但语法可以通过训练和验证数据来学习。
explain函数采用模型作为第一个形式，一个标签来区分后面的图形(这里我们区分训练和验证)，一个指针指向我们的data，以及一个响应变量SPAD(青少年生育率)的断言。
explainSGBt <- explain(sgbModel, label = "sgbt",
data = sgbDataTrain,
y = sgbDataTrain$SPAD)
explainSGBv <- explain(sgbModel, label = "sgbv",
data = sgbDataValidation,
y = sgbDataValidation$SPAD)
对于本文来说，explain对象本身并不特别重要。感兴趣的是能从他们身上找到的信息，函数model_performance是这些中的第一个。重用cowplot包中的plot_grid函数，绘制残差的分布和箱线图。图 8-9 中最左边的图表是残差绝对值的曲线图。记住正态分布残差的一般假设，然后想象绝对值变换将使正态分布看起来像两倍高度的半个分布，我们看到验证集看起来基本上仍然是正态的，尽管它比定型集更容易出错。第二张图显示了同样的事实，箱线图上的范围更大。
图 8-9
DALEX 残差可视化
performanceSGBt <- model_performance(explainSGBt)
performanceSGBv <- model_performance(explainSGBv)
plot_grid(
plot(performanceSGBt, performanceSGBv),
plot(performanceSGBt, performanceSGBv, geom = "boxplot"),
ncol = 2)
虽然我们已经看到了模型权重的相对影响，但是函数variable_importance计算了在没有来自各种变量的信息的情况下会出现的丢失损失。它对每个数据集都这样做，所以请注意，在图 8-10 中，训练数据的值不同于验证数据集的值。尽管如此，最顶端的关键变量确实与相对影响变量匹配(尽管顺序不同)。
图 8-10
DALEX 辍学损失
importanceSGBt <- variable_importance(explainSGBt)
importanceSGBv <- variable_importance(explainSGBv)
plot(importanceSGBt, importanceSGBv)
我们讨论的最后一个DALEX函数是variable_response函数，它接受一个单变量的附加形式参数。在这种情况下，我们选择部分相关图类型，并在图 8-11 中看到结果图。小学(SEPR)以外的女童人数与模型结果相对照。该图的价值在于，尽管模型本身可能不容易理解(不像简单的回归方程)，但青少年生育率和小学出勤率之间的关系是可以探究的。
图 8-11
DALEX 小学缺课计数与少女怀孕。
responseSGBprmt <- variable_response(explainSGBt, variable = "SEPR", type = "pdp")
responseSGBprmv <- variable_response(explainSGBv, variable = "SEPR", type = "pdp")
plot(responseSGBprmt, responseSGBprmv)
我们使用相同的设置运行代码，只是这次每 1000 人的死亡率(SPDY)如图 8-12 所示。这里，关系比较复杂。这也许说明了为什么这种类型的模型可能非常适合处理与预测值有复杂关系的数据。它还显示了为什么这个模型有过度拟合的风险。
图 8-12
DALEX 每 1，000 人死亡率与少女怀孕
responseSGBdynt <- variable_response(explainSGBt, variable = "SPDY", type = "pdp")
responseSGBdynv <- variable_response(explainSGBv, variable = "SPDY", type = "pdp")
plot(responseSGBdynt, responseSGBdynv)
下一步是什么？嗯，这个模型并没有显示出很高的精确度。扩展调整网格以允许更多的树、更深的树和其他可变模式，从而为算法提供更宽的搜索网格以获得最佳模型。我们把这样的扩展留给感兴趣的读者，因为结构已经构建好了，然后继续下一个模型。
多层感知器
多层感知器(MLP)是一种前馈人工神经网络。这实际上并不完全是最新和最伟大的深度学习或神经网络算法。尽管如此，这是一种人们倾向于从学习开始的神经网络类型。在它的核心，神经网络只是一个由算法优化的线性代数矩阵乘法。净效应是，不是每个预测变量有一个权重的标准回归方程，而是有一个相当大的权重数组。
神经网络的一个关键特征是使用多个方程(因此是矩阵)，而不是试图用单个方程来创建预测值。实际上，预测器的每一部分都将被单独处理、加权并映射到可能的输出。此外，在输入和输出之间可能有几层方程式，而不是使用单一的方程式。这些被恰当地命名为“隐藏层”，虽然它们可以给模型带来很大的灵活性，但大多数多层感知器作为一个模型很难理解。他们的可取之处是他们经常做出相当准确的预测。
事实上，经典的例子是通过图像视觉识别数字或字符。我们的例子要简单得多，可能不太适合这种方法(事实上不是)。这并不以任何方式否定该方法在其他情况下的功效。和以前一样，我们使用虚拟变量创建最后一轮训练和验证数据。
mlpDataTrain <- copy(trainData)
mlpDataTrain[,Year:=as.numeric(Year)]
mlpDataValidation <- copy(validationData)
mlpDataValidation[,Year:=as.numeric(Year)]
ddum <- dummyVars("˜.", data = mlpDataTrain)
mlpDataTrain <- data.table(predict(ddum, newdata = mlpDataTrain))
mlpDataValidation <- data.table(predict(ddum, newdata = mlpDataValidation))
rm(ddum)
模型和种子保持不变。这一次，使用了代表多层感知器的method = "mlpML"。尽管我们不使用主成分分析，但我们会继续扩展和集中我们的数据。第一次运行是通过默认方法完成的，该方法将第二层和第三层清空。这使得第一个隐藏层将被运行，并且默认情况下，模型在第一个层中模拟一些不同数量的节点。由于这本书的创建方式，警告被打印出来(事实上这通常是一个有用的特性)。对于这种模型，会发出许多警告，在这种情况下，这些警告仅显示超出本文范围的模型细微差别。函数supressWarnings用于其同名的结果。
set.seed(12345)
suppressWarnings(
mlpModel <- train(
SPAD ~ .,
data = mlpDataTrain,
method = "mlpML",
preProcess = c("scale", "center"),
metric = "RMSE",
trControl = trainControl(method = "cv",
number = 5)
)
)
mlpModel
## Multi-Layer Perceptron, with multiple layers
##
## 135 samples
## 20 predictor
##
## Pre-processing: scaled (20), centered (20)
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 108, 108, 108, 107, 109
## Resampling results across tuning parameters:
##
## layer1 RMSE Rsquared MAE
## 1 36 0.31 29
## 3 25 0.49 19
## 5 17 0.76 13
##
## Tuning parameter 'layer2' was held constant at a value of 0
##
## Tuning parameter 'layer3' was held constant at a value of 0
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were layer1 = 5, layer2 = 0
## and layer3 = 0.
从输出中可以看出，默认情况下，单个隐藏层被调整为五个节点。summary函数提供了更多的细节，显示了 20 个输入列、5 个隐藏节点和 1 个输出节点的激活函数权重和偏差。这也显示了我们网络的整体布局，20-5-1 布局。
summary(mlpModel)
## SNNS network definition file V1.4-3D
## generated at Fri Nov 02 19:23:07 2018
##
## network name : RSNNS_untitled
## source files :
## no. of units : 26
## no. of connections : 105
## no. of unit types : 0
## no. of site types : 0
##
##
## learning function : Std_Backpropagation
## update function : Topological_Order
##
##
## unit default section :
##
## act | bias | st | subnet | layer | act func | out func
## ---------|----------|----|--------|-------|--------------|-------------
## 0.00000 | 0.00000 | i | 0 | 1 | Act_Logistic | Out_Identity
## ---------|----------|----|--------|-------|--------------|-------------
##
##
## unit definition section :
##
## no. | typeName | unitName | act | bias | st | position | act func | out func | sites
## ----|----------|-------------------|----------|----------|----|----------|--------------|----------|-------
## 1 | | Input_CntNArbWr | -0.35224 | 0.25864 | i | 1, 0, 0 | Act_Identity | |
## 2 | | Input_CntNCEatB | -0.35224 | -0.07158 | i | 2, 0, 0 | Act_Identity | |
## 3 | | Input_`CntNEr&CA` | -0.35224 | 0.17340 | i | 3, 0, 0 | Act_Identity | |
## 4 | | Input_CntNErpnU | -0.35224 | 0.09913 | i | 4, 0, 0 | Act_Identity | |
## 5 | | Input_`CntNEsA&P` | -0.35224 | 0.02550 | i | 5, 0, 0 | Act_Identity | |
## 6 | | Input_`CntNLtA&C` | -0.35224 | -0.07856 | i | 6, 0, 0 | Act_Identity | |
## 7 | | Input_`CntNME&NA` | -0.35224 | 0.10749 | i | 7, 0, 0 | Act_Identity | |
## 8 | | Input_CntNNrthA | -0.35224 | -0.17845 | i | 8, 0, 0 | Act_Identity | |
## 9 | | Input_`CntNSb-SA` | 2.81793 | 0.20316 | i | 9, 0, 0 | Act_Identity | |
## 10 | | Input_Year | 1.51162 | 0.09500 | i | 10, 0, 0 | Act_Identity | |
## 11 | | Input_NYGD | -0.84754 | -0.12790 | i | 11, 0, 0 | Act_Identity | |
## 12 | | Input_NYGN | -0.83754 | -0.26720 | i | 12, 0, 0 | Act_Identity | |
## 13 | | Input_SEPR | 2.37160 | 0.23106 | i | 13, 0, 0 | Act_Identity | |
## 14 | | Input_ERMA | 2.28315 | -0.13108 | i | 14, 0, 0 | Act_Identity | |
## 15 | | Input_SESC | -1.27208 | -0.13525 | i | 15, 0, 0 | Act_Identity | |
## 16 | | Input_FEMA | -1.28836 | 0.21397 | i | 16, 0, 0 | Act_Identity | |
## 17 | | Input_SPDY | 0.55710 | 0.19976 | i | 17, 0, 0 | Act_Identity | |
## 18 | | Input_FEIN | -1.80461 | -0.15568 | i | 18, 0, 0 | Act_Identity | |
## 19 | | Input_MAIN | -1.77002 | 0.29051 | i | 19, 0, 0 | Act_Identity | |
## 20 | | Input_SPPO | 2.09958 | -0.29660 | i | 20, 0, 0 | Act_Identity | |
## 21 | | Hidden_2_1 | 1.00000 | 16.57569 | h | 1, 2, 0 | | |
## 22 | | Hidden_2_2 | 1.00000 | -42.96722| h | 2, 2, 0 | | |
## 23 | | Hidden_2_3 | 0.00000 | -49.55274| h | 3, 2, 0 | | |
## 24 | | Hidden_2_4 | 1.00000 | 39.86870 | h | 4, 2, 0 | | |
## 25 | | Hidden_2_5 | 1.00000 | 21.29272 | h | 5, 2, 0 | | |
## 26 | | Output_1 | 66.31960 |-1069.78821| o | 1, 4, 0 | Act_Identity | |
## ----|----------|-------------------|----------|----------|----|----------|--------------|----------|-------
##
##
## connection definition section :
##
## target | site | source:weight
## -------|------|---------------------------------------------------------
## 21 | | 20:48.96251, 19:-47.13240, 18:-46.89282, 17:29.71989, 16:-31.71886, 15:-30.48963, 14:58.16447, 13:58.57929, 12:-15.94363,
## 11:-15.49402, 10:21.26738, 9:66.89728, 8:-5.93687, 7:-24.55808, 6:-8.02125, 5:-5.53956, 4:-5.59685, 3:-5.55265,
## 2:-5.77468, 1:-5.90555
## 22 | | 20:19.36864, 19:-36.07077, 18:-30.54533, 17:30.76214, 16:-19.48156, 15:-10.87526, 14:40.41421, 13:38.18938, 12: 4.44169,
## 11: 0.46633, 10:46.79949, 9:49.01194, 8:15.15641, 7:-81.96992, 6:14.43461, 5:-6.70102, 4:-40.07672, 3:12.00658,
## 2: 7.56544, 1:31.19720
## 23 | | 20:-12.19030, 19: 2.85565, 18:31.55284, 17:84.06341, 16:77.33264, 15:103.85283, 14:-12.62139, 13:-9.16594, 12:81.62211,
## 11:81.21730, 10:-194.05856, 9:49.49640, 8:115.78496, 7:-4.44688, 6:179.36331, 5:-26.62095, 4:35.70350, 3:-54.66271,
## 2:33.46564, 1:-328.59045
## 24 | | 20:22.24673, 19:-17.17476, 18:-17.65513, 17:-30.85148, 16:-20.34034, 15:-17.87234, 14:19.58477, 13:16.31513, 12:-25.13864,
## 11:-24.56263, 10:-14.12056, 9:11.75429, 8:-14.03880, 7:-0.54804, 6:62.57944, 5:14.06488, 4:-13.82649, 3:-14.21823,
## 2:-31.15477, 1:-13.92656
## 25 | | 20:35.84281, 19:-31.83327, 18:-33.52740, 17: 7.88547, 16:-23.72048, 15:-24.41236, 14:36.01567, 13:37.40243, 12:-16.33773,
## 11:-15.88420, 10:19.37509, 9:42.81319, 8:-7.29130, 7:10.05342, 6:-4.69724, 5:-11.28413, 4:-7.28304, 3:-7.42553,
## 2:-7.40186, 1:-7.59370
## 26 | | 25: 5.30307, 24:31.84159, 23:30.69779, 22:36.21242, 21:-7.03748
## -------|------|--------------------------------------------------------
正如对相对简单的数据集使用过于复杂的方法所暗示的那样，该模型的表现并不那么好。它在验证数据上的表现甚至更差。注意，由于residuals函数被一个必需的包屏蔽了，我们明确地要求residuals的基础R版本。
mean(stats::residuals(mlpModel)ˆ2)
## [1] 462
mean((predict(mlpModel, mlpDataValidation) -
mlpDataValidation$SPAD)ˆ2)
## [1] 407
同样，这种方法的一个挑战是任何特定模型本身的精确机制都足够复杂，以至于很难理解是什么驱动了预测。我们再次使用包DALEX恰当命名的explain函数来显示预测器如何与响应交互。
回想一下，explain函数语法将模型用作第一种形式，一个用于区分后面图形的标签(这里我们区分训练和验证)，一个指向我们的data的指针，以及一个我们的响应变量SPAD的断言。
explainMLPt <- explain(mlpModel, label = "mlpt",
data = mlpDataTrain,
y = mlpDataTrain$SPAD)
explainMLPv <- explain(mlpModel, label = "mlpv",
data = mlpDataValidation,
y = mlpDataValidation$SPAD)
还是像以前一样，explain物体本身并不特别有趣。感兴趣的是能从他们身上找到的信息，函数model_performance是这些中的第一个。重用cowplot包中的plot_grid函数，绘制残差的分布和箱线图。图 8-13 中最左边的图表是残差绝对值的曲线图。记住正态分布残差的一般假设，然后想象绝对值变换将使正态分布看起来像两倍高度的半个分布，我们看到验证集看起来基本上仍然是正态的，尽管它比定型集更容易出错。我们还看到前面章节中的gbm方法不容易出错。第二张图表继续显示同样的事实，箱线图上的范围更大。在这两种情况下，图表也显示了过度训练的永久风险。模型数据看不到的验证数据具有更大的残差。
图 8-13
对比 SGB 和 MLP 方法的模型性能
performanceMLPt <- model_performance(explainMLPt)
performanceMLPv <- model_performance(explainMLPv)
plot_grid(
plot(performanceMLPt, performanceMLPv, performanceSGBt, performanceSGBv),
plot(performanceMLPt, performanceMLPv, performanceSGBt, performanceSGBv, geom = "boxplot"),
ncol = 2
)
我们使用variable_importance函数来计算在没有各种变量信息的情况下的辍学损失。有趣的是，在这些图中——以及之前的gbm模型中的图——哪些变量是最关键的。回想一下退学后损失最大的变量，如图 8-14 所示。这表明我们的多层感知器模型这些相同的数据非常不同。在这种情况下，这似乎是一个弱点。然而，从其他数据来看，这种差异很可能是一种优势。
图 8-14
确定主要变量
importanceMLPt <- variable_importance(explainMLPt)
importanceMLPv <- variable_importance(explainMLPv)
plot(importanceMLPt, importanceMLPv, importanceSGBt, importanceSGBv)
如前所述，我们以DALEX函数variable_response结束，它接受一个单变量的附加形式参数。在这种情况下，我们选择部分相关图类型，并在图 8-15 中看到结果图。小学以外的女生人数与模型结果相对照。该图的价值仍然在于，尽管模型本身可能不容易理解(不像简单的回归方程)，但青少年生育率和小学出勤率之间的关系是可以探索的。
图 8-15
了解小学入学率对少女怀孕的影响
responseMLPprmt <- variable_response(explainMLPt, variable = "SEPR", type = "pdp")
responseMLPprmv <- variable_response(explainMLPv, variable = "SEPR", type = "pdp")
plot(responseMLPprmt, responseMLPprmv, responseSGBprmt, responseSGBprmv)
我们使用如图 8-16 所示的每 1000 人死亡率来结束我们的DALEX探索。在这里，两种模式之间的对比更加鲜明。这再次证明了潜在的方法是完全不同的，因此mlpML可能适用于早期方法未能实现有效预测的其他数据。
图 8-16
DALEX 每 1，000 人死亡率与少女怀孕
responseMLPdynt <- variable_response(explainMLPt, variable = "SPDY", type = "pdp")
responseMLPdynv <- variable_response(explainMLPv, variable = "SPDY", type = "pdp")
plot(responseMLPdynt, responseMLPdynv, responseSGBdynt, responseSGBdynv)
在结束本节之前，我们考虑一下我们的模型是否调得不够好。由于其复杂性，多层感知器往往有许多调整参数。对于caret，尤其是mlpML，调整参数是隐藏层的数量，以及每个隐藏层中的节点数量。使用我们现在熟悉的tuneGride形式和expand.grid函数，我们让模型探索一千多个选项，以确定一个最佳模型。对于现实生活中的应用程序，感知器方法被认为是可能的最佳模型，拥有数百个节点并不罕见。这就产生了数以万计的变量，每个变量都被模型依次优化。这是对计算资源的重大投资，在这种情况下，将不会有回报。
set.seed(12345)
suppressWarnings(
mlpModelb <- train(
SPAD ~ .,
data = mlpDataTrain,
method = "mlpML",
preProcess = c("scale", "center"),
metric = "RMSE",
verbose = FALSE,
trControl = trainControl(method = "cv",
number = 5),
tuneGrid = expand.grid(
layer1 = 0:10,
layer2 = 0:10,
layer3 = 0:10
)
)
)
mlpModelb
## Multi-Layer Perceptron, with multiple layers
##
## 135 samples
## 20 predictor
##
## Pre-processing: scaled (20), centered (20)
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 108, 108, 108, 107, 109
## Resampling results across tuning parameters:
##
## layer1 layer2 layer3 RMSE Rsquared MAE
## 0 0 0 NaN NaN NaN
## 0 0 1 41 1.4e-01 30.4
## 0 0 2 36 3.0e-01 27.1
## 0 0 3 26 4.9e-01 17.5
## 0 0 4 18 7.4e-01 14.6
## 0 0 5 18 7.9e-01 14.2
## 0 0 6 16 7.7e-01 11.9
## 0 0 7 14 8.3e-01 10.7
## 0 0 8 16 8.0e-01 12.0
## 0 0 9 15 8.1e-01 9.8
## 0 0 10 11 8.8e-01 8.9
## 0 1 0 37 2.6e-01 29.9
## 0 1 1 33 3.0e-01 25.1
## 0 1 2 36 4.0e-01 28.5
## 0 1 3 38 3.2e-01 28.7
## 0 1 4 35 6.7e-02 24.7
## 0 1 5 33 2.3e-01 24.4
## 0 1 6 35 1.4e-01 26.1
## 0 1 7 34 2.0e-01 28.5
## 0 1 8 56 3.0e-01 50.3
## 0 1 9 41 3.7e-01 36.8
## 0 1 10 41 1.3e-01 32.9
## 0 2 0 28 3.5e-01 18.8
## 0 2 1 42 2.5e-01 32.2
## 0 2 2 40 1.7e-01 33.7
## 0 2 3 36 4.5e-01 27.7
## 0 2 4 34 2.5e-01 24.6
## 0 2 5 34 4.1e-01 28.1
## 0 2 6 32 3.0e-01 23.2
## 0 2 7 35 3.8e-01 27.6
## 0 2 8 39 2.6e-01 29.9
## 0 2 9 40 2.3e-01 31.9
## 0 2 10 249 3.3e-01 243.2
## 0 3 0 24 5.8e-01 19.1
## 0 3 1 40 2.3e-01 30.7
## 0 3 2 40 3.9e-04 33.2
## 0 3 3 38 2.7e-01 30.7
## 0 3 4 35 4.4e-01 26.2
## 0 3 5 36 8.9e-02 30.0
## 0 3 6 32 4.3e-01 21.5
## 0 3 7 32 2.1e-01 25.8
## 0 3 8 46 2.3e-01 37.9
## 0 3 9 30 3.0e-01 21.6
## 0 3 10 48 3.3e-01 42.4
## 0 4 0 22 5.7e-01 17.2
## 0 4 1 51 2.9e-01 42.1
## 0 4 2 38 2.3e-01 30.9
## 0 4 3 35 5.0e-02 27.3
## 0 4 4 36 2.2e-01 24.4
## 0 4 5 35 2.1e-01 29.9
## 0 4 6 31 2.6e-01 23.2
## 0 4 7 38 2.2e-01 28.8
## 0 4 8 31 4.1e-01 23.2
## 0 4 9 33 2.8e-01 27.1
## 0 4 10 35 3.1e-01 27.5
## 0 5 0 24 5.6e-01 16.5
## 0 5 1 37 3.7e-01 29.4
## 0 5 2 32 2.8e-02 25.2
## 0 5 3 34 2.6e-01 25.3
## 0 5 4 31 2.3e-01 22.6
## 0 5 5 40 2.0e-01 31.4
## 0 5 6 33 8.7e-02 25.3
## 0 5 7 30 2.7e-01 22.3
## 0 5 8 35 1.9e-01 25.3
## 0 5 9 28 3.3e-01 21.1
## 0 5 10 34 3.9e-01 27.1
## 0 6 0 18 7.4e-01 14.7
## 0 6 1 34 NaN 24.0
## 0 6 2 37 1.3e-01 26.5
## 0 6 3 32 3.9e-01 22.3
## 0 6 4 35 2.5e-01 24.4
## 0 6 5 33 2.1e-01 26.0
## 0 6 6 52 1.9e-01 42.8
## 0 6 7 30 2.7e-01 23.8
## 0 6 8 36 3.3e-01 27.8
## 0 6 9 54 1.5e-01 45.9
## 0 6 10 30 2.8e-01 23.7
## 0 7 0 22 7.0e-01 17.0
## 0 7 1 38 4.4e-02 30.6
## 0 7 2 42 1.5e-01 35.7
## 0 7 3 35 1.2e-01 24.6
## 0 7 4 33 2.4e-01 25.3
## 0 7 5 39 1.5e-01 32.4
## 0 7 6 32 1.7e-01 24.4
## 0 7 7 32 2.7e-01 26.3
## 0 7 8 41 1.3e-01 33.6
## 0 7 9 33 1.4e-01 25.6
## 0 7 10 31 2.5e-01 22.1
## 0 8 0 15 8.3e-01 11.9
## 0 8 1 34 NaN 25.8
## 0 8 2 33 6.2e-02 25.1
## 0 8 3 30 3.2e-01 22.5
## 0 8 4 43 1.8e-01 37.1
## 0 8 5 39 3.3e-01 31.3
## 0 8 6 34 1.1e-01 26.0
## 0 8 7 34 1.3e-01 24.7
## 0 8 8 45 1.2e-01 38.6
## 0 8 9 32 3.1e-01 23.3
## 0 8 10 50 4.4e-01 42.3
## 0 9 0 15 8.4e-01 11.9
## 0 9 1 34 NaN 24.7
## 0 9 2 34 1.6e-01 26.6
## 0 9 3 34 5.2e-01 25.5
## 0 9 4 43 2.4e-01 33.8
## 0 9 5 44 2.2e-01 37.4
## 0 9 6 37 3.6e-01 30.3
## 0 9 7 39 2.8e-01 27.3
## 0 9 8 35 1.3e-01 27.1
## 0 9 9 31 2.9e-01 22.3
## 0 9 10 51 2.2e-01 41.2
## 0 10 0 19 7.6e-01 15.1
## 0 10 1 44 NaN 35.5
## 0 10 2 33 NaN 26.3
## 0 10 3 37 4.8e-01 29.8
## 0 10 4 33 9.2e-02 24.2
## 0 10 5 35 1.4e-01 25.9
## 0 10 6 35 1.7e-01 27.7
## 0 10 7 32 1.6e-01 23.0
## 0 10 8 67 5.1e-02 62.2
## 0 10 9 31 1.6e-01 23.7
## 0 10 10 32 3.1e-01 25.1
## 1 0 0 39 4.1e-01 31.2
## 1 0 1 37 5.3e-01 29.1
## 1 0 2 40 3.5e-01 31.7
## 1 0 3 34 7.9e-02 24.5
## 1 0 4 36 3.8e-01 24.6
## 1 0 5 33 3.8e-01 24.7
## 1 0 6 35 3.4e-01 27.3
## 1 0 7 37 2.3e-01 27.4
## 1 0 8 30 3.8e-01 21.4
## 1 0 9 38 2.7e-01 30.3
## 1 0 10 40 2.8e-01 32.9
## 1 1 0 37 2.2e-01 29.0
## 1 1 1 42 7.3e-01 32.5
## 1 1 2 38 6.4e-01 31.6
## 1 1 3 33 NaN 25.0
## 1 1 4 36 4.5e-01 24.4
## 1 1 5 32 7.1e-01 26.1
## 1 1 6 34 2.7e-01 25.7
## 1 1 7 39 3.2e-01 31.9
## 1 1 8 37 5.4e-01 30.0
## 1 1 9 37 2.5e-01 28.4
## 1 1 10 35 9.0e-02 27.5
## 1 2 0 40 4.8e-01 34.1
## 1 2 1 34 NaN 25.7
## 1 2 2 33 5.0e-01 24.8
## 1 2 3 36 5.1e-01 26.6
## 1 2 4 40 4.6e-01 32.0
## 1 2 5 33 5.1e-01 23.7
## 1 2 6 41 5.5e-01 33.0
## 1 2 7 32 4.7e-01 23.3
## 1 2 8 33 3.5e-01 25.2
## 1 2 9 38 4.5e-01 29.5
## 1 2 10 39 4.5e-01 30.5
## 1 3 0 32 4.1e-01 24.1
## 1 3 1 34 NaN 26.3
## 1 3 2 34 2.8e-01 27.0
## 1 3 3 41 NaN 31.2
## 1 3 4 33 6.5e-01 23.4
## 1 3 5 31 5.4e-01 23.4
## 1 3 6 41 6.8e-01 33.8
## 1 3 7 44 5.0e-01 32.6
## 1 3 8 57 4.2e-01 47.5
## 1 3 9 41 4.7e-01 34.1
## 1 3 10 42 2.8e-01 32.6
## 1 4 0 38 4.3e-01 28.8
## [ reached getOption("max.print") -- omitted 1165 rows ]
##
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were layer1 = 0, layer2 = 0
## and layer3 = 10.
mean(stats::residuals(mlpModelb)²)
## [1] 668
mean((predict(mlpModelb, mlpDataValidation) -
mlpDataValidation$SPAD)²)
## [1] 552
summary(mlpModelb)
## SNNS network definition file V1.4-3D
## generated at Fri Nov 02 19:34:18 2018
##
## network name : RSNNS_untitled
## source files :
## no. of units : 31
## no. of connections : 210
## no. of unit types : 0
## no. of site types : 0
##
##
## learning function : Std_Backpropagation
## update function : Topological_Order
##
##
## unit default section :
##
## act | bias | st | subnet | layer | act func | out func
## ---------|----------|----|--------|-------|--------------|-------------
## 0.00000 | 0.00000 | i | 0 | 1 | Act_Logistic | Out_Identity
## ---------|----------|----|--------|-------|--------------|-------------
##
##
## unit definition section :
##
## no. | typeName | unitName | act | bias | st | position | act func | out func | sites
## ----|----------|-------------------|----------|----------|----|----------|--------------|----------|-------
## 1 | | Input_CntNArbWr | -0.35224 | 0.27237 | i | 1, 0, 0 | Act_Identity | |
## 2 | | Input_CntNCEatB | -0.35224 | 0.12640 | i | 2, 0, 0 | Act_Identity | |
## 3 | | Input_`CntNEr&CA` | -0.35224 | 0.14994 | i | 3, 0, 0 | Act_Identity | |
## 4 | | Input_CntNErpnU | -0.35224 | 0.05916 | i | 4, 0, 0 | Act_Identity | |
## 5 | | Input_`CntNEsA&P` | -0.35224 | -0.03508 | i | 5, 0, 0 | Act_Identity | |
## 6 | | Input_`CntNLtA&C` | -0.35224 | 0.20488 | i | 6, 0, 0 | Act_Identity | |
## 7 | | Input_`CntNME&NA` | -0.35224 | -0.18422 | i | 7, 0, 0 | Act_Identity | |
## 8 | | Input_CntNNrthA | -0.35224 | -0.24506 | i | 8, 0, 0 | Act_Identity | |
## 9 | | Input_`CntNSb-SA` | 2.81793 | -0.11938 | i | 9, 0, 0 | Act_Identity | |
## 10 | | Input_Year | 1.70822 | 0.19864 | i | 10, 0, 0 | Act_Identity | |
## 11 | | Input_NYGD | -0.84501 | 0.20640 | i | 11, 0, 0 | Act_Identity | |
## 12 | | Input_NYGN | -0.83288 | 0.20374 | i | 12, 0, 0 | Act_Identity | |
## 13 | | Input_SEPR | 2.38957 | 0.21969 | i | 13, 0, 0 | Act_Identity | |
## 14 | | Input_ERMA | 2.33470 | -0.14782 | i | 14, 0, 0 | Act_Identity | |
## 15 | | Input_SESC | -1.32473 | -0.08605 | i | 15, 0, 0 | Act_Identity | |
## 16 | | Input_FEMA | -1.37682 | 0.28955 | i | 16, 0, 0 | Act_Identity | |
## 17 | | Input_SPDY | 0.44706 | 0.01856 | i | 17, 0, 0 | Act_Identity | |
## 18 | | Input_FEIN | -1.72852 | -0.10129 | i | 18, 0, 0 | Act_Identity | |
## 19 | | Input_MAIN | -1.69223 | -0.26429 | i | 19, 0, 0 | Act_Identity | |
## 20 | | Input_SPPO | 2.06062 | -0.17451 | i | 20, 0, 0 | Act_Identity | |
## 21 | | Hidden_2_1 | 1.00000 | -8.20578 | h | 1, 2, 0 | | |
## 22 | | Hidden_2_2 | 0.00000 | -61.29456| h | 2, 2, 0 | | |
## 23 | | Hidden_2_3 | 1.00000 | -4.93293 | h | 3, 2, 0 | | |
## 24 | | Hidden_2_4 | 0.00000 | -1.26154 | h | 4, 2, 0 | | |
## 25 | | Hidden_2_5 | 1.00000 | -136.45082| h | 5, 2, 0 | | |
## 26 | | Hidden_2_6 | 1.00000 | -3.98742 | h | 6, 2, 0 | | |
## 27 | | Hidden_2_7 | 1.00000 | -38.50706| h | 7, 2, 0 | | |
## 28 | | Hidden_2_8 | 1.00000 | -59.23545| h | 8, 2, 0 | | |
## 29 | | Hidden_2_9 | 1.00000 | -13.07257| h | 9, 2, 0 | | |
## 30 | | Hidden_2_10 | 0.00000 | -3.80823 | h | 10, 2, 0 | | |
## 31 | | Output_1 | 80.51889 | 165.28026| o | 1, 4, 0 | Act_Identity | |
## ----|----------|-------------------|----------|----------|----|----------|--------------|----------|-------
##
##
## connection definition section :
##
## target | site | source:weight
## -------|------|--------------------------------------------------------
## 21 | | 20:14.99080, 19:-4.37922, 18:-5.64107, 17:-14.21278, 16:-7.40044, 15:-6.96982, 14: 2.04649, 13: 3.28454, 12:-1.09398,
## 11:-1.05998, 10:-10.35237, 9: 3.07458, 8: 2.81771, 7: 7.14501, 6:22.43020, 5:-17.89252, 4: 2.94234, 3:-3.41989,
## 2:-20.86507, 1: 2.99437
## 22 | | 20:-42.11395, 19:46.16928, 18:46.41003, 17: 5.73816, 16:50.23483, 15:47.74510, 14:-48.54859, 13:-49.08794, 12:49.66853,
## 11:51.15318, 10:-5.34949, 9:-55.80005, 8:25.68831, 7:-2.84146, 6:-12.38087, 5:-27.17203, 4:23.60743, 3:24.90088,
## 2: 1.37566, 1:22.42043
## 23 | | 20: 3.55224, 19:-3.82663, 18:-2.81083, 17:-8.62280, 16:-6.84271, 15:-6.68200, 14: 2.13178, 13: 1.19267, 12:-12.25461,
## 11:-12.00026, 10: 7.13982, 9:-0.63003, 8:-19.36786, 7: 1.56885, 6:18.87510, 5: 1.15033, 4: 1.62183, 3: 1.71615,
## 2:-6.58669, 1: 1.57223
## 24 | | 20:-3.03815, 19: 2.09218, 18: 5.45140, 17: 8.87992, 16: 8.22070, 15:10.87528, 14:-0.88254, 13:-2.19428, 12: 8.49697,
## 11: 8.60248, 10:-5.41809, 9: 0.47726, 8:10.84176, 7:-27.30257, 6:14.88415, 5:-4.98278, 4: 0.67941, 3: 1.06026,
## 2: 5.48360, 1:-0.97088
## 25 | | 0:115.12606, 19:-66.18334, 18:-81.21451, 17:59.74724, 16:-65.24387, 15:-91.23738, 14:48.93630, 13:63.63233, 12:17.89186,
## 11:13.69802, 10:-49.68809, 9:48.28717, 8:-4.37396, 7:36.57324, 6:-239.25703, 5:-116.16022, 4:45.37873, 3:47.11736,
## 2:-11.51299, 1:193.68028
## 26 | | 20: 1.62858, 19:-6.50968, 18:-10.44102, 17:-22.89452, 16:-21.50499, 15:-24.28332, 14:12.23182, 13: 8.47606, 12:-39.59196,
## 11:-38.13202, 10:14.99837, 9: 1.17484, 8:-42.70610, 7: 0.45940, 6:19.33571, 5:48.32096, 4: 2.15939, 3: 1.45809,
## 2:-34.17036, 1: 4.24288
## 27 | | 20:24.32948, 19:-21.69843, 18:-29.03406, 17:-39.08598, 16:-39.70897, 15:-33.39420, 14:22.44779, 13:21.58632, 12:-34.10506,
## 11:-30.81281, 10: 8.56867, 9:13.75277, 8:15.13938, 7: 4.94741, 6:29.65086, 5:20.66949, 4:-30.40738, 3:-89.47778,
## 2: 7.78840, 1:27.43818
## 28 | | 20:74.58337, 19:-73.31054, 18:-54.84992, 17:38.25743, 16:-26.25908, 15:-6.00562, 14:51.08411, 13:57.09294, 12:-2.95093,
## 11: 5.12013, 10:-112.35419, 9:81.15020, 8:53.54121, 7:63.02877, 6:169.16539, 5:-121.10607, 4:-34.08242, 3:-29.48601,
## 2:-13.01993, 1:-169.10889
## 29 | | 20:29.63905, 19:-10.05235, 18:-26.69251, 17:-51.13666, 16:-14.93525, 15:-15.09192, 14:20.62878, 13:20.32225, 12:-9.60404,
## 11:-9.00413, 10:72.53555, 9:19.84532, 8: 7.09636, 7:60.42725, 6:35.26074, 5:-6.58941, 4: 4.21469, 3:-123.86588,
## 2:-1.78124, 1: 5.16048
## 30 | | 20:-1.49977, 19: 0.78230, 18: 3.44946, 17: 9.52544, 16: 5.21607, 15: 7.46949, 14: 2.12706, 13: 0.74737, 12: 7.23410,
## 11: 7.26661, 10:-4.33287, 9: 1.54910, 8: 7.26920, 7:-28.82655, 6:10.88533, 5: 1.27544, 4: 1.88801, 3: 1.09380,
## 2: 3.79257, 1: 1.51486
## 31 | | 30: 9.60220, 29:-8.94521, 28:33.58144, 27: 1.55787, 26: 8.52473, 25:45.79420, 24:-4.24323, 23:-18.14175, 22: 5.18716,
## 21:18.14762
## -------|------|--------------------------------------------------------
不仅没有回报，使用几个多达十个节点的隐藏层的大胆尝试产生了一个最佳估计模型，该模型只包括一个只有八个节点的层。尽管在这种情况下感觉令人失望，但神经网络领域是机器学习的一个迷人和不断增长的领域。
多层感知器的介绍到此结束。它还总结了我们的监督机器学习方法的介绍。虽然在这一章中，我们只探讨了少数几个模型，我们已经探讨了几种不同的模型类型。最重要的是，通过使用caret，而不是每次练习一个新的结构，你有能力轻松地尝试几十个可能的模型。请记住，每种模型类型都倾向于拥有最适合的特定类型的数据，并且在处理或不处理分类或数字数据的能力方面可能有不同的标准。此外，一些模型对主成分分析有要求或限制。最后，记住模型可能如何处理缺失数据是很重要的。有些方法需要完整的数据，在这种情况下，关于缺失数据的章节可能会有所帮助。其他方法对缺失数据更稳健。
最后，我们委婉地提醒我们的读者，在确定哪种方法是最佳的之后，正确的最后一步是根据所有可用的数据对模型进行最后一次训练——尤其是在可用数据不多的情况下。对于模型调整，使用交叉验证来阻止过度拟合仍然是有帮助的。
8.3 总结
本章介绍并探索了作为多个模型的公共接口的caret包。通过这种权宜之计，检查最佳模型拟合的多种算法变得更加简单。此外，还研究了为建模做准备的数据预处理方法。最后，虽然没有明确编码，但是讨论了将原始数据分为训练集、验证集和测试集的原则。特别关注的概念是，所有模型都根据训练数据进行训练，使用验证数据选择最终的最佳模型，并通过原始的、未使用的测试数据对最终模型的真实性能进行评估。本章使用的功能汇总在表 8-2 中。
表 8-2
本章中描述的关键功能列表及其功能摘要
|功能
|
它的作用
|
| — | — |
| boxplot() | 生成箱线图，允许目视检查异常值和其他畸形数据。 |
| cbind() | 通过列将数据绑定在一起。 |
| colnames() | 列出数据集中列标题的名称。 |
| confusionMatrix() | 显示分类预测的实际输出与预测输出的实例。 |
| copy() | 创建数据的完整副本，而不是使用指针。这可以防止对副本的编辑级联到原始数据集。它确实能加倍记忆。 |
| cor() | 计算相关性。 |
| createDataPartition() | Caret 软件包功能，创建注意数据分段的训练/测试数据。在这种情况下，这一点很重要，因为数据是区域性的。 |
| dummyVars() | 一个热点编码回归分类数据。 |
| expand.grid() | 使用所有可能组合的所有输入来创建数据框。 |
| explain() | 在使用其他 DALEX 方法之前需要 DALEX 函数。 |
| fancyRpartPlot() | 更丰富多彩的树形图以及额外的有用数据。 |
| gbm | 插入符号的梯度增强机器形式。 |
| hist() | 绘制直方图。 |
| lapply() | 列表应用功能。 |
| mlpML | 多层感知器(多层)形式插入符号。 |
| model_performance() | 创建可绘制剩余对象的 DALEX 函数。 |
| par() | 可用于通过 R 的默认图形包设置图形参数(如用于绘制图形的列或行)。 |
| pca() | pcaMethods 软件包的主成分分析。 |
| plot_grid() | Cowplot 函数用于设置绘制多个图形的列和行—与 base R 的图形包函数 par()形成对比。 |
| prcomp() | 主成分分析从基地 R 统计软件包。 |
| predict() | 给定一个模型，从给定的输入预测输出。 |
| qqnorm() | 绘制 Q-Q 图。 |
| ranger | 插入符号的随机森林形式。 |
| rpart | 脱字符号的分类和回归树形式。 |
| sapply() | 简化应用功能。 |
| scale() | 基数 R 函数缩放和居中数据。 |
| shapiro.test() | 计算夏皮罗-威尔基的常态标准。 |
| stats::residuals() | 根据模型计算残差(base R stats 包版本)。 |
| summary() | 尽可能提供给定输入的摘要。 |
| svmLinear | 插入符号的支持向量机(线性)形式。 |
| svmPoly | 插入符号的支持向量机(多项式)形式。 |
| train() | 通过给定的训练数据训练指定模型的插入符号函数。 |
| trainControl() | 为所提供的培训方法提供输入的形式(和功能)(例如，交叉验证和折叠次数)。 |
| variable_importance() | 计算压差损失的 DALEX 函数。 |
| variable_response() | 计算部分相关性的 DALEX 函数。 |
| View() | 在新的查看窗口中打开数据框/表进行交互式查看。 |
扫一扫
160
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
