R-机器学习快速启动指南-二-R 机器学习快速启动指南（二）第三章：预测银行失败 - 描述性分析在本章中，我们将学习

R 机器学习快速启动指南（二）

原文：annas-archive.org/md5/13dee0bdb6445b090ed411f424dc82f4

译者：飞龙

协议：CC BY-NC-SA 4.0

第三章：预测银行失败 - 描述性分析

在本章中，我们将学习如何理解并准备我们的银行数据集以进行模型开发。我们将回答有关我们有多少变量以及它们的质量的问题。描述性分析对于理解我们的数据和分析信息质量可能存在的问题至关重要。我们将看到如何处理缺失值，将变量转换为不同的格式，以及如何分割我们的数据以训练和验证我们的预测模型。

具体来说，我们将涵盖以下主题：

数据概览
格式转换
抽样
处理缺失值和异常值
实施描述性分析

数据概览

首先，我们将分析数据集中变量的类型。为此，我们可以使用 class 函数，它告诉我们一个变量是数字、字符还是矩阵。例如，银行识别号 ID_RSSD 的类别可以这样获得：

class(Model_database$ID_RSSD)

 ## [1] "integer"

这个函数表明这个变量是一个没有小数的数字。

我们可以使用以下代码计算所有变量的相同信息并将其存储起来：

 classes<-as.data.frame(sapply(Model_database, class))
 classes<-cbind(colnames(Model_database),classes)
 colnames(classes)<-c("variable","class")

使用 sapply，迭代地对数据集上的 class 函数进行计算。然后，将变量的名称与仅在一个数据框中的类别结合起来，最后重命名生成的数据集：

head(classes)

 ##          variable   class
 ## ID_RSSD   ID_RSSD integer
 ## UBPR1795 UBPR1795 numeric
 ## UBPR4635 UBPR4635 numeric
 ## UBPRC233 UBPRC233 numeric
 ## UBPRD582 UBPRD582 numeric
 ## UBPRE386 UBPRE386 numeric

这个数据集包含四种不同类型的变量：

table(classes$class)

 ## character      Date   integer   numeric
 ##       462         1         4      1027

根据之前的步骤，我们知道只有具有 Date 格式的变量收集财务报表的日期。

熟悉我们的变量

一些变量的字符类型尚不明确。我们的数据集属于金融领域，并且只有金融比率作为数据，所以我们预计变量类型应该是整数或数值。让我们找出我们是否正确。

过滤 character 变量：

classes<-classes[classes$class=="character",]

head(classes)

 ##          variable     class
 ## UBPRE543 UBPRE543 character
 ## UBPRE586 UBPRE586 character
 ## UBPRE587 UBPRE587 character
 ## UBPRE594 UBPRE594 character
 ## UBPRFB64 UBPRFB64 character
 ## UBPRFB69 UBPRFB69 character

第一个变量，UBPRE543，衡量提供建筑贷款的银行的损失总额与授予的建筑贷款总额之比。正如我们所怀疑的，这个变量应该是数字、百分比或小数。

查找变量的缺失值

我们将使用以下代码计算这个变量，UBPRE543，随时间变化的缺失值数量，目的是了解这个变量的一些更多信息：

aggregate(UBPRE543 ~ Date, data=Model_database, function(x) {sum(is.na(x))}, na.action = NULL)

 ##          Date UBPRE543
 ## 1  2002-12-31     1127
 ## 2  2003-12-31      954
 ## 3  2004-12-31      772
 ## 4  2005-12-31      732
 ## 5  2006-12-31      639
 ## 6  2007-12-31      309
 ## 7  2008-12-31      110
 ## 8  2009-12-31       98
 ## 9  2010-12-31       91
 ## 10 2011-12-31       76
 ## 11 2012-12-31      132
 ## 12 2013-12-31       98
 ## 13 2014-12-31       85
 ## 14 2015-12-31       89
 ## 15 2016-12-31       68

如我们所见，这个比率在 2002 年到 2006 年之间显示了一些缺失值。

另一方面，我们可以使用 table 函数计算数据集中按年份的观测数：

table(Model_database$Date)

 ##
 ## 2002-12-31 2003-12-31 2004-12-31 2005-12-31 2006-12-31 2007-12-31
 ##       1127        954        772        732        639        652
 ## 2008-12-31 2009-12-31 2010-12-31 2011-12-31 2012-12-31 2013-12-31
 ##        686        671        587        533        664        615
 ## 2014-12-31 2015-12-31 2016-12-31
 ##        526        498        474

比较前两个表格，我们可以看到这个变量在最初几年没有提供信息。

在练习开始时，我们将 .txt 文件上传到 R 中时，由于这个变量在最初几年没有提供信息，R 自动将这个变量的格式分配为字符格式。

对于后来的年份，当变量有信息时，变量被读取为数值，但在将所有年份合并到数据框中时格式发生了变化，就在我们执行此代码时：

#database<-rbind(year2002,year2003,year2004,year2005,year2006,year2007,year2008, year2009,year2010,year2011,year2012,year2013,year2014,year2015,year2016)

在rbind函数中使用的第一个表格中变量的格式，属于 2002 年，固定并决定了合并后的结果表格的格式。

转换变量的格式

我们现在需要将这些变量全部转换为数值格式。从第二个变量（第一个是标识符）到数据框中的其余变量，变量将被显式转换为数值。最后两个变量也将被排除在此过程之外（Date和目标变量）。让我们使用以下代码将变量转换为数值格式：

for (k in 2:(ncol(Model_database)-2)) 
   {
    Model_database[,k]<-as.numeric(Model_database[,k])
   }

让我们看看这些更改是否已经应用：

table(sapply(Model_database, class))

 ##
 ##    Date integer numeric
 ##       1       1    1492

在继续开发之前，一旦我们解决了数据格式的问题，在下一节中，我们将指定样本的哪一部分用于开发，哪一部分用于验证模型。

采样

所有后续步骤和描述性分析都只考虑训练或开发样本。因此，我们的数据将被分为两个样本：

训练集：通常代表数据的 70%，用于训练模型（选择更适合模型的参数）。
验证集：通常代表数据的 30%，用于衡量模型在做出预测时的表现。

样本分区

尽管有众多方法可以实现数据分区，但caTools包是最有用的之一。此包包含一个名为sample.split的函数，该函数生成随机数以分割样本，但同时也保持原始数据集中不良和良好的比例在分离的样本中。

由于caTools包使用随机数，固定一个seed可以方便地保证结果的复现性：

set.seed(1234)

然后，使用sample.split函数：

library(caTools)
 index = sample.split(Model_database$Default, SplitRatio = .70)

此函数接受两个参数，目标变量和分区大小，在我们的案例中，是 70%。

它生成一个index，包含两个值，TRUE和FALSE，可以用来将数据集分割成两个所需的样本：

train<-subset(Model_database, index == TRUE)
test<-subset(Model_database, index == FALSE)

检查样本

让我们检查每个样本中观察值的数量以及失败银行的占比：

print("The development sample contains the following number of observations:")

 ## [1] "The development sample contains the following number of observations:"

nrow(train)

 ## [1] 7091

print("The average number of failed banks in the sample is:")

 ## [1] "The average number of failed banks in the sample is:"

(sum(train$Default)/nrow(train))

 ## [1] 0.04696094

print("The validation sample contains the following number of observations:")

 ## [1] "The validation sample contains the following number of observations:"

nrow(test)

 ## [1] 3039

print("The average number of failed banks in the sample is:")

 ## [1] "The average number of failed banks in the sample is:"

(sum(test$Default)/nrow(test))

 ## [1] 0.04705495

如所示，训练样本和测试样本分别占总样本的 70%和 30%。这两个样本中失败银行的比率保持大致相同，即 4.7%。

实施描述性分析

描述性统计分析有助于您正确理解您的数据。尽管 R 默认提供了一些函数来执行基本统计，但我们将使用两个更好的替代方案，即DataExplorer和fBasics包。

按照以下简单步骤进行：

由于数据集中变量的数量很高，我们将创建一个包含要用于描述性函数的变量名的列表：

Class<-as.data.frame(sapply(train, class))
 colnames(Class)<-"variable_class"
 Class$variable_name<-colnames(train)

numeric_vars<-Class[Class$variable_class=="numeric","variable_name"]

创建了一个包含 1,492 个变量的列表。将此列表传递给fBasics包中包含的basicStats函数：

library(fBasics)
 descriptives_num<-             as.data.frame(t(basicStats(train[,numeric_vars])))
 head(descriptives_num)

我们可以计算以下描述性统计量：

- 观察值数量（nobs）
- 缺失值的数量 (NAs)
- 最小值 (Minimum)
- 最大值 (Maximum)
- 第一和第三四分位数 (1. Quartile 和 3. Quartile)
- 中位数 (Median)
- 变量中值的总和 (Sum)
- 均值的标准误差 (SE Mean)
- 均值下置信限 (LCL Mean)
- 均值上置信限 (UCL Mean)
- 方差 (Variance)
- 标准差 (Stdev)
- 偏度 (Skewness)
- 峰度 (Kurtosis)

在这一步，我们将检测具有大量缺失值的变量，变量的范围和离散度，即使变量只有一个唯一值。

当变量的数量很高，如我们的案例，这项任务并不容易，我们需要一些时间来分析变量。变量的图形分析也很重要，并且是补充性的。

plot_histogram 函数对于可视化变量非常有用。此函数在 DataExplorer 包中可用：

library(DataExplorer)
plot_histogram(train[,1410:1441])

以下图表显示了前面代码的输出。这些图表显示了数据中一些变量的直方图。以下是输出结果的第一页：

这里是输出结果的第二页：

这种对变量分布的分析不仅需要理解变量的分布，还需要检测潜在的问题。

处理异常值

一个重要的问题是检测数据中的异常值。异常值是看起来与一组观察值不同的值。考虑一个正态分布的例子，其中分布尾部的值可以被认为是异常值。它们与样本中最近的值关系并不紧密。

有些算法对异常值非常敏感，因此其处理不是一个简单的问题。如果变量的数量较少，检测异常值会更容易。

Winsorization 过程

当异常值的数量很高时，我们需要使用自动程序来帮助自动检测它们。避免异常值问题的最有效方法之一是Winsorization 过程。

根据这种方法，异常值将被替换为固定值。如果一个变量的值小于特定的阈值，这个值将被替换为这个极限。对于变量中的高值，情况也是相同的。

理想情况下，这些极限或阈值基于百分位数。例如，在较低范围内选择 1、2.5 或 5 的百分位数，在较高范围内选择 95、97.5 和 99 的百分位数，可以用于 Winsorization 技术，尽管可以选择其他方法，例如使用四分位距。

实施 Winsorization

让我们将 Winsorization 方法付诸实践。首先，我们需要知道数据集中比率的位置：

head(colnames(train))

 ## [1] "ID_RSSD"  "UBPR1795" "UBPR4635" "UBPRC233" "UBPRD582" "UBPRE386"

tail(colnames(train))

 ## [1] "UBPRE541" "UBPRE542" "UBPRJ248" "UBPRK447" "Date"     "Default"

因此，我们需要将技术应用于所有变量，除了第一个和最后两个变量。

在训练集中完成的全部转换都应该在测试数据集中后续应用。测试样本中的修改将使用训练数据的限制来完成。我们将对两个数据集都进行 winsorization：

for (k in 2:(ncol(train)-2))
{
   variable<-as.character(colnames(train)[k])
   limits <- quantile(train[,k], probs=c(.01, .99), na.rm = TRUE)
   train[complete.cases(train[,k]) & train[,k] <         as.numeric(limits[1]),k] <-      as.numeric(limits[1])
   train[complete.cases(train[,k]) & train[,k] > as.numeric(limits[2]),k] <-      as.numeric(limits[2])
 test[complete.cases(test[,k]) & test[,k] < as.numeric(limits[1]),k]    <- as.numeric(limits[1])
 test[complete.cases(test[,k]) & test[,k] > as.numeric(limits[2]),k] <-as.numeric(limits[2])
 }

对于每个变量，这将计算训练集中的第一和第九十九百分位数。然后，替换超过第九十九百分位数值的异常值，或小于第一百分位数的值，用相应的百分位数值替换。这意味着它为第一和第九十九百分位数中固定的每个值建立了最大和最小值。这个程序对训练和测试样本都进行。

区分单值变量

现在，我们将计算一个变量所取的唯一值的数量。因此，如果一个变量只取一个单一值，它可以直接从数据集中删除。

sapply函数允许计算每个变量的n_distinct值。创建一个新的数据框，包含常量变量的名称：

library(dplyr)

unique_values<-as.data.frame(sapply(train, n_distinct))

在这个数据框中重命名变量的名称：

colnames(unique_values)<-"Unique_values"

在数据框中添加一个包含变量名称的列：

unique_values$variable_name<-colnames(train)

然后创建一个包含常量变量名称的列表：

variables_to_remove<-unique_values[unique_values$Unique_values==1,"variable_name"]
length(variables_to_remove)

 ## [1] 84

只有84个变量具有唯一的唯一值。这些变量将在train和test样本中删除：

 train<-train[, !colnames(train) %in% variables_to_remove]
 test<-test[, !colnames(test) %in% variables_to_remove]

winsorization 的一个问题是，如果一个变量显示的不同值数量很少，它可以只用一个值来替换所有值。这是因为一个变量可能在几个百分位数级别上取相同的值。了解每个程序及其对发展的影响的优缺点非常重要。

记得保存你的工作空间：

save.image("Data6.RData")

处理缺失信息

大多数算法在数据包含缺失值或自动处理这些值的预定动作时都会失败。在这种情况下，掌握控制权非常重要。

处理缺失信息的两种最常见的方法是：删除具有缺失值的观测值或用具体值（通常是中位数或平均值）替换它们。当进行值插补时，你可能会丢失重要信息。例如，变量的缺失值可能总是出现在目标变量的一个类别中。一个典型的例子是我们试图预测银行贷款的合格和不合格申请人的模型。

通常会有与过去某些支付问题相关的变量。有时，根据数据集的不同，这些变量显示缺失值仅仅是因为申请人没有先前的问题。这种情况通常发生在值插补可能导致我们丢失相关信息时。

当变量数量较低时，对缺失值的详细分析最为常见。如果变量数量较高，一些自动替代方法可能更有效。

在对缺失值采取行动之前，让我们通过分析列和行来找出缺失值的数量。然后您可以删除具有大量缺失值的变量和行（也称为观测值）。

在任何情况下，删除变量或观测值的阈值是主观的，并且取决于它们应用的特定案例。

使用DataExplorer包，我们可以找到我们数据中缺失值的百分比。让我们用少数几个变量试一试：

plot_missing(train[,c(6:8,1000:1020)])

上一行代码将打印出类似此图的图表：

上述图表仅是使用DataExplorer包表示某些变量缺失值的一个示例。此包还根据缺失值的数量提供有关使用变量的有用性建议。

我们还有另一种确定具有缺失值的变量的方法。这是一个替代方案，如果您无法访问DataExplorer包，或者只是不想使用它。这更是一种手动方式。让我们编写代码：

