R语言 - 基于[狗熊会]“百度搜索竞价广告创意研究”数据的定序回归、随机森林与xgboost - v0.5

继续挖坑

这次是基于狗熊会基础案例百度搜索竞价广告创意研究的定序回归与随机森林、xgboost。

————–更新记录—————

2017.12.17-v0.1 创建文档，代码占坑

2017.12.21-v0.2 更新为markdown，简单排版，增补内容，修改代码

2017.12.23-v0.3 增补随机森林oob error绘制、特征重要性、格搜索

2017.12.25-v0.4α 辛辛苦苦改了一晚上，没有保存，手贱乱点。。闷闷不乐ing。。

2017.12.25-v0.4β 耐着性子重新码了一遍，增补随机森林袋外误差细节，及随机森林其他内容

2017.12.26-v0.5 增补随机森林特征重要性与特征选择细节内容

欢迎转载，转载请注明出处http://blog.csdn.net/max_cola/article/details/78827475

简易特征工程
读取数据

高频词

摘要关键字是否飘红

摘要长度与标题长度

标题几处飘红

描述几处飘红

是否包含联系方式

定序回归
划分训练集和测试集

定序回归模型

模型筛选
此处应有AIC BIC挖坑

检验预测结果

改进方法
随机森林
随机森林建模与结果解读

混淆矩阵可视化

袋外误差 OOB ERROR

特征重要性与特征选择

使用格搜索优化超参数

xgboost

简易特征工程

读取数据

setwd("C:\\Users\\Pepsi\\Desktop\\狗熊案例：百度搜索竞价广告创意研究")
dat0 <- read.csv("百度广告创意.csv",stringsAsFactors = F)
str(dat0)

结果如下：

## 'data.frame':    4346 obs. of  5 variables:
##  $ keyword: chr  "口才培训" "口才培训" "口才培训" "口才培训" ...
##  $ ranking: int  1 2 3 4 5 1 2 3 4 5 ...
##  $ title  : chr  "北京卡耐基口才培训克服临场紧张无效全额退款." "卡耐基练好口才中国口才培训第一品牌!无效全退!" "北京口才培训-大钊训练是中国最专业的口才培训学校!" "成功卡耐基口才培训完全免费训练效果独特" ...
##  $ red    : chr  "北京 | 口才培训" "口才 | 口才培训 | 口才培训训练" "北京口才培训 | 训练" "口才培训 | 训练" ...
##  $ desc   : chr  "北京卡耐基的口才培训课程采用互动式教学克服临场紧张.费用3800地点大钟寺.北京卡耐基学校99成立领跑口才培训行业14年."| __truncated__ "卡耐基口才培训完全体验式授课多场景口才培训训练多主题口才培训缓解当众讲话紧张!口才培训教您规范体态教您重点突出老"| __truncated__ "大钊口才培训帮您克服讲话紧张提升自信心控制紧张讲话重点突出!最新的口才培训由著名大钊老师主讲机智风趣的课堂!听君"| __truncated__ "电话010-56077779卡耐基口才培训—独一无二特浸泡式口才培训体验教学!瞬间克服紧张各种场合站起能说;表达生动 说服力强!"| __truncated__ ...

同时根据第一行结果可以看到，原始数据由4346行，5列组成。

高频词

载入jieba分词

library(jiebaR)

wk <- worker(bylines = T) # 建立分词器
fw_l <- wk[dat0$desc] # 对desc进行分词
fw_l <- lapply(fw_l, function(x) x[nchar(x)>1]) # 取词长大于1
head(fw_l,2)

输出结果如下：

