前言

在通常的机器学习相关课程或资料中,我们学习的大多是对连续变量或者分类变量进行建模。但是现实世界中,等级变量(带有顺序的离散变量)也非常常见。2015年11月末开始的Kaggle数据分析比赛,Prudential Life Insurance Assessment,就是对等级变量进行建模的比赛,并使用不太常见的kappa值(其实就是类似于MSE的离散版本平方错误)进行评判指标。我非常幸运地参与了这场比赛,并在最后得到了第99名(总2619位有效参赛者),入围前10%,并学到了很多对等级变量建模的经验,希望有机会可以应用于之后的比赛和类似的工作。

数据

首先读取本次比赛的训练数据和测试数据,查看数据的基本结构和类型:

setwd("~/Desktop/SERN 151213")
train <- read.csv("train.csv")
test  <- read.csv("test.csv")
test$Response <- -1
dim(train)
## [1] 59381   128
dim(test)
## [1] 19765   128
str(train)
## 'data.frame':    59381 obs. of  128 variables:
##  $ Id                 : int  2 5 6 7 8 10 11 14 15 16 ...
##  $ Product_Info_1     : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Product_Info_2     : Factor w/ 19 levels "A1","A2","A3",..: 17 1 19 18 16 16 8 16 17 19 ...
##  $ Product_Info_3     : int  10 26 26 10 26 26 10 26 26 21 ...
##  $ Product_Info_4     : num  0.0769 0.0769 0.0769 0.4872 0.2308 ...
##  $ Product_Info_5     : int  2 2 2 2 2 3 2 2 2 2 ...
##  $ Product_Info_6     : int  1 3 3 3 3 1 3 3 3 3 ...
##  $ Product_Info_7     : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Ins_Age            : num  0.6418 0.0597 0.0299 0.1642 0.4179 ...
##  $ Ht                 : num  0.582 0.6 0.745 0.673 0.655 ...
##  $ Wt                 : num  0.149 0.132 0.289 0.205 0.234 ...
##  $ BMI                : num  0.323 0.272 0.429 0.352 0.424 ...
##  $ Employment_Info_1  : num  0.028 0 0.03 0.042 0.027 0.325 0.11 0.12 0.165 0.025 ...
##  $ Employment_Info_2  : int  12 1 9 9 9 15 1 12 9 1 ...
##  $ Employment_Info_3  : int  1 3 1 1 1 1 3 1 1 3 ...
##  $ Employment_Info_4  : num  0 0 0 0 0 0 NA 0 0 0 ...
##  $ Employment_Info_5  : int  3 2 2 3 2 2 3 2 2 3 ...
##  $ Employment_Info_6  : num  NA 0.0018 0.03 0.2 0.05 1 0.8 1 1 0.05 ...
##  $ InsuredInfo_1      : int  1 1 1 2 1 1 1 1 1 2 ...
##  $ InsuredInfo_2      : int  2 2 2 2 2 2 2 2 2 2 ...
##  $ InsuredInfo_3      : int  6 6 8 8 6 8 3 6 3 3 ...
##  $ InsuredInfo_4      : int  3 3 3 3 3 3 3 3 2 3 ...
##  $ InsuredInfo_5      : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ InsuredInfo_6      : int  2 2 1 2 2 1 2 1 1 2 ...
##  $ InsuredInfo_7      : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Insurance_History_1: int  1 2 2 2 2 2 1 1 1 2 ...
##  $ Insurance_History_2: int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Insurance_History_3: int  3 3 1 1 1 3 3 3 3 3 ...
##  $ Insurance_History_4: int  1 1 3 3 3 2 2 1 2 1 ...
##  $ Insurance_History_5: num  0.000667 0.000133 NA NA NA ...
##  $ Insurance_History_7: int  1 1 3 3 3 1 1 1 1 1 ...
##  $ Insurance_History_8: int  1 3 2 2 2 3 1 1 1 3 ...
##  $ Insurance_History_9: int  2 2 3 3 3 2 2 2 2 2 ...
##  $ Family_Hist_1      : int  2 2 3 3 2 2 3 2 3 3 ...
##  $ Family_Hist_2      : num  NA 0.188 0.304 0.42 0.464 ...
##  $ Family_Hist_3      : num  0.598 NA NA NA NA ...
##  $ Family_Hist_4      : num  NA 0.0845 0.2254 0.3521 0.4085 ...
##  $ Family_Hist_5      : num  0.527 NA NA NA NA ...
##  $ Medical_History_1  : int  4 5 10 0 NA 6 5 6 4 NA ...
##  $ Medical_History_2  : int  112 412 3 350 162 491 600 145 16 162 ...
##  $ Medical_History_3  : int  2 2 2 2 2 2 3 2 2 2 ...
##  $ Medical_History_4  : int  1 1 2 2 2 2 2 2 2 2 ...
##  $ Medical_History_5  : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Medical_History_6  : int  3 3 3 3 3 3 3 3 3 3 ...
##  $ Medical_History_7  : int  2 2 2 2 2 2 2 2 2 2 ...
##  $ Medical_History_8  : int  2 2 2 2 2 2 2 2 2 2 ...
##  $ Medical_History_9  : int  1 1 2 2 2 2 1 1 1 2 ...
##  $ Medical_History_10 : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ Medical_History_11 : int  3 3 3 3 3 3 3 3 3 3 ...
##  $ Medical_History_12 : int  2 2 2 2 2 2 2 2 2 2 ...
##  $ Medical_History_13 : int  3 3 3 3 3 3 3 3 3 3 ...
##  $ Medical_History_14 : int  3 3 3 3 3 3 3 3 3 3 ...
##  $ Medical_History_15 : int  240 0 NA NA NA NA NA NA NA NA ...
##  $ Medical_History_16 : int  3 1 1 1 1 1 1 1 1 3 ...
##  $ Medical_History_17 : int  3 3 3 3 3 3 3 3 3 3 ...
##  $ Medical_History_18 : int  1 1 1 1 1 2 1 1 1 1 ...
##  $ Medical_History_19 : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Medical_History_20 : int  2 2 2 2 2 2 2 2 2 2 ...
##  $ Medical_History_21 : int  1 1 1 2 1 2 1 1 1 1 ...
##  $ Medical_History_22 : int  2 2 2 2 2 2 2 2 2 2 ...
##  $ Medical_History_23 : int  3 3 3 3 3 3 3 3 3 1 ...
##  $ Medical_History_24 : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ Medical_History_25 : int  1 1 2 1 2 1 1 1 1 1 ...
##  $ Medical_History_26 : int  3 3 2 3 2 3 3 3 3 3 ...
##  $ Medical_History_27 : int  3 3 3 3 3 3 3 3 3 3 ...
##  $ Medical_History_28 : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Medical_History_29 : int  3 3 3 3 3 3 1 3 1 3 ...
##  $ Medical_History_30 : int  2 2 2 2 2 2 2 2 2 2 ...
##  $ Medical_History_31 : int  3 3 3 3 3 3 3 3 3 3 ...
##  $ Medical_History_32 : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ Medical_History_33 : int  1 3 3 3 3 3 3 3 3 3 ...
##  $ Medical_History_34 : int  3 1 3 3 3 1 3 3 3 3 ...
##  $ Medical_History_35 : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Medical_History_36 : int  2 2 3 2 3 2 2 2 2 2 ...
##  $ Medical_History_37 : int  2 2 2 2 2 2 2 2 2 2 ...
##  $ Medical_History_38 : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Medical_History_39 : int  3 3 3 3 3 3 3 3 3 3 ...
##  $ Medical_History_40 : int  3 3 3 3 3 3 3 3 3 3 ...
##  $ Medical_History_41 : int  3 1 1 1 1 3 3 1 3 1 ...
##  $ Medical_Keyword_1  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Medical_Keyword_2  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Medical_Keyword_3  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Medical_Keyword_4  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Medical_Keyword_5  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Medical_Keyword_6  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Medical_Keyword_7  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Medical_Keyword_8  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Medical_Keyword_9  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Medical_Keyword_10 : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Medical_Keyword_11 : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Medical_Keyword_12 : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Medical_Keyword_13 : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Medical_Keyword_14 : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Medical_Keyword_15 : int  0 0 0 0 0 0 0 0 0 1 ...
##  $ Medical_Keyword_16 : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Medical_Keyword_17 : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Medical_Keyword_18 : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Medical_Keyword_19 : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Medical_Keyword_20 : int  0 0 0 0 0 0 0 0 1 0 ...
##   [list output truncated]