ncol=rep(nrow(train) ,each=ncol(train))
 missingdata=as.data.frame(cbind(colnames=names(train),ncol,nmsg=as.integer(as.character(as.vector(apply(train, 2, function(x) length(which(is.na(x)))))))))
missingdata$nmsg=as.numeric(levels(missingdata$nmsg))[missingdata$nmsg]
missingdata=cbind(missingdata,percmissing=(missingdata$nmsg/ncol*100))

head(missingdata)

 ##   colnames ncol nmsg percmissing
 ## 1  ID_RSSD 7091    0           0
 ## 2 UBPR1795 7091    0           0
 ## 3 UBPR4635 7091    0           0
 ## 4 UBPRC233 7091    0           0
 ## 5 UBPRD582 7091    0           0
 ## 6 UBPRE386 7091    0           0

例如，我们可以检查缺失值占总观测值超过 99%的变量（这里只显示了少数几行）：

print(missingdata[missingdata$percmissing>=99,])

 ##      colnames ncol nmsg percmissing
 ## 19   UBPRE406 7091 7066    99.64744
 ## 26   UBPRE413 7091 7028    99.11155
 ## 35   UBPRFB69 7091 7038    99.25257
 ## 121  UBPRE137 7091 7048    99.39360
 ## 161  UBPRE184 7091 7046    99.36539
 ## 1347 UBPRE855 7091 7073    99.74616
 ## 1348 UBPRE856 7091 7047    99.37950
 ## 1356 UBPRE864 7091 7083    99.88718
 ## 1360 UBPRE868 7091 7056    99.50642

在这种情况下，我更倾向于不删除任何变量，考虑到其缺失值的数量。没有空变量：

print(missingdata[missingdata$percmissing==100,])

 ## [1] colnames    ncol        nmsg        percmissing
 ## <0 rows> (or 0-length row.names)

让我们通过分析行来计算缺失值的数量：

train$missingvalues<-rowSums(is.na(train[,2:1410]))/1409

现在绘制一个直方图来图形化描述缺失值的分布：

hist(train$missingvalues,main="Distribution of missing values",xlab="Percentage of missing values",border="blue", col="red",breaks=25)

上一段代码生成了以下图表：

可以使用以下代码获取银行缺失值百分比的摘要：

summary(train$missingvalues)

 ##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
 ## 0.06671 0.44003 0.47480 0.46779 0.50958 0.74663

虽然某些银行的缺失值数量很高，但建议现在不要删除任何观测值，而只删除最近创建的缺失变量：

train$missingvalues<-NULL

我们将使用一个有趣的包来可视化数据集中缺失值的数量，这个包是Amelia包。这是一个用于多重插补缺失数据的包，还包括一个图形界面。

让我们看看一些变量的例子：

library(Amelia)

missmap(train[,5:35], main = "Missing values vs observed",col=c("black", "grey"),,legend = FALSE)

上一段代码生成了以下图表：

虽然表示方式不太美观，但这个图表显示了x轴上的某些变量和数据集中的y轴上的观测值。黑色点表示数据集中存在缺失变量。

如所示，一些变量显示有大量的缺失值。

分析缺失值

我们已经提到，了解缺失值的来源以及它们是否可以提供一些信息是很重要的。在下面的例子中，我们将分析UBPRE628变量中呈现的缺失值。这个变量衡量的是一家银行的长期债务总额除以总银行股本资本。银行的资本很重要，因为在运营中面临损失的情况下，银行将使用这部分资本来吸收损失并避免未来的破产。资本越高，银行面对经济问题的缓冲能力就越强。

通常，与银行资本相关的债务比例越高，银行在未来如果发生新危机的情况下可能遇到的问题就越多。在新的金融危机发生时，银行可能无法偿还其债务，即使通过出售其资产。

根据我们的分析，这个变量显示了高比例的缺失值，具体来说，这个比率对于我们数据中的 23.97%的银行没有信息：

missingdata[missingdata$colnames=="UBPRE628",]

 ##     colnames ncol nmsg percmissing
 ## 281 UBPRE628 7091   17   0.2397405

理解结果

现在我们将创建一个辅助数据框来统计失败银行的数目并检查这个比率是否已告知：

missing_analysis<-train[,c("UBPRE628","Default")]

现在我们将创建一个标志来检查变量是否存在缺失：

missing_analysis$is_miss<-ifelse(is.na(missing_analysis$UBPRE628),"missing_ratio","complete_ratio")

最后，让我们总结一下数据集中两种情况下的现有违约数量：这个比率中缺失值的存在或缺失：

aggregate(missing_analysis$Default, by = list(missing_analysis$is_miss), sum)

 ##          Group.1   x
 ## 1 complete_ratio 319
 ## 2  missing_ratio  14

根据这个表格，只有14家失败的银行在这个比率中显示了缺失值。显然，我们可以从这个结论中得出，一家银行可能会故意不报告特定的比率，因为计算出的比率可能会让其他人意识到这家银行的糟糕经济状况。在这种情况下，如果观察到缺失值，我们不会观察到高比例的坏银行。

缺失值将通过计算训练数据集中非缺失观测值的比率平均值来估计。这意味着，如果验证数据集中存在缺失值，它们也可能存在于训练数据集中。让我们看一个例子：

train_nomiss<-train
test_nomiss<-test

 for(i in 2:(ncol(train_nomiss)-2))
   {
   train_nomiss[is.na(train_nomiss[,i]), i] <- mean(train_nomiss[,i],      na.rm =          TRUE)
   test_nomiss[is.na(test_nomiss[,i]), i] <- mean(train_nomiss[,i],      na.rm = TRUE) 
   }

我们可以使用Amelia包在训练和验证样本上检查这个过程是否成功（可能需要几分钟）。例如，你可以在过程执行后检查训练样本中是否存在缺失值：

missmap(train_nomiss[,2:(ncol(train_nomiss)-2)], main = "Missing values vs observed",col=c("black", "grey"),,legend = FALSE)

这里是前面代码的输出：

对测试样本执行相同的检查：

missmap(test_nomiss[,2:(ncol(train_nomiss)-2)], main = "Missing values vs observed",col=c("black", "grey"),,legend = FALSE)

再次，显示了新的输出：

两个地图以灰色绘制，表示没有缺失值。

现在我们将备份我们的工作空间并删除所有不必要的表格：

rm(list=setdiff(ls(), c("Model_database","train","test","train_nomiss","test_nomiss")))
save.image("Data7.RData")

摘要

在本章中，我们了解了一些准备和了解数据的基本重要步骤。在我们的数据集中有多少个变量可用？我们有什么样的信息？数据中是否有缺失值？我该如何处理缺失值和异常值？我希望你现在可以回答这些问题。

此外，在本章中，我们还学习了如何将我们的数据分割成训练集和验证集，以训练和验证我们即将到来的预测模型。在下一章，我们将更进一步，对这份数据进行单变量分析，这意味着分析变量是否对预测银行破产有用。

第四章：预测银行失败 - 单变量分析

近年来，大数据和机器学习在许多领域变得越来越受欢迎。人们普遍认为，变量越多，分类器就越准确。然而，这并不总是正确的。

在本章中，我们将通过分析每个变量的单独预测能力以及使用不同的替代方案来减少数据集中的变量数量。

在本章中，我们将涵盖以下主题：

特征选择算法
过滤方法
包装方法
内嵌方法
维度降低

特征选择算法

在预测银行失败的这一实际案例中，我们有许多变量或财务比率来训练分类器，因此我们预计会得到一个很好的预测模型。考虑到这一点，我们为什么要选择替代变量并减少它们的数量呢？

嗯，在某些情况下，通过添加新特征来增加问题的维度可能会降低我们模型的性能。这被称为维度诅咒问题。

根据这个问题，增加更多特征或增加特征空间的维度将需要收集更多数据。从这个意义上说，我们需要收集的新观察结果必须以指数速度快速增长，以维持学习过程并避免过拟合。

这个问题通常在变量数量与我们的数据中的观察数量之间的比率不是很高的情况下观察到。

特征选择对于从数据中识别和删除不必要的、不相关的和冗余变量，以及降低模型复杂性也是很有用的。

特征选择类别

在许多机器学习指南中可以找到三种特征选择算法的通用类别。这些包括以下内容：

过滤方法：在这个方法中，变量是根据与目标变量的相关性来选择的。因此，它是通过每个变量解释模型目标的能力来衡量的。这些方法在变量数量高时特别有用，并且有助于避免过拟合。作为缺点，值得提到的是，尽管这些变量在单独的情况下不是预测性的，并且以单变量的方式衡量，但当与其他变量结合时，它们可以成为预测性的。总之，这些方法不考虑变量之间的关系。
包装方法：包装方法通过评估变量的子集来检测变量之间的可能相互作用。在包装方法中，在预测模型中使用多个变量的组合，并根据模型的准确性为每个组合给出一个分数。因此，可以避免不相关的组合。然而，如果变量的数量很大，这些方法会非常耗时，并且始终存在过拟合的风险，尤其是在观察数量较低的情况下。
嵌入式方法：最后，通过嵌入式方法，在训练过程中学习并记住对提高模型准确性最有贡献的特征。正则化就是这样一种特征选择方法。它们也被称为惩罚方法，因为它们对优化参数施加约束，通常使模型变量数量减少。Lasso、弹性网络和岭回归是最常见的正则化方法。其他嵌入式特征选择算法的例子包括 Lasso、弹性网络和岭回归算法。我们将在稍后更详细地研究这些模型。

过滤方法

让我们从一种过滤方法开始，以减少第一步中的变量数量。为此，我们将测量一个变量的预测能力或其单独和正确分类目标变量的能力。

在这种情况下，我们试图找到能够正确区分清偿能力和非清偿能力银行的变量。为了衡量一个变量的预测能力，我们使用一个名为信息值（IV）的指标。

具体来说，给定一个分为n组的分组变量，每个组都有一定数量的良好银行和不良银行——或者在我们的案例中，清偿能力和非清偿能力银行——该预测器的信息值可以按以下方式计算：

IV 统计量通常根据其值进行解释：

< 0.02：分析变量不能准确区分目标变量的类别
0.02 到 0.1：变量与目标变量有较弱的关系
0.1 到 0.3：变量显示出中等强度的关系
> 0.3：变量是目标的好预测器

根据这个值，这个变量本身具有很强的预测性。因此，这个变量可以用于我们的模型。让我们看看一个变量IV的计算示例：

在先前的表中，我们可以看到UBPRE006变量的信息值计算，它代表贷款和损失预留总额除以银行的资产总额。

从广义上讲，当贷款发放时，其中一部分必须预留以备信用违约的情况；也就是说，银行在其损益表中做出两种类型的预留以覆盖所谓的信用风险：一种是贷款发放时做出的通用预留，另一种是针对未偿还信用的特定预留。

理论上，比率越高，银行破产的可能性就越大，因为如果预留水平高，这表明其贷款的信用质量将较低。

记住，在我们的样本中，失败银行的百分比为 4.70%。在这个例子中，UBPRE006 变量已被分为四个类别，以及一个额外的类别来衡量缺失值的水平。这可以在 BadRate 列中看到，作为失败银行的比率，其值低于 0.5487%。这非常低，仅代表该组中 0.80% 的银行。随着这个比率的增加，失败银行的比率也会更高。此外，在这个比率中没有银行有缺失值。

此表第一组中出现的值是根据此方程计算的：

在此表格的 IV 列中所有值的总和可以在第 6 行的 3.2803 列中找到。

根据这个值，这个变量本身具有很强的预测性。因此，在模型中使用这个变量可能是有用的。

另一方面，证据权重（WoE）是一个与信息值非常密切相关的指标。此指标也包括在先前的表格中。WoE 的计算方法如下：

实际上，WoE 方程是 IV 指标的一部分。如果 良好到不良 比率的优势为 1，则 WoE 的值将为 0。

如果一个组中不良的百分比大于良好的百分比，则优势比将小于 1，WoE 将是负数；如果组中良好的数量大于不良，则 WoE 值将是正数。

通常，WoE 的正值表示被放入该组的银行比样本中所有银行的平均状况更稳健。另一方面，负值越高，该组中的银行风险越大。

我们将计算训练集中每个变量的信息值。smbinning 包对此非常有用。

如已所见，一个重要步骤是分组变量，这个包会自动完成。

我们将进行两个不同的实验。因此，在缺失值插补前后，我们将计算训练集的信息值。我们将在本节后面讨论这些实验背后的原因。

此包假设如果一家银行稳健，则目标变量应取值为 1，否则为 0，这与我们在之前步骤中所做的正好相反。因此，第一步包括反转目标值：

aux_original<-train
aux_original$Defaultf<-as.numeric(as.character(aux_original$Default))
aux_original$Defaultf<-ifelse(aux_original$Default==1,0,1)
aux_nomiss<-train_nomiss
aux_nomiss$Defaultf<-as.numeric(as.character(aux_nomiss$Default))
aux_nomiss$Defaultf<-ifelse(aux_nomiss$Default==1,0,1)

接下来，我们运行以下代码（这是一个非常耗时的过程）：

library(smbinning)
table_iv<-matrix("NA",0,5)
table_iv<-data.frame(table_iv)
colnames(table_iv)<-c("Char","IV_original","Process_original","IV_nomiss","Process_nomiss")

 for (var in 1:length(aux_original[,2:1408]))
 {
 variable<-colnames(aux_original)[var+1]
 aux_original2<-aux_original[,c(variable,"Defaultf")]
 aux_nomiss2<-aux_nomiss[,c(variable,"Defaultf")]
 temp1<-smbinning.sumiv(aux_original2, "Defaultf")
 temp2<-smbinning.sumiv(aux_nomiss2, "Defaultf")
 colnames(temp1)<-c("Char","IV_original","Process_original")
 colnames(temp2)<-c("Char","IV_nomiss","Process_nomiss")
 temp2$Char<-NULL
 temp1<-cbind(temp1,temp2)
 table_iv<-rbind(table_iv,temp1)
 }

之前的代码创建了一个表格，其中将存储信息值（table_iv）。然后，对于 train 数据集中的每个变量，使用 smbinning.sumiv 函数计算信息值。

一旦过程完成，就会创建工作区备份：

save.image("Data8.RData")

让我们看看结果：

head(table_iv)
 ## Char IV_original      Process_original IV_nomiss
 ## 1 UBPR1795      2.6138    Numeric binning OK    2.6138
 ## 2 UBPR4635      2.5253    Numeric binning OK    2.5253
 ## 3 UBPRC233          NA No significant splits        NA
 ## 4 UBPRD582          NA    Uniques values < 5        NA
 ## 5 UBPRE386          NA No significant splits        NA
 ## 6 UBPRE388      0.5853    Numeric binning OK    0.5622
 ##          Process_nomiss
 ## 1    Numeric binning OK
 ## 2    Numeric binning OK
 ## 3 No significant splits
 ## 4    Uniques values < 5
 ## 5 No significant splits
 ## 6    Numeric binning OK

