Method functions for cutpoint estimation

The included methods for calculating cutpoints are:

• maximize_metric: Maximize the metric function
• minimize_metric: Minimize the metric function
• maximize_loess_metric: Maximize the metric function after LOESS smoothing
• minimize_loess_metric: Minimize the metric function after LOESS smoothing
• maximize_spline_metric: Maximize the metric function after spline smoothing
• minimize_spline_metric: Minimize the metric function after spline smoothing
• maximize_gam_metric: Maximize the metric function after smoothing via Generalized Additive Models
• minimize_gam_metric: Minimize the metric function after smoothing via Generalized Additive Models
• maximize_boot_metric: Bootstrap the optimal cutpoint when maximizing a metric
• minimize_boot_metric: Bootstrap the optimal cutpoint when minimizing a metric
• oc_manual: Specify the cutoff value manually
• oc_mean: Use the sample mean as the “optimal” cutpoint
• oc_median: Use the sample median as the “optimal” cutpoint
• oc_youden_kernel: Maximize the Youden-Index after kernel smoothing the distributions of the two classes
• oc_youden_normal: Maximize the Youden-Index parametrically assuming normally distributed data in both classes

Metric functions

The included metrics to be used with the minimization and maximization methods are:

• accuracy: Fraction correctly classified
• abs_d_sens_spec: The absolute difference of sensitivity and specificity
• abs_d_ppv_npv: The absolute difference between positive predictive value (PPV) and negative predictive value (NPV)
• roc01: Distance to the point (0,1) on ROC space
• cohens_kappa: Cohen’s Kappa
• sum_sens_spec: sensitivity + specificity
• sum_ppv_npv: The sum of positive predictive value (PPV) and negative predictive value (NPV)
• prod_sens_spec: sensitivity * specificity
• prod_ppv_npv: The product of positive predictive value (PPV) and negative predictive value (NPV)
• youden: Youden- or J-Index = sensitivity + specificity - 1
• odds_ratio: (Diagnostic) odds ratio
• risk_ratio: risk ratio (relative risk)
• p_chisquared: The p-value of a chi-squared test on the confusion matrix
• cost_misclassification: The sum of the misclassification cost of false positives and false negatives. Additional arguments: cost_fp, cost_fn
• total_utility: The total utility of true / false positives / negatives. Additional arguments: utility_tp, utility_tn, cost_fp, cost_fn
• F1_score: The F1-score (2 * TP) / (2 * TP + FP + FN)
• metric_constrain: Maximize a selected metric given a minimal value of another selected metric
• sens_constrain: Maximize sensitivity given a minimal value of specificity
• spec_constrain: Maximize specificity given a minimal value of sensitivity
• acc_constrain: Maximize accuracy given a minimal value of sensitivity

Furthermore, the following functions are included which can be used as metric functions but are more useful for plotting purposes, for example in plot_cutpointr, or for defining new metric functions: tp, fp, tn, fn, tpr, fpr, tnr, fnr, false_omission_rate, false_discovery_rate, ppv, npv, precision, recall, sensitivity, and specificity.

The inputs to the arguments method and metric are functions so that user-defined functions can easily be supplied instead of the built-in ones.

Separate subgroups and bootstrapping

Cutpoints can be separately estimated on subgroups that are defined by a third variable, gender in this case. Additionally, if boot_runs is larger zero, cutpointr will carry out the usual cutpoint calculation on the full sample, just as before, and additionally on boot_runs bootstrap samples. This offers a way of gauging the out-of-sample performance of the cutpoint estimation method. If a subgroup is given, the bootstrapping is carried out separately for every subgroup which is also reflected in the plots and output.

set.seed(12)
opt_cut_b <- cutpointr(suicide, dsi, suicide, boot_runs = 1000)

opt_cut_b