## [[1]]
##  [1] "北京"   "卡耐基" "口才"   "培训"   "课程"   "采用"   "互动式"
##  [8] "教学"   "克服"   "临场"   "紧张"   "费用"   "3800"   "地点"
## [15] "大钟寺" "北京"   "卡耐基" "学校"   "99"     "成立"   "领跑"
## [22] "口才"   "培训"   "行业"   "14"     "课程"   "采用"   "互动式"
## [29] "教学"   "无效"   "全额"   "退款"
##
## [[2]]
##  [1] "卡耐基"       "口才"         "培训"         "完全"
##  [5] "体验式"       "授课"         "场景"         "口才"
##  [9] "培训"         "训练"         "主题"         "口才"
## [13] "培训"         "缓解"         "当众"         "讲话"
## [17] "紧张"         "口才"         "培训"         "规范"
## [21] "体态"         "重点"         "突出"         "老师"
## [25] "众家"         "所长"         "口碑"         "超好"
## [29] "无效"         "全退"         "4000-777-069"

查看词频分布

w_freq <- freq(unlist(fw_l))
w_freq <- w_freq[order(-w_freq$freq),]
plot(w_freq$freq)

查看累计词频分布：

s_w <- sum(w_freq$freq)
for (i in 1:length(w_freq$freq)) {
w_freq$acc_freq[i] <- ifelse(i==1,w_freq$freq[1],w_freq$acc_freq[i-1]+w_freq$freq[i])
}
w_freq$acc_freq <- w_freq$acc_freq/s_w
plot(seq(0,1,length.out = nrow(w_freq)),w_freq$acc_freq,type = "l")

可以看到显著的二八现象

查看前20个关键词词频占比

w_freq$acc_freq[20]

前20个关键词占37.4%的比例，因此将 top20 作为高频词

## [1] 0.3740857

统计每条desc包含多少高频词

top20 <- w_freq[1:20,]
library(stringr)
dat_tmp <- dat0
dat_tmp$fw_num <- unlist(lapply(dat0$desc,function(x) sum(str_count(x,top20$char))))
summary(dat_tmp$fw_num)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##   0.000   2.000   6.000   6.514  10.000  23.000

hist(dat_tmp$fw_num)

摘要关键字是否飘红

for (i in 1:nrow(dat_tmp)) {
dat_tmp$kw_in[i] <- str_detect(dat_tmp$red[i],dat_tmp$keyword[i])
}
dat_tmp$kw_in <- 1*dat_tmp$kw_in
table(dat_tmp$kw_in)

##
##    0    1
## 3272 1074

摘要长度与标题长度

dat_tmp$des_l <- nchar(dat_tmp$desc)
summary(dat_tmp$des_l)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##   10.00   31.00   44.00   50.34   71.00   98.00

dat_tmp$tit_l <- nchar(dat_tmp$title)
summary(dat_tmp$tit_l)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##    4.00   15.00   19.00   18.41   22.00   31.00

标题几处飘红

tmp <- str_replace_all(dat_tmp$red," ","")
tmp <- strsplit(dat0$red,"\\|")

for (i in 1:nrow(dat_tmp)) {
dat_tmp$tit_red[i] <- ifelse(length(tmp[[i]])>0, str_count(dat_tmp$title[i],tmp[[i]]), 0)
}
table(dat_tmp$tit_red)

##
##    0    1    2    3
## 2329 1746  266    5

描述几处飘红

for (i in 1:nrow(dat_tmp)) {
dat_tmp$des_red[i] <- ifelse(length(tmp[[i]])>0, str_count(dat_tmp$desc[i],tmp[[i]]), 0)
}
table(dat_tmp$des_red)

##
##    0    1    2    3    4    5    6
## 3137  515  532   58   97    6    1

是否包含联系方式

dat_tmp$tel <- str_detect(dat_tmp$desc,"[\\d+(\\-\\d+)?]{7,}")*1 # 正则，七位以上数字，可以包含“-”
table(dat_tmp$tel)

##
##    0    1
## 2582 1764

修改数据格式，准备建模

dat_tmp$tel <- as.factor(dat_tmp$tel)
dat_tmp$tit_red <- as.integer(dat_tmp$tit_red)
dat_tmp$des_red <- as.integer(dat_tmp$des_red)
dat_tmp$kw_in <- as.factor(dat_tmp$kw_in)
dat_tmp$ranking <- as.factor(dat_tmp$ranking)
dat_tmp <- dat_tmp[,-c(1,3,4,5)]
str(dat_tmp)

