如何使用5种机器学习算法对罕见事件进行分类

机器学习是数据科学界的王冠，而监督学习是机器学习界这顶王冠上的宝石。

背景

几年前《哈佛商业评论》发表过一篇题为《数据科学家：21世纪最性感的工作》的文章。文章发表后，数据科学系或统计系备受大学生追捧，沉闷的数据科学家头回被认为很性感。

对一些行业而言，数据科学家已改变了公司结构，将许多决策交给了一线员工。能够从数据获得实用的业务洞察力从未如此容易。

据吴恩达称，监督学习算法为业界贡献了大部分价值。

监督学习为什么创造如此大的业务价值不容怀疑。银行用它来检测信用卡欺诈，交易员根据模型做出购买决定，工厂对生产线进行过滤以查找有缺陷的零部件。

这些业务场景有两个共同的特征：

二进制结果：欺诈vs不欺诈，购买vs不购买，有缺陷的vs没有缺陷。
不平均的数据分布：一个多数组vs一个少数组。

正如吴恩达最近指出，小数据、稳健性和人为因素是AI项目取得成功的三大障碍。在某种程度上，一个少数组方面的罕见事件问题也是一个小数据问题：机器学习算法从多数组学到更多信息，很容易对小数据组错误分类。

下面是几个事关重大的问题：

对于这些罕见事件，哪种机器学习方法性能更好?
什么度量指标?
有何美中不足?

本文试图通过运用5种机器学习方法处理实际数据集来回答上述问题，附有完整的R实现代码。

有关完整描述和原始数据集，请参阅原始数据集：https://archive.ics.uci.edu/ml/datasets/bank+marketing;有关完整的R代码，请查看我的Github：https://github.com/LeihuaYe/Machine-Learning-Classification-for-Imbalanced-Data。

业务问题

葡萄牙一家银行在实施一项新银行服务(定期存款)的营销策略，想知道哪些类型的客户已订购该服务，以便银行可以在将来调整营销策略，锁定特定人群。数据科学家与销售和营销团队合作，提出了统计解决方案，以识别未来订户。

R实现

以下面是模型选择流程和R实现。

1.导入、数据清理和探索性数据分析

不妨加载并清理原始数据集。

####load the dataset
banking=read.csv(“bank-additional-full.csv”,sep =”;”,header=T)##check for missing data and make sure no missing data
banking[!complete.cases(banking),]#re-code qualitative (factor) variables into numeric
banking$job= recode(banking$job, “‘admin.’=1;’blue-collar’=2;’entrepreneur’=3;’housemaid’=4;’management’=5;’retired’=6;’self-employed’=7;’services’=8;’student’=9;’technician’=10;’unemployed’=11;’unknown’=12”)#recode variable again
banking$marital = recode(banking$marital, “‘divorced’=1;’married’=2;’single’=3;’unknown’=4”)banking$education = recode(banking$education, “‘basic.4y’=1;’basic.6y’=2;’basic.9y’=3;’high.school’=4;’illiterate’=5;’professional.course’=6;’university.degree’=7;’unknown’=8”)banking$default = recode(banking$default, “‘no’=1;’yes’=2;’unknown’=3”)banking$housing = recode(banking$housing, “‘no’=1;’yes’=2;’unknown’=3”)banking$loan = recode(banking$loan, “‘no’=1;’yes’=2;’unknown’=3”)
banking$contact = recode(banking$loan, “‘cellular’=1;’telephone’=2;”)banking$month = recode(banking$month, “‘mar’=1;’apr’=2;’may’=3;’jun’=4;’jul’=5;’aug’=6;’sep’=7;’oct’=8;’nov’=9;’dec’=10”)banking$day_of_week = recode(banking$day_of_week, “‘mon’=1;’tue’=2;’wed’=3;’thu’=4;’fri’=5;”)banking$poutcome = recode(banking$poutcome, “‘failure’=1;’nonexistent’=2;’success’=3;”)#remove variable “pdays”, b/c it has no variation
banking$pdays=NULL #remove variable “pdays”, b/c itis collinear with the DV
banking$duration=NULL

