这篇文章介绍一个R包{cutpointr} [1]，它可以进一步提供鲁棒方法用于估计最佳切割点，包括基于Bootstrap的重采样技术和基于核估计的平滑方法，以及广义加性模型、样条平滑和局部回归等。这个R包的使用可以让分割点的结果更少的依赖当下的研究样本，并且更加稳定。
首先安装R包，并且载入：
install.packages("cutpointr")
library(cutpointr)
library(ggplot2)
R包get！
将使用这个R包自带的数据集，名为suicide，查看下数据集概况：
summary(suicide)
str(suicide)
其中dsi为一个连续变量，为评估自杀倾向的一个测验分数，分数越高代表自杀倾向越大。suicide为既往的自杀企图，no代表没有，yes代表至少一次。age和gender为人口学变量，其中的gender为一个分类变量，用于后续进行分组寻找切割点。
这个例子的目的在于寻找dsi的最优切割点用于区分受试者的自杀企图（suicide变量）的情况。
下面举一个最简单的例子用于说明这个R包的用法，代码如下：
mycut <- cutpointr(data = suicide,
x = dsi,
class = suicide,
silent = TRUE)
mycut
代码的前三行依次指定了数据集(data)，用于寻找切割点的连续变量（x），以及分组变量（class）。
从上面的输出结果可知，通过默认的函数找到的最优切割点（optimal_cutpoint这栏）为2，意味着：如果样本的dsi分数大于等于2，则被认为是阳性，即至少有一次自杀企图。上面的结果也包含了敏感度和特异度之和（sum_sens_spec这栏）为1.75，AUC为0.92，以及患病率等信息。
上面使用的代码省略了一部分信息，现在将省略掉的代码重新填上，这样可以更加清楚的了解这个函数是如何工作的，代码如下：
mycut1 <- cutpointr(data = suicide,
x = dsi,
class = suicide,
pos_class = "yes",
neg_class = "no",
direction = ">=",
method = maximize_metric,
metric = sum_sens_spec)
mycut1
上面代码的中：
从上述的结果可知，和之前的结果是一模一样的。
可以使用summary()函数展示更加详细的结果，代码如下：
summary(mycut1)
与之前相比，这个结果的排版看起来更加的有条理。
下面的例子为这个R包的核心内容，即采用bootstrap重采样法，它的输出结果将包含了“袋子内” （in-bag）以及“袋子外”（out-of-bag）的信息，从而让结果不太受当下样本数据的影响，让结果更加稳定。下面，跑一个试试手，代码如下：
set.seed(123)
mycut_boot <- cutpointr(data = suicide,
x = dsi,
class = suicide,
boot_runs = 1000,
silent = TRUE)
summary(mycut_boot)
从代码可知，上述代码重采样了1000次。
上述结果的上半部分与之前是一样的，后半部分(Bootstrap summary)增加了一段新的信息。比如，包含了最优切割点(optimal_cutpoint)在1000个样本中的分布情况，它的中位数为2。其他变量的后缀包含了"_b"以及"_oob"，分别指代“袋子内”以及“袋子外”。
还可以将上面的结果进行作图，更加直观，代码如下：
plot(mycut_boot)
左上的图片为自变量(dsi)在两个组别中的分布；右上的图片为ROC曲线；左下的图片为bootstrap后的切割点分布图；右下的图片为“袋子外”的评价指标（敏感度+特异度）的分布。
也可以使用下方代码将上面的图片依次画出来，代码如下：
plot_x(mycut_boot)
plot_roc(mycut_boot)
plot_cut_boot(mycut_boot)
plot_metric_boot(mycut_boot)
还可以画出切割点和评价指标（敏感度和特异度之和）的关系图，并且包含bootstrap置信区间，代码如下：
plot_metric(mycut_boot) +
ylab("Sensitivity + Specifity")
这个R包还有一个非常赞的功能，就是添加亚组分析，这也是实战中经常需要的。
这个数据集包含一个变量为gender，现在就按照gender分组，分别计算各自的分割点，代码如下：
set.seed(321)
mycut_boot_group <- cutpointr(data = suicide,
x = dsi,
class = suicide,
subgroup = gender,
boot_runs = 1000,
boot_stratify = TRUE,
silent = TRUE)
summary(mycut_boot_group)
从结果可知，输出的信息分成两部分，女性和男性。
从样本内数据的结果可知（男性和女性结果中的上半部分），女性组的模型表现似乎更优。女性组的评价指标（敏感度+特异度）为1.8，而男性组的为1.6。
同样，在Bootstrap的结果中（男性和女性结果中的下半部分），女性组占优的趋势保持不变。
同样，可以将上面的分组结果进行作图，代码如下：
plot(mycut_boot_group)
从图中也可以做出判断，比如，右上角的ROC曲线中女性的曲线下面积更大；并且在右下角的图片中，女性的评价指标的分布更加的“窄”，提示变异程度更小。
下面，使用其他方法估计切割点，这个R包提供了共11种方法用于选择。可以使用代码?cutpointr查看，如下图所示：
下面举一个例子，将使用method = oc_youden_kernel， 代码如下：
set.seed(3215)
mycut_boot_group_kernal <- cutpointr(data = suicide,
x = dsi,
class = suicide,
subgroup = gender,
method = oc_youden_kernel,
metric = sum_sens_spec,
boot_runs = 1000,
boot_stratify = TRUE,
silent = TRUE)
summary(mycut_boot_group)
各位读者可以与之前采用默认方法的结果进行比较。
此外，也可以选择其他评价指标，它提供了n个指标，见下表：
文献[1] 表2
假设现在希望通过将误分类的代价降到最低，需要改变方法（method）以及评价指标(metric)，并且假定一个假阴性的误分类比一个假阳性的误分类要付出更高的代价，它们的比值为10:1，代码如下：
set.seed(3215)
mycut_boot_group_cost <- cutpointr(data = suicide,
x = dsi,
class = suicide,
subgroup = gender,
method = minimize_metric,
metric = misclassification_cost,
cost_fp = 1,
cost_fn = 10,
boot_runs = 1000,
boot_stratify = TRUE,
silent = TRUE)
summary(mycut_boot_group_cost)
从结果可知，女性和男性组中找到的切割点与之前的结果相比是非常接近的。
好啦，今天的内容就到这里。如果有帮助，记得分享给需要的人！
参考文献
[1] cutpointr: Improved Estimation and Validation of Optimal Cutpoints in R
公众号的线上课程
1. 《R语言和统计新手课程》
2. 《回归：从入门到进阶》
【通过公众号菜单栏--线上课程--新手课程/回归课程】
统计咨询
公众号核心成员的成果发表
公众号核心成员担任SCI杂志Associate Editor!
《Frontiers in Neuroscience, Frontiers in Neurology and Frontiers in Psychiatry杂志的神经退行性病变板块》
▌本文由R语言和统计首发
▌课程相关咨询可添加R师妹微信: kefu_rstats
▌编辑：June
▌邮箱：contact@rstats.cn
▌网站：www.rstats.cn
▌我们致力于让R语言和统计变得简单！
