使用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))
不妨查看AUC值。
Area_Under_the_Curve
图5
此外，KNN拥有最大的AUC值(0.847)。