清理原始数据似乎很乏味，因为我们要为缺失的变量重新编码，并将定性变量转换成定量变量。清理实际数据要花更长的时间。有言道“数据科学家花80%的时间来清理数据、花20%的时间来构建模型。”

下一步，不妨探究结果变量的分布。

#EDA of the DV
plot(banking$y,main="Plot 1: Distribution of Dependent Variable")

图1

由此可见，相关变量(服务订购)并不均匀分布，“No”多过“Yes”。分布不平衡应该会发出一些警告信号，因为数据分布影响最终的统计模型。它很容易使用多数范例(majority case)开发的模型对少数范例(minority case)错误分类。

2. 数据分割

下一步，不妨将数据集分割成两部分：训练集和测试集。通常而言，我们坚持80–20分割：80%是训练集，20%是测试集。如果是时间序列数据，我们基于90%的数据训练模型，将剩余10%的数据作为测试数据集。

#split the dataset into training and test sets randomly
set.seed(1)#set seed so as to generate the same value each time we run the code#create an index to split the data: 80% training and 20% test
index = round(nrow(banking)*0.2,digits=0)#sample randomly throughout the dataset and keep the total number equal to the value of index
test.indices = sample(1:nrow(banking), index)#80% training set
banking.train=banking[-test.indices,] #20% test set
banking.test=banking[test.indices,] #Select the training set except the DV
YTrain = banking.train$y
XTrain = banking.train %>% select(-y)# Select the test set except the DV
YTest = banking.test$y
XTest = banking.test %>% select(-y)

这里，不妨创建一个空的跟踪记录。

records = matrix(NA, nrow=5, ncol=2)
colnames(records) <- c(“train.error”,”test.error”)
rownames(records) <- c(“Logistic”,”Tree”,”KNN”,”Random Forests”,”SVM”)

3. 训练模型

我们在这一节定义一个新的函数(calc_error_rate)，运用它计算每个机器学习模型的训练和测试误差。

calc_error_rate <- function(predicted.value, true.value)
{return(mean(true.value!=predicted.value))}

如果预测的标签与实际值不符，该函数就计算比率。

#1 逻辑回归模型