## # A tibble: 1 x 16
##   direction optimal_cutpoint method          sum_sens_spec      acc sensitivity
##   <chr>                <dbl> <chr>                   <dbl>    <dbl>       <dbl>
## 1 >=                       2 maximize_metric       1.75179 0.864662    0.888889
##   specificity      AUC pos_class neg_class prevalence outcome predictor
##         <dbl>    <dbl> <fct>     <fct>          <dbl> <chr>   <chr>    
## 1    0.862903 0.923779 yes       no         0.0676692 suicide dsi      
##   data               roc_curve                 boot                 
##   <list>             <list>                    <list>               
## 1 <tibble [532 x 2]> <roc_cutpointr [13 x 10]> <tibble [1,000 x 23]>

The returned object has the additional column boot which is a nested tibble that includes the cutpoints per bootstrap sample along with the metric calculated using the function in metric and various default metrics. The metrics are suffixed by _b to indicate in-bag results or _oob to indicate out-of-bag results:

opt_cut$boot

## [1] NA

The summary and plots include additional elements that summarize or display the bootstrap results:

summary(opt_cut)

## Method: maximize_metric 
## Predictor: dsi 
## Outcome: suicide 
## Direction: >= 
## 
##     AUC   n n_pos n_neg
##  0.9238 532    36   496
## 
##  optimal_cutpoint youden    acc sensitivity specificity tp fn fp  tn
##                 2 0.7518 0.8647      0.8889      0.8629 32  4 68 428
## 
## Predictor summary: 
##     Data Min.   5% 1st Qu. Median      Mean 3rd Qu.  95% Max.       SD NAs
##  Overall    0 0.00       0      0 0.9210526       1 5.00   11 1.852714   0
##       no    0 0.00       0      0 0.6330645       0 4.00   10 1.412225   0
##      yes    0 0.75       4      5 4.8888889       6 9.25   11 2.549821   0

plot(opt_cut)

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

Bootstrapped cutpoints

It has been shown that bagging can substantially improve performance of a wide range of types of models in regression as well as in classification tasks. This method is available for cutpoint estimation via the maximize_boot_metric and minimize_boot_metric functions. If one of these functions is used as method, boot_cut bootstrap samples are drawn, the cutpoint optimization is carried out in each one and a summary (e.g. the mean) of the resulting optimal cutpoints on the bootstrap samples is returned as the optimal cutpoint in cutpointr. Note that if bootstrap validation is run, i.e. if boot_runs is larger zero, an outer bootstrap will be executed. In the bootstrap validation routine boot_runs bootstrap samples are generated and each one is again bootstrapped boot_cut times. This may lead to long run times, so activating the built-in parallelization may be advisable.

The advantages of bootstrapping the optimal cutpoint are that the procedure doesn’t possess parameters that have to be tuned, unlike the LOESS smoothing, that it doesn’t rely on assumptions, unlike the Normal method, and that it is applicable to any metric that can be used with minimize_metric or maximize_metric, unlike the Kernel method. Furthermore, like Random Forests cannot be overfit by increasing the number of trees, the bootstrapped cutpoints cannot be overfit by running an excessive amount of boot_cut repetitions.

set.seed(100)
cutpointr(suicide, dsi, suicide, gender, 
          method = maximize_boot_metric,
          boot_cut = 200, summary_func = mean,
          metric = accuracy, silent = TRUE)

## # A tibble: 2 x 18
##   subgroup direction optimal_cutpoint method               accuracy      acc
##   <chr>    <chr>                <dbl> <chr>                   <dbl>    <dbl>
## 1 female   >=                 5.73246 maximize_boot_metric 0.956633 0.956633
## 2 male     >=                 8.41026 maximize_boot_metric 0.95     0.95    
##   sensitivity specificity      AUC pos_class neg_class prevalence outcome
##         <dbl>       <dbl>    <dbl> <fct>     <fct>          <dbl> <chr>  
## 1    0.444444    0.994521 0.944647 yes       no         0.0688776 suicide
## 2    0.222222    1        0.861747 yes       no         0.0642857 suicide
##   predictor grouping data               roc_curve                boot 
##   <chr>     <chr>    <list>             <list>                   <lgl>
## 1 dsi       gender   <tibble [392 x 2]> <roc_cutpointr [11 x 9]> NA   
## 2 dsi       gender   <tibble [140 x 2]> <roc_cutpointr [11 x 9]> NA

