In JianhuaHuang/streamlineR: Streamline Routine Modeling Work

流水化常规建模步骤：stremlineR包

互动演讲版本请点这里

# dir <- 'F:/Projects/Rpackage/streamlineR'
dir <- 'C:/Users/Jianhua/Dropbox/work_doc/Rpackage/streamlineR'
knitr::opts_chunk$set(echo = TRUE, root.dir = dir, warning=FALSE, message=FALSE)

这个streamlineR包的目的是把一些常规建模(尤其是逻辑回归跟生存模型）过程当中用得比较多的模块流水化。把纷繁的建模过程分解成一个个小步骤，每个步骤对应一到两个方程，从而使得整个建模过程更加清晰，效率更高。

安装R包

下面我们用一个简单的例子来说明这个包的使用。除了streamlineR这个包裹外，这个分析依赖于其他R包。在分析前，我们得先把这些包都安装上。

安装依赖包

# 如果某些镜像被屏蔽了，可以用下面的方程选择其他镜像进行安装
# chooseCRANmirror()
install.packages(c('dplyr', 'car', 'caret', 'e1071', 'knitr', 'reshape2',
  'corrplot','rpart', 'scales', 'shiny', 'survival', 'gridExtra', 'devtools', 
  'pec', 'MASS', 'pROC', 'manipulate'), dependencies = TRUE)

加载依赖包

安装完后，我们需要加载这些包：

# Load pacakges
sapply(c('dplyr', 'car', 'caret', 'e1071', 'knitr', 'reshape2', 'corrplot','rpart', 
  'scales', 'survival', 'gridExtra', 'devtools', 'pec', 'MASS', 'pROC', 
  'manipulate'), 
  require, character.only = TRUE)

从GitHub上安装streamlineR

安装了devtools包后，我们可以用 install_github 方程来安装streamlineR包

# 如果install_github不能用，你可以用下面的链接下载这个包到本地再安装
# https://api.github.com/repos/JianhuaHuang/streamlineR/zipball/master 
devtools::install_github('JianhuaHuang/streamlineR')
library(streamlineR)

数据准备(Survival PBC数据集)

在这个例子中，我们会使用survival包中的一个肝硬化数据集(PBC)来阐述整个streamlineR包的使用过程。关于这个数据集的详细信息可以看这里

加载数据

由于这个PBC数据集比较小，为了使后面的分析更加稳定，我们采用重复抽样的办法把数据量增加到10000条

dt <- survival::pbc %>%
  transmute(age = round(age), gender = sex, platelet, stage = as.character(stage), 
    time, status = as.numeric(status %in% c(1, 2))) %>%
  filter(!is.na(stage))

set.seed(1111)  # reproducible
dt <- dt[sample(nrow(dt), 10000, replace = T), ]
row.names(dt) <- NULL
dim(dt)
str(dt)
head(dt)

把数据分成训练数据跟测试数据

在做进一步分析前，整个数据集被分为两部分：训练数据跟测试数据。训练数据用于建模，而测试数据用于测试模型的可靠性。

set.seed(1111)
ind.train <- createDataPartition(dt$status, p = .7, list = FALSE)
dt.train <- dt[ind.train, ]
dt.test <- dt[-ind.train, ]
row.names(dt.train) <- NULL
row.names(dt.test) <- NULL
dim(dt.train)
dim(dt.test)

# preserve the original values
dt.train.bk <- dt.train
dt.test.bk <- dt.test

基于单变量回归的数据分箱: bin.knn

在建立回归模型之前，我们经常需要把一些与因变量存在非线性关系的数值型自变量进行一些转换，从而使得建模过程不受线性条件的约束。有很多方法可以实现这种转换，其中一种方法就是分箱。对数值型自变量分箱的目的是把邻近的，对因变量影响比较相似的值分成不同的组（箱），从而把数值型自变量转化成分类变量。

在streamlineR包中，主要有两个方程可以实现分箱操作：一个是基于回归系数的bin.knn，而另一个是基于决策树的bin.rpart。

bin.knn方法的整个过程可以分解为以下的步骤

1. 利用等量法（equal size），把待分箱的自变量的数值分成大概15-20个小组，每个小组有大概相等的数据量
2. 分成小组后，用这个自变量跟因变量建立回归模型
3. 获取各个小组的回归系数
4. 利用KNN算法，把相邻的且回归系数相近的小组进一步分成5个左右大组

bin.knn(status ~ platelet, data = dt.train, n.group = 5, min.bucket = 0.05)

结合manipulate包中的manipulate方程，通过调节n.group跟min.bucket两个参数，我们可以尝试把数据分成不同的组，从而得到满意的分箱效果。

manipulate(bin.knn(status ~ platelet, data = dt.train, n.group, min.bucket),
  n.group = slider(1, 10, step = 1, initial = 5, label = 'Number of Groups'),
  min.bucket = slider(0.01, .2, step = 0.01, initial = 0.05,
    label = 'Minimum Population'))

基于决策树的数据分箱: bin.rpart

bin.knn方程为数据分箱提供了直观的可视化过程，但是我们很难利用这个方法获得紧佳的分箱策略。而bin.rpart方程则弥补了这一缺陷。

决策树算法（递归分割）: rpart

rpart(recursivepartitioning，递归分割)是在决策树中广泛使用的一种算法，可以把数值变量一步步分成树状结构的节点。在rpart方程中，每个节点的划分都是基于最优算法，因此，我们最后得到的节点集也可以认为是最优的集合划分。基于rpart方程得到的分割点，我们可以进。基于rpart方程得到的分割点，bin.rpart方程进一步把数据分割成不同的箱。

rpart(formula = status ~ age, data = dt.train, 
  control = rpart.control(minbucket = .01 * nrow(dt.train)))

逻辑回归分箱

bin.rpart的使用跟rpart非常相似。两个主要的不同点就是rpart方程中的control参数被命名为rcontrol参数，以及加了一个n.group参数来控制箱子的个数。

下面的两个例子把年龄(age)跟血小板(platelet)进行了分箱操作，在rcontrol参数中，我们把minbucket设成了.01 * nrow(dt.train)，目的是保证每个箱子中包含至少1%的数据。

lg.bin.age <- bin.rpart(formula = status ~ age, data = dt.train, 
  rcontrol = rpart.control(minbucket = .01 * nrow(dt.train)))

str(lg.bin.age)

lg.bin.platelet <- bin.rpart(formula = status ~ platelet, data = dt.train, 
  rcontrol = rpart.control(minbucket = .01 * nrow(dt.train)))

生存模型分箱

相对于其他R包（例如smbinning跟woe），bin.rpart的一个优点是，它不仅仅可以对逻辑回归模型进行分箱，还可以对其他模型（例如生存模型）进行分箱操作。只要这个模型可以被 rpart识别，它就可以用于bin.rpart的分箱操作。

下面的例子展示了在生存模型中，对年龄进行分箱：

surv.bin.age <- bin.rpart(formula = Surv(time, status) ~ age, data = dt.train,
  rcontrol = rpart.control(minbucket = .05 * nrow(dt.train)))  ## cp = 0.05

自动寻找合适的分箱数

在rprat方程中，参数cp(complexity paramenter)可以用于控制决策树的层数，因此可以间接用来控制bin.rpart方程的分箱数目。在rpart中，cp的默认值是0.01。cp值越大，分箱数目越少，反之亦然。通过不断调整cp值，我们可以找到合适的分箱数目。为了自动化对cp的调整，我们引入了一个新的参数：n.group（分箱数目）。这个参数的值可以是一个数值（例如5）或值是一个向量（例如3：7，或者c(4, 6)）。输入这个参数以后，bin.rpart方程会自动调整cp值，直到分箱数目是n.group中的某一个值。

surv.bin.age2 <- bin.rpart(formula = Surv(time, status) ~ age, data = dt.train,
  rcontrol = rpart.control(minbucket = .05 * nrow(dt.train)， cp = .01), 
  n.group = 3:7)

替换数值型变量

bin.rpart方程会返会各个数值所对应的箱子，接着，我们可以用这些箱子代替掉原来的数值型变量。

# We don't need the time column anmore, delete it in both dt.train and dt.test
dt.train <- dplyr::select(dt.train, -time)  
dt.test <- dplyr::select(dt.test, -time)
head(dt.train)
dt.train$age <- lg.bin.age$bins
dt.train$platelet <- lg.bin.platelet$bins
head(dt.train)

分组统计： level.stat

对于因变量是二元变量的回归（例如逻辑回归跟生存模型），在建模前我们可以粗略了解一下各个组别间0/1分布的情况。level.stat方程可以用于统计各个组别的0/1频率，各个组的证据权重(Weight of Evidence (WOE)), 以及各个变量的信息值（Information Value (IV))。WOE跟IV值的计算可以参考这个例子。IV值越高，说明在这个变量里，组别间的0/1分布差别越大，因此这个变量在模型中就越重要。所以，在建模之前粗略统计一下各个变量的IV值有助于我们大概了解一下各个变量的重要性。