在这个表中，我们展示了在缺失值插补前后训练样本中每个变量的信息值。

在信息值之后，一列告诉我们计算状态。这里可能显示不同的消息如下：

数值分箱正常：信息值已正确计算，并且至少可以在变量中区分出两个不同的组，以区分好银行和坏银行。
无显著分割：由于变量无法区分银行的未来偿债能力，因此不计算信息值。未检测到银行的不同组。
唯一值少于 5：一个变量具有少于五个不同的值。在这种情况下，不计算信息值。

结果取决于缺失值之前是否已处理或是否包含在示例中：

table(table_iv$Process_original)
 ##
 ##    Numeric binning OK No significant splits    Uniques values < 5
 ##                   522                   807                    78
table(table_iv$Process_nomiss)
 ##
 ##    Numeric binning OK No significant splits    Uniques values < 5
 ##                   539                   790                    78

smbinning包将缺失值视为一个新类别，在某些情况下，这可能包含相关信息或与我们的目标变量中的一个类别更相关。

重要的是要意识到缺失变量的插补意味着对您的数据的一种限制。

让我们检查这个简单决策对预测能力的影响。计算如下，存在或不存在缺失值的IV差异：

diff_iv<-table_iv[complete.cases(table_iv) & table_iv$Process_original=="Numeric binning OK" &table_iv$Process_nomiss=="Numeric binning OK" ,]

diff_iv$diff<-(diff_iv$IV_nomiss - diff_iv$IV_original)
hist(diff_iv$diff, border=TRUE , col=rgb(0.8,0.2,0.8,0.7) , main="")

让我们看看输出结果：

这里是总结统计信息：

summary(diff_iv$diff)
 ##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
 ## -2.64110 -0.03052  0.00000 -0.05247  0.00000  0.22570

根据总结统计信息，平均而言，缺失值插补会降低变量的预测能力或信息值 5.247%。在某些情况下，变量甚至会增加预测能力。

我们可以根据其信息值对变量的预测能力进行分类。我们将考虑之前解释的阈值来定义变量是否显示强大的、中等的或弱的预测能力：

table_iv$IV_Category<-ifelse(table_iv$IV_nomiss >= 0.3, "1:Strong", ifelse(table_iv$IV_nomiss >= 0.1, "2:Medium","3:Weak"))

table(table_iv$IV_Category)
 ##
 ## 1:Strong 2:Medium   3:Weak
 ##      358      114       67

在这一步，我们可以移除预测能力低的变量。因此，从数据集中我们拥有的超过一千个变量中，我们只会选择被分类为强或中的变量：

table_iv<-table_iv[complete.cases(table_iv) & table_iv$IV_Category != "3:Weak",]

尽管存在限制，我们仍将使用之前已处理缺失值的数据库：

train<-train_nomiss
test<-test_nomiss

通常，对可能作为多元模型一部分的变量进行单变量转换可以提高其区分能力。

变量转换有助于减轻异常值对模型开发的影响，并捕捉变量与目标之间的非线性关系。

信用风险中最常见的做法之一是将变量转换为不同的类别，然后为每个类别分配其对应的 WoE 值。让我们看看这个例子。再次，使用smbinning包。

对于这一点，我们首先需要将目标值转换为由于包限制而具有的相反值：

train$Defaultf<-as.numeric(as.character(train$Default))
train$Defaultf<-ifelse(train$Default==1,0,1)

test$Defaultf<-as.numeric(as.character(test$Default))
test$Defaultf<-ifelse(test$Default==1,0,1)

我们将smbinning函数应用于UBPRD486变量或“一级杠杆资本”。从监管者的角度来看，“一级”比率代表银行财务实力的核心指标。比率越高，银行的偿付能力和实力就越强。

首先，我们分析这个变量在失败和非失败银行中的分布情况：

boxplot(train$UBPRD486~train$Default, horizontal=T, frame=F, col="lightgray",main="Tier One Leverage Ratio Distribution")

这是“一级”比率分布图：

通常情况下，根据之前的截图，失败银行在这个比率中显示的值较低。应用smbinning函数，创建了一个对象，并进行了变量分类。

可以进行以下一些图形分析：

smbinning.plot(result,option="dist")

以下截图很好地描述了分类情况：

在这个变量中，74.6%的银行显示的值高于 8.16。让我们看看按组划分的失败银行百分比：

smbinning.plot(result,option="badrate")

下面的结果图显示，在比率值低于或等于5.6的银行中观察到较高的不良率：

如注释所述，前面的截图显示，比率值越低，失败银行的数目就越高。这个变量非常有预测性，因为它很容易通过收集大量失败银行来找到组。因此，第一组中包含的 74.8%的银行破产了。也有可能通过运行以下代码来绘制每个组的WoE值：

smbinning.plot(result,option="WoE")

前面的代码提供了以下截图，其中显示了证据值的权重：

对于一些信用申请，随着评分模型的发展，将每个比率的原始值替换为对应证据值的权重，这是一种非常常见的做法，如图所示。例如，值低于或等于 5.6 的将被替换为-4.1。因此，WoE 变量随后用于训练模型，使用逻辑回归，这是最常见的方法。

smbinning包还帮助将原始变量转换为相应的组。根据我的经验，我没有发现很多证据表明 WoE 转换真的能提高模型的性能。因此，在这种情况下，我们不会转换我们的变量。

弱变量也被移除了：

relevant_vars<-as.vector(table_iv$Char)
relevant_vars<-c("ID_RSSD","Default","Defaultf", relevant_vars)
train<-train[,relevant_vars]
test<-test[,relevant_vars]

以下是将工作区保存的步骤：

save.image("Data9.RData")

我们可以继续过滤变量。到目前为止，我们数据集的维度如下：

dim(train)
 ## [1] 7091  465

当面对回归或分类问题时，如果移除了高度相关的属性，一些模型的表现会更好。相关性可以通过以下方式获得：

correlations <- cor(train[,4:ncol(train)])

属于 caret 包的 findCorrelation 函数在相关矩阵上执行搜索，并输出包含整数的向量。这些整数对应于如果删除，可以减少成对相关性的列。因此，此函数寻找具有更高相关性的属性。

它通过考虑成对相关列的绝对值来工作。它移除具有最大平均绝对值的变量，这是通过比较高度相关的变量来计算的：

## Loading required package: lattice
highlyCorrelated <- data.frame("Char"=findCorrelation(correlations,      cutoff=0.75,names = TRUE))

截断选项是选择高度相关变量的阈值：


correlated_vars<-as.vector(highlyCorrelated$Char)
non_correlated_vars<-!(colnames(train) %in% correlated_vars)

train<-train[,non_correlated_vars]
test<-test[,non_correlated_vars]

数据集中总共剩下 262 个变量：

ncol(train)
 #262

包装方法

如本节开头所述，包装方法评估变量子集以检测变量之间可能的相互作用，这比过滤方法提前一步。

在包装方法中，在预测模型中使用多个变量的组合，并根据模型精度对每个组合给出分数。

在包装方法中，分类器通过作为黑盒的多变量组合进行迭代训练，其唯一输出是重要特征的排名。

Boruta 包

R 中最知名的包装包之一称为 Boruta。此包主要基于 随机森林 算法。

虽然此算法将在本书的后续部分进行更详细的解释，但由 Breiman 于 2001 年提出的 Boruta 是一种挖掘数据和在样本上生成许多决策树并通过多数投票结合的工具。随机森林创建不同决策树的目的，是为了从不同类别的数据中获得最佳可能的分类。

随机森林的一个成功应用实例是在信用卡欺诈检测系统中。

在 Boruta 包中，使用数据集中其他变量的多个组合创建随机变量。

然后将新变量与原始变量结合，并训练不同的随机森林。通过比较随机变量与原始变量的重要性，获得不同特征的重要性。

只有比随机变量重要性更高的变量才被认为是重要的。如果变量数量很多，Boruta 包将非常耗时，尤其是因为算法会创建更多变量来对其特征进行排名。

让我们在 R 中启动 Boruta 算法。首先，建立一个 seed 以使练习可重复：

set.seed(123)

然后从训练数据集创建一个辅助表，并且没有删除相关变量：

aux<-train
aux$`ID_RSSD`<-NULL
aux$Defaultf<-NULL

最后，启动 Boruta 算法（这非常耗时，可能需要超过一小时）：

library(Boruta)
wrapper <- Boruta(Default ~. , data = aux, doTrace = 2,maxRuns = 100)

当打印 wrapper 对象时，它提供了数据集中特征的重要性。

Boruta 算法对我们数据库中的任何变量得出结论：

print(wrapper)
 ## Boruta performed 99 iterations in 1.15968 hours.
 ##  85 attributes confirmed important: UBPR2150, UBPR7402, UBPRA222,
 ## UBPRD488, UBPRD646 and 80 more;
 ##  139 attributes confirmed unimportant: UBPR0071, UBPR1590,
 ## UBPR1616, UBPR1658, UBPR1661 and 134 more;
 ##  35 tentative attributes left: UBPR2366, UBPR3816, UBPRE083,
 ## UBPRE085, UBPRE140 and 30 more;

许多变量被分类为重要或不重要，但在其他情况下，变量被分配到尝试性类别：

table(wrapper$finalDecision)
 ##
 ## Tentative Confirmed  Rejected
 ##        35        85       139

尝试性特征的重要性几乎与它们最好的随机特征相同。在这种情况下，Boruta无法就默认的随机森林迭代次数做出自信的决定。

在继续之前，让我们先备份工作空间：

save.image("Data10.RData")

通过TentativeRoughFix函数，可以对尝试性变量做出决定。为此，将中值特征 Z 分数与最重要的随机特征的中值 Z 分数进行比较，并做出决定：

wrapper <- TentativeRoughFix(wrapper)
print(wrapper)
 ## Boruta performed 99 iterations in 1.15968 hours.
 ## Tentatives roughfixed over the last 99 iterations.
 ##  108 attributes confirmed important: UBPR2150, UBPR3816, UBPR7402,
 ## UBPRA222, UBPRD488 and 103 more;
 ##  151 attributes confirmed unimportant: UBPR0071, UBPR1590,
 ## UBPR1616, UBPR1658, UBPR1661 and 146 more;

因此，根据此包，我们的训练样本将减少到只有99个变量。

Boruta不是唯一的包装方法。caret包还包括一个包装过滤器。在这种情况下，该算法被称为递归特征消除（RFE）。

在此算法中，首先使用所有独立变量训练一个模型，并计算特征的重要性。将不那么重要的变量（n）从样本中移除，并再次训练模型。这一步骤重复多次，直到所有变量都被使用。在每个迭代中，评估模型的性能。在性能最佳的模型中确定最佳预测变量。

在此算法中，除了随机森林（rfFuncs）之外，还有许多可用于训练的模型，例如以下这些：

线性回归，lmFuncs函数
朴素贝叶斯函数，nbFuncs
带袋的树函数，treebagFuncs

让我们看看这个算法如何在 R 中使用：

首先，固定一个seed以获得相同的结果：

library(caret)
set.seed(1234)

将目标变量转换为factor。因此，算法用于分类。如果不这样做，则假定是一个回归问题：

aux$Default<-as.factor(aux$Default)

最后，运行算法。随机森林被选为具有 10 折验证的分类器（此执行也很耗时）：

rfe_control <- rfeControl(functions=rfFuncs, method='cv', number=10)
recursive <- rfe(aux[,2:260], aux[,1], rfeControl=rfe_control)

如果打印recursive对象，将显示最重要的变量：

print(recursive, top=10)
 ##
 ## Recursive feature selection
 ##
 ## Outer resampling method: Cross-Validated (10 fold)
 ##
 ## Resampling performance over subset size:
 ##
 ##  Variables Accuracy  Kappa AccuracySD KappaSD Selected
 ##          4   0.9848 0.8224   0.005833 0.06490        
 ##          8   0.9866 0.8451   0.004475 0.04881        
 ##         16   0.9884 0.8685   0.005398 0.06002        
 ##        259   0.9886 0.8659   0.004019 0.04617        *
 ##
 ## The top 10 variables (out of 259):
 ##    UBPRD488, UBPRE626, UBPRE217, UBPRE170, UBPRE392, UBPRE636, UBPRE883, UBPRE394, UBPRE370, UBPRE074
plot(recursive, type=c("g", "o"), cex = 1.0)

将获得以下输出：

让我们获取所有变量的排名：

head(predictors(recursive))
 ##   [1] "UBPRD488" "UBPRE626" "UBPRE217" "UBPRE170" "UBPRE392"           "UBPRE636
head(recursive$resample, 10)
 ##    Variables  Accuracy     Kappa .cell1 .cell2 .cell3 .cell4 Resample
 ## 4        259 0.9915374 0.8988203    675      1      5     28   Fold01
 ## 8        259 0.9915374 0.8988203    675      1      5     28   Fold02
 ## 12       259 0.9830748 0.7976887    672      3      9     25   Fold03
 ## 16       259 0.9887165 0.8690976    673      3      5     28   Fold04
 ## 20       259 0.9929577 0.9169746    676      0      5     29   Fold05
 ## 24       259 0.9901269 0.8801237    675      1      6     27   Fold06
 ## 28       259 0.9873061 0.8590115    671      5      4     29   Fold07
 ## 32       259 0.9859155 0.8314180    674      2      8     26   Fold08
 ## 36       259 0.9929478 0.9169536    675      1      4     29   Fold09
 ## 40       259 0.9816384 0.7903799    669      6      7     26   Fold10

如您所见，我们成功获得了变量的排名。尽管如此，用户必须具体选择最终模型中将包含多少最终变量。

在这种情况下，只考虑Boruta包执行后的结果变量：

 predictors<-data.frame("decision"=wrapper$finalDecision)

 predictors<-cbind("variable"=row.names(predictors),predictors)

 predictors<-                    as.vector(predictors[predictors$decision=="Confirmed","variable"])

 train<-train[,c('ID_RSSD','Default',predictors)]

 test<-test[,c('ID_RSSD','Default',predictors)]

我们的样本已经减少：

ncol(train)
 ## [1] 110
save.image("Data11.RData")

嵌入式方法

过滤器和包装方法之间的主要区别在于，在过滤器方法中，例如嵌入式方法，你不能将学习和特征选择部分分开。

正则化方法是嵌入式特征选择方法中最常见的类型。

在此类分类问题中，逻辑回归方法无法处理当变量高度相关时的多重共线性问题。当观测数数量不比协变量数量 p 多很多时，可能会有很多变异性。因此，这种变异性甚至可以通过简单地添加更多参数来增加似然，从而导致过拟合。

如果变量高度相关或存在多重共线性，我们预计模型参数和方差会被夸大。高方差是因为错误指定的模型包含了冗余的预测变量。

为了解决这些局限性，一些方法已经出现：岭回归、Lasso 和弹性网络是最常见的方法。

岭回归

在 Ridge regression 中，回归系数的大小基于 L2 范数进行惩罚：

这里，L(B|y,x) 代表逻辑回归的似然函数，λ 是调整参数，用于控制这两个项对回归系数估计的相对影响。