LOESS smoothing for selecting a cutpoint

When using maximize_metric and minimize_metric the optimal cutpoint is selected by searching the maximum or minimum of the metric function. For example, we may want to minimize the misclassification cost. Since false negatives (a suicide attempt was not anticipated) can be regarded as much more severe than false positives we can set the cost of a false negative cost_fn for example to ten times the cost of a false positive.

opt_cut <- cutpointr(suicide, dsi, suicide, gender, method = minimize_metric,
                     metric = misclassification_cost, cost_fp = 1, cost_fn = 10)

## Assuming the positive class is yes

## Assuming the positive class has higher x values

plot_metric(opt_cut)

As this “optimal” cutpoint may depend on minor differences between the possible cutoffs, smoothing of the function of metric values by cutpoint value might be desirable, especially in small samples. The minimize_loess_metric and maximize_loess_metric functions can be used to smooth the function so that the optimal cutpoint is selected based on the smoothed metric values. Options to modify the smoothing, which is implemented using loess.as from the fANCOVA package, include:

• criterion: the criterion for automatic smoothing parameter selection: “aicc” denotes bias-corrected AIC criterion, “gcv” denotes generalized cross-validation.
• degree: the degree of the local polynomials to be used. It can be 0, 1 or 2.
• family: if “gaussian” fitting is by least-squares, and if “symmetric” a re-descending M estimator is used with Tukey’s biweight function.
• user.span: the user-defined parameter which controls the degree of smoothing.

Using parameters for the LOESS smoothing of criterion = "aicc", degree = 2, family = "symmetric", and user.span = 0.7 we get the following smoothed versions of the above metrics:

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)

The optimal cutpoint for the female subgroup changes to 3. Note, though, that there are no reliable rules for selecting the “best” smoothing parameters. Notably, the LOESS smoothing is sensitive to the number of unique cutpoints. A large number of unique cutpoints generally leads to a more volatile curve of metric values by cutpoint value, even after smoothing. Thus, the curve tends to be undersmoothed in that scenario. The unsmoothed metric values are returned in opt_cut$roc_curve in the column m_unsmoothed.

Smoothing via Generalized Additive Models for selecting a cutpoint

In a similar fashion, the function of metric values per cutpoint can be smoothed using Generalized Additive Models with smooth terms. Internally, mgcv::gam carries out the smoothing which can be customized via the arguments formula and optimizer, see help("gam", package = "mgcv"). Most importantly, the GAM can be specified by altering the default formula, for example the smoothing function could be configured to apply cubic regression splines ("cr") as the smooth term. As the suicide data has only very few unique cutpoints, it is not very suitable for showcasing the GAM smoothing, so we will use two classes of the iris data here. In this case, the purely empirical method and the GAM smoothing lead to identical cutpoints, but in practice the GAM smoothing tends to be more robust, especially with larger data. An attractive feature of the GAM smoothing is that the default values tend to work quite well and usually require no tuning, eliminating researcher degrees of freedom.

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)

## Assuming the positive class is virginica

## Assuming the positive class has higher x values

plot_metric(opt_cut)

Spline smoothing for selecting a cutpoint

Again in the same fashion the function of metric values per cutpoint can be smoothed using smoothing splines. By default, the number of knots is automatically chosen using the cutpoint_knots function. That function uses stats::.nknots.smspl, which is the default in stats::smooth.spline to pick the number of knots.

Alternatively, the number of knots can be set manually and also the other smoothing parameters of stats::smooth.spline can be set as desired. For details see ?maximize_spline_metric.

opt_cut <- cutpointr(suicide, dsi, suicide, gender, 
                     method = minimize_spline_metric, spar = 0.4,
                     metric = misclassification_cost, cost_fp = 1, cost_fn = 10)

## Assuming the positive class is yes