## 'data.frame':    4346 obs. of  8 variables:
##  $ ranking: Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5 1 2 3 4 5 ...
##  $ fw_num : int  17 13 11 8 6 21 17 4 10 14 ...
##  $ kw_in  : Factor w/ 2 levels "0","1": 2 2 2 2 2 1 1 2 2 2 ...
##  $ des_l  : int  77 84 65 83 62 80 77 55 73 69 ...
##  $ tit_l  : int  22 23 25 19 14 25 22 6 24 15 ...
##  $ tit_red: int  0 0 0 0 1 0 0 0 1 1 ...
##  $ des_red: int  0 0 0 0 2 0 0 0 2 4 ...
##  $ tel    : Factor w/ 2 levels "0","1": 1 2 1 2 2 1 1 2 1 1 ...

定序回归

划分训练集和测试集

set.seed(42)
index <- sample(1:nrow(dat_tmp),round(0.75*nrow(dat_tmp)))
d_train <- dat_tmp[index,]
d_test <- dat_tmp[-index,]
dim(d_train);dim(d_test)

## [1] 3260    8
## [1] 1086    8

定序回归模型

library(MASS)
dat_m1 <- polr(ranking~1, data = dat_tmp,method = "logistic",Hess = T)
dat_m2 <- polr(ranking~., data = dat_tmp,method = "logistic",Hess = T)
anova(dat_m1,dat_m2)

## Likelihood ratio tests of ordinal regression models
##
## Response: ranking
##                                                      Model Resid. df
## 1                                                        1      4342
## 2 fw_num + kw_in + des_l + tit_l + tit_red + des_red + tel      4335
##   Resid. Dev   Test    Df LR stat. Pr(Chi)
## 1   13903.71
## 2   13142.37 1 vs 2     7 761.3341       0

summary(dat_m2)

## Call:
## polr(formula = ranking ~ ., data = dat_tmp, Hess = T, method = "logistic")
##
## Coefficients:
##             Value Std. Error t value
## fw_num  -0.013140   0.008942  -1.470
## kw_in1   0.245292   0.066428   3.693
## des_l    0.007748   0.002037   3.803
## tit_l   -0.106180   0.006715 -15.812
## tit_red  0.016441   0.055915   0.294
## des_red  0.412252   0.038961  10.581
## tel1    -0.843273   0.067432 -12.505
##
## Intercepts:
##     Value    Std. Error t value
## 1|2  -3.0731   0.1315   -23.3779
## 2|3  -1.9475   0.1270   -15.3341
## 3|4  -0.9771   0.1237    -7.8985
## 4|5   0.1805   0.1234     1.4631
##
## Residual Deviance: 13142.37
## AIC: 13164.37

library(car)
Anova(dat_m2,type = "III")

## Analysis of Deviance Table (Type III tests)
##
## Response: ranking
##         LR Chisq Df Pr(>Chisq)
## fw_num     2.159  1  0.1417077
## kw_in     13.632  1  0.0002224 ***
## des_l     14.480  1  0.0001416 ***
## tit_l    261.242  1  < 2.2e-16 ***
## tit_red    0.086  1  0.7687373
## des_red  116.250  1  < 2.2e-16 ***
## tel      158.262  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

模型筛选

此处应有AIC | BIC，挖坑

检验预测结果

pred <- predict(dat_m2,d_test)
table(d_test$ranking,pred)

##    pred
##       1   2   3   4   5
##   1 124  52  17  31   9
##   2 119  58  34  32  11
##   3  83  32  51  33  27
##   4  30  19  65  47  46
##   5  23  10  61  36  36

准确率29%，比随机猜测20%好一些

zq <- sum(d_test$ranking==pred)
zq/nrow(d_test)

## [1] 0.2909761

偏差正负1以内67.5%，随机猜测准确率52%

pc <- abs(as.integer(d_test$ranking)-as.integer(pred))
sum(pc<2)/nrow(d_test)

## [1] 0.674954

定序回归useful，但是还不够准，一方面是模型因素，另一方面可以对自变量进一步挖掘，从模型因素考虑，可以使用其他更鲁棒的模型，比如随机森林和xgboost等。