想了解逻辑模型的简介，不妨看看这两篇文章：《机器学习101》(https://towardsdatascience.com/machine-learning-101-predicting-drug-use-using-logistic-regression-in-r-769be90eb03d)和《机器学习102》(https://towardsdatascience.com/machine-learning-102-logistic-regression-with-polynomial-features-98a208688c17)。

不妨添加一个逻辑模型，包括结果变量以外的所有其他变量。由于结果是二进制的，我们将模型设置为二项分布(“family-binomial”)。

glm.fit = glm(y ~ age+factor(job)+factor(marital)+factor(education)+factor(default)+factor(housing)+factor(loan)+factor(contact)+factor(month)+factor(day_of_week)+campaign+previous+factor(poutcome)+emp.var.rate+cons.price.idx+cons.conf.idx+euribor3m+nr.employed, data=banking.train, family=binomial)

下一步是获得训练误差。由于我们预测结果的类型并采用多数规则，于是将类型设置为响应式：如果先验概率超过或等于0.5，我们预测结果为yes，否则是no。

prob.training = predict(glm.fit,type=”response”)banking.train_glm = banking.train %>% #select all rows of the train
mutate(predicted.value=as.factor(ifelse(prob.training<=0.5, “no”, “yes”))) #create a new variable using mutate and set a majority rule using ifelse# get the training error
logit_traing_error <- calc_error_rate(predicted.value=banking.train_glm$predicted.value, true.value=YTrain)# get the test error of the logistic model
prob.test = predict(glm.fit,banking.test,type=”response”)banking.test_glm = banking.test %>% # select rows
mutate(predicted.value2=as.factor(ifelse(prob.test<=0.5, “no”, “yes”))) # set ruleslogit_test_error <- calc_error_rate(predicted.value=banking.test_glm$predicted.value2, true.value=YTest)# write down the training and test errors of the logistic model
records[1,] <- c(logit_traing_error,logit_test_error)#write into the first row

#2 决策树

若是决策树，我们遵循交叉验证，以识别最佳的分割节点。想大致了解决策树，请参阅此文：https://towardsdatascience.com/decision-trees-in-machine-learning-641b9c4e8052。

# finding the best nodes
# the total number of rows
nobs = nrow(banking.train)#build a DT model;
#please refer to this document (https://www.datacamp.com/community/tutorials/decision-trees-R) for constructing a DT model
bank_tree = tree(y~., data= banking.train,na.action = na.pass,
control = tree.control(nobs , mincut =2, minsize = 10, mindev = 1e-3))#cross validation to prune the tree
set.seed(3)
cv = cv.tree(bank_tree,FUN=prune.misclass, K=10)
cv#identify the best cv
best.size.cv = cv$size[which.min(cv$dev)]
best.size.cv#best = 3bank_tree.pruned<-prune.misclass(bank_tree, best=3)
summary(bank_tree.pruned)

交叉验证的最佳大小是3。

# Training and test errors of bank_tree.pruned
pred_train = predict(bank_tree.pruned, banking.train, type=”class”)
pred_test = predict(bank_tree.pruned, banking.test, type=”class”)# training error
DT_training_error <- calc_error_rate(predicted.value=pred_train, true.value=YTrain)# test error
DT_test_error <- calc_error_rate(predicted.value=pred_test, true.value=YTest)# write down the errors
records[2,] <- c(DT_training_error,DT_test_error)

#3 K最近邻(KNN)

作为一种非参数方法，KNN不需要任何分布的先验知识。简而言之，KNN将k个数量的最近邻分配给相关的单元。

想大致了解，不妨参阅这篇文章《R中的K最近邻入门指南：从菜鸟到高手》：https://towardsdatascience.com/beginners-guide-to-k-nearest-neighbors-in-r-from-zero-to-hero-d92cd4074bdb。想详细了解交叉验证和do.chunk函数，请参阅此文：https://towardsdatascience.com/beginners-guide-to-k-nearest-neighbors-in-r-from-zero-to-hero-d92cd4074bdb。

使用交叉验证，我们发现当k = 20时交叉验证误差最小。

nfold = 10
set.seed(1)# cut() divides the range into several intervals
folds = seq.int(nrow(banking.train)) %>%
cut(breaks = nfold, labels=FALSE) %>%
sampledo.chunk <- function(chunkid, folddef, Xdat, Ydat, k){
train = (folddef!=chunkid)# training indexXtr = Xdat[train,] # training set by the indexYtr = Ydat[train] # true label in training setXvl = Xdat[!train,] # test setYvl = Ydat[!train] # true label in test setpredYtr = knn(train = Xtr, test = Xtr, cl = Ytr, k = k) # predict training labelspredYvl = knn(train = Xtr, test = Xvl, cl = Ytr, k = k) # predict test labelsdata.frame(fold =chunkid, # k folds
train.error = calc_error_rate(predYtr, Ytr),#training error per fold
val.error = calc_error_rate(predYvl, Yvl)) # test error per fold
}# set error.folds to save validation errors
error.folds=NULL# create a sequence of data with an interval of 10
kvec = c(1, seq(10, 50, length.out=5))set.seed(1)for (j in kvec){
tmp = ldply(1:nfold, do.chunk, # apply do.function to each fold
folddef=folds, Xdat=XTrain, Ydat=YTrain, k=j) # required arguments
tmp$neighbors = j # track each value of neighbors
error.folds = rbind(error.folds, tmp) # combine the results
}#melt() in the package reshape2 melts wide-format data into long-format data
errors = melt(error.folds, id.vars=c(“fold”,”neighbors”), value.name= “error”)

随后，不妨找到尽量减少验证误差的最佳K数。

