互动演讲版本请点这里
# 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
包的目的是把一些常规建模(尤其是逻辑回归跟生存模型)过程当中用得比较多的模块流水化。把纷繁的建模过程分解成一个个小步骤,每个步骤对应一到两个方程,从而使得整个建模过程更加清晰,效率更高。
下面我们用一个简单的例子来说明这个包的使用。除了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)
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)来阐述整个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
方法的整个过程可以分解为以下的步骤
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)
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')
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')
把所有数据转换成WOE,以及检查完各个变量间的相关性后,我们可以开始建立逻辑回归模型
lg <- glm(status ~ ., dt.train, family=binomial(link='logit')) summary(lg)
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)
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
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)
coef2rate
对于非统计背景人士而言,另一个非常重要,但是比较难理解的概念是回归系数。在某些模型中,回归系数并不是跟因变量或自变量直接相关,而是跟转换后的变量相关(例如box-cox转换)。这给回归系数的理解带来一定的困难。为了让回归模型的结果更易于理解,我们可以利用得到的回归系数,来进一步计算各个组别的0/1预测概率。在计算预测概率的过程中,我们可以通过调节force.change参数来得到不同的结果:
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')
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.