改进方法

随机森林

讲道理，随机森林应该不需要划分训练集和测试集，这是根据随机森林的bootstrap抽样思想决定的，关于bootstrap思想，可以戳这里回顾。不过我之所以这么说，完全是因为Leo Breiman and Adele Cutler他俩是酱婶说的：

In random forests, there is no need for cross-validation or a separate test set to get an unbiased estimate of the test set error. It is estimated internally, during the run…

不过这里有一点还没弄清楚，既然不需要考虑测试集，那么为什么还会有过拟合？虽然RF很鲁棒，但我看资料时，依然时常见到有人提出遇到RF过拟合的问题。挖个坑，以后来填。如果哪位童鞋知道，请指教。

随机森林建模与结果解读

以下暂用训练集建模，有空填。

library(randomForest)

d_m3 <- randomForest(as.factor(ranking)~.,data = d_train)
d_m3

输出结果如下：

##
## Call:
##  randomForest(formula = as.factor(ranking) ~ ., data = d_train)
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 2
##
##         OOB estimate of  error rate: 56.26%
## Confusion matrix:
##     1   2   3   4   5 class.error
## 1 552 116  72  35   8   0.2950192
## 2 297 217 124  47  24   0.6939351
## 3 194 160 158 101  46   0.7602428
## 4  61  46  88 263 143   0.5623960
## 5  47  31  47 147 236   0.5354331

summary告诉我们，这是一个分类模型（而不是回归），生成了500棵树用于投票表决，No. of variables tried at each split: 2指的是

randomForest

参数中的

mtry

，分裂过程中随机抽取特征的数量，没有固定默认值，关于mtry细节敲

?randomForest

查看，还记得随机森林为什么叫“随机”么？一是生成树的过程随机抽取样本，二是生成树的过程随机抽取特征。（这里是个人理解，请指教）。然后是袋外估计错误率 56.26%，关于oob error，稍后会提。接下来就是混淆矩阵了。

混淆矩阵可视化

关于混淆矩阵，一个可视化方法如下：

library(reshape2)
library(ggplot2)
cf_df <- fit_rf$confusion[,1:5]
cf_df <- melt(cf_df)
ggplot(cf_df,aes(x=Var1,y=Var2,fill=value))+geom_tile()

结果输出如下：



当然，我们的预测精度只能算useful，绝对谈不上高，一个好看的混淆矩阵应该是对角线清晰的。更多混淆矩阵可视化方案可以戳这里。本文的目的在于对机器学习的方法进行探讨，因此没有深入讨论特征工程，同时竞价因素对排名有着极大的影响，而这点并没有在数据当中体现出来（也不可能体现出来，道理你知道的，囧），所以这个模型混淆矩阵看起来不是很完美，不过能做到useful就很有用了，当然，你可以更深入的挖掘特征工程，进一步提高预测精度。

接下来要对测试集进行预测，因为rf在过拟合问题上很鲁棒，所以这一部分计划在未来的修改中干掉了。

pred_rf <- predict(d_m3,d_test)
table(d_test$ranking,pred_rf)

##    pred_rf
##       1   2   3   4   5
##   1 148  40  28  13   4
##   2 114  72  50  15   3
##   3  72  44  66  36   8
##   4  20  15  24  86  62
##   5  12  11  17  64  62

zq_rf <- sum(d_test$ranking==pred_rf)
zq_rf/nrow(d_test)

随机森林准确率40%，比前文的定序回归模型有了很大提高

## [1] 0.3996317

pc_rf <- abs(as.integer(d_test$ranking)-as.integer(pred_rf))
sum(pc_rf<2)/nrow(d_test)

偏差1以内准确率80%，显著提高

## [1] 0.7992634

袋外误差 OOB ERROR

一方面，我们希望了解rf模型的准确率，另一方面，也希望了解以下生成树的数量，因为如果生成两百棵树和五百颗树的结果相差不多，干嘛浪费时间和算力去生成五百颗树呢？

plot(fit_rf)

这花花绿绿的画是个什么gui！！Σ(っ °Д °;)っ