## Assuming the positive class has higher x values

## nknots: 10
## nknots: 10

plot_metric(opt_cut)

cutpointr(suicide, dsi, suicide, gender, method = oc_youden_normal)

## Assuming the positive class is yes

## Assuming the positive class has higher x values

## # A tibble: 2 x 18
##   subgroup direction optimal_cutpoint method           sum_sens_spec      acc
##   <chr>    <chr>                <dbl> <chr>                    <dbl>    <dbl>
## 1 female   >=                 2.47775 oc_youden_normal       1.71618 0.895408
## 2 male     >=                 3.17226 oc_youden_normal       1.54453 0.864286
##   sensitivity specificity      AUC pos_class neg_class prevalence outcome
##         <dbl>       <dbl>    <dbl> <fct>     <fct>          <dbl> <chr>  
## 1    0.814815    0.901370 0.944647 yes       no         0.0688776 suicide
## 2    0.666667    0.877863 0.861747 yes       no         0.0642857 suicide
##   predictor grouping data               roc_curve                boot 
##   <chr>     <chr>    <list>             <list>                   <lgl>
## 1 dsi       gender   <tibble [392 x 2]> <roc_cutpointr [11 x 9]> NA   
## 2 dsi       gender   <tibble [140 x 2]> <roc_cutpointr [11 x 9]> NA

cutpointr(suicide, dsi, suicide, gender, method = oc_youden_kernel)

## Assuming the positive class is yes

## Assuming the positive class has higher x values

## # A tibble: 2 x 18
##   subgroup direction optimal_cutpoint method           sum_sens_spec      acc
##   <chr>    <chr>                <dbl> <chr>                    <dbl>    <dbl>
## 1 female   >=                 1.18128 oc_youden_kernel       1.80812 0.885204
## 2 male     >=                 1.31636 oc_youden_kernel       1.58694 0.807143
##   sensitivity specificity      AUC pos_class neg_class prevalence outcome
##         <dbl>       <dbl>    <dbl> <fct>     <fct>          <dbl> <chr>  
## 1    0.925926    0.882192 0.944647 yes       no         0.0688776 suicide
## 2    0.777778    0.809160 0.861747 yes       no         0.0642857 suicide
##   predictor grouping data               roc_curve                boot 
##   <chr>     <chr>    <list>             <list>                   <lgl>
## 1 dsi       gender   <tibble [392 x 2]> <roc_cutpointr [11 x 9]> NA   
## 2 dsi       gender   <tibble [140 x 2]> <roc_cutpointr [11 x 9]> NA

Calculating only the ROC curve

When running cutpointr, a ROC curve is by default returned in the column roc_curve. This ROC curve can be plotted using plot_roc. Alternatively, if only the ROC curve is desired and no cutpoint needs to be calculated, the ROC curve can be created using roc() and plotted using plot_cutpointr. The roc function, unlike cutpointr, does not determine direction, pos_class or neg_class automatically.

roc_curve <- roc(data = suicide, x = dsi, class = suicide,
    pos_class = "yes", neg_class = "no", direction = ">=")
auc(roc_curve)

## [1] 0.9237791

head(roc_curve)

## # A tibble: 6 x 9
##   x.sorted    tp    fp    tn    fn    tpr   tnr     fpr   fnr
##      <dbl> <dbl> <dbl> <int> <int>  <dbl> <dbl>   <dbl> <dbl>
## 1      Inf     0     0   496    36 0      1     0       1    
## 2       11     1     0   496    35 0.0278 1     0       0.972
## 3       10     2     1   495    34 0.0556 0.998 0.00202 0.944
## 4        9     3     1   495    33 0.0833 0.998 0.00202 0.917
## 5        8     4     1   495    32 0.111  0.998 0.00202 0.889
## 6        7     7     1   495    29 0.194  0.998 0.00202 0.806

plot_roc(roc_curve)

Midpoints