如上所示,原始数据一共127个特征,大部分的特征名只提供了对应特征非常粗略的意义,并且所有的连续型变量是经过“标准化”调整以后的数值(具体调整方法不明,大部分被调整至0-1的范围内),导致本次比赛的人工特征设计变得异常困难。根据主办方提供的数据类型的信息,我将读入的数据按提供的数据类型信息进行一一转换,得益于Rfactor的特性,我不需要对分类数据做特别的one-hot-encoding操作。不过在数据转换之前,利用读入的原始数据类型,我构建了三个新特征,分别是总NA出现的次数,前缀为“Medical_Keyword”的dummy variable出现1的次数以及BMI和年龄交互项。在完成上述处理后,我使用均数填充连续变量中的NA值,使用众数(自写mode函数)填充离散变量中的NA值。在这里插几句题外话,其实h2o那些基于树的模型也好Xgboost也好,是支持自动处理NA值的(具体方式待考,可能是赋一个极端值,而且在Xgboost中不同版本的处理方式据说略有差异)。之所以我在这里对NA进行人工赋值,是因为在接下来的建模过程中,训练数据略大以及很多特征为binary类型,我会将原始的数据转换成稀疏矩阵以减少内存消耗提高计算速度,而那些算法包无法对稀疏矩阵的NA自动处理(目前为止的版本)。