嗯，这里的纵坐标是指oob error，就是前文提到的袋外误差，具体内容可以戳这里，这是翻译自stackoverflow的一个回答，原文参见链接博客的引用。

实际上这张图最关键的就是黑色那条线，这代表了模型整体的袋外误差，也可以理解为1-准确率。

从图中可以看到，当生成树的数量大于100时，我们的RF模型误差就已经不再明显下降了，因此我们可以在模型中设置

ntree=200

来改进我们的算法，精简算力，你可以自己动手试试。

好了，黑线我们清楚了，那么其他花花绿绿的线又是个什么gui？

实际上，这张图画的是

fit_rf$err.rate

中的数据，你可以尝试画一下

plot(fit_rf$err.rate[,1],type="l")

来深入体会一下。

fit_rf$err.rate[,1]

是模型整体误差，剩下五条线分别对应了五个分类，即每一个分类的分错概率。如果你像我一样忍不住想要弄清这些东西到底代表啥的话（虽然这有点没啥用），可以执行以下命令：

plot(fit_rf)
legend("top", cex =0.5, legend = colnames(fit_rf$err.rate), lty=c(1,2,3,4,5,6), col=c(1,2,3,4,5,6), horiz=T)

输出结果如下：

特征重要性与特征选择

做线性回归时，模型summary可以直接看到p-value，后面的几个小星号很是方便，哪个特征显著，哪个特征不显著，一眼就可以看出来，类似的，我们也希望了解到哪些特征对随机森林效果显著，这一点对高维数据集尤其重要，对于动不动就成百上千个维度的数据，显然并不是每个特征都对RF模型有显著贡献，因此我们需要引入“特征重要性”这个概念，可以通过

varImpPlot

函数刻画：

varImpPlot(fit_rf)

很直观的，我们可以看到前三个维度非常重要，后面四个维度相对不那么重要。

可问题是，不那么重要到底有多重要？有没有可以量化的标准指导我们如何筛选特征？幸运的是，借助caret包，我们可以轻松实现这个目的：

library(caret)
control <- rfeControl(functions=rfFuncs, method="cv", number=10)
set.seed(24)
results <- rfe(d_train[,-1], d_train[,1], sizes=c(1:8), rfeControl=control)
results

输出结果如下：

Recursive feature selection

Outer resampling method: Cross-Validated (10 fold)

Resampling performance over subset size:

Variables Accuracy  Kappa AccuracySD KappaSD Selected
1   0.3702 0.2047    0.02382 0.02909
2   0.4147 0.2620    0.02541 0.03270
3   0.4325 0.2856    0.02805 0.03544        *
4   0.4298 0.2813    0.03301 0.04148
5   0.4289 0.2797    0.02500 0.03190
6   0.4270 0.2768    0.01929 0.02446
7   0.4304 0.2808    0.01916 0.02459

The top 3 variables (out of 3):
des_l, tit_l, fw_num

results结果选择了前3个特征，这与我们之前观察到的结果是一致的，但是并不绝对，因为我们设置了随机种子，你可以修改随机种子查看是否会出现其他情况。

接下来可以使用

predictors

查看建议的特征选择：

predictors(results)

输出caret建议特征如下：

[1] "des_l"  "tit_l"  "fw_num"

将results结果可视化：

plot(results, type=c("g", "o"))

可以清晰的看到为什么选择这三个特征。

但是这里我并没有针对这三个特征重新建模，一是因为这个结果是通过设置随机种子实现的，如果你去掉随机种子的话，结果可能会发生变化，说明选3个特征还是7个特征差别并不是特别大；二是因为我们的数据维度很小，没必要对特征进行筛选，但是如果真的遇到成百上千的特征，特征选择会很有帮助。

使用格搜索优化超参数

ctrl <- trainControl(method = "repeatedcv", number = 5, repeats = 5)
grid_rf <- expand.grid(.mtry = c(2, 4, 8, 16))
set.seed(300)
m_rf <- train(as.factor(ranking)~ ., data = d_train, method = "rf",
metric = "Kappa", trControl = ctrl,
tuneGrid = grid_rf)
m_rf