So far - which is the default in cutpointr - we have considered all unique values of the predictor as possible cutpoints. An alternative could be to use a sequence of equidistant values instead, for example in the case of the suicide data all integers in \([0, 10]\). However, with very sparse data and small intervals between the candidate cutpoints (i.e. a ‘dense’ sequence like seq(0, 10, by = 0.01)) this leads to the uninformative evaluation of large ranges of cutpoints that all result in the same metric value. A more elegant alternative, not only for the case of sparse data, that is supported by cutpointr is the use of a mean value of the optimal cutpoint and the next highest (if direction = ">=") or the next lowest (if direction = "<=") predictor value in the data. The result is an optimal cutpoint that is equal to the cutpoint that would be obtained using an infinitely dense sequence of candidate cutpoints and is thus usually more efficient computationally. This behavior can be activated by setting use_midpoints = TRUE, which is the default. If we use this setting, we obtain an optimal cutpoint of 1.5 for the complete sample on the suicide data instead of 2 when maximizing the sum of sensitivity and specificity.

Assume the following small data set:

dat <- data.frame(outcome = c("neg", "neg", "neg", "pos", "pos", "pos", "pos"),
                  pred    = c(1, 2, 3, 8, 11, 11, 12))

Since the distance of the optimal cutpoint (8) to the next lowest observation (3) is rather large we arrive at a range of possible cutpoints that all maximize the metric. In the case of this kind of sparseness it might for example be desirable to classify a new observation with a predictor value of 4 as belonging to the negative class. If use_midpoints is set to TRUE, the mean of the optimal cutpoint and the next lowest observation is returned as the optimal cutpoint, if direction is >=. The mean of the optimal cutpoint and the next highest observation is returned as the optimal cutpoint, if direction = "<=".

opt_cut <- cutpointr(dat, x = pred, class = outcome, use_midpoints = TRUE)

## Assuming the positive class is pos

## Assuming the positive class has higher x values

plot_x(opt_cut)

A simulation demonstrates more clearly that setting use_midpoints = TRUE avoids biasing the cutpoints. To simulate the bias of the metric functions, the predictor values of both classes were drawn from normal distributions with constant standard deviations of 10, a constant mean of the negative class of 100 and higher mean values of the positive class that are selected in such a way that optimal Youden-Index values of 0.2, 0.4, 0.6, and 0.8 result in the population. Samples of 9 different sizes were drawn and the cutpoints that maximize the Youden-Index were estimated. The simulation was repeated 10000 times. As can be seen by the mean error, use_midpoints = TRUE eliminates the bias that is introduced by otherwise selecting the value of an observation as the optimal cutpoint. If direction = ">=", as in this case, the observation that represents the optimal cutpoint is the highest possible cutpoint that leads to the optimal metric value and thus the biases are positive. The methods oc_youden_normal and oc_youden_kernel are always unbiased, as they don’t select a cutpoint based on the ROC-curve or the function of metric values per cutpoint.

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

Finding all cutpoints with acceptable performance

By default, most packages only return the “best” cutpoint and disregard other cutpoints with quite similar performance, even if the performance differences are minuscule. cutpointr makes this process more explicit via the tol_metric argument. For example, if all cutpoints are of interest that achieve at least an accuracy within 0.05 of the optimally achievable accuracy, tol_metric can be set to 0.05 and also those cutpoints will be returned.

In the case of the suicide data and when maximizing the sum of sensitivity and specificity, empirically the cutpoints 2 and 3 lead to quite similar performances. If tol_metric is set to 0.05, both will be returned.

opt_cut <- cutpointr(suicide, dsi, suicide, metric = sum_sens_spec, 
                     tol_metric = 0.05, break_ties = c)

## Assuming the positive class is yes

## Assuming the positive class has higher x values

## Multiple optimal cutpoints found, applying break_ties.

library(tidyr)
opt_cut %>% 
    select(optimal_cutpoint, sum_sens_spec) %>% 
    unnest(cols = c(optimal_cutpoint, sum_sens_spec))