train$NA_count <- rowSums(is.na(train))
test$NA_count <- rowSums(is.na(test))
train$Medical_Keyword_count <- rowSums(train[, paste("Medical_Keyword_", as.character(1:48), sep = "")])
test$Medical_Keyword_count <- rowSums(test[, paste("Medical_Keyword_", as.character(1:48), sep = "")])
train$BMI_Age <- train$BMI * train$Ins_Age
test$BMI_Age <- test$BMI * test$Ins_Age
nominal_features <- c("Product_Info_1", "Product_Info_2", "Product_Info_3", 
                      "Product_Info_5", "Product_Info_6", "Product_Info_7", 
                      "Employment_Info_2", "Employment_Info_3", "Employment_Info_5", 
                      "InsuredInfo_1", "InsuredInfo_2", "InsuredInfo_3", 
                      "InsuredInfo_4", "InsuredInfo_5", "InsuredInfo_6", "InsuredInfo_7",
                      "Insurance_History_1", "Insurance_History_2", "Insurance_History_3", 
                      "Insurance_History_4", "Insurance_History_7", "Insurance_History_8", 
                      "Insurance_History_9", "Family_Hist_1", "Medical_History_2", 
                      "Medical_History_3", "Medical_History_4", "Medical_History_5", 
                      "Medical_History_6", "Medical_History_7", "Medical_History_8", 
                      "Medical_History_9", "Medical_History_11", "Medical_History_12", 
                      "Medical_History_13", "Medical_History_14", "Medical_History_16", 
                      "Medical_History_17", "Medical_History_18", "Medical_History_19", 
                      "Medical_History_20", "Medical_History_21", "Medical_History_22", 
                      "Medical_History_23", "Medical_History_25", "Medical_History_26", 
                      "Medical_History_27", "Medical_History_28", "Medical_History_29", 
                      "Medical_History_30", "Medical_History_31", "Medical_History_33", 
                      "Medical_History_34", "Medical_History_35", "Medical_History_36", 
                      "Medical_History_37", "Medical_History_38", "Medical_History_39", 
                      "Medical_History_40", "Medical_History_41")