岭回归的局限性

岭回归将所有预测变量包含在最终模型中。然而，当变量数量 p 很大时，它通常在模型解释上显示出问题。

Lasso

Lasso 代表正则化的另一种选择，并且它克服了岭回归的缺点，减少了最终模型中的预测变量数量。这次，它使用 L1 惩罚来惩罚回归系数的大小：

当 λ 足够大时，它迫使一些系数估计值恰好等于零，从而获得更简洁的模型。

Lasso 的局限性

有时，Lasso 也显示出重要的弱点：如果协变量的数量 p 远远大于观测数的数量，选定的变量数量将受到观测数数量的限制。

弹性网络

Elastic net 尝试克服岭回归和 Lasso 模型的局限性，并在变量高度相关时表现良好。

弹性网络使用所有变量来训练模型，但它也试图结合两种先前使用的方法（岭回归和 Lasso 回归）的优点。因此，弹性网络根据 L1 范数和 L2 范数对回归系数的大小进行惩罚，如下所示：

弹性网络的缺点

弹性网络（Elastic net）涉及选择 λ[1] 和 λ[2] 作为模型良好性能的关键值。这些参数通常通过交叉验证技术获得。从这些方法中，Lasso 和弹性网络通常用于特征选择。目前，我们的数据集中有 96 个变量；我们决定不减少变量的数量。

维度降低

维度投影，或特征投影，包括将高维空间中的数据转换为低维空间。

高维数大大增加了计算复杂度，甚至可能增加过拟合的风险。

降维技术对特征选择也很有用。在这种情况下，变量通过不同的组合转换为其他新变量。这些组合通过较少的变量从复杂的数据库中提取和总结相关信息。

存在着不同的算法，以下是最重要的：

主成分分析（PCA）
Sammon 映射
奇异值分解（SVD）
Isomap
局部线性嵌入（LLE）
拉普拉斯特征映射
t 分布随机邻域嵌入（t-SNE）

尽管在诸如故障预测模型或信用风险等情况下，降维并不常见，但我们将看到我们数据中的一个例子。

我们还将看到 PCA 和 t-SNE 的应用，它们是最常用的算法。

PCA 是一种通过变量的线性变换来提取数据集上重要变量的方法。因此，我们可以将主成分定义为原始变量的归一化线性组合。

第一主成分是变量的线性组合，它捕捉了数据集中最大的方差。第一成分中捕捉到的方差越大，该成分捕捉到的信息就越多。第一成分只用一行就能最好地总结我们数据中的最大信息。第二和后续的主成分也是原始变量的线性组合，它们捕捉了数据中剩余的方差。

当变量高度相关时，PCA 也被使用。这种方法的主要特性之一是不同组件之间的相关性为零。

让我们看看在 R 中的实现。为此，我们使用rstat包中包含的prcomp函数：

pca <- prcomp(train[,3:ncol(train)], retx=TRUE, center=TRUE, scale=TRUE)

在实现 PCA 方法之前，变量应该被标准化。这意味着我们应该确保变量具有等于零的均值和等于 1 的标准差。

这可以通过使用同一函数中的scale和center选项作为参数来完成：

names(pca)
## [1] "sdev"     "rotation" "center"   "scale"    "x"

center和scale向量包含我们所使用的变量的均值和标准差。

旋转测量返回主成分。我们获得与样本中变量相同数量的主成分。

让我们打印出这些组件的外观。例如，前四个组件的第一行如下所示：

pca$rotation[1:10,1:4]
 ##                  PC1          PC2         PC3          PC4
 ## UBPRE395 -0.05140105  0.027212743  0.01091903 -0.029884263
 ## UBPRE543  0.13068409 -0.002667109  0.03250766 -0.010948699
 ## UBPRE586  0.13347952 -0.013729338  0.02583513 -0.030875234
 ## UBPRFB60  0.17390861 -0.042970061  0.02813868  0.016505787
 ## UBPRE389  0.07980840  0.069097429  0.08331793  0.064870471
 ## UBPRE393  0.08976446  0.115336263  0.02076018 -0.012963786
 ## UBPRE394  0.16230020  0.119853462  0.07177180  0.009503902
 ## UBPRE396  0.06572403  0.033857693  0.07952204 -0.005602078
 ## UBPRE417 -0.06109615 -0.060368186 -0.01204455 -0.155802734
 ## UBPRE419  0.08178735  0.074713474  0.11134947  0.069892907

每个组件解释了总方差的一部分。每个组件解释的方差比例可以按以下方式计算：

让我们先计算每个组件的方差：

pca_variances =pca$sdev²

然后将每个方差除以成分方差的和：

prop_var_explained <- pca_variances/sum(pca_variances)

head(prop_var_explained,10)

 ##  [1] 0.10254590 0.06510543 0.04688792 0.04055387 0.03637036          0.03576523
 ##  [7] 0.02628578 0.02409343 0.02305206 0.02091978

第一主成分解释了大约 10%的方差。第二成分解释了 6%的方差，以此类推。

我们可以使用此代码图形化地观察总方差及其贡献：

plot(pca, type = "l",main = " Variance of Principal components")

上述代码生成了以下内容：

让我们运行代码来绘制方差图：

plot(prop_var_explained, xlab = "Principal Component",
              ylab = "Proportion of Variance Explained",
              type = "b")

上述代码生成了以下图形：

前面的截图有助于确定解释总方差重要部分的变量数量或主成分数量。

因此，这些成分可以用来建模，而不是使用完整的变量列表。绘制累积解释方差很有趣：

plot(cumsum(prop_var_explained), xlab = "Principal Component",
 ylab = "Cumulative Proportion of Variance Explained",
 type = "b")

上述代码生成了以下图形：

根据前面的截图，前 20 个成分解释了我们数据集中大约 60%的总方差。

我们可以选择使用这 20 个成分来创建我们的模型。这种方法在信用风险模型中并不常见，所以我们不会使用这些转换。

然而，评估我们的数据集的外观很重要。在下面的截图中，我们使用前两个成分展示了数据的图形表示。

此外，我们根据相应的目标变量在图中对每个银行进行分类。这次，我们使用ggfortify包：

library(ggfortify)