col.x <- c('age', 'gender', 'platelet', 'stage')
stat.train <- level.stat(dt.train, x = col.x, y = 'status', flag.0 = 0, flag.1 = 1)
head(stat.train)

可视化分组统计结果: ggstat

作为level.stat的一个附属方程，ggstat可以用于对分组统计的结果进行可视化。ggstat方程是基于ggplot作图系统，因此ggstat所作的图可以直接用ggsave方程保存。

利用默认值作图

ggstat方程提供了多个不同的参数。除了第一个参数以外，其他参数都有默认值。因此，要想ggstat正确运行，我们只需要提供第一个参数(data)的值。通常情况下，这个data参数的值就是我们用level.stat方程时所输出的结果。

ggstat(data = stat.train, var = "Variable.IV", x = "Group", y = "Rate.1", 
    y.label = "Perc.1", y.label.col = "red", y.title = NULL, 
    bar.col = "cornflowerblue", width = "Rate.group", width.label = "Perc.group", 
    width.label.col = "black", ncol = NULL, theme = "classic", 
    background = "white")

统一条形宽度

ggstat(stat.train, width = .2)

用WOE作图

ggstat(stat.train, y = 'WOE', y.label = 'WOE.round', width = .2, 
  width.label = NULL, ncol = 4)

改变字体跟图形颜色