continuous_features <- c("Product_Info_4", "Ins_Age", "Ht", "Wt", "BMI", "Employment_Info_1", 
                         "Employment_Info_4", "Employment_Info_6", "Insurance_History_5", 
                         "Family_Hist_2", "Family_Hist_3", "Family_Hist_4", "Family_Hist_5")
discrete_features <- c("Medical_History_1", "Medical_History_10", "Medical_History_15", 
                       "Medical_History_24", "Medical_History_32")
InsuredInfo_nominal <- c("InsuredInfo_1", "InsuredInfo_2", "InsuredInfo_4", 
                         "InsuredInfo_5", "InsuredInfo_6", "InsuredInfo_7")
Insurance_History_nominal <- c("Insurance_History_1", "Insurance_History_2", "Insurance_History_3", 
                               "Insurance_History_4", "Insurance_History_7", "Insurance_History_8", 
                               "Insurance_History_9")
Medical_History_nominal <- c("Medical_History_3", "Medical_History_4", "Medical_History_5", 
                             "Medical_History_6", "Medical_History_7", "Medical_History_8", 
                             "Medical_History_9", "Medical_History_11", "Medical_History_12", 
                             "Medical_History_13", "Medical_History_14", "Medical_History_16", 
                             "Medical_History_17", "Medical_History_18", "Medical_History_19", 
                             "Medical_History_20", "Medical_History_21", "Medical_History_22", 
                             "Medical_History_23", "Medical_History_25", "Medical_History_26", 
                             "Medical_History_27", "Medical_History_28", "Medical_History_29", 
                             "Medical_History_30", "Medical_History_31", "Medical_History_33", 
                             "Medical_History_34", "Medical_History_35", "Medical_History_36", 
                             "Medical_History_37", "Medical_History_38", "Medical_History_39", 
                             "Medical_History_40", "Medical_History_41")
for (feature in nominal_features) {
  train[, feature] <- as.factor(train[, feature])
  test[, feature] <- as.factor(test[, feature])
}
for (feature in continuous_features) {
  train[is.na(train[, feature]), feature] <- median(train[, feature], na.rm = TRUE)
  train[, feature] <- as.numeric(train[, feature])
  test[is.na(test[, feature]), feature] <- median(test[, feature], na.rm = TRUE)
  test[, feature] <- as.numeric(test[, feature])
}
Mode <- function(x) {
  x <- x[!is.na(x)]
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}
for (feature in discrete_features) {
  train[is.na(train[, feature]), feature] <- Mode(train[, feature])
  train[, feature] <- as.numeric(train[, feature])
  test[is.na(test[, feature]), feature] <- Mode(test[, feature])
  test[, feature] <- as.numeric(test[, feature])
}

接着我们看看这次的预测目标Response是长什么样子的。

table(train$Response)
## 
##     1     2     3     4     5     6     7     8 
##  6207  6552  1013  1428  5432 11233  8027 19489
hist(train$Response)

根据主办方的说明,Response代表的是某一种risk的程度。从直方图上看,呈现两头略高中间低陷的分布。我曾经花了很多时间探索为什么会产生这样的一个分布形态,因为训练数据接近5万,在这样一个数量级之下,根据中心极限定理,会多少接近一个钟形的分布。有人根据Response的均数与方差接近(均数5.64,方差6.04)来判断它是一个Poisson分布,我并不认同,无论从分布的形态还是Response背后的直观意义。我的理解是Response对应于某种风险率倍数的关系,当然可能性多种多样,我也仅仅停留于猜测。重要的是,后面使用到的方法并不会受到其分布形态的影响。

关键点

在进一步的工作之前,我想说Kaggle的比赛或者机器学习的工作,其实大部分也是step-by-step程序化的东西,就像我作为麻醉科大夫上个全身麻醉,也不过是镇静镇痛加肌松,然后根据不同的患者情况和手术需要进行微调。在这个比赛中,我觉得关键点有三:可视化、验证(cv)和ensemble或者说stacking。