val.error.means = errors %>%
filter(variable== “val.error” ) %>%
group_by(neighbors, variable) %>%
summarise_each(funs(mean), error) %>%
ungroup() %>%
filter(error==min(error))#the best number of neighbors =20
numneighbor = max(val.error.means$neighbors)
numneighbor## [20]

遵循同一步，我们查找训练误差和测试误差。

#training error
set.seed(20)
pred.YTtrain = knn(train=XTrain, test=XTrain, cl=YTrain, k=20)
knn_traing_error <- calc_error_rate(predicted.value=pred.YTtrain, true.value=YTrain)#test error =0.095set.seed(20)
pred.YTest = knn(train=XTrain, test=XTest, cl=YTrain, k=20)
knn_test_error <- calc_error_rate(predicted.value=pred.YTest, true.value=YTest)records[3,] <- c(knn_traing_error,knn_test_error)

#4 随机森林

我们遵循构建随机森林模型的标准步骤。想大致了解随机森林，参阅此文：https://towardsdatascience.com/understanding-random-forest-58381e0602d2。

# build a RF model with default settings
set.seed(1)
RF_banking_train = randomForest(y ~ ., data=banking.train, importance=TRUE)# predicting outcome classes using training and test sets
pred_train_RF = predict(RF_banking_train, banking.train, type=”class”)pred_test_RF = predict(RF_banking_train, banking.test, type=”class”)# training error
RF_training_error <- calc_error_rate(predicted.value=pred_train_RF, true.value=YTrain)# test error
RF_test_error <- calc_error_rate(predicted.value=pred_test_RF, true.value=YTest)records[4,] <- c(RF_training_error,RF_test_error)

#5 支持向量机

同样，我们遵循构建支持向量机的标准步骤。想大致了解该方法，请参阅此文：https://towardsdatascience.com/support-vector-machine-introduction-to-machine-learning-algorithms-934a444fca47。

set.seed(1)
tune.out=tune(svm, y ~., data=banking.train,
kernel=”radial”,ranges=list(cost=c(0.1,1,10)))# find the best parameters
summary(tune.out)$best.parameters# the best model
best_model = tune.out$best.modelsvm_fit=svm(y~., data=banking.train,kernel=”radial”,gamma=0.05555556,cost=1,probability=TRUE)# using training/test sets to predict outcome classes
svm_best_train = predict(svm_fit,banking.train,type=”class”)
svm_best_test = predict(svm_fit,banking.test,type=”class”)# training error
svm_training_error <- calc_error_rate(predicted.value=svm_best_train, true.value=YTrain)# test error
svm_test_error <- calc_error_rate(predicted.value=svm_best_test, true.value=YTest)records[5,] <- c(svm_training_error,svm_test_error)

4. 模型度量指标

我们已构建了遵循模型选择过程的所有机器学习模型，并获得了训练误差和测试误差。这一节将使用一些模型的度量指标选择最佳模型。

4.1 训练/测试误差

可以使用训练/测试误差找到最佳模型吗?

现在不妨看看结果。

records

图2

这里，随机森林的训练误差最小，不过其他方法有类似的测试误差。你可能注意到，训练误差和测试误差很接近，很难说清楚哪个明显胜出。

此外，分类精度(无论是训练误差还是测试误差)都不应该是高度不平衡数据集的度量指标。这是由于数据集以多数范例为主，即使随机猜测也会得出50%的准确性。更糟糕的是，高度精确的模型可能严重“处罚”少数范例。因此，不妨查看另一个度量指标：ROC曲线。

4.2受试者工作特征(ROC)曲线

ROC是一种图形表示，显示分类模型在所有分类阈值下有怎样的表现。我们更喜欢比其他分类器更快逼近1的分类器。

ROC曲线在同一个图中绘制不同阈值下的两个参数：真阳率(True Positive Rate)和假阳率(False Positive Rate)。

TPR (Recall) = TP/(TP+FN)

FPR = FP/(TN+FP)

图3

