https://github.com/cran/cutpointr
Tip revision: 2900dc24d2c5a7d8fdb3f1abb1540fb704e51742 authored by Christian Thiele on 15 February 2021, 13:40:03 UTC
version 1.1.0
version 1.1.0
Tip revision: 2900dc2
cutpointr.R
## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(fig.width = 6, fig.height = 5, fig.align = "center")
options(rmarkdown.html_vignette.check_title = FALSE)
load("vignettedata/vignettedata.Rdata")
## ----CRAN, eval = FALSE-------------------------------------------------------
# install.packages("cutpointr")
## -----------------------------------------------------------------------------
library(cutpointr)
data(suicide)
head(suicide)
cp <- cutpointr(suicide, dsi, suicide,
method = maximize_metric, metric = sum_sens_spec)
## -----------------------------------------------------------------------------
summary(cp)
## -----------------------------------------------------------------------------
plot(cp)
## -----------------------------------------------------------------------------
opt_cut <- cutpointr(suicide, dsi, suicide, direction = ">=", pos_class = "yes",
neg_class = "no", method = maximize_metric, metric = youden)
## -----------------------------------------------------------------------------
plot_metric(opt_cut)
## -----------------------------------------------------------------------------
predict(opt_cut, newdata = data.frame(dsi = 0:5))
## ----separate subgroups and bootstrapping, eval = FALSE-----------------------
# set.seed(12)
# opt_cut_b <- cutpointr(suicide, dsi, suicide, boot_runs = 1000)
## -----------------------------------------------------------------------------
opt_cut_b
## -----------------------------------------------------------------------------
opt_cut$boot
## -----------------------------------------------------------------------------
summary(opt_cut)
plot(opt_cut)
## ---- eval = FALSE------------------------------------------------------------
# library(doParallel)
# cl <- makeCluster(2) # 2 cores
# registerDoParallel(cl)
# registerDoRNG(12) # Reproducible parallel loops using doRNG
# opt_cut <- cutpointr(suicide, dsi, suicide, gender, pos_class = "yes",
# direction = ">=", boot_runs = 1000, allowParallel = TRUE)
# stopCluster(cl)
# opt_cut
## ---- cache=TRUE--------------------------------------------------------------
set.seed(100)
cutpointr(suicide, dsi, suicide, gender,
method = maximize_boot_metric,
boot_cut = 200, summary_func = mean,
metric = accuracy, silent = TRUE)
## -----------------------------------------------------------------------------
opt_cut <- cutpointr(suicide, dsi, suicide, gender, method = minimize_metric,
metric = misclassification_cost, cost_fp = 1, cost_fn = 10)
## -----------------------------------------------------------------------------
plot_metric(opt_cut)
## ---- message = FALSE---------------------------------------------------------
opt_cut <- cutpointr(suicide, dsi, suicide, gender,
method = minimize_loess_metric,
criterion = "aicc", family = "symmetric",
degree = 2, user.span = 0.7,
metric = misclassification_cost, cost_fp = 1, cost_fn = 10)
## -----------------------------------------------------------------------------
plot_metric(opt_cut)
## -----------------------------------------------------------------------------
library(ggplot2)
exdat <- iris
exdat <- exdat[exdat$Species != "setosa", ]
opt_cut <- cutpointr(exdat, Petal.Length, Species,
method = minimize_gam_metric,
formula = m ~ s(x.sorted, bs = "cr"),
metric = abs_d_sens_spec)
plot_metric(opt_cut)
## -----------------------------------------------------------------------------
opt_cut <- cutpointr(suicide, dsi, suicide, gender,
method = minimize_spline_metric, spar = 0.4,
metric = misclassification_cost, cost_fp = 1, cost_fn = 10)
plot_metric(opt_cut)
## -----------------------------------------------------------------------------
cutpointr(suicide, dsi, suicide, gender, method = oc_youden_normal)
## -----------------------------------------------------------------------------
cutpointr(suicide, dsi, suicide, gender, method = oc_youden_kernel)
## ---- fig.width=4, fig.height=3-----------------------------------------------
roc_curve <- roc(data = suicide, x = dsi, class = suicide,
pos_class = "yes", neg_class = "no", direction = ">=")
auc(roc_curve)
head(roc_curve)
plot_roc(roc_curve)
## -----------------------------------------------------------------------------
dat <- data.frame(outcome = c("neg", "neg", "neg", "pos", "pos", "pos", "pos"),
pred = c(1, 2, 3, 8, 11, 11, 12))
## -----------------------------------------------------------------------------
opt_cut <- cutpointr(dat, x = pred, class = outcome, use_midpoints = TRUE)
plot_x(opt_cut)
## ---- echo = FALSE------------------------------------------------------------
plotdat_nomidpoints <- structure(list(sim_nr = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L,
6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 8L,
8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 10L,
10L, 10L, 10L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 11L, 11L,
11L, 11L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 13L, 13L, 13L,
13L, 13L, 13L, 13L, 13L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L,
15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 16L, 16L, 16L, 16L, 16L,
16L, 16L, 16L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 18L, 18L,
18L, 18L, 18L, 18L, 18L, 18L, 19L, 19L, 19L, 19L, 19L, 19L, 19L,
19L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 21L, 21L, 21L, 21L,
21L, 21L, 21L, 21L, 22L, 22L, 22L, 22L, 22L, 22L, 22L, 22L, 23L,
23L, 23L, 23L, 23L, 23L, 23L, 23L, 24L, 24L, 24L, 24L, 24L, 24L,
24L, 24L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 26L, 26L, 26L,
26L, 26L, 26L, 26L, 26L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L,
28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 29L, 29L, 29L, 29L, 29L,
29L, 29L, 29L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 31L, 31L,
31L, 31L, 31L, 31L, 31L, 31L, 32L, 32L, 32L, 32L, 32L, 32L, 32L,
32L, 33L, 33L, 33L, 33L, 33L, 33L, 33L, 33L, 34L, 34L, 34L, 34L,
34L, 34L, 34L, 34L, 35L, 35L, 35L, 35L, 35L, 35L, 35L, 35L, 36L,
36L, 36L, 36L, 36L, 36L, 36L, 36L), method = structure(c(1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L), .Label = c("emp",
"normal", "loess", "boot", "spline", "spline_20", "kernel", "gam"
), class = "factor"), n = c(30, 30, 30, 30, 30, 30, 30, 30, 50,
50, 50, 50, 50, 50, 50, 50, 75, 75, 75, 75, 75, 75, 75, 75, 100,
100, 100, 100, 100, 100, 100, 100, 150, 150, 150, 150, 150, 150,
150, 150, 250, 250, 250, 250, 250, 250, 250, 250, 500, 500, 500,
500, 500, 500, 500, 500, 750, 750, 750, 750, 750, 750, 750, 750,
1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 30, 30, 30, 30,
30, 30, 30, 30, 50, 50, 50, 50, 50, 50, 50, 50, 75, 75, 75, 75,
75, 75, 75, 75, 100, 100, 100, 100, 100, 100, 100, 100, 150,
150, 150, 150, 150, 150, 150, 150, 250, 250, 250, 250, 250, 250,
250, 250, 500, 500, 500, 500, 500, 500, 500, 500, 750, 750, 750,
750, 750, 750, 750, 750, 1000, 1000, 1000, 1000, 1000, 1000,
1000, 1000, 30, 30, 30, 30, 30, 30, 30, 30, 50, 50, 50, 50, 50,
50, 50, 50, 75, 75, 75, 75, 75, 75, 75, 75, 100, 100, 100, 100,
100, 100, 100, 100, 150, 150, 150, 150, 150, 150, 150, 150, 250,
250, 250, 250, 250, 250, 250, 250, 500, 500, 500, 500, 500, 500,
500, 500, 750, 750, 750, 750, 750, 750, 750, 750, 1000, 1000,
1000, 1000, 1000, 1000, 1000, 1000, 30, 30, 30, 30, 30, 30, 30,
30, 50, 50, 50, 50, 50, 50, 50, 50, 75, 75, 75, 75, 75, 75, 75,
75, 100, 100, 100, 100, 100, 100, 100, 100, 150, 150, 150, 150,
150, 150, 150, 150, 250, 250, 250, 250, 250, 250, 250, 250, 500,
500, 500, 500, 500, 500, 500, 500, 750, 750, 750, 750, 750, 750,
750, 750, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000), mean_err = c(0.532157164015659,
0.0344907054484091, 1.09430750651166, 0.847845162156675, 1.72337372126503,
0.893756658507988, 0.0430309247027736, 0.785821459035346, 0.368063404388512,
0.0256197760404459, 0.54480529648463, 0.54385597929651, 0.657325657699579,
0.578611116865437, 0.0400491342691897, 0.515688005217413, 0.256713912589642,
0.0444582875885996, 0.326975493112402, 0.371128780921122, 0.473515115741104,
0.389519558405289, 0.105044360789378, 0.301924717299333, 0.207750921776918,
-0.00318128936770314, 0.215170156089776, 0.27218780048926, 0.260519564021842,
0.236792923882582, 0.0209319074923902, 0.232957055204834, 0.0726605917614469,
-0.00282823355849125, 0.0753216783313991, 0.147121931849656,
0.0986417724955371, 0.10048009778446, -0.0117861260923649, 0.0650845904350442,
0.0985144485083747, 0.00601003227249322, 0.107439979908118, 0.120777421732797,
0.098470489820427, 0.0940946984227826, 0.0340166854141625, 0.107851118082414,
0.0249685210781582, -0.00275219600614378, 0.0258069390207201,
0.0303381972274654, -0.000994602151869198, 0.00196854833683764,
-0.0172319489562159, 0.0230957871932473, 0.00787486424680835,
-0.018438041997315, -0.000808033567394628, -0.00151904153864496,
-0.0258118523805697, -0.020984156892953, -0.0411584927473141,
-0.0462075435919094, 0.0481149217843661, -0.0115241085997692,
0.0278045708419731, 0.0358588316426625, 0.0424473909450939, 0.0379773233328197,
0.0298772985879321, 0.0494939492036379, 0.561286668337778, 0.0210874502384648,
0.607711822769155, 0.944733906256477, 1.32069801051061, 0.623604782428556,
0.0138075109769806, 0.640859854412358, 0.284873604303057, -0.0170357985365701,
0.273426417633118, 0.524432737895336, 0.355003110979807, 0.312837607951434,
-0.0316296929553873, 0.270109834098986, 0.174110335581819, 0.0253101719615279,
0.199956222742702, 0.375485416120771, 0.278956745944806, 0.244245525945888,
0.0325126314233263, 0.253352659868514, 0.154840760004461, 0.00231589709639472,
0.154412480165179, 0.264847742842386, 0.196572744185608, 0.182934783774992,
0.0207139021497755, 0.17351041376412, 0.14190910348156, -0.000766834484010096,
0.159975205214477, 0.191222128926019, 0.119768252112669, 0.12033372914036,
-0.00429047209392067, 0.120982527821078, 0.0756304869484406,
-0.00890884219048113, 0.0727782693168392, 0.118690444738942,
0.0814898789647033, 0.0799724348001957, 0.0182926240912726, 0.0887155007804252,
0.00799604720502299, 0.000599148388616836, -0.00567769035990384,
0.0358412514670032, 0.0308474979074875, 0.0341668723768997, -0.000180318451026095,
0.0180733341290925, 0.00456876236626807, 0.00150574966485876,
0.011152095953916, 0.0176039119729626, 0.00608274255434991, 0.0146257828313115,
-0.0108877417404102, 0.00341198000323035, -0.00198459880370283,
0.0026551895445694, 0.00199040664534129, 0.0150165794544221,
0.00646287144368147, 0.00999205240904708, -0.00850278571195971,
0.000833666619266177, 0.714730067273087, -0.00916546079360956,
0.662799490366986, 1.18552468844156, 1.25901933062308, 0.672701515532179,
0.0311066140197676, 0.699068058809396, 0.451908043813962, 0.0131716226592205,
0.429551369887738, 0.697928133235757, 0.445367768988423, 0.408463448982185,
0.0318707241721211, 0.406284953982951, 0.257736307754364, -0.00525924423719458,
0.236977000055322, 0.429144726596141, 0.291381752107184, 0.267557606428613,
0.0103657879176852, 0.254728590646094, 0.187783398487578, -0.00216064381479362,
0.209025103860707, 0.318293592390017, 0.216751610346408, 0.195630579126633,
-0.00355723971644246, 0.174111826408428, 0.151010324964235, 0.0152409223899092,
0.159002511320467, 0.214643583389694, 0.136211731513269, 0.138948149207635,
0.00736196817594524, 0.115637867729083, 0.0491348055596302, -0.00133957946235943,
0.0507437758212659, 0.103956325245849, 0.0641182216839426, 0.0721933081297794,
-0.0124376134651938, 0.0632317888879588, 0.0322195712438111,
0.00170122889182022, 0.0287526766624194, 0.0589662164030242,
0.0348535721709848, 0.039527944642463, -0.00617539706415593,
0.0274246010641889, 0.0325877909680824, 0.00530528253248245,
0.0221776555499961, 0.0389702052631117, 0.0221602091288215, 0.0254478639695596,
-0.0016189234058987, 0.0197144417326668, -0.00632485262604172,
-0.00364979854195596, -0.00276076468388984, 0.0126267527301874,
0.0123592498266038, 0.0154921632247644, -0.00591512196680815,
0.0098685016547149, 1.19276916750486, -0.0296831640401583, 0.99406393593888,
1.95758669445116, 1.45842790978446, 0.916899913902239, -0.0240222410217233,
1.00771193034927, 0.748151091428865, 0.001671855025917, 0.665180535306263,
1.1777049634557, 0.578603609273264, 0.546625362141714, 0.0292152981607387,
0.615230814912951, 0.417886753756131, 0.00324593885807739, 0.406076310942717,
0.732191741449251, 0.352684769616612, 0.326901376027897, -0.000759357576337989,
0.350075431324921, 0.310927617707656, -0.0107255472998434, 0.28102101085112,
0.514683023356017, 0.24913510139508, 0.235155452507568, -0.0220885572014814,
0.243370611433649, 0.209652330609093, -0.00502865663759991, 0.2172246261346,
0.356540958804122, 0.172121720418057, 0.17487914828986, 0.00365942442127361,
0.176594681455494, 0.126956927057327, -0.00270525933073803, 0.120116234221594,
0.210827536708082, 0.101520409193932, 0.101379097920023, 0.00238043252144371,
0.113027315928011, 0.0598624378953727, -0.00538838415690431,
0.0568400730102315, 0.0978115258288965, 0.0454207684906316, 0.0473140143579152,
-0.00165813015281622, 0.0521772135812508, 0.0530224090961669,
-0.000416993405198653, 0.0353236911458531, 0.0605493601241619,
0.0316204159297213, 0.0344789374555544, -0.00446984887315054,
0.0328807595695966, 0.0396438546423947, -0.00331466369719113,
0.0379029847219126, 0.0572435100638761, 0.0253269328104989, 0.0235663211070417,
0.00220241478536399, 0.0307132312422208), youden = c(0.2, 0.2,
0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2,
0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2,
0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2,
0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2,
0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2,
0.2, 0.2, 0.2, 0.2, 0.2, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4,
0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4,
0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4,
0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4,
0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4,
0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.6,
0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6,
0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6,
0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6,
0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6,
0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6,
0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8,
0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8,
0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8,
0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8,
0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8,
0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8
)), row.names = c(NA, -288L), group_sizes = c(8L, 8L, 8L, 8L,
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L,
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L
), biggest_group_size = 8L, class = c("grouped_df", "tbl_df",
"tbl", "data.frame"), groups = structure(list(sim_nr = 1:36,
.rows = list(1:8, 9:16, 17:24, 25:32, 33:40, 41:48, 49:56,
57:64, 65:72, 73:80, 81:88, 89:96, 97:104, 105:112, 113:120,
121:128, 129:136, 137:144, 145:152, 153:160, 161:168,
169:176, 177:184, 185:192, 193:200, 201:208, 209:216,
217:224, 225:232, 233:240, 241:248, 249:256, 257:264,
265:272, 273:280, 281:288)), row.names = c(NA, -36L), class = c("tbl_df",
"tbl", "data.frame"), .drop = TRUE))
## ---- echo = FALSE------------------------------------------------------------
library(dplyr)
ggplot(plotdat_nomidpoints %>% filter(!(method %in% c("spline_20"))),
aes(x = n, y = mean_err, color = method, shape = method)) +
geom_line() + geom_point() +
facet_wrap(~ youden, scales = "fixed") +
scale_shape_manual(values = 1:nlevels(plotdat_nomidpoints$method)) +
scale_x_log10(breaks = c(30, 50, 75, 100, 150, 250, 500, 750, 1000)) +
ggtitle("Bias of all methods when use_midpoints = FALSE",
"normally distributed data, 10000 repetitions of simulation")
## -----------------------------------------------------------------------------
opt_cut <- cutpointr(suicide, dsi, suicide, metric = sum_sens_spec,
tol_metric = 0.05, break_ties = c)
library(tidyr)
opt_cut %>%
select(optimal_cutpoint, sum_sens_spec) %>%
unnest(cols = c(optimal_cutpoint, sum_sens_spec))
## ---- eval = FALSE------------------------------------------------------------
# set.seed(100)
# opt_cut_manual <- cutpointr(suicide, dsi, suicide, method = oc_manual,
# cutpoint = mean(suicide$dsi), boot_runs = 1000)
# set.seed(100)
# opt_cut_mean <- cutpointr(suicide, dsi, suicide, method = oc_mean, boot_runs = 1000)
## ---- eval = FALSE------------------------------------------------------------
# myvar <- "dsi"
# cutpointr(suicide, !!myvar, suicide)
## -----------------------------------------------------------------------------
mcp <- multi_cutpointr(suicide, class = suicide, pos_class = "yes",
use_midpoints = TRUE, silent = TRUE)
summary(mcp)
## ---- eval = FALSE, message = FALSE-------------------------------------------
# set.seed(123)
# opt_cut_b_g <- cutpointr(suicide, dsi, suicide, gender, boot_runs = 1000)
## ---- message = FALSE---------------------------------------------------------
# Using dplyr and tidyr
library(tidyr)
opt_cut_b_g %>%
group_by(subgroup) %>%
select(subgroup, boot) %>%
unnest(cols = boot) %>%
summarise(sd_oc_boot = sd(optimal_cutpoint),
m_oc_boot = mean(optimal_cutpoint),
m_acc_oob = mean(acc_oob))
## -----------------------------------------------------------------------------
cutpointr(suicide, dsi, suicide, gender, metric = youden, silent = TRUE) %>%
add_metric(list(ppv, npv)) %>%
select(subgroup, optimal_cutpoint, youden, ppv, npv)
## -----------------------------------------------------------------------------
roc(data = suicide, x = dsi, class = suicide, pos_class = "yes",
neg_class = "no", direction = ">=") %>%
add_metric(list(cohens_kappa, F1_score)) %>%
select(x.sorted, tp, fp, tn, fn, cohens_kappa, F1_score) %>%
head()
## ---- eval = FALSE------------------------------------------------------------
# mean_cut <- function(data, x, ...) {
# oc <- mean(data[[x]])
# return(data.frame(optimal_cutpoint = oc))
# }
## -----------------------------------------------------------------------------
misclassification_cost
## ---- fig.width=4, fig.height=3-----------------------------------------------
plot_cut_boot(opt_cut_b_g)
plot_metric(opt_cut_b_g, conf_lvl = 0.9)
plot_metric_boot(opt_cut_b_g)
plot_precision_recall(opt_cut_b_g)
plot_sensitivity_specificity(opt_cut_b_g)
plot_roc(opt_cut_b_g)
## ---- fig.width=4, fig.height=3-----------------------------------------------
p <- plot_x(opt_cut_b_g)
p + ggtitle("Distribution of dsi") + theme_minimal() + xlab("Depression score")
## ---- fig.width=4, fig.height=3, cache=FALSE----------------------------------
plot_cutpointr(opt_cut_b, xvar = cutpoints, yvar = sum_sens_spec, conf_lvl = 0.9)
plot_cutpointr(opt_cut_b, xvar = fpr, yvar = tpr, aspect_ratio = 1, conf_lvl = 0)
plot_cutpointr(opt_cut_b, xvar = cutpoint, yvar = tp, conf_lvl = 0.9) + geom_point()
## ---- fig.width=4, fig.height=3-----------------------------------------------
opt_cut_b_g %>%
select(data, subgroup) %>%
unnest(cols = data) %>%
ggplot(aes(x = suicide, y = dsi)) +
geom_boxplot(alpha = 0.3) + facet_grid(~subgroup)
## ---- eval = FALSE------------------------------------------------------------
# # Return cutpoint that maximizes the sum of sensitivity and specificiy
# # ROCR package
# rocr_sensspec <- function(x, class) {
# pred <- ROCR::prediction(x, class)
# perf <- ROCR::performance(pred, "sens", "spec")
# sens <- slot(perf, "y.values")[[1]]
# spec <- slot(perf, "x.values")[[1]]
# cut <- slot(perf, "alpha.values")[[1]]
# cut[which.max(sens + spec)]
# }
#
# # pROC package
# proc_sensspec <- function(x, class) {
# r <- pROC::roc(class, x, algorithm = 2, levels = c(0, 1), direction = "<")
# pROC::coords(r, "best", ret="threshold", transpose = FALSE)[1]
# }
## ---- eval = FALSE, echo = FALSE----------------------------------------------
# library(OptimalCutpoints)
# library(ThresholdROC)
# n <- 100
# set.seed(123)
# dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
# x_pos <- dat$x[dat$y == 1]
# x_neg <- dat$x[dat$y == 0]
# bench_100 <- microbenchmark::microbenchmark(
# cutpointr(dat, x, y, pos_class = 1, neg_class = 0,
# direction = ">=", metric = youden, break_ties = mean),
# rocr_sensspec(dat$x, dat$y),
# proc_sensspec(dat$x, dat$y),
# optimal.cutpoints(X = "x", status = "y", tag.healthy = 0, methods = "Youden",
# data = dat),
# thres2(k1 = x_neg, k2 = x_pos, rho = 0.5,
# method = "empirical", ci = FALSE),
# times = 100, unit = "ms"
# )
#
# n <- 1000
# set.seed(123)
# dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
# x_pos <- dat$x[dat$y == 1]
# x_neg <- dat$x[dat$y == 0]
# bench_1000 <- microbenchmark::microbenchmark(
# cutpointr(dat, x, y, pos_class = 1, neg_class = 0,
# direction = ">=", metric = youden, break_ties = mean),
# rocr_sensspec(dat$x, dat$y),
# proc_sensspec(dat$x, dat$y),
# optimal.cutpoints(X = "x", status = "y", tag.healthy = 0, methods = "Youden",
# data = dat),
# thres2(k1 = x_neg, k2 = x_pos, rho = 0.5,
# method = "empirical", ci = FALSE),
# times = 100, unit = "ms"
# )
#
# n <- 10000
# set.seed(123)
# dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
# x_pos <- dat$x[dat$y == 1]
# x_neg <- dat$x[dat$y == 0]
# bench_10000 <- microbenchmark::microbenchmark(
# cutpointr(dat, x, y, pos_class = 1, neg_class = 0,
# direction = ">=", metric = youden, break_ties = mean, silent = TRUE),
# rocr_sensspec(dat$x, dat$y),
# optimal.cutpoints(X = "x", status = "y", tag.healthy = 0, methods = "Youden",
# data = dat),
# proc_sensspec(dat$x, dat$y),
# thres2(k1 = x_neg, k2 = x_pos, rho = 0.5,
# method = "empirical", ci = FALSE),
# times = 100
# )
#
# n <- 1e5
# set.seed(123)
# dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
# bench_1e5 <- microbenchmark::microbenchmark(
# cutpointr(dat, x, y, pos_class = 1, neg_class = 0,
# direction = ">=", metric = youden, break_ties = mean),
# rocr_sensspec(dat$x, dat$y),
# proc_sensspec(dat$x, dat$y),
# times = 100, unit = "ms"
# )
#
# n <- 1e6
# set.seed(123)
# dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
# bench_1e6 <- microbenchmark::microbenchmark(
# cutpointr(dat, x, y, pos_class = 1, neg_class = 0,
# direction = ">=", metric = youden, break_ties = mean),
# rocr_sensspec(dat$x, dat$y),
# proc_sensspec(dat$x, dat$y),
# times = 30, unit = "ms"
# )
#
# n <- 1e7
# set.seed(123)
# dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
# bench_1e7 <- microbenchmark::microbenchmark(
# cutpointr(dat, x, y, pos_class = 1, neg_class = 0,
# direction = ">=", metric = youden, break_ties = mean),
# rocr_sensspec(dat$x, dat$y),
# proc_sensspec(dat$x, dat$y),
# times = 30, unit = "ms"
# )
#
# results <- rbind(
# data.frame(time = summary(bench_100)$median,
# Solution = summary(bench_100)$expr,
# n = 100),
# data.frame(time = summary(bench_1000)$median,
# Solution = summary(bench_1000)$expr,
# n = 1000),
# data.frame(time = summary(bench_10000)$median,
# Solution = summary(bench_10000)$expr,
# n = 10000),
# data.frame(time = summary(bench_1e5)$median,
# Solution = summary(bench_1e5)$expr,
# n = 1e5),
# data.frame(time = summary(bench_1e6)$median,
# Solution = summary(bench_1e6)$expr,
# n = 1e6),
# data.frame(time = summary(bench_1e7)$median,
# Solution = summary(bench_1e7)$expr,
# n = 1e7)
# )
# results$Solution <- as.character(results$Solution)
# results$Solution[grep(pattern = "cutpointr", x = results$Solution)] <- "cutpointr"
# results$Solution[grep(pattern = "rocr", x = results$Solution)] <- "ROCR"
# results$Solution[grep(pattern = "optimal", x = results$Solution)] <- "OptimalCutpoints"
# results$Solution[grep(pattern = "proc", x = results$Solution)] <- "pROC"
# results$Solution[grep(pattern = "thres", x = results$Solution)] <- "ThresholdROC"
#
# results$task <- "Cutpoint Estimation"
## ---- echo = FALSE------------------------------------------------------------
# These are the original results on our system
# dput(results)
results <- structure(list(time = c(4.5018015, 1.812802, 0.662101, 2.2887015,
1.194301, 4.839401, 2.1764015, 0.981001, 45.0568005, 36.2398515,
8.5662515, 5.667101, 2538.612001, 4.031701, 2503.8012505, 45.384501,
43.118751, 37.150151, 465.003201, 607.023851, 583.0950005, 5467.332801,
7850.2587, 7339.356101), Solution = c("cutpointr", "ROCR", "pROC",
"OptimalCutpoints", "ThresholdROC", "cutpointr", "ROCR", "pROC",
"OptimalCutpoints", "ThresholdROC", "cutpointr", "ROCR", "OptimalCutpoints",
"pROC", "ThresholdROC", "cutpointr", "ROCR", "pROC", "cutpointr",
"ROCR", "pROC", "cutpointr", "ROCR", "pROC"), n = c(100, 100,
100, 100, 100, 1000, 1000, 1000, 1000, 1000, 10000, 10000, 10000,
10000, 10000, 1e+05, 1e+05, 1e+05, 1e+06, 1e+06, 1e+06, 1e+07,
1e+07, 1e+07), task = c("Cutpoint Estimation", "Cutpoint Estimation",
"Cutpoint Estimation", "Cutpoint Estimation", "Cutpoint Estimation",
"Cutpoint Estimation", "Cutpoint Estimation", "Cutpoint Estimation",
"Cutpoint Estimation", "Cutpoint Estimation", "Cutpoint Estimation",
"Cutpoint Estimation", "Cutpoint Estimation", "Cutpoint Estimation",
"Cutpoint Estimation", "Cutpoint Estimation", "Cutpoint Estimation",
"Cutpoint Estimation", "Cutpoint Estimation", "Cutpoint Estimation",
"Cutpoint Estimation", "Cutpoint Estimation", "Cutpoint Estimation",
"Cutpoint Estimation")), row.names = c(NA, -24L), class = "data.frame")
## ---- eval = FALSE------------------------------------------------------------
# # ROCR package
# rocr_roc <- function(x, class) {
# pred <- ROCR::prediction(x, class)
# perf <- ROCR::performance(pred, "sens", "spec")
# return(NULL)
# }
#
# # pROC package
# proc_roc <- function(x, class) {
# r <- pROC::roc(class, x, algorithm = 2, levels = c(0, 1), direction = "<")
# return(NULL)
# }
## ---- eval = FALSE, echo = FALSE----------------------------------------------
# n <- 100
# set.seed(123)
# dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
# bench_100 <- microbenchmark::microbenchmark(
# cutpointr::roc(dat, "x", "y", pos_class = 1,
# neg_class = 0, direction = ">="),
# rocr_roc(dat$x, dat$y),
# proc_roc(dat$x, dat$y),
# times = 100, unit = "ms"
# )
# n <- 1000
# set.seed(123)
# dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
# bench_1000 <- microbenchmark::microbenchmark(
# cutpointr::roc(dat, "x", "y", pos_class = 1, neg_class = 0,
# direction = ">="),
# rocr_roc(dat$x, dat$y),
# proc_roc(dat$x, dat$y),
# times = 100, unit = "ms"
# )
# n <- 10000
# set.seed(123)
# dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
# bench_10000 <- microbenchmark::microbenchmark(
# cutpointr::roc(dat, "x", "y", pos_class = 1, neg_class = 0,
# direction = ">="),
# rocr_roc(dat$x, dat$y),
# proc_roc(dat$x, dat$y),
# times = 100, unit = "ms"
# )
# n <- 1e5
# set.seed(123)
# dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
# bench_1e5 <- microbenchmark::microbenchmark(
# cutpointr::roc(dat, "x", "y", pos_class = 1, neg_class = 0,
# direction = ">="),
# rocr_roc(dat$x, dat$y),
# proc_roc(dat$x, dat$y),
# times = 100, unit = "ms"
# )
# n <- 1e6
# set.seed(123)
# dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
# bench_1e6 <- microbenchmark::microbenchmark(
# cutpointr::roc(dat, "x", "y", pos_class = 1, neg_class = 0,
# direction = ">="),
# rocr_roc(dat$x, dat$y),
# proc_roc(dat$x, dat$y),
# times = 30, unit = "ms"
# )
# n <- 1e7
# set.seed(123)
# dat <- data.frame(x = rnorm(n), y = sample(c(0:1), size = n, replace = TRUE))
# bench_1e7 <- microbenchmark::microbenchmark(
# cutpointr::roc(dat, "x", "y", pos_class = 1, neg_class = 0,
# direction = ">="),
# rocr_roc(dat$x, dat$y),
# proc_roc(dat$x, dat$y),
# times = 30, unit = "ms"
# )
#
# results_roc <- rbind(
# data.frame(time = summary(bench_100)$median,
# Solution = summary(bench_100)$expr,
# n = 100),
# data.frame(time = summary(bench_1000)$median,
# Solution = summary(bench_1000)$expr,
# n = 1000),
# data.frame(time = summary(bench_10000)$median,
# Solution = summary(bench_10000)$expr,
# n = 10000),
# data.frame(time = summary(bench_1e5)$median,
# Solution = summary(bench_1e5)$expr,
# n = 1e5),
# data.frame(time = summary(bench_1e6)$median,
# Solution = summary(bench_1e6)$expr,
# n = 1e6),
# data.frame(time = summary(bench_1e7)$median,
# Solution = summary(bench_1e7)$expr,
# n = 1e7)
# )
# results_roc$Solution <- as.character(results_roc$Solution)
# results_roc$Solution[grep(pattern = "cutpointr", x = results_roc$Solution)] <- "cutpointr"
# results_roc$Solution[grep(pattern = "rocr", x = results_roc$Solution)] <- "ROCR"
# results_roc$Solution[grep(pattern = "proc", x = results_roc$Solution)] <- "pROC"
# results_roc$task <- "ROC curve calculation"
## ---- echo = FALSE------------------------------------------------------------
# Our results
results_roc <- structure(list(time = c(0.7973505, 1.732651, 0.447701, 0.859301,
2.0358515, 0.694802, 1.878151, 5.662151, 3.6580505, 11.099251,
42.8208515, 35.3293005, 159.8100505, 612.471901, 610.4337005,
2032.693551, 7806.3854515, 7081.897251), Solution = c("cutpointr",
"ROCR", "pROC", "cutpointr", "ROCR", "pROC", "cutpointr", "ROCR",
"pROC", "cutpointr", "ROCR", "pROC", "cutpointr", "ROCR", "pROC",
"cutpointr", "ROCR", "pROC"), n = c(100, 100, 100, 1000, 1000,
1000, 10000, 10000, 10000, 1e+05, 1e+05, 1e+05, 1e+06, 1e+06,
1e+06, 1e+07, 1e+07, 1e+07), task = c("ROC curve calculation",
"ROC curve calculation", "ROC curve calculation", "ROC curve calculation",
"ROC curve calculation", "ROC curve calculation", "ROC curve calculation",
"ROC curve calculation", "ROC curve calculation", "ROC curve calculation",
"ROC curve calculation", "ROC curve calculation", "ROC curve calculation",
"ROC curve calculation", "ROC curve calculation", "ROC curve calculation",
"ROC curve calculation", "ROC curve calculation")), row.names = c(NA,
-18L), class = "data.frame")
## ---- echo = FALSE------------------------------------------------------------
results_all <- dplyr::bind_rows(results, results_roc)
ggplot(results_all, aes(x = n, y = time, col = Solution, shape = Solution)) +
geom_point(size = 3) + geom_line() +
scale_y_log10(breaks = c(0.5, 1, 2, 3, 5, 10, 25, 100, 250, 1000, 5000, 1e4, 15000)) +
scale_x_log10(breaks = c(100, 1000, 1e4, 1e5, 1e6, 1e7)) +
ylab("Median Time (milliseconds, log scale)") + xlab("Sample Size (log scale)") +
theme_bw() +
theme(legend.position = "bottom",
legend.key.width = unit(0.8, "cm"),
panel.spacing = unit(1, "lines")) +
facet_grid(~task)
## ---- echo = FALSE------------------------------------------------------------
res_table <- tidyr::spread(results_all, Solution, time) %>%
arrange(task)
knitr::kable(res_table)