首先是可视化。可视化的意义不仅仅是进行变量的筛选,或者说在特征数目非常大的情况下,用可视化方法筛选变量,效率上可能远远不及使用相关系数或者使用RandomForest中的变量importance来筛选。可视化的意义更多的在于让你知道数据的分布,特征间的关系,以让你判断是不是需要对数据进行一些数学形式上变换,是不是需要对离群点(outlier)进行一些调整和处理,是不是需要对数据进行粗粒化(离散化或者减少小数点后有效位数)以减少noise的干扰。在本次比赛中,由于之后的特征设计涉及无监督聚类,可视化也是让你判断聚类效果最简便直观的方法。所以在进一步的数据处理之前,我写了一些可视化函数,可以输入特征名(分类特征)来得到该特征的每个类别对应Response的箱盒图(由于数据偏态分布,我使用箱盒图)。同时在title的位置会显示每个类别所占比重,使用parallel包多线程计算加速计算效率。当然该函数也可以略作修改呈现其他的图例,在此不一一例举。

library(ggplot2)
library(gridExtra)
library(parallel)
plot_box <- function(data, col_name) {
  data <- data.frame(x = data[, col_name], Response = data$Response)
  prop <- round(table(data$x)/nrow(data), digits = 2)
  prop_str <- paste("class ", names(prop), " : ", as.vector(prop), sep = "", collapse = "; ")
  plt <- ggplot(data, aes(x = factor(x), y = Response)) + geom_boxplot() + xlab(col_name) + ggtitle(prop_str) + theme_light()
  return (plt)
}
do_plots <- function(data, fun, col_names_vector, ncol = 3) {
  plt_list <- mclapply(col_names_vector, mc.preschedule = TRUE, mc.cores = 8, function(col_name) {
    plt <- fun(data = data, col_name = col_name)
  })
  do.call("grid.arrange", c(plt_list, ncol = ncol))
}
do_plots(data = train, fun = plot_box, col_names_vector = paste("Product_Info_", c(1:3, 5:7), sep = ""), ncol = 3)

其次是验证,不仅仅是交叉验证(cv),更要使用独立的数据集进行验证。我觉得只有当自己的cv结果,独立的验证集的结果和Kaggle的public leaderboard的测试数据的结果在建模效果的提升上步调接近一致,才比较靠谱。当然这个仅仅是个人观点,但仅仅依靠public leaderboard的结果判断建模效果,是绝对不可取的。 本比赛中我使用训练数据30%的数据作为独立的验证集,仅参与建模完成以后的测试,不以任何形式参与建模。分割方法如下:

set.seed(1126)
val_size <- round(nrow(train) * 0.3)
val_idx <- sample(1:nrow(train), val_size)
train_idx <- setdiff(1:nrow(train), val_idx)
full_train_idx <- 1:nrow(train)

最后是ensemble的部分。其实本次比赛最后ensemble的部分完全依赖的是独特的思路,本身的ensemble的方法相当的简单。Ensemble数据集的基本构成参考了之前Kaggle冠军的solution,使用stacking的方法,类似于cv,将数据等分为n份,底层的模型用“n-1份”的数据训练。对训练集,使用“1份”数据进行预测,一共n个迭代,构成总观察数与原始训练集观测数相等的ensemble训练集。对于测试集,使用“n-1份”训练的模型进行预测,一共n个迭代,然后取预测值的均值构成总观察数与原始测试集观测数相等的ensemble测试集。这个方法相较传统的将数据分割成专为底层模型训练的数据集和专为ensemble训练的数据集的方法,提供了更多的数据而不失稳定性。在训练数据和测试数据来自于同一整体的情况下,ensemble训练集和ensemble测试集在个meta-feature(底层模型的预测值)的分布是相当接近的,接近程度随着分割数n(n-fold stacking)的增大而增大。在本次比赛中,考虑到数据并没有非常大,我使用了10-fold stacking,一般情况下5-fold已经可以满足需要。

特征设计

尽管本比赛特征设计很难,但本着特征为王,调参为辅,ensemble为辅的理念,这一步过程还是非常必要的。在这次比赛中特征设计的思路基本分为两路:利用无监督方法(选择可以有KMeans,LCA,t-SNE或autoenconder)构建meta-feature和利用有限的特征信息(如之前已经完成的年龄和BMI的交互)进行特征的组合。