为了快速找到可用的颜色，我们可以用display.col方程来显示常用的颜色名字。

display.col()  # reference to http://sape.inf.usi.ch/quick-reference/ggplot2/colour 

在ggstat中，我们可以直接用这些颜色来改变图的样式：

ggstat(stat.train, width = .2, y.label.col = 'white', bar.col = 'black',
  width.label.col = 'blue', background = 'green')

利用WOE代替分类变量: replace.woe

逻辑回归的一个重要性质是WOE值与log-odd存在完美的线性关系（补充资料）。因此，我们可以把分类变量的各个值代替成与之对应的WOE值，从而把分类变量进一步转换回数值变量。把分类变量代替成WOE可以降低回归方程的自由度，并且易于检测各个自变量间的相关性。replace.woe利用level.stat所得到的WOE值，来代替掉分类变量里的各个类别。

replace.woe(data = dt.train, stat = stat.train, replace = FALSE) %>%
  head

dt.train <- replace.woe(data = dt.train, stat = stat.train, 
  replace = TRUE)
head(dt.train)

自变量间的相关性：ggcorr

在建模前，我们通常也要粗略看一下各个自变量间的相关性，来去掉一部分高度相关的自变量。在这个包中我们主要用来了ggcorr函数。ggplot函数相当于corrplot包中的corrplot函数的简化版。通过调节几十个参数，corrplot函数可以提供多样的可视化效果；而ggcorr函数只提供了一些比较重要的调节参数，大大简化了 corrplot的功能。corrplot函数是基于R基础作图系统，而ggcorr是基于ggplot的作图系统，因此可以直接用ggsave保存结果。与corrplot相比，ggcorr的另一特点是去掉了对角线的点，因为这些点过于抢眼。

下面两幅图分别展示了corrplot跟ggcorr的默认效果：

cor.mat <- cor(dt.train[, col.x])
corrplot(cor.mat)

ggcorr(cor.mat)

利用随机生成的一些数据，我们可以看一下ggcorr生成的各种样式的图：

set.seed(1111)
data.random <- matrix(runif(100), 10)
colnames(data.random) <- paste('Variable', 1:10)
cor.random <- cor(data.random)

ggcorr(cor.random)  # default output
ggcorr(cor.random, var.position = 'diagonal', add.legend = FALSE)
ggcorr(cor.random, lower = TRUE)
ggcorr(cor.random, lower = TRUE, var.position = 'diagonal', high = 'blue', 
  low = 'green')

建立逻辑(logistic)回归模型

把所有数据转换成WOE，以及检查完各个变量间的相关性后，我们可以开始建立逻辑回归模型