Random Forest

3260 samples
7 predictor
5 classes: '1', '2', '3', '4', '5'

No pre-processing
Resampling: Cross-Validated (5 fold, repeated 5 times)
Summary of sample sizes: 2608, 2609, 2607, 2607, 2609, 2606, ...
Resampling results across tuning parameters:

mtry  Accuracy   Kappa
2    0.4282253  0.2781285
4    0.4341736  0.2877875
8    0.4287764  0.2808918
16    0.4266310  0.2781918

Kappa was used to select the optimal model using the largest value.
The final value used for the model was mtry = 4.

通常应该选取kappa值最大的参数为模型参数，但本例kappa值差距不大，优化意义不大。

超参数优化后的建模如下：

fit_rf <- randomForest(as.factor(ranking)~.,data=d_train,mtry=4)
fit_rf

xgboost

library(xgboost)

dat_xg <- dat_tmp
str(dat_xg)

## 'data.frame':    4346 obs. of  8 variables:
##  $ ranking: Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5 1 2 3 4 5 ...
##  $ fw_num : int  17 13 11 8 6 21 17 4 10 14 ...
##  $ kw_in  : Factor w/ 2 levels "0","1": 2 2 2 2 2 1 1 2 2 2 ...
##  $ des_l  : int  77 84 65 83 62 80 77 55 73 69 ...
##  $ tit_l  : int  22 23 25 19 14 25 22 6 24 15 ...
##  $ tit_red: int  0 0 0 0 1 0 0 0 1 1 ...
##  $ des_red: int  0 0 0 0 2 0 0 0 2 4 ...
##  $ tel    : Factor w/ 2 levels "0","1": 1 2 1 2 2 1 1 2 1 1 ...

重新修改数据格式，xgboost要求数据类型不能为int，修改为num

for (i in c(2,4:7)) {
dat_xg[,i] <- as.numeric(dat_xg[,i])
}
str(dat_xg)

## 'data.frame':    4346 obs. of  8 variables:
##  $ ranking: Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5 1 2 3 4 5 ...
##  $ fw_num : num  17 13 11 8 6 21 17 4 10 14 ...
##  $ kw_in  : Factor w/ 2 levels "0","1": 2 2 2 2 2 1 1 2 2 2 ...
##  $ des_l  : num  77 84 65 83 62 80 77 55 73 69 ...
##  $ tit_l  : num  22 23 25 19 14 25 22 6 24 15 ...
##  $ tit_red: num  0 0 0 0 1 0 0 0 1 1 ...
##  $ des_red: num  0 0 0 0 2 0 0 0 2 4 ...
##  $ tel    : Factor w/ 2 levels "0","1": 1 2 1 2 2 1 1 2 1 1 ...

重新划分训练集和测试集

set.seed(42)
index <- sample(1:nrow(dat_xg),round(0.75*nrow(dat_xg)))
d_train <- dat_xg[index,]
d_test <- dat_xg[-index,]

xgb <- xgboost(data = data.matrix(d_train[,-1]), label = d_train$ranking, max.depth = 6, eta = 0.1,nround = 2, nthread=4, num_class=6, objective = "multi:softmax")

## [1]  train-merror:0.515644
## [2]  train-merror:0.515951

查看预测结果

y_pred <- predict(xgb, data.matrix(d_test[,-1]))
table(d_test$ranking,y_pred)

##    y_pred
##       1   2   3   4   5
##   1 142  51  27  10   3
##   2 109  70  49  19   7
##   3  69  65  47  34  11
##   4  31  20  22  94  40
##   5  14  16  22  64  50

zq_xg <- sum(d_test$ranking==y_pred)
zq_xg/nrow(d_test)

正确率37%，比随机森林略低，实际上xgboost参数众多，修改参数也许有效，有空回来填坑。

## [1] 0.3710866

偏差1以内准确率77%

pc_xg <- abs(as.integer(d_test$ranking)-as.integer(y_pred))
sum(pc_xg<2)/nrow(d_test)

## [1] 0.7707182