Add Best PC and update vignette (#4)

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: stats4phc
 Title: Performance Evaluation for the Prognostic Value of Predictive Models Intended to Support 
   Personalized Healthcare Through Predictiveness Curves and Positive / Negative Predictive Values
-Version: 0.1
+Version: 0.1.1
 Authors@R: c(
     person("Ondrej", "Slama", email = "[email protected]", role = c("aut", "cre")),
     person("Darrick", "Shen", email = "[email protected]", role = "aut"),

R/PV.R

@@ -116,12 +116,13 @@ predictionColours <- function(x, show.best) {
     "NPV" = "#53B400",
     "Best PPV" = "gray65",
     "PPV" = "red",
+    "Best PC" = "black",
     "PC" = "plum",
     "1-NPV" = "royalblue2",
     "Best 1-NPV" = "gray50"
   )
   if (show.best) {
-    x <- c(x, paste("Best", setdiff(x, "PC")))
+    x <- c(x, paste("Best", x))
   }
   return(clrs[names(clrs) %in% x]) # keep the ordering above
 }

R/riskProfile.R

@@ -386,12 +386,13 @@ riskProfile <- function(outcome,
       distinct() %>%
       arrange(.data$percentile, .data$score) %>%
       mutate(
+        PC = ifelse(.data$percentile <= 1 - prev, 0, 1),
         PPV = bestPPV(perc = .data$percentile, prev = prev),
         `1-NPV` = bestMNPV(perc = .data$percentile, prev = prev),
         NPV = 1 - .data$`1-NPV`
       ) %>%
       tidyr::pivot_longer(
-        cols = c("PPV", "1-NPV", "NPV"), names_to = "pv", values_to = "pvValue"
+        cols = c("PC", "PPV", "1-NPV", "NPV"), names_to = "pv", values_to = "pvValue"
       ) %>%
       filter(.data$pv %in% include) %>%
       mutate(
@@ -404,10 +405,10 @@ riskProfile <- function(outcome,
         geom_line(
           data = best,
           aes(x = .data[[xvar]], y = .data$pvValue, linewidth = .data$pv),
-          colour = "gray70"
+          colour = "gray60"
         ) +
         scale_linewidth_manual(
-          values = c("Best 1-NPV" = 0.3, "Best PPV" = 0.3, "Best NPV" = 0.3),
+          values = c("Best 1-NPV" = 0.3, "Best PPV" = 0.3, "Best PC" = 0.3, "Best NPV" = 0.3),
           name = "Best PVs"
         )
     } else {

README.md

@@ -30,7 +30,7 @@ remotes::install_github(repo = "genentech/stats4phc")
 For reproducibility, refer to a specific version tag, for example
 
 ``` r
-remotes::install_github(repo = "genentech/stats4phc", ref = "v0.1")
+remotes::install_github(repo = "genentech/stats4phc", ref = "v0.1.1")
 ```
 
 

vignettes/stats4phc.Rmd

@@ -21,77 +21,105 @@ This package provides functions for performance evaluation for the prognostic va
 -   Calibration
 -   Sensitivity and specificity
 
-To begin with, let's align on terminology.
+The intention is not to replace standard metrics evaluation (like Brier score, log-loss, or AUC). 
+On the contrary, the mentioned quantities should be checked in addition in order to get the overall sense of model behaviour.
+
 
 
 ## Terminology
 
+To begin with, let's align on terminology.
 Below we define the terms that will be used across the article:
 
 -   outcome: the true observation of the quantity of interest;
 
--   score:
+-   score / risk:
 
     -   either a raw value (e.g. a biomarker) for the purpose of measuring (or approximating) the outcome,
 
-    -   or a prediction score given by a predictive model, where the outcome was modeled as the response;
+    -   or a probability given by a predictive model, where the outcome was modeled as the response;
 
 -   estimate: output of a statistical methodology, where score is used as independent variable and outcome as a dependent variable.
 
 
+Let's now load the package and the example data.
+
+```{r}
+library(stats4phc)
+
+# Read in example data
+auroc <- read.csv(system.file("extdata", "sample.csv", package = "stats4phc"))
+rscore <- auroc$predicted_calibrated  # vector of already calibrated model probabilities
+truth <- as.numeric(auroc$actual)     # vector of 0s or 1s
+```
+
+
 # Predictiveness curves
 
 Predictiveness curves are an insightful visualization to assess the inherent ability of prognostic models to provide predictions to individual patients. Cumulative versions of predictiveness curves represent positive predictive values (PPV) and 1 - negative predictive values (1 - NPV) and are also informative if the eventual goal is to use a cutoff for clinical decision making.
 
 You can use `riskProfile` function to visualize and assess all these quantities.
 
-Here is an example:
+Note that method "asis" (below on the graphs) means that the score (or model probabilities in our case) are taken as is, 
+i.e. there is no estimation or smoothing. 
 
-```{r}
-library(stats4phc)
 
-# Read in example data
-auroc <- read.csv(system.file("extdata", "sample.csv", package = "stats4phc"))
-rscore <- auroc$predicted_calibrated
-truth <- as.numeric(auroc$actual)
-```
+## Predictiveness curve
+
+Let's start with predictiveness curve:
 
 ```{r fig.height=5, fig.width=8}
-p <- riskProfile(outcome = truth, score = rscore)
+p <- riskProfile(outcome = truth, score = rscore, include = "PC")
 p$plot
 ```
 
-```{r}
-head(p$data)
+Ideally, all subjects in the population that have the condition (=> prevalence) are marked as having the condition (predicted risk = 1) and all subjects without the condition (=> 1 - prevalence) are marked as not having the condition (predicted risk = 0).
+This implies that the ideal predictiveness curve is 0 for all subjects not having the condition, 
+and then it steps (jumps) at `1 - prevalence` to 1 for all the subjects having the condition (see the gray line). 
+
+In reality, the curves are not step functions. The more flat the curves get, the less discrimination, and therefore utility, there is in the model.
+
+One can also investigate the tails of predictiveness curve. The model is more useful if these regions have very low or very high predicted risks relatively to the rest of the data.
+
+
+## Positive / Negative predictive values
+
+Now let's plot PPV and 1-NPV:
+
+```{r fig.height=5, fig.width=8}
+p <- riskProfile(outcome = truth, score = rscore, include = c("PPV", "1-NPV"))
+p$plot
 ```
 
-There is an extensive documentation of this function with examples if you run `?riskProfile`.
+Again, in an ideal case, they both are as close to the gray lines as possible.
 
-To briefly highlight the functionalities:
+In an ideal scenario:
 
--   Use the `methods` argument to specify the desired estimation method (see the last section) or use `"asis"` for no estimation.
+- in terms of PPV: If all the subjects with the condition are predicted perfectly, then PPV = TP / PP = 1 (TP = true positive, PP = predicted positive). 
+Hence, all the subjects with the condition must be higher than `1 - prevalence` on the prediction percentile for PPV = 1. 
 
--   You can adjust the prevalence by setting `prev.adj` to the desired amount.
+- in terms of 1-NPV: If all the subjects without the condition are predicted perfectly, then NPV = TN / PN = 1 (TN = true negative, PN = predicted negative). 
+Hence, all the subjects without the condition must be lower than `1 - prevalence` on the risk percentile for 1-NPV = 0.
 
--   `show.nonparam.pv` controls whether to show/hide the non-parametric estimation of PPV, 1-NPV, and NPV.
 
--   `show.best.pv` controls whether to show/hide the theoretically best possible PPV, 1-NPV, NPV.
+## Output settings
 
--   Use `include` argument to specify what quantities to show:
+Note that:
 
-    -   PC = predictiveness curve
+- most importantly, in case of a biomarker or if the model probabilities are not calibrated well, you can use a smoother, see `methods` argument and the last section of the vignette.
 
-    -   PPV = positive predictive value
+- the prevalence can be adjusted by setting it in `prev.adj`. 
 
-    -   NPV = negative predictive value
+- you can also plot "NPV" by adjusting the `include` parameter.
 
-    -   1-NPV = 1 - negative predictive value
+- you can also access the underlying data:
 
--   `plot.raw` sets whether to plot raw values or percentiles.
+```{r}
+p <- riskProfile(outcome = truth, score = rscore, include = c("PPV", "1-NPV"))
+head(p$data)
+```
 
--   `rev.order` sets whether to reverse the order of the score (useful if higher score refers to lower outcome).
 
--   The output is the plot itself and the underlying data.
 
 # Calibration
 
@@ -104,37 +132,41 @@ p <- calibrationProfile(outcome = truth, score = rscore)
 p$plot
 ```
 
-```{r}
-head(p$data)
-```
+In an ideal scenario, the fitted curves should be identical with the identity line. 
+
+In reality, the closer they are to the identity line, the better. 
+
+Note that you can also quantify calibration through discrimination and miscalibration index, 
+see this [blog post](https://stats4datascience.com/posts/three_metrics_binary/) 
+and [modsculpt](https://github.com/Genentech/modsculpt)
+R package (metrics functions).
+
 
-There is an extensive documentation of this function with examples if you run `?calibrationProfile`.
+## Output settings
 
-To briefly highlight the functionalities:
+Note that:
 
--   Use the `methods` argument to specify the desired estimation method (see the last section). In this case, `"asis"` is not allowed.
+- in case of a biomarker or if the model probabilities are not calibrated well, you can use a smoother, see `methods` argument and the last section of the vignette. In this case, `"asis"` is not allowed.
 
--   Use `include` argument to specify what additional quantities to show:
+- use `include` argument to specify what additional quantities to show:
 
-    -   `"loess"`: Adds non-parametric Loess fit.
+  - `"loess"`: Adds non-parametric Loess fit.
 
-    -   `"citl"`: Adds "Calibration in the Large", an overall mean of outcome and score.
+  - `"citl"`: Adds "Calibration in the Large", an overall mean of outcome and score.
 
-    -   `"rug"`: Adds "rug", i.e. ticks on x-axis showing the individual data points (top axis shows score for outcome == 1, bottom axis shows score for outcome == 0).
+  - `"rug"`: Adds "rug", i.e. ticks on x-axis showing the individual data points (top axis shows score for outcome == 1, bottom axis shows score for outcome == 0).
 
-    -   `"datapoints"`: Similar to rug, just shows jittered points instead of ticks.
+  - `"datapoints"`: Similar to rug, just shows jittered points instead of ticks.
 
--   `plot.raw` sets whether to plot raw values or percentiles.
+- use `margin.type` to add a marginal plot through `ggExtra::ggMarginal`. You can select one of `c("density", "histogram", "boxplot", "violin", "densigram")`. It adds the selected 1d graph on top of the calibrtion plot and is suitable for investigating the score.
 
--   `rev.order` sets whether to reverse the order of the score (useful if higher score refers to lower outcome).
+- you can again access the underlying data with `p$data`.
 
--   Use `margin.type` to add a marginal plot through `ggExtra::ggMarginal`. You can select one of `c("density", "histogram", "boxplot", "violin", "densigram")`. It adds the selected 1d graph on top of the calibrtion plot and is suitable for investigating the score.
 
--   The output is the plot itself and the underlying data.
 
 # Sensitivity and specificity
 
-Ultimately, we provide a sensitivity and specificity plot. For these quantities you need to define a cutoff with which you can trasnform the numeric score to binary. We use data-driven cutoffs, meaning that every single value of score is taken as the cutoff, allowing us to visualize the sensitivity and specificity as a function of score. This graph may inform you of the best suitable cutoff for your model, although we usually recommend to output the whole score range, not just the binary decisions.
+Ultimately, we provide a sensitivity and specificity plot as a function of score (threshold is data-driven). This graph may inform you of the best suitable cutoff for your model, although we usually recommend to output the whole score range, not just the binary decisions.
 
 You can use `sensSpec` function to visualize and assess sensitivity and specificity.
 
@@ -143,21 +175,19 @@ p <- sensSpec(outcome = truth, score = rscore)
 p$plot
 ```
 
-```{r}
-head(p$data)
-```
+Again, the ideal scenario would be having a model following the gray lines.
+Since there is a trade-off between sensitivity and specificity, the graph may guide you which threshold (or thresholds) to choose, depending if one is more important than the other.
 
-There is an extensive documentation of this function with examples if you run `?sensSpec`.
 
-To briefly highlight the functionalities:
+## Output settings
 
--   Use the `methods` argument to specify the desired estimation method (see the last section) or use `"asis"` for no estimation.
+Note that:
 
--   `show.best` controls whether to show/hide the theoretically best possible sensitivity and specificity.
+- in case of a biomarker or if the model probabilities are not calibrated well, you can use a smoother, see `methods` argument and the last section of the vignette. 
+
+- you can again access the underlying data with `p$data`.
 
--   `plot.raw` sets whether to plot raw values or percentiles.
 
--   `rev.order` sets whether to reverse the order of the score (useful if higher score refers to lower outcome).
 
 # Adjusting the graphs
 
@@ -188,6 +218,7 @@ p$plot +
 
 Otherwise, you can use the `$data` element to construct your own graph as well.
 
+
 # Estimations in stats4phc
 
 For all the plotting functions from this package, there is a possibility to define an estimation function, which will be applied on the given score. In `calibrationProfile`, this serves as a calibration curve. In `riskProfile`, this smooths the given score. All of this is always driven by the `methods` argument, which is available in each of the plotting functions.