## # A tibble: 2 x 2
##   optimal_cutpoint sum_sens_spec
##              <dbl>         <dbl>
## 1                2          1.75
## 2                1          1.70

Manual and mean / median cutpoints

Using the oc_manual function the optimal cutpoint will not be determined based on, for example, a metric but is instead set manually using the cutpoint argument. This is useful for supplying and evaluating cutpoints that were found in the literature or in other external sources.

The oc_manual function could also be used to set the cutpoint to the sample mean using cutpoint = mean(data$x). However, this may introduce bias into the bootstrap validation procedure, since the actual mean of the population is not known and thus the mean to be used as the cutpoint should be automatically determined in every resample. To do so, the oc_mean and oc_median functions can be used.

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)

Nonstandard evaluation via tidyeval

The arguments to cutpointr do not need to be enclosed in quotes. This is possible thanks to nonstandard evaluation of the arguments, which are evaluated on data.

Functions that use nonstandard evaluation are often not suitable for programming with. The use of nonstandard evaluation may lead to scoping problems and subsequent obvious as well as possibly subtle errors. cutpointr uses tidyeval internally and accordingly the same rules as for programming with dplyr apply. Arguments can be unquoted with !!:

myvar <- "dsi"
cutpointr(suicide, !!myvar, suicide)

ROC curve and optimal cutpoint for multiple variables

Alternatively, we can map the standard evaluation version cutpointr to the column names. If direction and / or pos_class and neg_class are unspecified, these parameters will automatically be determined by cutpointr so that the AUC values for all variables will be \(> 0.5\).

We could do this manually, e.g. using purrr::map, but to make this task more convenient multi_cutpointr can be used to achieve the same result. It maps multiple predictor columns to cutpointr, by default all numeric columns except for the class column.

mcp <- multi_cutpointr(suicide, class = suicide, pos_class = "yes", 
                use_midpoints = TRUE, silent = TRUE) 
summary(mcp)

## Method: maximize_metric 
## Predictor: age, dsi 
## Outcome: suicide 
## 
## Predictor: age 
## -------------------------------------------------------------------------------- 
##  direction    AUC   n n_pos n_neg
##         <= 0.5257 532    36   496
## 
##  optimal_cutpoint sum_sens_spec    acc sensitivity specificity tp fn  fp tn
##              55.5        1.1154 0.1992      0.9722      0.1431 35  1 425 71
## 
## Predictor summary: 
##     Data Min. 5% 1st Qu. Median    Mean 3rd Qu.   95% Max.      SD NAs
##  Overall   18 19      24   28.0 34.1259   41.25 65.00   83 15.0542   0
##       no   18 19      24   28.0 34.2218   41.25 65.50   83 15.1857   0
##      yes   18 18      22   27.5 32.8056   41.25 54.25   69 13.2273   0
## 
## Predictor: dsi 
## -------------------------------------------------------------------------------- 
##  direction    AUC   n n_pos n_neg
##         >= 0.9238 532    36   496
## 
##  optimal_cutpoint sum_sens_spec    acc sensitivity specificity tp fn fp  tn
##               1.5        1.7518 0.8647      0.8889      0.8629 32  4 68 428
## 
## Predictor summary: 
##     Data Min.   5% 1st Qu. Median   Mean 3rd Qu.  95% Max.     SD NAs
##  Overall    0 0.00       0      0 0.9211       1 5.00   11 1.8527   0
##       no    0 0.00       0      0 0.6331       0 4.00   10 1.4122   0
##      yes    0 0.75       4      5 4.8889       6 9.25   11 2.5498   0

Accessing data, roc_curve, and boot