首先我尝试使用KMeans进行聚类,不过考虑到大量的分类变量以及前缀为“Medical_Keyword”的变量类别为0的在数量上占有绝对的优势,使用所有的特征会使各观测之间的相似度过高导致聚类效果不佳,我先使用RandomForest对特征进行筛选,提取占总体重要性95%重要性的变量进行聚类。

library(h2o)
## Loading required package: statmod
## 
## ----------------------------------------------------------------------
## 
## Your next step is to start H2O:
##     > h2o.init()
## 
## For H2O package documentation, ask for help:
##     > ??h2o
## 
## After starting H2O, you can use the Web UI at http://localhost:54321
## For more information visit http://docs.h2o.ai
## 
## ----------------------------------------------------------------------
## 
## 
## Attaching package: 'h2o'
## 
## The following objects are masked from 'package:stats':
## 
##     sd, var
## 
## The following objects are masked from 'package:base':
## 
##     %*%, %in%, apply, as.factor, as.numeric, colnames, colnames<-,
##     ifelse, is.factor, is.numeric, log, trunc
h2o.init(nthreads = -1)
trainHex <- as.h2o(train)
features <- names(trainHex)[-c(1, 128)]
rfHex <- h2o.randomForest(x = features,
                          y = "Response", 
                          model_id = "rfHex_features",
                          training_frame = trainHex,
                          ntrees = 300,
                          seed = 1126
)
rfHex_importance <- as.data.frame(rfHex@model$variable_importances)
which(cumsum(rfHex_importance$percentage) > 0.95)[1]
important_features <- rfHex_importance$variable[1:which(cumsum(rfHex_importance$percentage) > 0.95)[1]]
important_features
##  [1] "Medical_History_2"     "BMI"                  
##  [3] "Medical_History_15"    "Medical_History_23"   
##  [5] "Wt"                    "BMI_Age"              
##  [7] "Product_Info_2"        "Medical_Keyword_3"    
##  [9] "Medical_Keyword_count" "Product_Info_4"       
## [11] "Medical_History_4"     "InsuredInfo_3"        
## [13] "Ins_Age"               "Employment_Info_2"    
## [15] "Medical_Keyword_15"    "Ht"                   
## [17] "Medical_History_39"    "Family_Hist_3"        
## [19] "Employment_Info_6"     "Family_Hist_5"        
## [21] "Medical_History_40"    "Medical_History_1"    
## [23] "Employment_Info_1"     "Family_Hist_4"        
## [25] "NA_count"              "InsuredInfo_6"        
## [27] "Family_Hist_2"         "Insurance_History_5"  
## [29] "Medical_History_30"    "Medical_History_5"    
## [31] "InsuredInfo_5"         "Medical_History_13"   
## [33] "Product_Info_3"        "Insurance_History_2"  
## [35] "Medical_History_20"    "Medical_History_28"   
## [37] "Medical_History_18"    "InsuredInfo_7"        
## [39] "Insurance_History_8"   "Employment_Info_4"    
## [41] "InsuredInfo_2"         "Medical_Keyword_38"   
## [43] "Insurance_History_4"   "Medical_History_11"   
## [45] "Family_Hist_1"         "Insurance_History_7"  
## [47] "Medical_History_41"    "InsuredInfo_1"        
## [49] "Medical_History_17"    "Medical_Keyword_23"   
## [51] "Medical_History_27"    "Medical_History_33"   
## [53] "Medical_History_36"    "Medical_History_9"    
## [55] "Medical_History_7"     "Medical_History_8"    
## [57] "Medical_History_16"    "Medical_Keyword_37"

使用无监督算法的技巧之一是需要将训练数据和测试数据一同训练。所以我先将训练数据和测试数据进行合并。KMeans算法简单直观,h2o的KMeans在速度上也非常理想,但是KMeans非常依赖于初始值和分类数,所以通常需要一个个分类数一个个random seed进行尝试和挑选。在这里我就完全依据之前写的可视化函数判断聚类的效果是不是理想。其基于的假设是,如果是理想的分类,那么每个类别对应的Response分布之间应该或多或少呈现差异。满意的分类数和seed需要记录好以便之后的数据处理可以重现结果。这里运行的结果仅仅是我在若干次尝试中感觉还比较满意的分类。

