R/nomogrammer.r

In sbalci/ClinicoPathJamoviModule: Analysis for Clinicopathological Research

#' @import ggplot2
#' @import scales
#' @importFrom stats line


## Create simple Fagan nomograms as ggplot objects
##   Based on Perl web-implementation (https://araw.mede.uic.edu/cgi-bin/testcalc.pl)
##   Authors: AM. Chekroud* & A. Schwartz (* adam dot chekroud at yale . edu)
##   December 2016

nomogrammer <- function(Prevalence,
                        Sens = NULL,
                        Spec = NULL,
                        Plr = NULL,
                        Nlr = NULL,
                        Detail = FALSE,
                        NullLine = FALSE,
                        LabelSize = (14/5),
                        Verbose = FALSE){

## Function inputs:
    # Prevalence (prior probability) as a number between 0 and 1
    # Either
        # Sens & Spec
        # model sensitivity and specificity as a number between 0 and 1
    # Or
        # Likelihood ratios
        # Positive and Negative LRs (numeric)

## Function options:
    # Detail: If true, will overlay key statistics onto the plot
    # NullLine: If true, will add a line from prior prob through LR = 1
    # LabelSize: Tweak this number to change the label sizes
    # Verbose: Print out relevant metrics in the console

## Function returns:
    # ggplot object



######################################
########## Libraries & Functions #####
######################################

## Libraries
# require(ggplot2)
# require(scales)






## Helper functions
##   (defined inside nomogrammer, so remain local only & wont clutter user env)
odds         <- function(p){
    # Function converts probability into odds
    o <- p/(1-p)
    return(o)
}

logodds      <- function(p){
    # Function returns logodds for a probability
    lo <- log10(p/(1-p))
    return(lo)
}

logodds_to_p <- function(lo){
    # Function goes from logodds back to a probability
    o <- 10^lo
    p <- o/(1+o)
    return(p)
}

p2percent <- function(p){
    # Function turns numeric probability into string percentage
    # e.g. 0.6346111 -> 63.5%
    scales::percent(signif(p, digits = 3))}


######################################
########## Calculations     ##########
######################################

## Checking inputs

## Prevalence
# needs to exist
if(missing(Prevalence)){
    stop("Prevalence is missing")
}
# needs to be numeric
if(!is.numeric(Prevalence)){stop("Prevalence should be numeric")}
# needs to be a prob not a percent
if((Prevalence > 1) | (Prevalence <= 0)){stop("Prevalence should be a probability (did you give a %?)")}

# Did user give sens & spec?
if(missing(Sens) | missing(Spec)){
    sensspec <- FALSE
} else{ sensspec <- TRUE}
# if yes, make sure they are numbers
if(sensspec == TRUE){
    if(!is.numeric(Sens)){stop("Sensitivity should be numeric")}
    if(!is.numeric(Spec)){stop("Specificity should be numeric")}
    # numbers that are probabilities not percentages
    if((Sens > 1) | (Sens <= 0)){stop("Sensitivity should be a probability (did you give a %?)")}
    if((Spec > 1) | (Spec <= 0)){stop("Specificity should be a probability (did you give a %?)")}
}


# Did user give PLR & NLR?
if(missing(Plr) | missing(Nlr)){
    plrnlr <- FALSE
} else{plrnlr <- TRUE}
# if yes, make sure they are numbers
if(plrnlr == TRUE){
    if(!is.numeric(Plr)){stop("PLR should be numeric")}
    if(!is.numeric(Nlr)){stop("NLR should be numeric")}
    # numbers that vaguely make sense
    if(Plr < 1){stop("PLR shouldn't be less than 1")}
    if(Nlr < 0){stop("NLR shouldn't be below zero")}
    if(Nlr > 1){stop("NLR shouldn't be more than 1")}
}

# Did they give a valid sensspec and plrnlr? If yes, ignore the LRs and tell them
if((sensspec == TRUE) && (plrnlr == TRUE) ){
    warning("You provided sens/spec as well as likelihood ratios-- I ignored the LRs!")
}


## If sens/spec provided, we calculate posterior probabilities & odds using sens & spec
##  otherwise, if plr and nlr provided, we calculate posteriors using them
##  if neither exist, then return an error
if(sensspec == TRUE){
    prior_prob  <- Prevalence
    prior_odds  <- odds(prior_prob)
    sensitivity <- Sens
    specificity <- Spec
    PLR <- sensitivity/(1-specificity)
    NLR <- (1-sensitivity)/specificity
    post_odds_pos  <- prior_odds * PLR
    post_odds_neg  <- prior_odds * NLR
    post_prob_pos  <- post_odds_pos/(1+post_odds_pos)
    post_prob_neg  <- post_odds_neg/(1+post_odds_neg)
} else if(plrnlr == TRUE){
    prior_prob  <- Prevalence
    prior_odds  <- odds(prior_prob)
    PLR <- Plr
    NLR <- Nlr
    sensitivity <- (PLR*(1-NLR))/(PLR-NLR)    ## TODO: check Adam's math!
    specificity <- (1-PLR)/(NLR-PLR)          ## TODO: check Adam's math!
    post_odds_pos  <- prior_odds * PLR
    post_odds_neg  <- prior_odds * NLR
    post_prob_pos  <- post_odds_pos/(1+post_odds_pos)
    post_prob_neg  <- post_odds_neg/(1+post_odds_neg)
} else{
        stop("Couldn't find sens & spec, or positive & negative likelihood ratios")
    }



######################################
########## Plotting (prep)  ##########
######################################


## Set common theme preferences up front
theme_set(theme_bw() +
              theme(axis.text.x = element_blank(),
                    axis.ticks.x = element_blank(),
                    axis.title.x = element_blank(),
                    axis.title.y = element_text(angle = 0),
                    axis.title.y.right = element_text(angle = 0),
                    axis.line = element_blank(),
                    panel.grid = element_blank(),
                    legend.position = "none"
                    )
          )

## Setting up the points of interest along the y-axes

# Select probabilities of interest (nb as percentages)
ticks_prob    <- c(0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 30,
                   40, 50, 60, 70, 80, 90, 95, 99)
# Convert % to odds
ticks_odds    <- odds(ticks_prob/100)
# Convert % to logodds
ticks_logodds <- logodds(ticks_prob/100)

# Select the likelihood ratios of interest (for the middle y-axis)
ticks_lrs     <- sort(c(10^(-3:3), 2*(10^(-3:2)), 5*(10^(-3:2))))
# Log10 them since plot is in logodds space
ticks_log_lrs <- log10(ticks_lrs)




## Fixing particular x-coordinates
left     <- 0
right    <- 1
middle   <- 0.5
midright <- 0.75

## Lay out the four key plot points
##  (the start and finish of the positive and negative lines)

# Initially these are expressed as probabilities
df <- data.frame(x=c(left, right, left, right),
                 y=c(prior_prob, post_prob_pos, prior_prob, post_prob_neg),
                 line = c("pos", "pos", "neg", "neg"))

adj_min      <- range(ticks_logodds)[1]
adj_max      <- range(ticks_logodds)[2]
adj_diff     <- adj_max - adj_min
scale_factor <- abs(adj_min) - adj_diff/2
#df$lo_y <- ifelse(df$x==left,(10/adj_diff)*logodds(1-df$y)-1,logodds(df$y))

# Convert probabilities to logodds for plotting
df$lo_y  <- ifelse(df$x==left,logodds(1-df$y)-scale_factor,logodds(df$y))
# zero         <- data.frame(x = c(left,right),
#                            y = c(0,0),
#                            line = c('pos','pos'),
#                            lo_y = c(-scale_factor,0))





rescale   <- range(ticks_logodds) + abs(adj_min) - adj_diff/2
rescale_x_breaks  <- ticks_logodds + abs(adj_min) - adj_diff/2



######################################
########## Plot             ##########
######################################


p <- ggplot(df) +
        geom_line(aes(x = x, y = lo_y, color = line), size = 1) +
        geom_vline(xintercept = middle) +
        annotate(geom = "text",
                 x = rep(middle+.075, length(ticks_log_lrs)),
                 y = (ticks_log_lrs-scale_factor)/2,
                 label = ticks_lrs,
                 size = rel(LabelSize)) +
        annotate(geom="point",
                 x = rep(middle, length(ticks_log_lrs)),
                 y = (ticks_log_lrs-scale_factor)/2,
                 size = 1) +
        scale_x_continuous(expand = c(0,0)) +
        scale_y_continuous(expand = c(0,0),
                           limits = rescale,
                           breaks = -rescale_x_breaks,
                           labels = ticks_prob,
                           name = "prior \n prob.",
                           sec.axis = sec_axis(trans = ~.,
                                               name = "posterior \n prob.",
                                               labels = ticks_prob,
                                               breaks = ticks_logodds))

## Optional overlay text: prevalence, PLR/NLR, and posterior probabilities
detailedAnnotation <- paste(
    paste0("prevalence = ", p2percent(prior_prob)),
    paste("PLR =", signif(PLR, 3),", NLR =", signif(NLR, 3)),
    paste("post. pos =", p2percent(post_prob_pos),
          ", neg =", p2percent(post_prob_neg)),
    sep = "\n")


## Optional amendments to the plot

## Do we add the null line i.e. LR = 1, illustrating an uninformative model
if(NullLine == TRUE){
    ## If yes, first calculate the start and end points
    uninformative <- data.frame(
        x = c(left,right),
        lo_y = c( (logodds(1-prior_prob) - scale_factor) , logodds(prior_prob))
    )

    p <- p + geom_line(aes(x = x, y = lo_y), data = uninformative,
                       color = "gray",
                       lty = 2,
                       inherit.aes = FALSE)
}


## Do we add the detailed stats to the top right?
if(Detail == TRUE){
    p <- p + annotate(geom = "text",
                      x = midright,
                      y = 2,
                      label = detailedAnnotation,
                      size = rel(LabelSize))
    }

if(Verbose == TRUE){
    writeLines(
        text = c(
            paste0("prevalence = ", p2percent(prior_prob)),
            paste("PLR =", signif(PLR, 3)),
            paste("NLR =", signif(NLR, 3)),
            paste("posterior probability (positive) =", p2percent(post_prob_pos)),
            paste("posterior probability (negative) =", p2percent(post_prob_neg)),
            paste("sensitivity =", p2percent(sensitivity)),
            paste("specificity =", p2percent(specificity))
            # sep = "\n"
        )
    )
    }


return(p)

}






### TODO:
# Allow ppl to input the confusion matrix
# input_TP
# input_FP
# input_FN
# input_TN
#                            Obs
#                   present       absent
#            +----------------+----------------+
#      pos   | True Positive  | False Positive |
# Pred       +----------------+----------------+
#      neg   | False Negative | True Negative  |
#            +----------------+----------------+