The object returned by cutpointr is of the classes cutpointr, tbl_df, tbl, and data.frame. Thus, it can be handled like a usual data frame. The columns data, roc_curve, and boot consist of nested data frames, which means that these are list columns whose elements are data frames. They can either be accessed using [ or by using functions from the tidyverse. If subgroups were given, the output contains one row per subgroup and the function that accesses the data should be mapped to every row or the data should be grouped by subgroup.

set.seed(123)
opt_cut_b_g <- cutpointr(suicide, dsi, suicide, gender, boot_runs = 1000)

# 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))

## # A tibble: 2 x 4
##   subgroup sd_oc_boot m_oc_boot m_acc_oob
##   <chr>         <dbl>     <dbl>     <dbl>
## 1 female        0.766      2.17     0.880
## 2 male          1.51       2.92     0.806

Adding metrics to the result of cutpointr() or roc()

By default, the output of cutpointr includes the optimized metric and several other metrics. The add_metric function adds further metrics. Here, we’re adding the negative predictive value (NPV) and the positive predictive value (PPV) at the optimal cutpoint per subgroup:

cutpointr(suicide, dsi, suicide, gender, metric = youden, silent = TRUE) %>% 
    add_metric(list(ppv, npv)) %>% 
    select(subgroup, optimal_cutpoint, youden, ppv, npv)

## # A tibble: 2 x 5
##   subgroup optimal_cutpoint   youden      ppv      npv
##   <chr>               <dbl>    <dbl>    <dbl>    <dbl>
## 1 female                  2 0.808118 0.367647 0.993827
## 2 male                    3 0.625106 0.259259 0.982301

In the same fashion, additional metric columns can be added to a roc_cutpointr object:

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()

## # A tibble: 6 x 7
##   x.sorted    tp    fp    tn    fn cohens_kappa F1_score
##      <dbl> <dbl> <dbl> <int> <int>        <dbl>    <dbl>
## 1      Inf     0     0   496    36       0        0     
## 2       11     1     0   496    35       0.0506   0.0541
## 3       10     2     1   495    34       0.0931   0.103 
## 4        9     3     1   495    33       0.138    0.15  
## 5        8     4     1   495    32       0.182    0.195 
## 6        7     7     1   495    29       0.301    0.318

User-defined functions

mean_cut <- function(data, x, ...) {
    oc <- mean(data[[x]])
    return(data.frame(optimal_cutpoint = oc))
}

misclassification_cost

## function (tp, fp, tn, fn, cost_fp = 1, cost_fn = 1, ...) 
## {
##     misclassification_cost <- cost_fp * fp + cost_fn * fn
##     misclassification_cost <- matrix(misclassification_cost, 
##         ncol = 1)
##     colnames(misclassification_cost) <- "misclassification_cost"
##     return(misclassification_cost)
## }
## <bytecode: 0x00000000247b08c0>
## <environment: namespace:cutpointr>

Flexible plotting function

Using plot_cutpointr any metric can be chosen to be plotted on the x- or y-axis and results of cutpointr() as well as roc() can be plotted. If a cutpointr object is to be plotted, it is thus irrelevant which metric function was chosen for cutpoint estimation. Any metric that can be calculated based on the ROC curve can be subsequently plotted as only the true / false positives / negatives over all cutpoints are needed. That way, not only the above plots can be produced, but also any combination of two metrics (or metric functions) and / or cutpoints. The built-in metric functions as well as user-defined functions or anonymous functions can be supplied to xvar and yvar. If bootstrapping was run, confidence intervals can be plotted around the y-variable. This is especially useful if the cutpoints, available in the cutpoints function, are placed on the x-axis. Note that confidence intervals can only be correctly plotted if the values of xvar are constant across bootstrap samples. For example, confidence intervals for TPR by FPR (a ROC curve) cannot be plotted easily, as the values of the false positive rate vary per bootstrap sample.

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()

Manual plotting

Since cutpointr returns a data.frame with the original data, bootstrap results, and the ROC curve in nested tibbles, these data can be conveniently extracted and plotted manually. The relevant nested tibbles are in the columns data, roc_curve and boot. The following is an example of accessing and plotting the grouped data.

opt_cut_b_g %>% 
    select(data, subgroup) %>% 
    unnest(cols = data) %>% 
    ggplot(aes(x = suicide, y = dsi)) + 
    geom_boxplot(alpha = 0.3) + facet_grid(~subgroup)