在很大程度上，ROC曲线不仅衡量分类准确度，还在TPR和FPR之间达到了很好的平衡。这是罕见事件所需要的，因为我们还想在多数范例和少数范例之间达到平衡。

# load the library
library(ROCR)#creating a tracking record
Area_Under_the_Curve = matrix(NA, nrow=5, ncol=1)
colnames(Area_Under_the_Curve) <- c(“AUC”)
rownames(Area_Under_the_Curve) <- c(“Logistic”,”Tree”,”KNN”,”Random Forests”,”SVM”)########### logistic regression ###########
# ROC
prob_test <- predict(glm.fit,banking.test,type=”response”)
pred_logit<- prediction(prob_test,banking.test$y)
performance_logit <- performance(pred_logit,measure = “tpr”, x.measure=”fpr”)########### Decision Tree ###########
# ROC
pred_DT<-predict(bank_tree.pruned, banking.test,type=”vector”)
pred_DT <- prediction(pred_DT[,2],banking.test$y)
performance_DT <- performance(pred_DT,measure = “tpr”,x.measure= “fpr”)########### KNN ###########
# ROC
knn_model = knn(train=XTrain, test=XTrain, cl=YTrain, k=20,prob=TRUE)prob <- attr(knn_model, “prob”)
prob <- 2*ifelse(knn_model == “-1”, prob,1-prob) — 1
pred_knn <- prediction(prob, YTrain)
performance_knn <- performance(pred_knn, “tpr”, “fpr”)########### Random Forests ###########
# ROC
pred_RF<-predict(RF_banking_train, banking.test,type=”prob”)
pred_class_RF <- prediction(pred_RF[,2],banking.test$y)
performance_RF <- performance(pred_class_RF,measure = “tpr”,x.measure= “fpr”)########### SVM ###########
# ROC
svm_fit_prob = predict(svm_fit,type=”prob”,newdata=banking.test,probability=TRUE)
svm_fit_prob_ROCR = prediction(attr(svm_fit_prob,”probabilities”)[,2],banking.test$y==”yes”)
performance_svm <- performance(svm_fit_prob_ROCR, “tpr”,”fpr”)

不妨绘制ROC曲线。

我们添加一条直线，以显示随机分配的概率。我们的分类器其表现胜过随机猜测，是不是?

#logit
plot(performance_logit,col=2,lwd=2,main=”ROC Curves for These Five Classification Methods”)legend(0.6, 0.6, c(‘logistic’, ‘Decision Tree’, ‘KNN’,’Random Forests’,’SVM’), 2:6)#decision tree
plot(performance_DT,col=3,lwd=2,add=TRUE)#knn
plot(performance_knn,col=4,lwd=2,add=TRUE)#RF
plot(performance_RF,col=5,lwd=2,add=TRUE)# SVM
plot(performance_svm,col=6,lwd=2,add=TRUE)abline(0,1)

图4

这里已分出胜负。

据ROC曲线显示，KNN(蓝色线)高于其他所有方法。

4.3 曲线下面积(AUC)

顾名思义，AUC是ROC曲线下的面积。它是直观的AUC曲线的数学表示。AUC给出了分类器在可能的分类阈值下性能如何的合并结果。

########### Logit ###########
auc_logit = performance(pred_logit, “auc”)@y.values
Area_Under_the_Curve[1,] <-c(as.numeric(auc_logit))########### Decision Tree ###########
auc_dt = performance(pred_DT,”auc”)@y.values
Area_Under_the_Curve[2,] <- c(as.numeric(auc_dt))########### KNN ###########
auc_knn <- performance(pred_knn,”auc”)@y.values
Area_Under_the_Curve[3,] <- c(as.numeric(auc_knn))########### Random Forests ###########
auc_RF = performance(pred_class_RF,”auc”)@y.values
Area_Under_the_Curve[4,] <- c(as.numeric(auc_RF))########### SVM ###########
auc_svm<-performance(svm_fit_prob_ROCR,”auc”)@y.values[[1]]
Area_Under_the_Curve[5,] <- c(as.numeric(auc_svm))