data <- rbind(train, test)
dataHex <- as.h2o(data)
KM <- h2o.kmeans(training_frame = dataHex,
                 x = important_features,
                 k = 8,
                 model_id = "KM",
                 max_iterations = 1000,
                 seed = 1126)
KM
## Model Details:
## ==============
## 
## H2OClusteringModel: kmeans
## Model ID:  KM 
## Model Summary: 
##   number_of_rows number_of_clusters number_of_categorical_columns
## 1          79146                  8                            35
##   number_of_iterations within_cluster_sum_of_squares total_sum_of_squares
## 1                   49                 1764859.90665       84106098.00000
##   between_cluster_sum_of_squares
## 1                 82341238.09335
## 
## 
## H2OClusteringMetrics: kmeans
## ** Reported on training data. **
## 
## 
## Total Within SS:  1764860
## Between SS:  82341238
## Total SS:  84106098 
## Centroid Statistics: 
##   centroid        size within_cluster_sum_of_squares
## 1        1 14476.00000                  271241.61785
## 2        2   544.00000                   15689.38722
## 3        3  7373.00000                  289241.93919
## 4        4  4767.00000                  142625.87094
## 5        5 11536.00000                  305740.46524
## 6        6  6335.00000                  150969.37010
## 7        7 20227.00000                  314917.33476
## 8        8 13888.00000                  274433.92134
data$KM <- as.factor(as.data.frame(predict(KM, dataHex))[, 1])

看一下聚类的效果。

plot_box(data = data[1:nrow(train), ], col_name = "KM")

我还尝试了使用latent classify analysis(LCA)进行聚类。LCA并不是一个常用的无监督方法,因为它只能处理分类变量,但这次的比赛分类变量占有主导地位,所以LCA看似是一个很好的选择。但因为我使用LCA工具包poLCA只能够支持单线程的计算,所以LCA的计算相当的耗时,导致参数的调整不是非常的方便,试了几个结果加入模型后性能也并未提高,可能的原因是尽管LCA的机制和KMeans不同(LCA基于概率分布而不是空间距离),在结果上还是有很高的相关性,最后所幸弃之不用了。在这里仅贴出示例代码, formulacbind的部分包括了需要聚类的feature,也可以在~ 1的部分加入其它项或者交互项,你可以把他当成虚拟一个变量对需要聚类的feature进行建模,然后得到这个变量:

library(poLCA)
formula <- cbind(Product_Info_1, Product_Info_2, Product_Info_3, 
                 Product_Info_5, Product_Info_6, Product_Info_7) ~ 1
set.seed(1126)
LCA_Product_Info <- poLCA(formula, data_Product_Info, nclass = 8, nrep = 1, maxiter = 50000)
data$Product_Info_LCA <- as.factor(LCA_Product_Info$predclass)

我也使用h2o的autoenconder进行了尝试,不过一方面对autoenconder了解不深,一方面数据比较奇葩,效果非常糟糕。t-SNE由于时间关系没有尝试。

在人工特征方面其实进行了相当多的尝试,包括将一些连续变量分割成二分类,按Product_Info_2的字母前缀重新构成分类,或者将BMI调整至围绕0点分布然后取绝对值这些帮助模型识别pattern的操作(BMI过低或过高都不好吧),但是结果往往都是增加了noise,令人沮丧。可能这些特征对单纯线性模型会有帮助,而Xgboost这类复杂的模型已经足以识别这些特征。所以最后仅仅增加了先前的三个特征而已,simple is good。

建模