全模型

lg <- glm(status ~ ., dt.train, family=binomial(link='logit'))
summary(lg)

分步(stepwise)变量筛选

lg.aic <- stepAIC(lg, k =  qchisq(0.05, 1, lower.tail=F))   # p to enter: 0.05
summary(lg.aic)

data.frame(vif(lg.aic))  # check the multicollinearity between predictors

测试集数据准备: bin.custom & replace.woe

以上的所有分析都是基于训练集的数据。为了测试模型的预测性能，我们可以用测试集的数据来进行检验。

数据分箱: bin.custom

这里，我们主要用到了bin.custom方程。bin.custom方程可以根据给定的分割点，对数值变量进行分箱。因为我们已经利用训练集的数据找到了最佳的分割点，我们可以直接用这些分割点来对训练集的数据进行分箱操作。

head(dt.test)
dt.test$age <-  bin.custom(dt.test$age, cut.p = lg.bin.age$cut.points)
dt.test$platelet <- bin.custom(dt.test$platelet, cut.p = lg.bin.platelet$cut.points)
head(dt.test)

利用WOE代替分类变量： replace.woe

接着，我们可样可以用WOE来代替分箱后的结果

dt.test <- replace.woe(dt.test, stat = stat.train, replace = TRUE)
head(dt.test)

检测模型性能: perf.auc & perf.decile

把分类变量转换成WOE后，我们可以用其测试模型的性能。在这个包里，我们主要提供两种不同的方法来对检测模型的性能：perf.auc以及perf.decile

基于AUC的性能检测: perf.auc

对于逻辑回归，perf.auc方程会画出训练集跟测试集的ROC(Receiver Operating Characteristic)曲线，并标注各自的AUC(Area Under Curve)。对于生存模型，perf.auc会画出各个生存时间的AUC，并给出整个综合的AUC值(iAUC)

perf.auc(model = lg.aic, dt.train, dt.test)

基于十分位的性能检测: perf.decile

尽管AUC在统计中广泛应用，在非统计专业人士眼中，它仍然是比较难以理解的概念。perf.decile方程提供了一种更简单且直观的方法来描述模型的性能。perf.decile方程通过以下几个步骤来比较预测值跟真实值的差别： 1. 对所有预测的概率进行排序，把排序后的数据等量分成十组(decile) 2. 计算每一组的平均预测值跟真实值 3. 利用预测值跟真实值作为X跟Y轴来作图，并加上对角线 在这个图中，我们可以非常直观的看到预测值跟真实值是否吻合。对于没有统计背景的人，这种图也是非常浅显易懂。

pred.test <- predict(lg.aic, newdata = dt.test, type = 'response')
perf.decile(actual = dt.test$status, pred = pred.test, add.legend = TRUE)

回归系数转换成原始0/1率: coef2rate

对于非统计背景人士而言，另一个非常重要，但是比较难理解的概念是回归系数。在某些模型中，回归系数并不是跟因变量或自变量直接相关，而是跟转换后的变量相关（例如box-cox转换）。这给回归系数的理解带来一定的困难。为了让回归模型的结果更易于理解，我们可以利用得到的回归系数，来进一步计算各个组别的0/1预测概率。在计算预测概率的过程中，我们可以通过调节force.change参数来得到不同的结果：

• 如果force.change是FALSE: 给定的数据会被直接用来计算概率，然后进一步计算各组的平均值。
• 如果force.change是TRUE: 这个概念源于对回归系数的解释：在所有其他变量保持不变的情况下，如果某个变量由A1变成A2，因变量的值会怎么变化。因此，当force.change设成TRUE后，这个coef2rate方程逐个改变自变量的值（其他自变量保持不变），用所有数据来预测各个组别之间预测值的差别。

pred.stat <- coef2rate(model = lg.aic, data = dt.test, 
  stat = stat.train, force.change = TRUE)
head(pred.stat)

计算到预测值后，我们同样可以用ggstat方程来看预测值的区别：

ggstat(pred.stat, y = 'Pred.Rate.1',  y.label = 'Pred.Perc.1')

参考(reference)

• streamlineR package information: https://github.com/JianhuaHuang/streamlineR
• Submit bug: http://github.com/JianhuaHuang/streamlineR/issues

JianhuaHuang/streamlineR documentation built on May 7, 2019, 10:40 a.m.