本文分析数据是基于这个数据集：https://www.kaggle.com/uciml/student-alcohol-consumption。
1、分析目的和数据集处理
（1）目的
本文希望以两个学生酒精消费程度变量为起点，利用问卷中获得的其他信息，建立对于学生酒精消费状况的分类预测模型。最终实现，一方面是学习和对比几个常见的机器学习分类模型，呈现几个算法的特点；另一方面是对模型的可用性进行一些讨论，例如结合本数据获得的预测模型，就是否能够能让学校管理者或学校社工基于一定的信息，能够预判哪些学生会有较高的酒精消费程度，从而据此做出干预。
本文更侧重于对于数据分析方法的讨论。因为本数据集所涉及的的问题，也是一个在社会学和流行病学领域很值得研究的问题，但囿于能力有限，本文在此方面不做详细展开。因此，本文对于数据的信效度也不做展开讨论。
（2）数据集处理
该数据集是针对一个中学的选上葡萄牙语课和数学课两个课程的学生，进行问卷调查获得的数据。两个课程分别获得数据葡萄牙语课（649条）和数学课（395条）数据。
模型建立希望基于尽可能多的数据量，希望将两个数据集合并。要达到这一目的，需要排除两门课程的学生在两个目标变量上有显著差异的可能。在此可借助单因素方差分析来确定课程（course）这一变量的不同分组，是否在两个目标变量中有明显的组间差异。
d3withrepeat <- rbind(mat, por)
d3 <- d3withrepeat %>% distinct(school,sex,age,address,famsize,Pstatus,
Medu,Fedu,Mjob,Fjob,reason,
guardian,traveltime,studytime,failures,
schoolsup,famsup,activities,
nursery,higher,internet,
romantic,famrel,freetime,goout,
Dalc,Walc,health,absences,
.keep_all = TRUE)合并两个数据集后获得1044条数据，去重后获得959条数据，即两门课程都上了的学生有85名。
接下来就是检验两门课程学生的目标变量有没有显著差异，考虑到目标变量不是一个连续性的变量，遂采用Wilcoxon秩检验：
> wilcox.test(Dalc ~ course, d3)
Wilcoxon rank sum test with continuity correction
data: Dalc by course
W = 1e+05, p-value = 0.6
alternative hypothesis: true location shift is not equal to 0
> wilcox.test(Walc ~ course, d3)
Wilcoxon rank sum test with continuity correction
data: Walc by course
W = 1e+05, p-value = 0.7
alternative hypothesis: true location shift is not equal to 0 两个目标变量的p值均大于0.05，故可以确认两门课程的学生在两个目标变量间没有显著差异，两个数据集可以合并。
2、问题背景理解与目标变量处理
数据分析开始前，需要对分析对象的背景和大致的成果、假设做一个简单描述。本数据集所关注的高酒精消费程度或酗酒行为形成原因可大致归纳为：
- （1）生理基础（此处不讨论）
- （2）社会结构：性别、年龄、种族、职业、收入、健康状况等；
- （3）社会网络状况：亲子关系状况、亲密关系状况、日常娱乐方式、社会网络支持状况等；
- （4）社会政策与文化（此处不讨论）
关于酒精消费和酗酒问题深入的研究，可以参考这两篇文章：
http://www.robinroom.net/quagmire.pdf
http://www.robinroom.net/socscien.pdf
可以看到，该数据集对于社会结构和社会网络状况两方面的覆盖已经很详细，所以本文所建立的预测模型，主要基于学生的社会结构和社会网络状况方面的信息。
为达到本文了解酒精消费程度较高学生状况的目的，首先需要明确什么程度的酒精消费是值得关注的。我们可以先看看两个目标变量工作日酒精消费程度（Dalc）和周末酒精消费程度(Walc)的分布情况。
D.b <- d3 %>% ggplot(aes(x=Dalc))+
geom_bar(stat = "count")
W.b <- d3 %>% ggplot(aes(x=Walc))+
geom_bar(stat = "count")
grid.arrange(D.b, W.b, nrow = 2)
从上图可看出，学生在工作日酒精消费程度大都处于较低的程度，而在周末则是呈现一个由最低程度向最高程度梯度较缓的下降。我们可以据此分别设定一个新变量：酒精消费预警程度，具体是将工作日酒精消费程度为1的设定为『Safe』，2为『General』，3~5为『Risk』；将周末酒精消费程度为1~2的设定为『Safe』，3~4的为『General』，5为『Risk』。
d3$dr.dalc <- mapvalues(d3$Dalc,
from = 1:5,
to = c("Safe", "General", "Risk", "Risk", "Risk"))
d3$dr.walc <- mapvalues(d3$Walc,
from = 1:5,
to = c("Safe", "Safe", "General", "General", "Risk"))
接下来是将工作日和周末的酒精消费预警程度进行合并，形成总体酒精消费预警程度。具体是将工作日和周末中都是『Safe』的，取值为『Safe』；将二者中有一个是『General』，同时没有『Risk』的，取值为『General』；将二者中任意一个有『Risk』的，取值为『Risk』。
d3$All.alc[d3$dr.dalc == "Safe"& d3$dr.walc == "Safe"] <- "Safe"
d3$All.alc[d3$dr.dalc != "Risk" & d3$dr.walc != "Risk" &
d3$dr.dalc == "General"|d3$dr.walc == "General"] <- "General"
d3$All.alc[d3$dr.dalc == "Risk"| d3$dr.walc == "Risk"] <- "Risk"
d3$fAll.alc <- factor(d3$All.alc)我们最终获得的总体酒精消费预警程度如下，将136人标记为『Risk』的程度。
> summary(d3$fAll.alc)
Safe General Risk
514 309 136 3、数据探索
响应变量和预测变量间的交叉分布图，是理解数据集最为直观的方法，同时也因为其信息表达较为模糊的精确性，也使得需要进一步的求证。
本文的数据探索，将对响应变量的不同类别在每个预测变量值下的对应关系，进行逐一地筛查，并展示响应变量不同分类在其不同变量值中，所占比重有明显差别的变量，从而说明这些预测变量在不同预警类型下的区别。
构图所用的自编函数：
pile.f.bar <- function(x, y, xn, yn,"Freq", title = NULL){
z <- as.data.frame(table(x, y))
colnames(z) <- c(xn, yn, Freq)
bar.p <- ggplot(z, aes_string(x= xn, y = Freq, fill = yn))+
geom_bar(stat = "identity",position="stack",width=0.9)+
ggtitle(title)+
scale_fill_brewer()
return(bar.p)
}
#reliable on ggplot2
#example：
bar.sex <-pile.f.bar(d3$sex, d3$fAll.alc,"sex","All.alc","Freq")
bar.age <-pile.f.bar(d3$age, d3$fAll.alc,"age","All.alc","Freq")
bar.x <- pile.f.bar(d3$x, d3$fAll.alc,"x","All.alc","Freq")
......（1）个人特征
从上图可看出，达到『Risk』程度的男性学生，其占总体比例明显高于女性；总体年龄段峰值集中在16，17岁两个年龄段，而『Risk』程度的学生峰值集中在17，18岁两个年龄段；不同预警程度的学生的健康状况没有明显区别；『Risk』程度学生在对更高教育的计划上，所在比例更高，但总体趋势相近。
（2）家庭关系
父母同居状况、监护人角色和家庭关系3个变量中，『Risk』程度的学生分布趋势和总体分布趋势区别较小；少于3人的家庭的达到『Risk』的程度，比之总体要高；母亲受教育程度和父亲受教育状况中达到『Risk』程度的学生分布状况与总体分布状况关系不大。
（3）家庭经济状况
不同预警程度学生父母亲的工作状况的总体趋势较为相近，不过『Risk』程度学生父母是服务行业的比例略高；额外课程支付、额外网络活动、课外辅导和网络接入状况方面，『Risk』程度的学生所占比率与总体趋势相近；家庭支持中『Risk』程度学生缺乏家庭支出比例明显高于其他类型。
（4）学习状况
在到学校的时间、学校支持力度和学生课业挂科的比例中达到『Risk』程度的学生分布比例与总体相近；达到『Risk』程度的学生学习时长分布状况明显与总体区别，每周小于2小时总体占比明显大于其他类型；课后*时间也稍多于其他程度的占比。
（5）同辈关系
在是否有亲密关系上，达到『Risk』程度的学生分布比例与总体相近；与朋友外出的分布上，『Risk』程度学生的分布与总体明显区别，显示出更频繁的与朋友外出。
（6）学习成绩
获得三个时段成绩，并构建箱线图。
d3$avgG <- rowMeans(cbind(d3$G1,d3$G2,d3$G3))
d3 %>% ggplot(aes(x= All.alc, y = avgG, group = All.alc))+
geom_boxplot(fill = "darkblue", colour = "gray",
outlier.colour = "red", outlier.shape = 1)
相比于其他程度的学生，达到『Risk』的学生成绩分布更偏低，该程度内部差异性较小。
（7）描述性分析汇总
通过比较达到『Risk』程度的学生分布情况与总体分布状况，对比该程度的学生与其他程度学生占总比的状况，可以将变量分为以下三类：
- a、『Risk』程度的学生分布情况与总体分布状况有明显区别，且呈现特定的分布趋势，其包括：sex、age、Fjob、famsup、studytime、freetime、goout。
- b、『Risk』程度的学生分布情况与总体分布状况区别较小，呈现大致类似的分布趋势，其包括：high、health、Mjob、paid、activities、nursery、internet、traveltime、failsure、schoolsup、romantic。
- c、『Risk』程度的学生分布情况与总体分布状况有区别，但没有特定分布趋势，其包括：Medu、Fedu。
4、模型建立与评估
(1)数据准备和建立交叉验证
#build the model database
d4 <- d3[,1:34]
d4$fAll.alc <- d3$fAll.alc
dwalc.v <- names(d4) %in% c("Dalc","Walc")
d4 <- d4[!dwalc.v]
is.factor(d4$fAll.alc)
set.seed(100)
train <- sample(nrow(d4),0.7*nrow(d4))
d4.train <- d4[train,]
d4.test <- d4[-train,]
#cross-validation
t.control<- trainControl(method="repeatedcv", number=10,
repeats=3)(2)模型建立
为便于不同模型的对比，本文尽可能的使用一致的建模函数：caret::train。对于机器学习的入门者来说，caret确实是个很棒的包，里面大部分R所涉及的算法都可以用万能的train()来开始对这个算法的认识。关于caret支持哪些算法，可以具体看这http://topepo.github.io/caret/available-models.html。
a、递归分割树（rpart）
首先是较为传统的分类决策树。rpart作为一种无参算法，对于数据的线性可分没有要求，也是是我在考虑算法时，首先想到的经典算法。
#rpart
rp.model <- train(fAll.alc~., d4.train, trControl=t.control, metric="Kappa",
method="rpart")
predi.rp<- predict(rp.model,d4.test)
cM.rp<- confusionMatrix(predi.rp,d4.test$fAll.alc)
混淆矩阵结果：
> cM.rp
Confusion Matrix and Statistics
Reference
Prediction General Risk Safe
General 20 22 8
Risk 16 9 3
Safe 58 18 134
Overall Statistics
Accuracy : 0.5659722
95% CI : (0.5065642, 0.6240087)
No Information Rate : 0.5034722
P-Value [Acc > NIR] : 0.01946895
Kappa : 0.2245056
Mcnemar's Test P-Value : 0.0000000001000841
Statistics by Class:
Class: General Class: Risk Class: Safe
Sensitivity 0.21276596 0.18367347 0.9241379
Specificity 0.84536082 0.92050209 0.4685315
Pos Pred Value 0.40000000 0.32142857 0.6380952
Neg Pred Value 0.68907563 0.84615385 0.8589744
Prevalence 0.32638889 0.17013889 0.5034722
Detection Rate 0.06944444 0.03125000 0.4652778
Detection Prevalence 0.17361111 0.09722222 0.7291667
Balanced Accuracy 0.52906339 0.55208778 0.6963347b、朴素贝叶斯（Naive Bayes）
作为一种基于概率的分类器，对较小的样本量有好的表现。另一方面，它也要求所有特征变量是独立的，本案例在此方面有满足程度较低，但作为一类较为独特的算法，仍然可以试试。
#naive Bayes
nb.model <- train(fAll.alc~., d4.train, trControl = t.control, metric="Kappa", method="nb")
predi.nb<- predict(nb.model,d4.test)
cM.nb<- confusionMatrix(predi.nb,d4.test$fAll.alc)混淆矩阵结果：
> cM.nb
Confusion Matrix and Statistics
Reference
Prediction General Risk Safe
General 40 25 42
Risk 17 17 18
Safe 37 7 85
Overall Statistics
Accuracy : 0.4930556
95% CI : (0.4339105, 0.5523449)
No Information Rate : 0.5034722
P-Value [Acc > NIR] : 0.66001449
Kappa : 0.1856371
Mcnemar's Test P-Value : 0.08281805
Statistics by Class:
Class: General Class: Risk Class: Safe
Sensitivity 0.4255319 0.34693878 0.5862069
Specificity 0.6546392 0.85355649 0.6923077
Pos Pred Value 0.3738318 0.32692308 0.6589147
Neg Pred Value 0.7016575 0.86440678 0.6226415
Prevalence 0.3263889 0.17013889 0.5034722
Detection Rate 0.1388889 0.05902778 0.2951389
Detection Prevalence 0.3715278 0.18055556 0.4479167
Balanced Accuracy 0.5400855 0.60024763 0.6392573c、支持向量机（SVM）
SVM利用了面向工程问题的核函数，能够有较高的预测率的同时，也能较好的解决模型在过度适应、局部最优和多重共线方面的问题。
#SVM
svm.model <- train(fAll.alc~., d4.train, trControl=t.control, metric="Kappa",
method="svmRadialWeights")
predi.svm<- predict(svm.model,d4.test)
cM.svm<- confusionMatrix(predi.svm,d4.test$fAll.alc)混淆矩阵结果：
> cM.svm
Confusion Matrix and Statistics
Reference
Prediction General Risk Safe
General 38 34 18
Risk 0 4 1
Safe 56 11 126
Overall Statistics
Accuracy : 0.5833333
95% CI : (0.5240449, 0.640889)
No Information Rate : 0.5034722
P-Value [Acc > NIR] : 0.003939441
Kappa : 0.2528214
Mcnemar's Test P-Value : 0.0000000000002369114
Statistics by Class:
Class: General Class: Risk Class: Safe
Sensitivity 0.4042553 0.08163265 0.8689655
Specificity 0.7319588 0.99581590 0.5314685
Pos Pred Value 0.4222222 0.80000000 0.6528497
Neg Pred Value 0.7171717 0.84098940 0.8000000
Prevalence 0.3263889 0.17013889 0.5034722
Detection Rate 0.1319444 0.01388889 0.4375000
Detection Prevalence 0.3125000 0.01736111 0.6701389
Balanced Accuracy 0.5681070 0.53872428 0.7002170d、神经网络（Neural Networks）
作为一种深度学习类的算法，神经网络能够检测因变量和自变量之间的非线性关系，同时也是一种无参算法，能够避免参数估计过程中的失误。
#neural network
nnt.model <- train(fAll.alc~., d4.train, trControl = t.control,metric="Kappa",
method = "nnet")
predi.nnt<- predict(nnt.model,d4.test)
cM.nnt<- confusionMatrix(predi.nnt,d4.test$fAll.alc)混淆矩阵结果：
> cM.nnt
Confusion Matrix and Statistics
Reference
Prediction General Risk Safe
General 32 16 27
Risk 16 22 5
Safe 46 11 113
Overall Statistics
Accuracy : 0.5798611
95% CI : (0.5205431, 0.6375187)
No Information Rate : 0.5034722
P-Value [Acc > NIR] : 0.005564622
Kappa : 0.2907992
Mcnemar's Test P-Value : 0.065929434
Statistics by Class:
Class: General Class: Risk Class: Safe
Sensitivity 0.3404255 0.44897959 0.7793103
Specificity 0.7783505 0.91213389 0.6013986
Pos Pred Value 0.4266667 0.51162791 0.6647059
Neg Pred Value 0.7089202 0.88979592 0.7288136
Prevalence 0.3263889 0.17013889 0.5034722
Detection Rate 0.1111111 0.07638889 0.3923611
Detection Prevalence 0.2604167 0.14930556 0.5902778
Balanced Accuracy 0.5593880 0.68055674 0.6903545e、随机森林（Random Forest）
RF作为一种集成算法，其本身包含多个决策树的分类器，其结果理论上也会比其他单一算法要好。
#rf
rf.model <- train(fAll.alc~., d4.train, trControl = t.control,metric="Kappa",
method = "rf")
predi.rf<- predict(rf.model,d4.test)
cM.rf<- confusionMatrix(predi.rf,d4.test$fAll.alc)混淆矩阵结果：
> cM.rf
Confusion Matrix and Statistics
Reference
Prediction General Risk Safe
General 51 14 16
Risk 5 24 1
Safe 38 11 128
Overall Statistics
Accuracy : 0.7048611
95% CI : (0.6485267, 0.7569231)
No Information Rate : 0.5034722
P-Value [Acc > NIR] : 0.000000000003057811
Kappa : 0.4920635
Mcnemar's Test P-Value : 0.000080553103510669
Statistics by Class:
Class: General Class: Risk Class: Safe
Sensitivity 0.5425532 0.48979592 0.8827586
Specificity 0.8453608 0.97489540 0.6573427
Pos Pred Value 0.6296296 0.80000000 0.7231638
Neg Pred Value 0.7922705 0.90310078 0.8468468
Prevalence 0.3263889 0.17013889 0.5034722
Detection Rate 0.1770833 0.08333333 0.4444444
Detection Prevalence 0.2812500 0.10416667 0.6145833
Balanced Accuracy 0.6939570 0.73234566 0.7700506（3）模型评估
作为一个多类分类案例，评价模型的好坏，单元格正确划分了概率准确率外（即测试样本中『Risk』、『General』和『Normal』被预测正确的样本所占比重）是最为直观和显要的。
同时，结合模型运用的目的，即模型能够尽可能地准确确定哪个学生有处于『Risk』状态的可能，还可以从正类样本被正确预测的概率敏感度（即处于『Risk』状态的学生被正确识别的概率），以及被预测为正类的样本中，正确的概率精确率（被预测为『Risk』状态的学生确实是的概率）另外两个指标入手。
a、准确率
> Accuracy
Accuracy
rpart 0.57
Naive Bayes 0.49
SVM 0.58
Neural Net 0.58
RF 0.70可以看到，除了朴素贝叶斯外，其他模型的准确率都处于居中的水平，尤以随机森林准确率最高（毕竟是个分类器集）。
其实现代码为：
cM.rpacc <- cM.rp$overall['Accuracy']
cM.nbacc <- cM.nb$overall['Accuracy']
cM.svmacc <- cM.svm$overall['Accuracy']
cM.nntacc <- cM.nnt$overall['Accuracy']
cM.rfacc <- cM.rf$overall['Accuracy']
Accuracy <- as.data.frame(round(cbind(cM.rpacc,cM.nbacc,cM.svmacc,
cM.nntacc,cM.rfacc),2))
colnames(Accuracy) <- c("rpart","Naive Bayes","SVM","Neural Net","RF")
Accuracy <- t(Accuracy)b、敏感度
> Sensitive
rpart Naive Bayes SVM Neural Net RF
Class: General 0.21 0.43 0.40 0.34 0.54
Class: Risk 0.18 0.35 0.08 0.45 0.49
Class: Safe 0.92 0.59 0.87 0.78 0.88各个模型对于Risk类别的敏感度不是很高，也只有随机森林和神经网络处于可以接受的水平。说明用当前的模型寻找到处于Risk状态的学生只有不到一半的概率。
其实现代码为：
rp.sen <- cM.rp$byClass[,1]
nb.sen <- cM.nb$byClass[,1]
svm.sen <- cM.svm$byClass[,1]
nnt.sen <- cM.nnt$byClass[,1]
rf.sen <- cM.rf$byClass[,1]
Sensitive <- as.data.frame(round(cbind(rp.sen,nb.sen,svm.sen,
nnt.sen,rf.sen),2))
colnames(Sensitive) <- c("rpart","Naive Bayes","SVM","Neural Net","RF")c、精确率
> Precision
rpart Naive Bayes SVM Neural Net RF
Class: General 0.40 0.37 0.42 0.43 0.63
Class: Risk 0.32 0.33 0.80 0.51 0.80
Class: Safe 0.64 0.66 0.65 0.66 0.72随机森林和SVM预测『Risk』类别的精确率比较高，说明这两个模型预测为『Risk』的学生，确实处于『Risk』的概率较大。
其实现代码为：
rp.prec <- cM.rp$byClass[,5]
nb.prec <- cM.nb$byClass[,5]
svm.prec <- cM.svm$byClass[,5]
nnt.prec <- cM.nnt$byClass[,5]
rf.prec <- cM.rf$byClass[,5]
Precision <- as.data.frame(round(cbind(rp.prec,nb.prec,svm.prec,
nnt.prec,rf.prec),2))
colnames(Precision) <- c("rpart","Naive Bayes","SVM","Neural Net","RF")(4)基于模型的解读
a、基于预测变量重要性的解读
plot(varImp(rf.model))
从这张预测变量的重要性图可以看出，与朋友外出时间、缺勤次数、与家人关系、考试成绩的影响力、空余时间、性别为男等几个变量，是值得关注的重要变量。
b、基于决策树的解读
rp.model2 <- rpart(fAll.alc~., d4.train)
#Prune the classication tree
##Find the minimum cross-validation error
min(rp.model2$cptable[,"xerror"])
##Locate the record with the minimum cross-validation errors
which.min(rp.model2$cptable[,"xerror"])
##Get the cost complexity parameter of the record with the minimum cross-validation errors
rp.cp = rp.model2$cptable[6,"CP"]
##Prune the tree
prune.tree <- prune(rp.model2, cp= churn.cp)
draw.tree(prune.tree, cex=0.6, nodeinfo=TRUE, col=gray(0:8 / 8))
根据以上剪枝处理后的分类树，以及之前的变量重要性表，可据此大致描绘出具有高酗酒风险的学生的特征：
- 和朋友外出大于2.5是一个重要零界点，这一变量作为首要影响变量，说明这一数据集所调查的中学生中，酗酒风险还是同辈群体的压力推动的；
- 男学生酗酒的风险明显高于女学生；
- 期末成绩低于中位值（13.5），酗酒风险明显；
- 每周缺勤超过7.5次，酗酒风险明显。
最后，结合常识可以知道，根据模型所聚焦的对象，也必然是其他方面的问题学生，模型的价值并没有明确展现。此外，模型的准确率不是很高，一方面是模型构建本身有待改进，另一方面也说明这一数据集并不是很好的建模材料。