在建模之前,我们再次回顾一下Response的特点以及比赛的评判metric。Response为1~8的整数,metric为kappa值,也是一种平方错误,在kappa的判定下,将1判断为8的错误是将1判断为2的错误的49倍。因此,使用multi:softmax类似的方式进行多分类建模是不合适的,因为各类别之间并不是平等的关系,保留Response的等级关系是建模的关键。换句话说,一个理想模型,g(x),最低要求是,对于实际Response比较小的x,输出的值需要比Response比较大的x来的小,反之亦然。此外,kappa值可以通过Metrics包中的ScoreQuadraticWeightedKappa函数进行计算,不需要自写函数。

在模型的选择方面,我分别尝试了h2o的RandomForest,glm和Xgboost。因为数据量的关系只能在高运算效率的工具包中进行选择了,抱歉我无法忍受单线程包的耗时,不禁感叹caret的那些东西已经落后于时代了。Xgboost的性能远超前两个,甚至在使用三者进行ensemble的尝试中,整体性能也会因为RandomForest和glm性能的低下而不及单一Xgboost。其实这次的比赛reg:linear的Xgboost模型已经能够取得不错的效果,因为即便我们不能从底层改写算法针对kappa进行优化,对MSE的优化也可以是kappa优化的一个很好的近似。所以我先拿reg:linear的Xgboost的举例。

首先Xgboost只接受矩阵的输入,所以我们把先前的数据转换成稀疏矩阵。在这里说明一点,无论你使用普通的model.matrix也好,还是sparse.model.matrix也好,在矩阵转换的时候请将测试集和训练集一同转换(rbind在一起或者其他方法),因为对于factor类型特征,分开转换会导致encode的方式不同,从而在建模中产生错误,这个也是我之前吃了亏的。

library(Matrix)
data_mat <- sparse.model.matrix(Response ~. -1 - Id, data)
train_idx <- 1:nrow(train)
test_idx <- (nrow(train) + 1):nrow(data)
train_mat <- data_mat[train_idx, ]
test_mat <- data_mat[test_idx, ]
library(xgboost)
train_DM <- xgb.DMatrix(data = train_mat, label = train$Response)

如下的代码可以做一次简单的cv:

param_linear <- list(
  objective = "reg:linear",
  booster = "gbtree",
  eta = 0.1, # Control the learning rate
  max.depth = 10, # Maximum depth of the tree
  min_child_weight = 35,
  subsample = 0.8, # subsample ratio of the training instance
  colsample_bytree = 0.5, # subsample ratio of columns when constructing each tree
  gamma = 0.01
)
set.seed(1126)
history <- xgb.cv(
  data = train_DM,
  params = param_linear,
  early.stop.round = 500, # training with a validation set will stop if the performance keeps getting worse consecutively for k rounds
  nround = 15000, # number of trees
  verbose = TRUE, 
  print.every.n = 10,
  nfold = 4, # number of CV folds
  maximize = FALSE # the lower the evaluation score the better
)

熟悉Xgboost的朋友应该都很清楚,参数verboseprint.every.nearly.stop.round对避免模型过拟合和模型性能评判的帮助很大。在调参方面,max.depthmin_child_weightsubsamplecolsample_bytree对模型性能很重要,eta理论上越小越好,但eta越小nround就需要越大,这个需要在性能和计算效率之间trade-off。我个人的感觉是eta最后调,初始eta可以略高,这样nround需求数就少,cv的效率高,grid search也好,一个个调也好(我发现一个个调的结果和grid search的结果其实差不多,难道Xgboost参数的local optima即为global optima?)。等确定其他参数,再降低eta提高nround,让它慢慢gradient descent。 当参数确定后,就可以建模了(这里的参数仅为示例,并非实际最优参数,但基本可以反映出真实比赛的情况):

set.seed(1126)
xgb_reg <- xgb.train(data = train_DM,
                     params = param_linear,
                     nround = 120 # number of trees
)
train_pred <- predict(xgb_reg, train_DM)
test_pred <- predict(xgb_reg, test_mat)

由于reg:linear的Xgboost仅仅做了回归,我们当然不能期望它输入等级或者整数的形式,在进一步处理之前,首先看一下训练集输入模型后输出的样子:

summary(train_pred)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.4499  4.5580  5.9850  5.6380  7.0190  8.4090
hist(train_pred)