train$Default<-as.factor(train$Default)
autoplot(pca, data = train, colour = 'Default'

此截图显示了失败和非失败银行的分类图：

只看两个成分非常有趣。尽管这些成分只解释了大约 17%的总方差，但失败和非失败的银行在某种程度上是区分开的。

降维技术

您应该考虑，主成分假设变量的线性变换，但还有其他非线性降维技术。

对我来说，最有趣的技巧之一是 Laurens van der Maaten 开发的 t-SNE，他说：

“作为一个理智的检查，尝试对你的数据进行 PCA，将其降低到二维。如果这也给出糟糕的结果，那么可能你的数据一开始就没有太多好的结构。如果 PCA 运行良好但 t-SNE 不行，我相当确信你做错了什么。”

让我们看看 t-SNE 在我们数据集上的应用示例。通常，建议您设置一个seed：

set.seed(1234)

我们需要使用Rtsne包。此包包含执行算法的Rtsne函数。最重要的参数如下：

pca：这确定在运行 t-SNE 之前是否执行主成分分析。
perplexity：这是信息的一个度量（定义为香农熵的 2 次方）。perplexity参数确定了每个观察中最近邻的数量。这个参数对算法很有用，因为它使它能够在你的数据观察中找到局部和全局关系之间的平衡。

运行算法的代码如下：

library(Rtsne)

tsne= Rtsne(as.matrix(train[,3:ncol(train)]), check_duplicates=TRUE, pca=TRUE, perplexity=75, theta=0.5, dims=2,max_iter = 2000,verbose=TRUE)

这个过程需要几分钟才能完成。有关算法工作原理的更多信息也包含在包文档及其参考文献中。

通常，完整的数据集被简化为只有两个向量：

tsne_vectors = as.data.frame(tsne$Y)

 head(tsne_vectors)
 ##           V1          V2
 ## 1  -4.300888 -14.9082526
 ## 2   4.618766  44.8443129
 ## 3  21.554283   3.2569812
 ## 4  45.518532   0.7150365
 ## 5  12.098218   4.9833460
 ## 6 -14.510530  31.7903585

让我们根据其向量绘制我们的训练数据集：

ggplot(tsne_vectors, aes(x=V1, y=V2)) +
   geom_point(size=0.25) +
   guides(colour=guide_legend(override.aes=list(size=6))) +
   xlab("") + ylab("") +
   ggtitle("t-SNE") +
   theme_light(base_size=20) +
   theme(axis.text.x=element_blank(),
         axis.text.y=element_blank()) +
   scale_colour_brewer(palette = "Set2")

前面的代码生成了以下图表：

现在让我们再次绘制它，为每个目标值分配一个颜色，以及失败和未失败的银行：

plot(tsne$Y, t='n', main="tsne",xlab="Vector X",ylab="Vector y")
 text(tsne$Y, labels=as.vector(train$Default), col=c('red', 'blue')[as.numeric(train$Default)])

前面的代码生成了以下图表：

我们可以看到许多失败的银行被放置在结果的双变量地图的同一部分。然而，t-SNE 的一个主要弱点是算法的黑盒性质。基于结果无法对额外数据进行推断，这在使用 PCA 时是不会发生的。

t-SNE 主要用于探索性数据分析，它也被用作聚类算法的输入。

在这个实际案例中，我们试图在信用风险分析过程中保持准确，我们将忽略 PCA 和 t-SNE 的结果，并继续使用我们的原始维度。

一旦我们选定了最具有预测性的变量，我们将尝试使用不同的算法将它们结合起来。目标是开发一个具有最高准确性的模型来预测银行未来的破产。

在继续之前，让我们保存工作空间：

rm(list=setdiff(ls(), c("Model_database","train","test","table_iv")))

save.image("~/Data12.RData")

摘要

在本章中，我们看到了如何通过单变量分析减少我们问题数据的空间样本，并分析了数据。因此，在下一章中，我们将看到这些变量如何结合以获得一个准确的模型，其中将测试多个算法。

第五章：预测银行失败 - 多元分析

在本章中，我们将应用不同的算法，目的是通过我们的预测因子的组合来获得一个好的模型。在信用风险应用中，如信用评分和评级，最常用的算法是逻辑回归。在本章中，我们将看到其他算法如何应用于解决逻辑回归的一些弱点。

在本章中，我们将介绍以下主题：

逻辑回归
正则化方法
测试随机森林模型
梯度提升
神经网络中的深度学习
支持向量机
集成方法
自动机器学习

逻辑回归

从数学上讲，二元逻辑模型有一个具有两个分类值的因变量。在我们的例子中，这些值与银行是否有偿付能力相关。

在逻辑模型中，对数几率指的是一个类别的对数几率，它是一个或多个独立变量的线性组合，如下所示：

逻辑回归算法的系数（beta 值，β）必须使用最大似然估计来估计。最大似然估计涉及获取回归系数的值，以最小化模型预测的概率与实际观察案例之间的误差。

逻辑回归对异常值的存在非常敏感，因此应避免变量之间的高相关性。在 R 中应用逻辑回归的方法如下：

set.seed(1234)
LogisticRegression=glm(train$Default~.,data=train[,2:ncol(train)],family=binomial())
 ## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

代码运行没有问题，但出现了一个警告信息。如果变量高度相关或存在多重共线性，模型参数和方差膨胀是预期的。

高方差不是由于准确的或好的预测因子，而是由于一个未指定且有冗余预测因子的模型。因此，通过简单地添加更多参数来增加最大似然，这会导致过拟合。

我们可以使用summary()函数观察模型的参数：

summary(LogisticRegression)
 ## 
 ## Call:
 ## glm(formula = train$Default ~ ., family = binomial(), data =         train[, 
 ##     2:ncol(train)])
 ## 
 ## Deviance Residuals: 
 ##     Min       1Q   Median       3Q      Max  
 ## -3.9330  -0.0210  -0.0066  -0.0013   4.8724  
 ## 
 ## Coefficients:
 ##                   Estimate     Std. Error z value Pr(>|z|)  
 ## (Intercept) -11.7599825009   6.9560247460  -1.691   0.0909 .
 ## UBPRE395     -0.0575725641   0.0561441397  -1.025   0.3052  
 ## UBPRE543      0.0014008963   0.0294470630   0.048   0.9621  
 ##         ....                             .....                           ....                            ....             ....
 ## UBPRE021     -0.0114148389   0.0057016025  -2.002   0.0453 *
 ## UBPRE023      0.4950212919   0.2459506994   2.013   0.0441 *
 ## UBPRK447     -0.0210028916   0.0192296299  -1.092   0.2747  
 ## ---
 ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 ## 
 ## (Dispersion parameter for binomial family taken to be 1)
 ## 
 ##     Null deviance: 2687.03  on 7090  degrees of freedom
 ## Residual deviance:  284.23  on 6982  degrees of freedom
 ## AIC: 502.23
 ## 
 ## Number of Fisher Scoring iterations: 13

我们可以看到，前一个表的最后一列中的大多数变量都是不显著的。在这种情况下，应该减少回归中的变量数量，或者遵循另一种方法，例如惩罚或正则化方法。

正则化方法

使用正则化方法有三种常见的方法：

Lasso
岭回归
弹性网络

在本节中，我们将看到这些方法如何在 R 中实现。对于这些模型，我们将使用h2o包。这个包提供了一个开源的预测分析平台，用于机器学习，基于内存参数，分布式、快速且可扩展。它有助于创建基于大数据的模型，并且由于其提高了生产质量，因此最适合企业应用。

有关h2o软件包的更多信息，请访问其文档cran.r-project.org/web/packages/h2o/index.html。

此软件包非常有用，因为它在一个软件包中总结了几个常见的机器学习算法。此外，这些算法可以在我们的计算机上并行执行，因为它非常快。该软件包包括广义线性朴素贝叶斯、分布式随机森林、梯度提升和深度学习等。

不需要具备高级编程知识，因为该软件包自带用户界面。

让我们看看这个软件包是如何工作的。首先，应该加载软件包：

library(h2o)

使用h2o.init方法初始化 H2O。此方法接受可在软件包文档中找到的其他选项：

h2o.init()

建立我们的模型的第一步是将我们的数据放入 H2O 集群/Java 进程。在此步骤之前，我们将确保我们的目标是作为factor变量考虑：

train$Default<-as.factor(train$Default)

test$Default<-as.factor(test$Default)

现在，让我们将我们的数据上传到h2o集群：

as.h2o(train[,2:ncol(train)],destination_frame="train")

as.h2o(test[,2:ncol(test)],destination_frame="test")

如果您关闭 R 并在稍后重新启动它，您将需要再次上传数据集，如前面的代码所示。

我们可以使用以下命令检查数据是否已正确上传：

h2o.ls()

 ##     key
 ## 1  test
 ## 2 train

该软件包包含一个易于使用的界面，允许我们在浏览器中运行时创建不同的模型。通常，可以通过在网页浏览器中写入以下地址来启动界面，http://localhost:54321/flow/index.html。您将面对一个类似于以下截图的页面。在模型选项卡中，我们可以看到该软件包中实现的所有可用模型的列表：

首先，我们将开发正则化模型。为此，必须选择广义线性建模…。本模块包括以下内容：

高斯回归
泊松回归
二项式回归（分类）
多项式分类
伽马回归
有序回归

如以下截图所示，我们应该填写必要的参数来训练我们的模型：

我们将填写以下字段：

model_id：在这里，我们可以指定模型可以引用的名称。
training_frame：我们希望用于构建和训练模型的数据库可以在这里提及，因为这将是我们的训练数据集。
validation_frame：在这里，提到的数据集将用于检查模型的准确性。
nfolds：为了验证，我们需要在此处提及一定数量的折数。在我们的案例中，nfolds 的值是5。
seed：这指定了算法将使用的种子。我们将使用随机数生成器（RNG）为算法中需要随机数的组件提供随机数。
response_column：这是用作因变量的列。在我们的案例中，该列命名为 Default。
ignored_columns：在本节中，可以在训练过程中忽略变量。在我们的情况下，所有变量都被认为是相关的。
ignore_const_cols：这是一个标志，表示包应避免常数变量。
family：这指定了模型类型。在我们的情况下，我们想要训练一个回归模型，因此 family 应该固定为二项式，因为我们的目标变量有两个可能的值。
solver：这指定了要使用的求解器。我们不会更改此值，因为无论选择哪个求解器，都没有观察到显著差异。因此，我们将保持默认值。
alpha：在这里，您必须从 L1 到 L2 选择正则化分布的值。如果您选择 1，它将是一个 Lasso 回归。如果您选择 0，它将是一个 Ridge 回归。如果您选择 0 和 1 之间的任何值，您将得到 Lasso 和 Ridge 的混合。在我们的情况下，我们将选择 1。Lasso 模型的主要优点之一是减少变量的数量，因为训练模型将非相关变量的系数设为零，从而产生简单但同时又准确的模型。
lambda_search：此参数启动对正则化强度的搜索。
标准化：如果此标志被标记，则表示数值列将被转换，以具有零均值和零单位方差。

最后，构建模型按钮训练模型。尽管可以选择其他选项，但前面的规格已经足够：

我们可以看到模型训练得很快。查看按钮为我们提供了有关模型的一些有趣细节：

模型参数
得分历史
训练和验证样本的接收者操作特征（ROC）曲线
标准化系数幅度
交叉验证、训练和验证样本的增益/提升表
交叉验证模型
指标
系数

让我们看看一些主要结果：

如我们所见，我们的 Lasso 模型是用 108 个不同的变量训练的，但只有 56 个变量导致系数大于零的模型。

该模型提供了几乎完美的分类。在训练样本中，曲线下面积（AUC）达到 99.51%。在验证样本中，此值略低，为 98.65%。标准化变量也相关：

如果一个变量以蓝色显示，这表示系数是正的。如果是负的，颜色是橙色。

如我们所见，UBPRE626 看起来是一个重要的变量。它显示总贷款和租赁融资应收账款超过实际股本总额的次数。这里的正号意味着更高的比率，这也意味着银行在其运营中失败的概率更高。

根据这个图，前五个相关变量如下：

UBPRE626：净贷款和租赁融资应收账款超过总股本次数
UBPRE545：应计和未实现贷款及租赁总额，除以贷款和租赁损失准备金
UBPRE170：总股本
UBPRE394：其他建筑和土地开发贷款，除以平均总贷款和租赁
UBPRE672：证券年度化实现收益（或损失）的四分之一，除以平均资产

在评估信用风险时，了解哪些变量最为重要以及这些变量的经济相关性非常重要。例如，如果非正常贷款或有问题贷款越多，银行的偿债能力就越高，那就没有意义了。我们并不关心模型中变量的经济意义，但这对于金融机构开发的一些模型来说是一个关键问题。如果变量没有预期的符号，它们必须从模型中移除。

在某些情况下，有必要测试不同的参数组合，直到获得最佳模型。例如，在最近训练的正则化模型中，我们可以尝试不同的 alpha 参数值。要同时测试不同的参数，需要使用代码执行算法。让我们看看如何做。这次我们将再次训练正则化模型，但这次使用一些代码。首先，我们从 h2o 系统中删除所有对象，包括最近创建的模型：

h2o.removeAll()
 ## [1] 0

然后，我们再次上传我们的训练和验证样本：

as.h2o(train[,2:ncol(train)],destination_frame="train")
as.h2o(test[,2:ncol(test)],destination_frame="test")

让我们编写我们的模型。创建一个空参数网格如下：

grid_id <- 'glm_grid'

然后，我们在网格中分配不同的参数进行测试：

hyper_parameters <- list( alpha = c(0, .5, 1) )
stopping_metric <- 'auc'
glm_grid <- h2o.grid(
     algorithm = "glm",
     grid_id = grid_id,
     hyper_params = hyper_parameters,
     training_frame = training,
     nfolds=5,
     x=2:110,
     y=1,
     lambda_search = TRUE,
     family = "binomial", seed=1234)

如我们所见，参数与我们用来训练先前模型的参数完全相同。唯一的区别是我们现在同时使用不同的 alpha 值，这对应于岭回归、Lasso 和弹性网络。模型使用以下代码进行训练：

results_glm <- h2o.getGrid(
     grid_id = grid_id,
     sort_by = stopping_metric,
     decreasing = TRUE)

根据之前的代码，网格中的不同模型应该按照 AUC 指标排序。因此，我们对第一个模型感兴趣：

best_GLM <- h2o.getModel(results_glm@model_ids[[1]])

让我们来看看这个模型的一些细节：

best_GLM@model$model_summary$regularization
 ## [1] "Ridge ( lambda = 0.006918 )"

表现最好的模型是一个岭回归模型。模型的性能可以通过以下方式获得：

perf_train<-h2o.performance(model = best_GLM,newdata = training)
perf_train
 ## H2OBinomialMetrics: glm
 ##
 ## MSE:  0.006359316
 ## RMSE:  0.07974532
 ## LogLoss:  0.02561085
 ## Mean Per-Class Error:  0.06116986
 ## AUC:  0.9953735
 ## Gini:  0.990747
 ## R²:  0.8579102
 ## Residual Deviance:  363.213
 ## AIC:  581.213
 ##
 ## Confusion Matrix (vertical: actual; across: predicted) for F1-              optimal threshold:
 ##           0   1    Error      Rate
 ## 0      6743  15 0.002220  =15/6758
 ## 1        40 293 0.120120   =40/333
 ## Totals 6783 308 0.007756  =55/7091
 ##
 ## Maximum Metrics: Maximum metrics at their respective thresholds
 ##                         metric threshold    value idx
 ## 1                       max f1  0.540987 0.914197 144
 ## 2                       max f2  0.157131 0.931659 206
 ## 3                 max f0point5  0.617239 0.941021 132
 ## 4                 max accuracy  0.547359 0.992244 143
 ## 5                max precision  0.999897 1.000000   0
 ## 6                   max recall  0.001351 1.000000 383
 ## 7              max specificity  0.999897 1.000000   0
 ## 8             max absolute_mcc  0.540987 0.910901 144
 ## 9   max min_per_class_accuracy  0.056411 0.972973 265
 ## 10 max mean_per_class_accuracy  0.087402 0.977216 239
 ##

AUC 和 Gini 指数，作为性能的主要指标，仅略高于我们最初训练的 Lasso 模型——至少在训练样本中是这样。

模型在测试样本中的性能也很高：

perf_test<-h2o.performance(model = best_GLM,newdata = as.h2o(test))
perf_test
 ## H2OBinomialMetrics: glm
 ##
 ## MSE:  0.01070733
 ## RMSE:  0.1034762
 ## LogLoss:  0.04052454
 ## Mean Per-Class Error:  0.0467923
 ## AUC:  0.9875425
 ## Gini:  0.975085
 ## R²:  0.7612146
 ## Residual Deviance:  246.3081
 ## AIC:  464.3081
 ##
 ## Confusion Matrix (vertical: actual; across: predicted) for F1-            optimal threshold:
 ##           0   1    Error      Rate
 ## 0      2868  28 0.009669  =28/2896
 ## 1        12 131 0.083916   =12/143
 ## Totals 2880 159 0.013162  =40/3039
 ##
 ## Maximum Metrics: Maximum metrics at their respective thresholds
 ##                         metric threshold    value idx
 ## 1                       max f1  0.174545 0.867550 125
 ## 2                       max f2  0.102341 0.904826 138
 ## 3                 max f0point5  0.586261 0.885167  89
 ## 4                 max accuracy  0.309187 0.987167 107
 ## 5                max precision  0.999961 1.000000   0
 ## 6                   max recall  0.000386 1.000000 388
 ## 7              max specificity  0.999961 1.000000   0
 ## 8             max absolute_mcc  0.174545 0.861985 125
 ## 9   max min_per_class_accuracy  0.027830 0.955456 210
 ## 10 max mean_per_class_accuracy  0.102341 0.965295 138

与 Lasso 模型相比，结果没有显著差异。尽管如此，Lasso 模型的系数数量较少，这使得它更容易解释且更简洁。

岭回归中的系数总数等于数据集中的变量数量和模型的截距：

head(best_GLM@model$coefficients)
 ##    Intercept     UBPRE395     UBPRE543     UBPRE586     UBPRFB60
 ## -8.448270911 -0.004167366 -0.003376142 -0.001531582  0.027969152
 ##     UBPRE389
 ## -0.004031844

现在，我们将每个模型的预测结果存储在一个新的数据框中。我们可以将不同模型的成果结合起来，得到一个额外的模型。最初，我们的数据框将只包含每个银行的 ID 和目标变量：

summary_models_train<-train[,c("ID_RSSD","Default")]
summary_models_test<-test[,c("ID_RSSD","Default")]

让我们计算模型预测并将它们存储在汇总数据框中：

summary_models_train$GLM<-as.vector(h2o.predict(best_GLM,training)[3])
summary_models_test$GLM<-as.vector(h2o.predict(best_GLM,validation)[3])

当我们运行前面的代码来计算模型的性能时，我们还获得了一个混淆矩阵。例如，在测试样本中，我们得到以下结果：

perf_test@metrics$cm$table
 ## Confusion Matrix: Row labels: Actual class; Column labels: Predicted class
 ##           0   1  Error         Rate
 ## 0      2868  28 0.0097 = 28 / 2,896
 ## 1        12 131 0.0839 =   12 / 143
 ## Totals 2880 159 0.0132 = 40 / 3,039

此软件包将违约概率高于 0.5 的银行分类为失败银行，否则为成功银行。

根据这个假设，40家银行被错误分类（28+12）。然而，0.5 的截止点实际上并不正确，因为在样本中失败银行与非失败银行的比率是不同的。

失败银行的比率实际上仅为 4.696%，如下代码所示：

mean(as.numeric(as.character(train$Default)))
 ## [1] 0.04696094

因此，如果一个银行的违约概率高于这个比例，将其视为失败银行更为合适：

aux<-summary_models_test
aux$pred<-ifelse(summary_models_test$GLM>0.04696094,1,0)

因此，测试样本的新混淆表如下：

table(aux$Default,aux$pred)
 ##   
 ##        0    1
 ##   0 2818   78
 ##   1    8  135

根据这张表，模型错误分类了 86 家银行（78+8）。几乎所有失败银行都被正确分类。要获得比这个更好的算法将非常困难。

模型可以使用h2o.saveModel本地保存：

h2o.saveModel(object= best_GLM, path=getwd(), force=TRUE)

我们从工作区中移除无关的对象，并按以下方式保存：

rm(list=setdiff(ls(), c("Model_database","train","test","summary_models_train","summary_models_test","training","validation")))

save.image("Data13.RData")

记住，如果你关闭 R 并再次加载此工作区，你应该再次将你的训练和测试样本转换为h2o格式：

training<-as.h2o(train[,2:ncol(train)],destination_frame=“train”)
validation<-as.h2o(test[,2:ncol(test)],destination_frame=“test”)

测试随机森林模型

随机森林是一组决策树。在决策树中，基于独立变量的训练样本将被分割成两个或更多同质集合。此算法处理分类变量和连续变量。使用递归选择方法选择最佳属性，并将其分割成叶节点。这会一直持续到满足停止循环的准则。通过叶节点的扩展创建的每个树都被视为一个弱学习器。这个弱学习器建立在子集的行和列之上。树的数量越多，方差越低。分类和回归随机森林都会计算所有树的平均预测，以做出最终预测。

当训练随机森林时，可以设置一些不同的参数。其中最常见的参数包括树的数量、最大变量数、终端节点的大小以及每棵树的深度。应该进行多次测试，以在性能和过拟合之间找到平衡。例如，树的数量和深度越高，训练集上的准确率越好，但这增加了过拟合的风险。为了获得这种平衡，应该在验证集上测试几个参数及其参数组合，然后在训练过程中进行交叉验证。

再次，这个算法在 h2o 包中很容易实现，可以使用浏览器中的可视化指南。参数网格应该通过编码实现。代码几乎与先前的模型相同。然而，这次过程更耗时：

grid_space <- list()
 grid_space$ntrees <- c(25, 50, 75)
 grid_space$max_depth <- c(4, 10, 20)
 grid_space$mtries <- c(10, 14, 20)
 grid_space$seed <- c(1234)

 grid <- h2o.grid("randomForest", grid_id="RF_grid", x=2:110,y=1,training_frame=training, nfolds=5, hyper_params=grid_space)

results_grid <- h2o.getGrid(grid_id = "RF_grid",
                              sort_by = "auc",
                              decreasing = TRUE)

print(results_grid)

 ## H2O Grid Details
 ## ================
 ## 
 ## Grid ID: RF_grid 
 ## Used hyper parameters: 
 ##   -  max_depth 
 ##   -  mtries 
 ##   -  ntrees 
 ##   -  seed 
 ## Number of models: 27 
 ## Number of failed models: 0 
 ## 
 ## Hyper-Parameter Search Summary: ordered by decreasing auc
 ##   max_depth mtries ntrees seed        model_ids                auc
 ## 1        20     20     75 1234 RF_grid_model_26 0.9928546480780869
 ## 2        10     10     75 1234 RF_grid_model_19 0.9922021014799943
 ## 3        10     10     50 1234 RF_grid_model_10 0.9921534437663471
 ## 4        10     20     75 1234 RF_grid_model_25 0.9920343545676484
 ## 5        10     20     50 1234 RF_grid_model_16 0.9919039341205663
 ## 
 ## ---
 ##    max_depth mtries ntrees seed       model_ids                auc
 ## 22        20     20     25 1234 RF_grid_model_8 0.9879017816277361
 ## 23        20     10     25 1234 RF_grid_model_2 0.9876307203918924
 ## 24        10     20     25 1234 RF_grid_model_7 0.9873765449379537
 ## 25        10     14     25 1234 RF_grid_model_4  0.986949956763511
 ## 26         4     10     25 1234 RF_grid_model_0  0.984477522802471
 ## 27        20     14     25 1234 RF_grid_model_5  0.980687331308817

如我们所见，树的数量（ntrees）、深度（max_depth）和每棵树中要考虑的变量数量（mtries）的组合被测试。使用 AUC 指标对生成的模型进行排序。

根据 preceding 规范，已经训练了 27 个不同的模型。选择第一个模型，或者准确率最高的模型：

best_RF <- h2o.getModel(results_grid@model_ids[[1]])

该模型在训练和测试样本上的性能如下：

h2o.performance(model = best_RF,newdata = training)
 ## H2OBinomialMetrics: drf
 ## 
 ## MSE:  0.001317125
 ## RMSE:  0.03629222
 ## LogLoss:  0.009026859
 ## Mean Per-Class Error:  0
 ## AUC:  1
 ## Gini:  1
 ## 
 ## Confusion Matrix (vertical: actual; across: predicted) for F1-            optimal          threshold:
 ##           0   1    Error     Rate
 ## 0      6758   0 0.000000  =0/6758
 ## 1         0 333 0.000000   =0/333
 ## Totals 6758 333 0.000000  =0/7091
 ## 
 ## Maximum Metrics: Maximum metrics at their respective thresholds
 ##                         metric threshold    value idx
 ## 1                       max f1  0.586667 1.000000  29
 ## 2                       max f2  0.586667 1.000000  29
 ## 3                 max f0point5  0.586667 1.000000  29
 ## 4                 max accuracy  0.586667 1.000000  29
 ## 5                max precision  1.000000 1.000000   0
 ## 6                   max recall  0.586667 1.000000  29
 ## 7              max specificity  1.000000 1.000000   0
 ## 8             max absolute_mcc  0.586667 1.000000  29
 ## 9   max min_per_class_accuracy  0.586667 1.000000  29
 ## 10 max mean_per_class_accuracy  0.586667 1.000000  29

如您所见，代码与先前的模型完全相同。

现在，我们使用以下代码找到测试或验证样本的性能：

h2o.performance(model = best_RF,newdata = validation)
 ## H2OBinomialMetrics: drf
 ## 
 ## MSE:  0.00940672
 ## RMSE:  0.09698825
 ## LogLoss:  0.05488315
 ## Mean Per-Class Error:  0.06220299
 ## AUC:  0.9882138
 ## Gini:  0.9764276
 ## 
 ## Confusion Matrix (vertical: actual; across: predicted) for F1-            optimal              threshold:
 ##           0   1    Error      Rate
 ## 0      2880  16 0.005525  =16/2896
 ## 1        17 126 0.118881   =17/143
 ## Totals 2897 142 0.010859  =33/3039
 ## 
 ## Maximum Metrics: Maximum metrics at their respective thresholds
 ##                         metric threshold    value idx
 ## 1                       max f1  0.346667 0.884211  44
 ## 2                       max f2  0.280000 0.897790  49
 ## 3                 max f0point5  0.760000 0.897196  18
 ## 4                 max accuracy  0.346667 0.989141  44
 ## 5                max precision  1.000000 1.000000   0
 ## 6                   max recall  0.000000 1.000000  70
 ## 7              max specificity  1.000000 1.000000   0
 ## 8             max absolute_mcc  0.346667 0.878520  44
 ## 9   max min_per_class_accuracy  0.106667 0.965035  62
 ## 10 max mean_per_class_accuracy  0.106667 0.968878  62

两个样本的结果几乎完美。可以获取变量的重要性：

var_importance<-data.frame(best_RF@model$variable_importances)
 h2o.varimp_plot(best_RF,20)

上述代码生成以下输出：

就像岭回归一样，破产的概率将存储在训练和验证样本中：

summary_models_train$RF<-as.vector(h2o.predict(best_RF,training)[3])
summary_models_test$RF<-as.vector(h2o.predict(best_RF,validation)[3])

最后，我们可以计算混淆矩阵。记住，根据观察到的坏银行在总样本中的比例，确定将银行分类为破产概率的截止值：

aux<-summary_models_test
 aux$pred<-ifelse(summary_models_test$RF>0.04696094,1,0)
 table(aux$Default,aux$pred)
 ##        0    1
 ##   0 2753  143
 ##   1    5  138

如果将随机森林和岭回归模型进行比较，我们可以看到随机森林只错误分类了五家失败的银行，而岭回归中有 12 家银行被错误分类。尽管如此，随机森林将更多有偿付能力的银行错误分类为失败的银行，这意味着它有较高的假阳性率。

无关对象再次从工作区中移除。此外，我们备份了我们的工作区：

rm(list=setdiff(ls(), c("Model_database","train","test","summary_models_train","summary_models_test","training","validation")))
save.image("Data14.RData")

梯度提升

梯度提升意味着将弱预测器和平均预测器结合以获得一个强预测器。这确保了鲁棒性。它与随机森林类似，主要基于决策树。区别在于样本在树之间没有修改；只有不同观察的权重被修改。

提升通过使用先前训练的树的信息按顺序训练树。为此，我们首先需要使用训练数据集创建决策树。然后，我们需要创建另一个模型，它除了纠正训练模型中发生的错误之外什么都不做。这个过程会按顺序重复，直到达到指定的树的数量或达到某个停止规则。

关于该算法的更具体细节可以在h2o包的文档中找到。在训练算法时，我们需要定义诸如我们将要组合的树的数量和每个节点中的最小观察值等参数，就像我们在随机森林中做的那样。

收缩参数，或者说提升学习的学习速率，可以改变模型的性能。我们需要考虑许多实验的结果，以确定最佳参数，确保高精度。

我们的参数网格收集了树的数量和max_depth参数的不同组合：

grid_space <- list()
grid_space$ntrees <- c(25,75,100)
grid_space$max_depth = c(4,6,8,12,16,20)

通过执行以下代码将训练不同的模型：

gbm_grid <- h2o.grid(hyper_params = grid_space,
   algorithm = "gbm",
   grid_id ="Grid1", 
   x=2:110,
   y=1,
   training_frame = training,seed=1234)

网格是按照AUC排序的。结果如下：

results_gbm <- h2o.getGrid("Grid1", sort_by = "AUC", decreasing = TRUE)    
results_gbm
 ## H2O Grid Details
 ## ================
 ## 
 ## Grid ID: Grid1 
 ## Used hyper parameters: 
 ##   -  max_depth 
 ##   -  ntrees 
 ## Number of models: 18 
 ## Number of failed models: 0 
 ## 
 ## Hyper-Parameter Search Summary: ordered by decreasing AUC
 ##    max_depth ntrees      model_ids                auc
 ## 1         16    100 Grid1_model_16                1.0
 ## 2          4    100 Grid1_model_12                1.0
 ## 3         16     25  Grid1_model_4                1.0
 ## 4         20     75 Grid1_model_11                1.0
 ## 5          6     75  Grid1_model_7                1.0
 ## 6         20    100 Grid1_model_17                1.0
 ## 7          8     75  Grid1_model_8                1.0
 ## 8         20     25  Grid1_model_5                1.0
 ## 9         12     75  Grid1_model_9                1.0
 ## 10        16     75 Grid1_model_10                1.0
 ## 11         6    100 Grid1_model_13                1.0
 ## 12        12    100 Grid1_model_15                1.0
 ## 13         8    100 Grid1_model_14                1.0
 ## 14         4     75  Grid1_model_6 0.9999986669119549
 ## 15        12     25  Grid1_model_3 0.9999986669119549
 ## 16         8     25  Grid1_model_2 0.9999922236530701
 ## 17         6     25  Grid1_model_1 0.9998680242835318
 ## 18         4     25  Grid1_model_0 0.9977795196794901

大多数模型获得了完美的分类。这可能是一个过拟合的迹象。让我们看看第一个模型在验证样本上的性能：

best_GBM <- h2o.getModel(results_gbm@model_ids[[1]])
h2o.performance(model = best_GBM,newdata = as.h2o(test))
 ## H2OBinomialMetrics: gbm
 ## 
 ## MSE:  0.01053012
 ## RMSE:  0.1026164
 ## LogLoss:  0.06001792
 ## Mean Per-Class Error:  0.05905179
 ## AUC:  0.9876222
 ## Gini:  0.9752444
 ## 
 ## Confusion Matrix (vertical: actual; across: predicted) for F1-            optimal          threshold:
 ##           0   1    Error      Rate
 ## 0      2878  18 0.006215  =18/2896
 ## 1        16 127 0.111888   =16/143
 ## Totals 2894 145 0.011188  =34/3039
 ## 
 ## Maximum Metrics: Maximum metrics at their respective thresholds
 ##                         metric threshold    value idx
 ## 1                       max f1  0.076792 0.881944 143
 ## 2                       max f2  0.010250 0.892857 154
 ## 3                 max f0point5  0.852630 0.906902 118
 ## 4                 max accuracy  0.076792 0.988812 143
 ## 5                max precision  0.999962 1.000000   0
 ## 6                   max recall  0.000006 1.000000 392
 ## 7              max specificity  0.999962 1.000000   0
 ## 8             max absolute_mcc  0.076792 0.876096 143
 ## 9   max min_per_class_accuracy  0.000181 0.958042 246
 ## 10 max mean_per_class_accuracy  0.000816 0.963611 203

结果非常好，甚至在测试样本中也是如此。选择了第一个模型，并将预测存储起来，就像之前的模型一样：

summary_models_train$GBM<-as.vector(h2o.predict(best_GBM,training)[3])
summary_models_test$GBM<-as.vector(h2o.predict(best_GBM,validation)[3])

最后，测试样本中的混淆表被计算出来：

aux<-summary_models_test
aux$pred<-ifelse(summary_models_test$GBM>0.04696094,1,0)
table(aux$Default,aux$pred)
 ##    
 ##        0    1
 ##   0 2876   20
 ##   1   15  128

总共有 14 家破产银行和 25 家非破产银行被错误分类。让我们保存结果：

rm(list=setdiff(ls(), c("Model_database","train","test","summary_models_train","summary_models_test","training","validation")))
save.image("Data15.RData")

神经网络中的深度学习

对于机器学习，我们需要能够处理非线性且无关数据集的系统。这对于我们预测破产问题非常重要，因为违约变量和解释变量之间的关系很少是线性的。因此，使用神经网络是最佳解决方案。

人工神经网络（ANNs）长期以来一直被用来解决破产问题。ANN 是一个具有多个相互连接处理器的计算机系统。这些处理器通过处理信息和动态响应提供的输入来提供输出。ANN 的一个突出和基本示例是多层感知器（MLP）。MLP 可以表示如下：

除了输入节点外，每个节点都是一个使用非线性激活函数的神经元，该函数被发送进来。

从其图表中可以看出，多层感知器（MLP）不过是一个前馈神经网络。这意味着提供的输入信息将只向前移动。这种类型的网络通常由一个输入层、一个隐藏层和一个输出层组成。输入层代表模型的输入数据，或变量。在我们的案例中，这些是金融变量。在这个层中不进行任何计算。隐藏层是进行中间处理或计算的地方。它们执行计算，然后将权重（信号或信息）从输入层传递到下一层。最后，输出层从隐藏层接收输入并计算网络的输出。输入节点使用非线性激活函数将信息从一层传递到下一层。激活函数的目的是将输入信号转换为输出信号，以模拟复杂的非线性模式。

感知器网络通过在处理完每一组数据后修改权重来学习。这些权重指定了处理输入时发生的错误数量，这是通过比较期望的输出获得的。

深度学习与 MLP 有何不同？MLP 只是深度学习算法的一种。在许多情况下，深度学习与 MLP 网络不同，但这仅仅是因为计算复杂性和隐藏层数量的不同。深度学习可以被视为具有两个或更多隐藏层的 MLP。当包含两个或更多隐藏层时，学习过程也应该有所不同，因为 MLP 中使用的反向传播学习规则将失效。感知器更新规则容易产生消失和爆炸梯度，这使得训练多于一层或两层的网络变得困难。

设计神经网络

在设计多层网络时，确保你确定适当的层数以获得更好的准确性和精度。通常，对于许多模型来说，仅仅一个隐藏层就足以解决分类问题。然而，使用多个隐藏层已经在语音识别或物体检测等领域证明了其有用性。另一个需要考虑的是每个隐藏层中的神经元数量。这是一个非常重要的方面。估计这些值时的错误可能导致过度拟合（当添加太多神经元时）和欠拟合（当添加的神经元不足时）等问题。

训练神经网络

h2o包帮助我们训练神经网络。深度学习模型有许多输入参数。在这个练习中，以下参数将被测试：

hyper_params <- list(
   hidden=list(c(5),c(80,80,80),c(75,75)),
   input_dropout_ratio=c(0.05,0.1,0.15,0.2,0.25),
   rate=c(0.01,0.02,0.10))

在这里，我们将测试三种结构：首先，一个只包含一个包含 25 个神经元的隐藏层的网络，然后是一个包含三个隐藏层，每个层有 32 个神经元的网络，最后是一个包含两个隐藏层，每个层有 64 个神经元的网络。

神经网络学习，神经元逐渐在特定变量的值上实现专业化。如果神经元在训练集中过于专业化，就有很高的过拟合风险。为了避免过拟合，包含了input_dropout_ratio命令。dropout 技术是神经网络模型的一种正则化方法，用于提高神经网络的泛化能力。

在训练过程中，dropout 方法会随机选择神经元并在训练中忽略它们。在实践中，在每一步训练中，都会创建一个不同的网络，因为一些随机单元被移除，并使用反向传播进行训练，就像通常一样。这迫使网络学习具有相同输入和输出的多个独立表示，从而提高泛化能力。

要获取有关 dropout 方法的更多信息，我建议阅读原始论文，Dropout: *一种简单防止神经网络过拟合的方法，作者为 Nitish Srivastava 等人。该论文可在www.cs.toronto.edu/~hinton/absps/JMLRdropout.pdf找到。

建议的输入层 dropout 比值为 0.1 或 0.2。最后，使用rate命令，我们可以指定学习率。记住，如果学习率设置得太高，模型可能变得不稳定；如果设置得太低，则收敛将非常缓慢。

让我们编写一些训练代码：

deep_grid <- h2o.grid(
   algorithm="deeplearning",
   grid_id="dl_grid", 
   training_frame=training,
   validation_frame=as.h2o(test),
   x=2:110,
   y=1,
   epochs=2,
   stopping_metric="AUC",
   stopping_tolerance=1e-2,
   stopping_rounds=2,
   score_duty_cycle=0.01,  
   l1=1e-5,
   l2=1e-5,
   activation=c("Rectifier"),
   nfolds=5,
   hyper_params=hyper_params,standardize=TRUE,seed=1234)

上述代码中的参数可以描述如下：

epochs：这里指定的值决定了在学习过程中数据集需要流式传输的次数。
stopping_metric：这指定了用于早期停止的指标，在我们的案例中是 AUC。
stopping_tolerance和stopping_rounds：这些参数分别确定在模型停止学习前的容忍值和停止值，当stopping_metric在指定轮次后没有改善时，防止模型继续学习。当指定交叉验证（如我们案例中所示）时，此选项将应用于所有交叉验证模型。在我们的案例中，我们设置了stopping_tolerance=1e-2和stopping_rounds = 2，这意味着模型将在 2 轮迭代后停止训练，或者如果没有至少 2%（1e-2）的改善。
score_duty_cycle：这表示在评分和训练之间花费多少时间。这些值是介于 0 到 1 之间的百分比。较低的值表示更多的训练。此选项的默认值为 0.1，表示应该花费 10%的时间进行评分，剩余的 90%用于训练。
l1和l2：这里添加的值是正则化指数，它确保了更好的泛化能力和稳定性。
activation：可以在此处提及如tanh、tanh with dropout、Maxout等激活函数。
nfolds：这表示交叉验证的折数。由于要测试多个配置，因此训练过程非常耗时。可以通过运行以下代码来获得不同配置的性能：

results_deep <- h2o.getGrid("dl_grid",sort_by="auc",decreasing=TRUE)
results_deep
 ## H2O Grid Details
 ## ================
 ## 
 ## Grid ID: dl_grid 
 ## Used hyper parameters: 
 ##   -  hidden 
 ##   -  input_dropout_ratio 
 ##   -  rate 
 ## Number of models: 45 
 ## Number of failed models: 0 
 ## 
 ## Hyper-Parameter Search Summary: ordered by decreasing auc
 ##         hidden input_dropout_ratio rate        model_ids
 ## 1     [75, 75]                0.25 0.01 dl_grid_model_14
 ## 2     [75, 75]                0.25  0.1 dl_grid_model_44
 ## 3     [75, 75]                 0.2 0.01 dl_grid_model_11
 ## 4 [80, 80, 80]                0.25 0.02 dl_grid_model_28
 ## 5     [75, 75]                 0.1 0.01  dl_grid_model_5
 ##                  auc
 ## 1 0.9844357527103902
 ## 2 0.9841366966255987
 ## 3 0.9831344365969994
 ## 4 0.9830902225101693
 ## 5 0.9830724480029008
 ## 
 ## ---
 ##    hidden input_dropout_ratio rate        model_ids                auc
 ## 40    [5]                 0.1  0.1 dl_grid_model_33 0.9603608491593103
 ## 41    [5]                 0.1 0.01  dl_grid_model_3 0.9599749201702442
 ## 42    [5]                 0.2 0.01  dl_grid_model_9 0.9599749201702442
 ## 43    [5]                 0.2 0.02 dl_grid_model_24 0.9591890647676383
 ## 44    [5]                0.05 0.02 dl_grid_model_15 0.9587149297862527
 ## 45    [5]                0.15  0.1 dl_grid_model_36 0.9575646969846437

使用三个具有 32 个单位的隐藏层、0.25的 dropout 比率和0.02的学习率获得了最佳模型。

模型选择如下：

best_deep <- h2o.getModel(results_deep@model_ids[[1]])

测试样本的性能如下：

h2o.performance(model = best_deep,newdata = validation)
 ## H2OBinomialMetrics: deeplearning
 ## 
 ## MSE:  0.02464987
 ## RMSE:  0.1570028
 ## LogLoss:  0.1674725
 ## Mean Per-Class Error:  0.1162044
 ## AUC:  0.9794568
 ## Gini:  0.9589137
 ## 
 ## Confusion Matrix (vertical: actual; across: predicted) for F1-         optimal threshold:
 ##           0   1    Error      Rate
 ## 0      2871  25 0.008633  =25/2896
 ## 1        32 111 0.223776   =32/143
 ## Totals 2903 136 0.018756  =57/3039
 ## 
 ## Maximum Metrics: Maximum metrics at their respective thresholds
 ##                         metric threshold    value idx
 ## 1                       max f1  0.001538 0.795699 135
 ## 2                       max f2  0.000682 0.812672 153
 ## 3                 max f0point5  0.011028 0.831904 109
 ## 4                 max accuracy  0.001538 0.981244 135
 ## 5                max precision  0.999998 1.000000   0
 ## 6                   max recall  0.000000 1.000000 398
 ## 7              max specificity  0.999998 1.000000   0
 ## 8             max absolute_mcc  0.001538 0.786148 135
 ## 9   max min_per_class_accuracy  0.000017 0.937063 285
 ## 10 max mean_per_class_accuracy  0.000009 0.943153 314

与先前的模型一样，训练和验证样本的预测被存储：

summary_models_train$deep<-as.vector(h2o.predict(best_deep,training)[3])
summary_models_test$deep<- as.vector(h2o.predict(best_deep,validation)[3])

验证样本的混淆矩阵也获得了：

aux<-summary_models_test
 aux$pred<-ifelse(summary_models_test$deep>0.04696094,1,0)
 table(aux$Default,aux$pred)
 ##    
 ##        0    1
 ##   0 2886   10
 ##   1   61   82

rm(list=setdiff(ls(), c("Model_database","train","test","summary_models_train","summary_models_test","training","validation")))
save.image("Data16.RData")

支持向量机

支持向量机（SVM）算法是一种监督学习技术。为了理解这个算法，请查看以下图中关于最优超平面和最大边界的图示：

在这个分类问题中，我们只有两个类别，它们对应于许多可能的解决方案。如图所示，支持向量机通过计算最优超平面并最大化类别间的边界来对这些对象进行分类。这两者都将最大限度地区分类别。位于边界最近的样本被称为支持向量。然后，问题被处理为一个优化问题，可以通过优化技术来解决，其中最常见的是使用拉格朗日乘数法。

即使在可分线性问题中，如图所示，有时也不总是能够获得完美的分离。在这些情况下，支持向量机模型是最大化边界同时最小化误分类数量的模型。在现实世界中，问题之间的距离太远，无法进行线性分离，至少在没有先前的数据预处理或转换的情况下。以下图中显示了线性可分问题和非线性可分问题之间的区别：

为了处理非线性问题，核函数将数据映射到不同的空间。这意味着数据被转换到更高维的空间。这种技术被称为核技巧，因为有时可以在数据中执行类之间的线性分离，从而进行转换。

支持向量机算法的优点如下：

支持向量机很简单
支持向量机是统计和机器学习技术的结合
支持向量机可以用于解决像我们问题陈述中的金融问题

选择支持向量机参数

让我们讨论一些我们可能需要的参数，这样我们就可以使用支持向量机。

支持向量机核参数

支持向量机的主要困难之一是选择将数据转换的核函数。以下是最常用的转换：

线性
多项式
径向基
Sigmoid

成本参数

控制清算和训练误差与模型复杂度之间的交易将由成本参数（C）负责。如果你为C设置一个相对较小的数值，就会有更多的训练误差。如果C是一个较大的数值，你可能会得到一个过拟合的模型，这意味着你的训练模型已经学会了所有的训练数据，但这个模型可能无法在任何一个其他样本上正常工作。你可以将成本设置为 1.001 到 100 之间的任何值。

Gamma 参数

在使用高斯核时需要 gamma 参数。此参数将计算每个训练样本可以产生的影响水平。在这里，你可能认为较低的值是远的，而较高的值是近的。

Gamma 实际上是与我们之前看到的支持向量相反。因此，在 SVM 中，应该测试所有三个参数的不同值。

训练 SVM 模型

h2o包中没有 SVM 算法。为了训练 SVM 分类器，我们将使用caret包。记住，我们的目标值有两个不同的值：

levels(train$Default)
 ## [1] "0" "1"

尽管这个变量的不同值（0和1）在其他算法中不会显示问题，但在这种情况下，我们需要在这里进行一点转换。目标变量的类别只能取X0或X1这样的值，因此我们需要对它们进行转换。让我们为这个任务编写一些代码：

levels(train$Default) <- make.names(levels(factor(train$Default)))
levels(train$Default)
 ## [1] "X0" "X1"

这些值也在测试样本中进行了转换：

test$Default<-as.factor(test$Default)
levels(test$Default) <- make.names(levels(factor(test$Default)))
 levels(test$Default)
 ## [1] "X0" "X1"

我们将以与h2o包类似的方式创建一个网格，其中包含成本和 gamma 参数的不同值：

svmGrid <- expand.grid(sigma= 2^c(-20, -15,-10, -5, 0), C= 2^c(2:5))
print(svmGrid)
 ##              sigma  C
 ## 1  0.0000009536743  4
 ## 2  0.0000305175781  4
 ## 3  0.0009765625000  4
 ## 4  0.0312500000000  4
 ## 5  1.0000000000000  4
 ## 6  0.0000009536743  8
 ## 7  0.0000305175781  8
 ## 8  0.0009765625000  8
 ## 9  0.0312500000000  8
 ## 10 1.0000000000000  8
 ## 11 0.0000009536743 16
 ## 12 0.0000305175781 16
 ## 13 0.0009765625000 16
 ## 14 0.0312500000000 16
 ## 15 1.0000000000000 16
 ## 16 0.0000009536743 32
 ## 17 0.0000305175781 32
 ## 18 0.0009765625000 32
 ## 19 0.0312500000000 32
 ## 20 1.0000000000000 32

然后，我们将运行以下代码来训练不同的模型：

library(caret)
set.seed(1234)

SVM <- train(Default ~ ., data = train[,2:ncol(train)], 
 method = "svmRadial",
 standardize=TRUE,
 tuneGrid = svmGrid,
 metric = "ROC",
 allowParallel=TRUE,
 trControl = trainControl(method = "cv", 5, classProbs = TRUE, 
 summaryFunction=twoClassSummary))

为了训练一个SVM分类器，应该将train()方法与method参数作为svmRadial传递，这是我们选择的核。TuneGrid代表成本和 gamma 参数的不同组合。模型的准确性使用ROC指标来衡量。使用 5 折交叉验证。

模型训练完成后，我们可以如下查看结果：

print(SVM)
 ## Support Vector Machines with Radial Basis Function Kernel 
 ## 
 ## 7091 samples
 ##  108 predictor
 ##    2 classes: 'X0', 'X1' 
 ## 
 ## No pre-processing
 ## Resampling: Cross-Validated (5 fold) 
 ## Summary of sample sizes: 5674, 5673, 5672, 5672, 5673 
 ## Resampling results across tuning parameters:
 ## 
 ##   sigma            C   ROC        Sens       Spec     
 ##   0.0000009536743   4  0.9879069  0.9899383  0.8710086
 ##   0.0000009536743   8  0.9879135  0.9903822  0.8710086
 ##   0.0000009536743  16  0.9879092  0.9900863  0.8710086
 ##   0.0000009536743  32  0.9880736  0.9909741  0.8679783
 ##   0.0000305175781   4  0.9894669  0.9943777  0.8380371
 ##   0.0000305175781   8  0.9903574  0.9957094  0.8439168
 ##   0.0000305175781  16  0.9903018  0.9958573  0.8499774
 ##   0.0000305175781  32  0.9903865  0.9958572  0.8619629
 ##   0.0009765625000   4  0.9917597  0.9960052  0.8739937
 ##   0.0009765625000   8  0.9913792  0.9963011  0.8590231
 ##   0.0009765625000  16  0.9900214  0.9960050  0.8379919
 ##   0.0009765625000  32  0.9883768  0.9961529  0.8410222
 ##   0.0312500000000   4  0.9824358  0.9789899  0.9159656
 ##   0.0312500000000   8  0.9824358  0.9767682  0.8735414
 ##   0.0312500000000  16  0.9824358  0.9783977  0.8622343
 ##   0.0312500000000  32  0.9824358  0.9755850  0.9189959
 ##   1.0000000000000   4  0.4348777  1.0000000  0.0000000
 ##   1.0000000000000   8  0.4336278  1.0000000  0.0000000
 ##   1.0000000000000  16  0.4273365  1.0000000  0.0000000
 ##   1.0000000000000  32  0.4325194  1.0000000  0.0000000
 ## 
 ## ROC was used to select the optimal model using the largest value.
 ## The final values used for the model were sigma = 0.0009765625 and C     = 4.

总结来说，最佳参数如下：

SVM$bestTune
 ##          sigma C
 ## 9 0.0009765625 4

此外，我们可以通过以下方式访问具有最佳参数的模型：

SVM$finalModel
 ## Support Vector Machine object of class "ksvm" 
 ## 
 ## SV type: C-svc  (classification) 
 ##  parameter : cost C = 4 
 ## 
 ## Gaussian Radial Basis kernel function. 
 ##  Hyperparameter : sigma =  0.0009765625 
 ## 
 ## Number of Support Vectors : 252 
 ## 
 ## Objective Function Value : -619.9088 
 ## Training error : 0.007333 
 ## Probability model included.

模型的性能不是直接获得的，就像在h2o包中一样。这并不难做，但我们需要使用ROCR包：

library(ROCR)
SVM_pred<-as.numeric(unlist(predict(SVM, newdata =test, type = "prob")[2]))
pred2 <- prediction(SVM_pred,test$Default)
pred3 <- performance(pred2,"tpr","fpr")
plot(pred3, lwd=1, colorize=FALSE)
lines(x=c(0, 1), y=c(0, 1), col="red", lwd=1, lty=3);

上述代码生成了以下图表：

Gini 指数可以计算为2*ROC -1。我们可以使用Hmisc包来计算 ROC，然后计算 Gini 指数，如下所示：

library(Hmisc)
print("Gini indicator of SVM in the test sample is:")
 ## [1] "Gini indicator of SVM in the test sample is:"

print(abs(as.numeric(2*rcorr.cens(SVM_pred,test[,'Default'])[1]-1)))
 ## [1] 0.9766884

在测试样本中，Gini 指数达到 0.9766。与之前的模型一样，混淆矩阵是使用验证或测试样本计算的。为此，首先，存储训练和测试样本的概率：

summary_models_train$SVM<-as.numeric(unlist(predict(SVM, newdata =train, type = "prob")[2]))
summary_models_test$SVM<- as.numeric(unlist(predict(SVM, newdata =test, type = "prob")[2]))

现在在测试样本上计算混淆表：

aux<-summary_models_test
aux$pred<-ifelse(summary_models_test$SVM>0.04696094,1,0)
table(aux$Default,aux$pred)
 ##    
 ##        0    1
 ##   0 2828   68
 ##   1    8  135

SVM 在对银行进行分类方面做得很好。在测试样本中只有 76 家银行（68+8）被错误分类。现在，创建一个新的工作区备份：

rm(list=setdiff(ls(), c("Model_database","train","test","summary_models_train","summary_models_test","train_woe","test_woe")))
save.image("~/Data17.RData")

集成

到目前为止，我们已经训练了五个不同的模型。预测结果存储在两个数据框中，一个用于训练样本，另一个用于验证样本：

head(summary_models_train)
 ##    ID_RSSD Default          GLM RF            GBM              deep
 ## 4       37       0 0.0013554364  0 0.000005755001 0.000000018217172
 ## 21     242       0 0.0006967876  0 0.000005755001 0.000000002088871
 ## 38     279       0 0.0028306028  0 0.000005240935 0.000003555978680
 ## 52     354       0 0.0013898732  0 0.000005707480 0.000000782777042
 ## 78     457       0 0.0021731695  0 0.000005755001 0.000000012535539
 ## 81     505       0 0.0011344433  0 0.000005461855 0.000000012267744
 ##             SVM
 ## 4  0.0006227083
 ## 21 0.0002813123
 ## 38 0.0010763298
 ## 52 0.0009740568
 ## 78 0.0021555739
 ## 81 0.0005557417

让我们总结一下先前训练的模型的准确性。首先，将使用基尼指数计算每个分类器的预测能力。以下代码计算了训练样本和验证样本的基尼指数：

gini_models<-as.data.frame(names(summary_models_train[,3:ncol(summary_models_train)]))
colnames(gini_models)<-"Char"

for (i in 3:ncol(summary_models_train))
{

   gini_models$Gini_train[i-2]<-(abs(as.numeric(2*rcorr.cens(summary_models_train[,i],summary_models_train$Default)[1]-1)))

   gini_models$Gini_test[i-2]<-(abs(as.numeric(2*rcorr.cens(summary_models_test[,i],summary_models_test$Default)[1]-1)))

}

结果存储在一个名为 gini_models 的数据框中。训练样本和测试样本之间预测能力的差异也被计算：

gini_models$var_train_test<-(gini_models$Gini_train-gini_models$Gini_test)/gini_models$Gini_train
print(gini_models)

 ##   Char Gini_train Gini_test var_train_test
 ## 1  GLM  0.9906977 0.9748967     0.01594943
 ## 2   RF  1.0000000 0.9764276     0.02357242
 ## 3  GBM  1.0000000 0.9754665     0.02453348
 ## 4 deep  0.9855324 0.9589837     0.02693848
 ## 5  SVM  0.9920815 0.9766884     0.01551595

模型之间并没有太多显著的差异。在测试样本中，SVM 是预测能力最高的模型。另一方面，深度学习模型获得了最差的结果。

这些结果表明，从当前的财务报表中找到将在不到一年内失败的银行并不非常困难，这正是我们定义的目标变量。

我们还可以看到每个模型的预测能力，这取决于正确分类的银行数量：

decisions_train <- summary_models_train

decisions_test <- summary_models_test

现在，让我们创建新的数据框，其中银行根据预测概率被分类为清偿能力银行或非清偿能力银行，就像我们对每个模型所做的那样：

for (m in 3:ncol(decisions_train))
{

   decisions_train[,m]<-ifelse(decisions_train[,m]>0.04696094,1,0)

   decisions_test[,m]<-ifelse(decisions_test[,m]>0.04696094,1,0)

 }

现在，创建了一个函数，用于计算银行被正确和非正确分类的数量：

accuracy_function <- function(dataframe, observed, predicted)
{
 bads<-sum(as.numeric(as.character(dataframe[,observed])))
  goods<-nrow(dataframe)-bads
   y <- as.vector(table(dataframe[,predicted], dataframe[,observed]))
   names(y) <- c("TN", "FP", "FN", "TP")
  return(y)
 }

通过运行前面的函数，我们将看到每个模型的性能摘要。首先，该函数应用于训练样本：

print("Accuracy GLM model:")
 ## [1] "Accuracy GLM model:"
accuracy_function(decisions_train,"Default","GLM")
 ##   TN   FP   FN   TP 
 ## 6584  174    9  324

print("Accuracy RF model:")
 ## [1] "Accuracy RF model:"
accuracy_function(decisions_train,"Default","RF")
 ##   TN   FP   FN   TP 
 ## 6608  150    0  333

print("Accuracy GBM model:")
 ## [1] "Accuracy GBM model:"
accuracy_function(decisions_train,"Default","GBM")
 ##   TN   FP   FN   TP 
 ## 6758    0    0  333

print("Accuracy deep model:")
 ## [1] "Accuracy deep model:"
accuracy_function(decisions_train,"Default","deep")
 ##   TN   FP   FN   TP 
 ## 6747   11  104  229

print("Accuracy SVM model:")
 ## [1] "Accuracy SVM model:"
accuracy_function(decisions_train,"Default","SVM")
 ##   TN   FP   FN   TP 
 ## 6614  144    7  326

然后，我们可以看到测试样本中不同模型的结果：

print("Accuracy GLM model:")
 ## [1] "Accuracy GLM model:"
accuracy_function(decisions_test,"Default","GLM")
 ##   TN   FP   FN   TP
 ## 2818   78    8  135

print("Accuracy RF model:")
 ## [1] "Accuracy RF model:"
accuracy_function(decisions_test,"Default","RF")
 ##   TN   FP   FN   TP
 ## 2753  143    5  138

print("Accuracy GBM model:")
 ## [1] "Accuracy GBM model:"
accuracy_function(decisions_test,"Default","GBM")
 ##   TN   FP   FN   TP
 ## 2876   20   15  128

print("Accuracy deep model:")
 ## [1] "Accuracy deep model:"
accuracy_function(decisions_test,"Default","deep")
 ##   TN   FP   FN   TP
 ## 2886   10   61   82

print("Accuracy SVM model:")
 ## [1] "Accuracy SVM model:"
accuracy_function(decisions_test,"Default","SVM")
 ##   TN   FP   FN   TP
 ## 2828   68    8  135

根据在测试样本上测量的表格，RF 是最准确的失败银行分类器，但它也将 138 家清偿能力银行错误分类为失败，提供了错误的警报。

不同模型的结果是相关的：

correlations<-cor(summary_models_train[,3:ncol(summary_models_train)], use="pairwise", method="pearson")

print(correlations)
 ##            GLM        RF       GBM      deep       SVM
 ## GLM  1.0000000 0.9616688 0.9270350 0.8010252 0.9910695
 ## RF   0.9616688 1.0000000 0.9876728 0.7603979 0.9719735
 ## GBM  0.9270350 0.9876728 1.0000000 0.7283464 0.9457436
 ## deep 0.8010252 0.7603979 0.7283464 1.0000000 0.7879191
 ## SVM  0.9910695 0.9719735 0.9457436 0.7879191 1.0000000

将不同模型的结果结合起来以获得更好的模型可能很有趣。在这里，集成（Ensemble）的概念派上用场。集成是一种将不同的算法组合起来以创建更稳健模型的技巧。这个组合模型包含了所有基学习器的预测。结果模型将比单独运行模型时的准确性更高。实际上，我们之前开发的一些模型是集成模型，例如；随机森林或梯度提升机（GBM）。创建集成时有很多选项。在本节中，我们将探讨从最简单到更复杂的不同替代方案。

平均模型

这简单定义为取模型预测的平均值：

summary_models_test$avg<-(summary_models_test$GLM + summary_models_test$RF + summary_models_test$GBM + summary_models_test$deep + summary_models_test$SVM)/5

因此，银行最终失败的概率将被计算为前五个模型失败概率的简单平均值。

这个简单集成模型的预测能力如下：

abs(as.numeric(2*rcorr.cens(summary_models_test[,"avg"],summary_models_test$Default)[1]-1))
 ## [1] 0.9771665

我们可以创建以下混淆矩阵：

aux<-summary_models_test
aux$pred<-ifelse(summary_models_test$avg>0.04696094,1,0)
table(aux$Default,aux$pred) 
 ##        0    1
 ##   0 2834   62
 ##   1    7  136

这个组合模型只错误地将7家破产银行和62家非破产银行分类。显然，这个平均模型比所有单个模型的表现都要好。

为了增加一些保守性，我们可能会认为更好的方法是从不同的模型中分配最高的失败概率。然而，这种方法不太可能成功，因为我们之前观察到随机森林对某些银行产生了误报。

多数投票

这定义为在预测分类问题的结果时，选择具有最大投票的预测。首先，我们需要为每个模型分配一个投票。这一步已经在decisions_test数据框中完成。如果五个模型中有三个将银行分类为非清偿，则该银行将被分类为非清偿。让我们看看这种方法的结果：

decisions_test$votes<-rowSums(decisions_test[,3:7])
decisions_test$majority_vote<-ifelse(decisions_test$votes>2,1,0) 

table(decisions_test$Default,decisions_test$majority_vote)
 ##        0    1
 ##   0 2844   52
 ##   1    8  135

结果似乎不如单个模型或考虑平均概率的集成模型好：

rm(list=setdiff(ls(), c("Model_database","train","test","summary_models_train","summary_models_test","train_woe","test_woe","decisions_train","decisions_test")))
save.image("~/Data18.RData")
rm(list=ls())

模型模型

这涉及到使用另一个机器学习模型（如 Lasso、GBM 或随机森林）结合模型的单个输出（如随机森林或 SVM）。集成模型的顶层可以是任何模型，即使底层使用了相同的技巧（如随机森林）。最复杂的算法（如随机森林、梯度提升、SVM 等）并不总是比简单的算法（如树或逻辑回归）表现更好。

在这种情况下，不会训练额外的示例和算法，因为之前的结果已经足够。

自动机器学习

现在我们已经学会了如何开发一个强大的模型来预测银行破产，我们将测试一个最终选项来开发不同的模型。具体来说，我们将尝试自动机器学习（autoML），它包含在h2o包中。我们执行的过程是自动完成的，无需任何先验知识，通过autoML函数构建许多模型并找到最佳模型。此函数通过尝试不同的参数网格来训练不同的模型。此外，堆叠集成或基于先前训练模型的模型也被训练，以找到更准确或更具预测性的模型。

在我看来，在启动任何模型之前使用这个函数强烈推荐，以获得一个参考起点的基本概念。使用自动方法，我们可以评估最可靠的算法、最重要的潜在变量或我们可能获得的准确度参考。

为了测试这个函数，我们将加载一个先前的工作空间：

load("~/Data12.RData")

Data12.RData 包含在启动任何模型之前的训练和测试样本。

我们还需要加载h2o包。此外，在h2o空间中创建的所有对象都将被删除：

library(h2o)
h2o.init()
h2o.removeAll()

标准化变量

在先前的模型中，我们固定了一个参数来标准化数据。然而，这个选项在autoML函数中不可用。因此，变量将首先进行标准化。列将具有零均值和单位方差。我们需要标准化变量，因为否则结果将会有主导变量，这些变量似乎比其他属性具有更高的方差，这是由于它们的规模造成的。

标准化使用caret包完成。首先，我们选择要标准化的数值列的名称：

library(caret)
features <- setdiff(names(train), c("ID_RSSD","Default"))

变量使用preProcess函数进行转换：

pre_process <- preProcess(x = train[, features], 
                method = c( "center", "scale"))

之前的功能存储了在任意数据集上进行标准化所需的参数。使用predict函数，我们可以实际应用这种转换。训练和测试样本都必须进行转换：

# apply to both training & test
train <- cbind(train[,"Default"],predict(pre_process, train[, features]))
test <- cbind(test[,"Default"],predict(pre_process, test[, features]))

colnames(train)[1]<-"Default"
colnames(test)[1]<-"Default"

现在，我们准备创建不同的模型。训练和测试样本被转换为h2o表：

train <- as.h2o(train)
test <- as.h2o(test)

我们需要目标和预测变量的名称：

y <- "Default"
x <- setdiff(names(train), y)

在所有这些基本预处理步骤之后，是时候创建一个模型了。h2o.automl函数实现了 autoML 方法。以下参数是必需的：

x：预测变量的名称
y：目标列名称
training_frame：用于创建模型的训练数据集
leaderboard_frame：h2o用于确保模型不会过度拟合数据的验证数据集

有更多参数，但前面的列表包含了最小要求。也可以排除一些算法，例如。在我们的案例中，我们将固定要训练的模型的最大数量和 AUC 作为停止指标标准。

让我们训练一些模型：

AML_models <- h2o.automl(y = y, x = x,
                   training_frame = train,
                   max_models = 10,stopping_metric ="AUC",
                   seed = 1234,sort_metric ="AUC")

我们可以如下访问训练模型的Leaderboard：

Leaderboard <- AML_models@leaderboard
print(Leaderboard)
 ##                                                model_id       auc
 ## 1             GBM_grid_0_AutoML_20190105_000223_model_4 0.9945125
 ## 2 StackedEnsemble_BestOfFamily_0_AutoML_20190105_000223 0.9943324
 ## 3    StackedEnsemble_AllModels_0_AutoML_20190105_000223 0.9942727
 ## 4             GLM_grid_0_AutoML_20190105_000223_model_0 0.9941941
 ## 5             GBM_grid_0_AutoML_20190105_000223_model_1 0.9930208
 ## 6             GBM_grid_0_AutoML_20190105_000223_model_5 0.9926648
 ##      logloss mean_per_class_error       rmse         mse
 ## 1 0.03801166           0.04984862 0.09966934 0.009933978
 ## 2 0.03566530           0.03747844 0.09228175 0.008515921
 ## 3 0.03589846           0.03929486 0.09251204 0.008558478
 ## 4 0.03026294           0.05200978 0.08904775 0.007929502
 ## 5 0.03664414           0.06546054 0.09659713 0.009331005
 ## 6 0.13078645           0.08500441 0.18747430 0.035146615
 ## 
 ## [12 rows x 6 columns]

根据Leaderboard，如果考虑准确度，则梯度提升模型是最佳的。让我们获取这个提升模型的预测：

leader_model <- AML_models@leaderpred_test <- as.data.frame(h2o.predict(object = leader_model, newdata = test))

可以使用以下代码打印出模型的完整细节（由于长度原因，这些结果没有打印在这本书中）：

print(leader_model)

同样，也可以分析堆叠模型中单个模型的重要性。让我们看看最佳模型在测试样本上的准确度：

head(pred_test)
 ##   predict        p0          p1
 ## 1       0 0.9977300 0.002270014
 ## 2       0 0.9977240 0.002275971
 ## 3       0 0.9819248 0.018075249
 ## 4       0 0.9975793 0.002420683
 ## 5       0 0.9977238 0.002276235
 ## 6       0 0.9977240 0.002276009

重要的是要记住，预测列是根据模型预测的类别，但我们需要考虑破产预测概率中的 50%作为阈值。

与先前的算法一样，我们应该在我们的样本中定义观察到的违约率：

pred_test$predict<-ifelse(pred_test$p1>0.04696094,1,0)

现在，我们将添加观察到的类别，然后计算准确度表：

pred_test<-cbind(as.data.frame(test[,"Default"]),pred_test)
table(pred_test$Default,pred_test$predict)
 ##        0    1
 ##   0 2810   86
 ##   1    6  137

自动模型获得了非常好的性能，但略逊于 SVM 模型。类似于h2o包中的先前模型，该模型可以保存以供将来使用：

h2o.saveModel(leader_model, path = "AML_model")

摘要

在本章中，我们使用了不同的模型和算法来尝试优化我们的模型。所有算法都获得了良好的结果。在其他问题中可能不会是这样。你可以在你的问题中尝试使用不同的算法，并测试最佳参数组合以解决你的特定问题。结合不同的算法或集成也可能是一个不错的选择。

在下一章中，我们将继续探讨其他实际问题——特别是，欧洲国家经济失衡的数据可视化。