R语言

继续挖坑
这次是基于狗熊会基础案例百度搜索竞价广告创意研究的定序回归与随机森林、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β 耐着性子重新码了一遍，增补随机森林袋外误差细节，及随机森林其他内容

---

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

简易特征工程

读取数据

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__ "电话[phone]卡耐基口才培训—独一无二特浸泡式口才培训体验教学!瞬间克服紧张各种场合站起能说;表达生动 说服力强!"| __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]) # 取词长大于1head(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_wplot(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 <- dat0dat_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_intable(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      [phone]  < 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

plot(fit_rf)

这花花绿绿的画是个什么gui！！Σ(っ °Д °;)っ
嗯，这里的纵坐标是指oob error，就是前文提到的袋外误差，具体内容可以戳这里，这是翻译自stackoverflow的一个回答，原文参见链接博客的引用。
实际上这张图最关键的就是黑色那条线，这代表了模型整体的袋外误差，也可以理解为1-准确率。
好了，黑线我们清楚了，那么其他花花绿绿的线又是个什么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)

输出结果如下：

varImpPlot(fit_rf)

使用格搜索优化超参数

library(caret)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-processingResampling: 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.2781918Kappa 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_tmpstr(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