R/riskyr_class.R

Defines functions print.summary.riskyr summary.riskyr plot.riskyr riskyr

Documented in plot.riskyr print.summary.riskyr riskyr summary.riskyr

## riskyr_class.R | riskyr
## 2022 08 09
## Define riskyr class and corresponding methods
## -----------------------------------------------

## - Re-define df_scenarios as a list of riskyr objects.

## (A) Create riskyr objects: ---------------

## Get some exemplary scenarios with "riskyr" class attribute -----

# scenario1 <- df_scenarios[1, ]  # get scenario 3 of df_scenarios
# class(scenario1) <- "riskyr"

# scenario2 <- df_scenarios[2, ]  # get scenario 2 of df_scenarios
# class(scenario2) <- "riskyr"


## (1) Function to create diagnostic riskyr scenarios: ------

## riskyr Documentation: ------

#' Create a riskyr scenario.
#'
#' \code{riskyr} creates a scenario of class "riskyr",
#' which can be visualized by the \code{plot} method \code{\link{plot.riskyr}}
#' and summarized by the \code{summary} method \code{\link{summary.riskyr}}.
#'
#' Beyond basic scenario information (i.e., text elements describing a scenario)
#' only the population size \code{\link{N}} and the essential probabilities
#' \code{\link{prev}}, \code{\link{sens}}, \code{\link{spec}}, and \code{\link{fart}}
#' are used and returned.
#'
#' Note:
#'
#' \itemize{
#'
#'   \item Basic text information and some numeric parameters
#'   (see \code{\link{num}} and \code{\link{init_num}})
#'   are integral parts of a \code{riskyr} scenario.
#'
#'   \item By contrast, basic \emph{color} information
#'   (see \code{\link{pal}} and \code{\link{init_pal}})
#'   is not an integral part, but independently defined.
#'
#'   \item The names of \emph{probabilities}
#'   (see \code{\link{prob}}) are currently
#'   not an integral part of \code{txt} and \code{riskyr} scenarios
#'   (but defined in \code{prob_lbl_def} and \code{label_prob}).
#' }
#'
#' @format An object of class "riskyr" with textual and numeric information
#' describing a risk-related scenario.
#'
#' @return An object of class "riskyr" describing a risk-related scenario.
#'
#' Scenario-specific titles and text labels (see \code{\link{txt}}).
#'
#' @param scen_lbl The current scenario title (sometimes in Title Caps).
#'
#' @param popu_lbl A brief description of the current population or sample.
#' @param N_lbl A label for the current population \code{\link{popu}} or sample.
#'
#' @param cond_lbl A label for the \emph{condition} or feature (e.g., some disease) currently considered.
#' @param cond_true_lbl A label for the \emph{presence} of the current condition
#' or \code{\link{cond_true}} cases (the condition's true state of TRUE).
#' @param cond_false_lbl A label for the \emph{absence} of the current condition
#' or \code{\link{cond_false}} cases (the condition's true state of FALSE).
#'
#' @param dec_lbl A label for the \emph{decision} or judgment (e.g., some diagnostic test) currently made.
#' @param dec_pos_lbl A label for \emph{positive} decisions
#' or \code{\link{dec_pos}} cases (e.g., predicting the presence of the condition).
#' @param dec_neg_lbl A label for \emph{negative} decisions
#' or \code{\link{dec_neg}} cases (e.g., predicting the absence of the condition).
#'
#' @param acc_lbl A label for \emph{accuracy} (i.e., correspondence between condition and decision or judgment).
#' @param dec_cor_lbl A label for \emph{correct} (or accurate) decisions or judgments.
#' @param dec_err_lbl A label for \emph{incorrect} (or erroneous) decisions or judgments.
#'
#' @param sdt_lbl A label for the combination of \emph{condition} and \emph{decision} currently made.
#' @param hi_lbl A label for \emph{hits} or \emph{true positives} \code{\link{hi}}
#' (i.e., correct decisions of the presence of the condition, when the condition is actually present).
#' @param mi_lbl A label for \emph{misses} or \emph{false negatives} \code{\link{mi}}
#' (i.e., incorrect decisions of the absence of the condition when the condition is actually present).
#' @param fa_lbl A label for \emph{false alarms} or \emph{false positives} \code{\link{fa}}
#' (i.e., incorrect decisions of the presence of the condition when the condition is actually absent).
#' @param cr_lbl A label for \emph{correct rejections} or \emph{true negatives} \code{\link{cr}}
#' (i.e., a correct decision of the absence of the condition, when the condition is actually absent).
#'
#' Essential probabilities:
#'
#' @param prev The condition's prevalence \code{\link{prev}}
#' (i.e., the probability of condition being \code{TRUE}).
#'
#' @param sens The decision's sensitivity \code{\link{sens}}
#' (i.e., the conditional probability of a positive decision
#' provided that the condition is \code{TRUE}).
#' \code{sens} is optional when its complement \code{mirt} is provided.
#'
#' @param spec The decision's specificity value \code{\link{spec}}
#' (i.e., the conditional probability
#' of a negative decision provided that the condition is \code{FALSE}).
#' \code{spec} is optional when its complement \code{fart} is provided.
#'
#' @param fart The decision's false alarm rate \code{\link{fart}}
#' (i.e., the conditional probability
#' of a positive decision provided that the condition is \code{FALSE}).
#' \code{fart} is optional when its complement \code{spec} is provided.
#'
#' Essential frequencies:
#'
#' @param N The number of individuals in the scenario's population.
#' A suitable value of \code{\link{N}} is computed, if not provided.
#'
#' @param hi The number of hits \code{\link{hi}} (or true positives).
#' @param mi The number of misses \code{\link{mi}} (or false negatives).
#' @param fa The number of false alarms \code{\link{fa}} (or false positives).
#' @param cr The number of correct rejections \code{\link{cr}} (or true negatives).
#'
#' Details and source information:
#'
#' @param scen_lng Language of the current scenario (as character code).
#' Options: \code{"en"} for English, \code{"de"} for German.
#'
#' @param scen_txt A longer text description of the current scenario
#' (which may extend over several lines).
#'
#' @param scen_src Source information for the current scenario.
#'
#' @param scen_apa Source information for the current scenario
#' according to the American Psychological Association (APA style).
#'
#' @param round  Boolean value that determines whether frequency values
#' are rounded to the nearest integer.
#' Default: \code{round = TRUE}.
#'
#' Note: Only rounding when using \code{\link{comp_freq_prob}}
#' (i.e., computing \code{\link{freq}} from \code{\link{prob}} description).
#'
#' @param sample  Boolean value that determines whether frequency values
#' are sampled from \code{N}, given the probability values of
#' \code{prev}, \code{sens}, and \code{spec}.
#' Default: \code{sample = FALSE}.
#'
#' Note: Only sampling when using \code{\link{comp_freq_prob}}
#' (i.e., computing \code{\link{freq}} from \code{\link{prob}} description).
#'
#' @examples
#' # Defining scenarios: -----
#' # (a) minimal information:
#' hustosis <- riskyr(scen_lbl = "Screening for hustosis",
#'                    N = 1000, prev = .04, sens = .80, spec = .95)
#'
#' # (2) detailed information:
#' scen_reoffend <- riskyr(scen_lbl = "Identify reoffenders",
#'                         cond_lbl = "being a reoffender",
#'                         popu_lbl = "Prisoners",
#'                         cond_true_lbl = "has reoffended",
#'                         cond_false_lbl = "has not reoffended",
#'                         dec_lbl = "test result",
#'                         dec_pos_lbl = "will reoffend",
#'                         dec_neg_lbl = "will not reoffend",
#'                         sdt_lbl = "combination",
#'                         hi_lbl = "reoffender found", mi_lbl = "reoffender missed",
#'                         fa_lbl = "false accusation", cr_lbl = "correct release",
#'                         prev = .45,  # prevalence of being a reoffender.
#'                         sens = .98,
#'                         spec = .46, fart = NA,  # (provide 1 of 2)
#'                         N = 753,
#'                         scen_src = "Example scenario")
#'
#' # Using scenarios: -----
#' summary(hustosis)
#' plot(hustosis)
#'
#' summary(scen_reoffend)
#' plot(scen_reoffend)
#'
#' # 2 ways of defining the same scenario:
#' s1 <- riskyr(prev = .5, sens = .5, spec = .5, N = 100)  # s1: define by 3 prob & N
#' s2 <- riskyr(hi = 25, mi = 25, fa = 25, cr = 25)        # s2: same scenario by 4 freq
#' all.equal(s1, s2)  # should be TRUE
#'
#' # Rounding and sampling:
#' s3 <- riskyr(prev = 1/3, sens = 2/3, spec = 6/7, N = 100, round = FALSE)  # s3: w/o rounding
#' s4 <- riskyr(prev = 1/3, sens = 2/3, spec = 6/7, N = 100, sample = TRUE)  # s4: with sampling
#'
#' # Note:
#' riskyr(prev = .5, sens = .5, spec = .5, hi = 25, mi = 25, fa = 25, cr = 25)  # works (consistent)
#' riskyr(prev = .5, sens = .5, spec = .5, hi = 25, mi = 25, fa = 25)           # works (ignores freq)
#'
#' ## Watch out for:
#' # riskyr(hi = 25, mi = 25, fa = 25, cr = 25, N = 101)  # warns, uses actual sum of freq
#' # riskyr(prev = .4, sens = .5, spec = .5, hi = 25, mi = 25, fa = 25, cr = 25)  # warns, uses freq
#'
#' @family riskyr scenario functions
#' @family functions initializing scenario information
#'
#' @seealso
#' \code{\link{init_num}} and \code{\link{num}} for basic numeric parameters;
#' \code{\link{init_txt}} and \code{\link{txt}} for current text settings;
#' \code{\link{init_pal}} and \code{\link{pal}} for current color settings.
#'
#' @export

## riskyr Definition: ------

riskyr <- function(#
  # (1) Scenario label:
  scen_lbl = txt$scen_lbl,  # OR: "" (to hide header/plot titles)
  # Population:
  popu_lbl = txt$popu_lbl,
  N_lbl    = txt$N_lbl,
  # (2) Three perspectives:
  # a. by condition:
  cond_lbl = txt$cond_lbl,
  cond_true_lbl = txt$cond_true_lbl, cond_false_lbl = txt$cond_false_lbl,
  # b. by decision:
  dec_lbl = txt$dec_lbl,
  dec_pos_lbl = txt$dec_pos_lbl, dec_neg_lbl = txt$dec_neg_lbl,
  # c. by accuracy:
  acc_lbl = txt$acc_lbl,
  dec_cor_lbl = txt$dec_cor_lbl, dec_err_lbl = txt$dec_err_lbl,
  # (3) 4 SDT cases/cells/combinations:
  sdt_lbl = txt$sdt_lbl,
  hi_lbl = txt$hi_lbl, mi_lbl = txt$mi_lbl,
  fa_lbl = txt$fa_lbl, cr_lbl = txt$cr_lbl,
  # (4) Essential probabilities:
  prev = NA,
  sens = NA,
  spec = NA, fart = NA,
  # (5) Essential frequencies:
  N = NA,  # WAS: freq$N,
  hi = NA, mi = NA,
  fa = NA, cr = NA,
  # (6) Scenario details:
  scen_lng = txt$scen_lng,
  scen_txt = txt$scen_txt,
  scen_src = txt$scen_src,
  scen_apa = txt$scen_apa,
  # (7) Creating freq from prob (by description):
  round = TRUE,   # round freq values to integers?
  sample = FALSE  # sample freq values from probabilities?
) {

  ## (0): Initialize some stuff: ------
  freqs <- NA
  probs <- NA
  prob_quintet <- NA

  ## Case_1: Using 4 frequencies: -------

  if (!any(is.na(c(hi, mi, fa, cr)))) {  # all four frequencies are provided:

    # (a) Check consistency of N: -----

    if (is.na(N)) { # check, whether N is NA:

      N <- sum(c(hi, mi, fa, cr))  # set N to sum of frequencies.

    } else { # N is provided:

      # check, whether N matches the sum of frequencies:
      N.sum <- sum(c(hi, mi, fa, cr))

      if (N != N.sum) {

        msg <- paste0("N = ", N, ", but the sum of frequencies = ", N.sum, ". Using the latter...")
        warning(msg)

        N <- N.sum  # use sum of frequencies as N

      }
    }

    # (b) Compute probabilities: -----
    probs_calc <- comp_prob_freq(hi, mi, fa, cr)
    probs      <- c(probs_calc$prev, probs_calc$sens, probs_calc$mirt, probs_calc$spec, probs_calc$fart)
    need_probs <- FALSE  # flag that probs are no longer needed

    # (c) Calculate key frequencies from 4 essential frequencies:
    freqs <- comp_freq_freq(hi, mi, fa, cr)

  } else {  # if not all 4 essential frequencies are provided:

    probs <- NA         # set probs to NA (and use inputs below).
    need_probs <- TRUE  # set flag that probs are still needed

  } # if (!any(is.na(c(hi, mi, fa, cr))))...


  ## Case_2: Using 3 essential probabilities: ------

  if (is_valid_prob_set(prev = prev,
                        sens = sens, mirt = NA,
                        spec = spec, fart = fart,
                        tol = .01)) {  # a valid set of probabilities is provided:

    # (a) Compute the complete quintet of probabilities:
    prob_quintet <- comp_complete_prob_set(prev, sens, mirt = NA, spec, fart)
    # sens <- prob_quintet[2] # gets sens (if not provided)
    # mirt <- prob_quintet[3] # gets mirt (if not provided)
    # spec <- prob_quintet[4] # gets spec (if not provided)
    # fart <- prob_quintet[5] # gets fart (if not provided)

    # (b) If frequencies have been provided, test whether the probabilities match:
    if (!any(is.na(probs))) { # if probs has been calculated:

      if (!is.null(prob_quintet)) {  # check, whether prob_quintet has been calculated.

        ## Do they equal the provided probabilities?
        if (any(abs(probs - prob_quintet) > .01)) {

          warning("Probabilities provided differ from probabilities calculated from frequencies.\nUsing probabilities from frequencies...")

        }
      }  # end check prob_quintet.

    } else {  # if no frequencies have been provided (probs is NA): ------

      # (c) Set probs to the computed prob quintet:
      probs <- prob_quintet

      # (d) if no N has been provided:
      if (is.na(N)) {

        N <- comp_min_N(prev = probs[1], sens = probs[2], spec = probs[4],
                        min_freq = 1)  # calculate a suitable N.
      }

      # (e) Calculate frequencies from probabilities (by description):
      freqs <- comp_freq_prob(prev = probs[1], sens = probs[2], mirt = probs[3],
                              spec = probs[4], probs[5], N = N,
                              round = round,   # round freq values to integers?
                              sample = sample  # sample freq values from probabilities?
      )

    }
  } # if (is_valid_prob_set(...

  else { # no valid set of probabilities was provided:

    if (need_probs) {  # not all 4 essential frequencies were provided above:

      ## Case_3: Not all frequencies OR a valid set of probabilities were provided:
      warning("Neither a full set of frequencies nor a valid set of probabilities was provided.")

    }
  }


  ## Case_4: Something is missing: ------

  # if (is.na(freqs)) {
  #   warning("Frequencies were not provided or could not be computed.")
  # }
  #
  # if (is.na(probs)) {
  #   warning("Probabilities were not provided or could not be computed.")
  # }
  #
  # prob_quintet <- probs  # both should be the same by now (not needed?).


  ## Define riskyr scenario object (as a list): ------

  object <- list(#
    # (1) Scenario text:
    # Title label:
    scen_lbl = scen_lbl,
    # population:
    popu_lbl = popu_lbl,
    N_lbl = N_lbl,
    # a. by condition:
    cond_lbl = cond_lbl,
    cond_true_lbl = cond_true_lbl, cond_false_lbl = cond_false_lbl,
    # b. by decision:
    dec_lbl = dec_lbl,
    dec_pos_lbl = dec_pos_lbl, dec_neg_lbl = dec_neg_lbl,
    # c. by accuracy:
    acc_lbl = acc_lbl,
    dec_cor_lbl = dec_cor_lbl, dec_err_lbl = dec_err_lbl,
    # 4 SDT cases:
    sdt_lbl = sdt_lbl,
    hi_lbl = hi_lbl, mi_lbl = mi_lbl, fa_lbl = fa_lbl, cr_lbl = cr_lbl,
    # (2) Numeric info:
    # Probabilities:
    prev = probs[1],
    sens = probs[2],
    spec = probs[4], fart = probs[5],
    # Frequencies:
    N = N,
    hi = freqs$hi, mi = freqs$mi,
    fa = freqs$fa, cr = freqs$cr,
    # (+) Scenario details:
    scen_lng = scen_lng, scen_txt = scen_txt,
    scen_src = scen_src, scen_apa = scen_apa
  )

  ## Add class riskyr:
  class(object) <- "riskyr"

  return(object)

} # riskyr().


## Check: -------
# test.obj <- riskyr()  # initialize with default parameters
# names(test.obj)

## 2 ways to define the same scenario:
# s1 <- riskyr(prev = .5, sens = .5, spec = .5, N = 100)  # s1
# s2 <- riskyr(hi = 25, mi = 25, fa = 25, cr = 25)        # s2: same in terms of freq
# all.equal(s1, s2)  # should be TRUE
#
## Rounding and sampling:
# s3 <- riskyr(prev = 1/3, sens = 2/3, spec = 6/7, N = 100, round = FALSE)  # s3: w/o rounding
# s4 <- riskyr(prev = 1/3, sens = 2/3, spec = 6/7, N = 100, sample = TRUE)  # s4: with sampling
#
## Note:
# riskyr(prev = .5, sens = .5, spec = .5, hi = 25, mi = 25, fa = 25, cr = 25)  # works (consistent)
# riskyr(prev = .5, sens = .5, spec = .5, hi = 25, mi = 25, fa = 25)           # works (ignores freq)
#
## Watch out for:
# riskyr(hi = 25, mi = 25, fa = 25, cr = 25, N = 101)  # warn: use sum of freq
# riskyr(prev = .4, sens = .5, spec = .5, hi = 25, mi = 25, fa = 25, cr = 25)  # warn: use freq

## Check incomplete or inconsistent entries:
## riskyr(prev = NA, hi = NA)

## Compare with df_scenarios:
# names(df_scenarios)
# all.equal(names(test.obj), names(df_scenarios))

# # cat(
# #   paste0(
# #     paste0(names(scenarios$scen1), " = ", names(scenarios$scen1)),
# #     collapse = ", "))



## (B) Create scenarios: Define scenarios as a list of riskyr objects -----------

## Note: Convert the data frame df_scenarios into a list "scenarios"
##       of riskyr objects:

scenarios <- NULL # initialize

## Note helper stuff:
# cat(paste0("#'   \\item ", df_scenarios$scen_lbl[-1], "\n#'\n"))

## scenarios: Documentation ------

#' A collection of riskyr scenarios from various sources (as list).
#'
#' \code{scenarios} is a list of scenarios of class \code{riskyr}
#' collected from the scientific literature and other sources
#' and to be used by visualization and summary functions.
#'
#' \code{scenarios} currently contains the following scenarios
#' (n1 to n12 in English language, n13 to n25 in German language):
#'
#' \enumerate{
#'
#'   \item Bowel cancer screening
#'   \item Cab problem
#'   \item Hemoccult test
#'   \item Mammography screening
#'   \item Mammography (freq)
#'   \item Mammography (prob)
#'   \item Mushrooms
#'   \item Musical town
#'   \item PSA test (baseline)
#'   \item PSA test (patients)
#'   \item Psylicraptis screening
#'   \item Sepsis
#'
#'   \item Amniozentese (in German language)
#'   \item HIV-Test 1
#'   \item HIV-Test 2
#'   \item HIV-Test 3
#'   \item HIV-Test 4
#'   \item Mammografie 1
#'   \item Mammografie 2
#'   \item Mammografie 3
#'   \item Mammografie 4
#'   \item Nackenfaltentest (NFT) 1
#'   \item Nackenfaltentest (NFT) 2
#'   \item Sigmoidoskopie 1
#'   \item Sigmoidoskopie 2
#' }
#'
#' Variables describing a scenario:
#'
#' \enumerate{
#'
#'   \item \code{scen_lbl}: Text label for current scenario.
#'   \item \code{scen_lng}: Language of current scenario (en/de).
#'   \item \code{scen_txt}: Description text of current scenario.
#'
#'   \item \code{popu_lbl}: Text label for current population.
#'
#'   \item \code{cond_lbl}: Text label for current condition.
#'   \item \code{cond_true_lbl}: Text label for \code{\link{cond_true}} cases.
#'   \item \code{cond_false_lbl}: Text label for \code{\link{cond_false}} cases.
#'
#'   \item \code{dec_lbl}: Text label for current decision.
#'   \item \code{dec_pos_lbl}: Text label for \code{\link{dec_pos}} cases.
#'   \item \code{dec_neg_lbl}: Text label for \code{\link{dec_neg}} cases.
#'
#'   \item \code{hi_lbl}: Text label for cases of hits \code{\link{hi}}.
#'   \item \code{mi_lbl}: Text label for cases of misses \code{\link{mi}}.
#'   \item \code{fa_lbl}: Text label for cases of false alarms \code{\link{fa}}.
#'   \item \code{cr_lbl}: Text label for cases of correct rejections \code{\link{cr}}.
#'
#'   \item \code{prev}: Value of current prevalence \code{\link{prev}}.
#'   \item \code{sens}: Value of current sensitivity \code{\link{sens}}.
#'   \item \code{spec}: Value of current specificity \code{\link{spec}}.
#'   \item \code{fart}: Value of current false alarm rate \code{\link{fart}}.
#'
#'   \item \code{N}: Current population size \code{\link{N}}.
#'
#'   \item \code{scen_src}: Source information for current scenario.
#'   \item \code{scen_apa}: Source information in APA format.
#'
#' }
#'
#' Note that names of variables (columns)
#' correspond to a subset of \code{\link{init_txt}} (to initialize \code{\link{txt}})
#' and \code{\link{init_num}} (to initialize \code{\link{num}}).
#'
#' The variables \code{scen_src} and \code{scen_apa}
#' provide a scenario's source information.
#'
#' The information of \code{scenarios} is also contained in an
#' R data frame \code{\link{df_scenarios}} (and generated from
#' the corresponding \code{.rda} file in \code{/data/}).
#'
#' @format A list with currently 25 scenarios of class \code{riskyr}
#' which are each described by 21 variables.
#'
#' @seealso
#' \code{\link{riskyr}} initializes a \code{riskyr} scenario.
#'
#' @export

## scenarios: Definition ------

scenarios <- vector("list", nrow(df_scenarios))  # initialize scenarios as a list (from df_scenarios)
names(scenarios) <- paste0("n", 1:nrow(df_scenarios))

for (i in 1:nrow(df_scenarios)) {  # for each scenario i in df_scenarios:

  ## (1) define scenario s:
  s <- df_scenarios[i, ]

  ## (2) pass scenario s to riskyr function:
  cur_scen <- riskyr(#
    # Scenario label:
    scen_lbl = s$scen_lbl,
    # Population:
    popu_lbl = s$popu_lbl,
    N_lbl    = txt$N_lbl,  # use txt default (as currently not set in data)
    # a. by condition:
    cond_lbl = s$cond_lbl,
    cond_true_lbl = s$cond_true_lbl,
    cond_false_lbl = s$cond_false_lbl,
    # b. by decision:
    dec_lbl = s$dec_lbl,
    dec_pos_lbl = s$dec_pos_lbl,
    dec_neg_lbl = s$dec_neg_lbl,
    # c. by accuracy:
    acc_lbl = txt$acc_lbl,  # use txt default (as currently not set in data)
    dec_cor_lbl = txt$dec_cor_lbl,
    dec_err_lbl = txt$dec_err_lbl,
    # 4 SDT cases:
    sdt_lbl = txt$sdt_lbl,  # use txt default (as currently not set in data)
    hi_lbl = s$hi_lbl,
    mi_lbl = s$mi_lbl,
    fa_lbl = s$fa_lbl,
    cr_lbl = s$cr_lbl,
    # Probabilities:
    prev = s$prev,
    sens = s$sens,
    spec = s$spec,
    fart = s$fart,
    # Frequencies:
    N = s$N,
    # Scenario details:
    scen_lng = s$scen_lng,
    scen_txt = s$scen_txt,
    scen_src = s$scen_src,
    scen_apa = s$scen_apa
  )

  # (3) Add cur_scen (riskyr object) as i-th element of scenarios
  scenarios[[i]] <- cur_scen

} # end for ...


## Check: --------
# length(scenarios)
# scenarios$n25  # => shows elements of a scenario



## (C) Handle riskyr objects: ------------------


## 1. plot.riskyr function: -------

## Testing dots:
# test_fun <- function(...) {
#   plot_icons(...)
# }
#
# test_fun(N = 100, blubb = 5, prev = .7)
#
## ok.


## plot.riskyr Documentation: ------

#' Plot a riskyr scenario.
#'
#' \code{plot.riskyr} is a method that allows to generate
#' different plot types from a \code{"riskyr"} object.
#'
#' \code{plot.riskyr} also uses the text settings
#' specified in the "riskyr" object.
#'
#' @param x An object of class "riskyr", usually a result of a call to \code{\link{riskyr}}.
#' Pre-defined \code{\link{scenarios}} are also of type "riskyr".
#'
#' @param type The type of plot to be generated.
#'
#' @param main Text label for main plot title.
#' Default: \code{main = NULL} (using \code{x$scen_lbl} per default).
#'
#' @param sub Text label for plot subtitle (on 2nd line).
#' Default: \code{sub = NULL} (using \code{sub = "type"} shows plot type).
#'
#' The following plot types are currently available:
#'
#' \enumerate{
#'
#'   \item \code{type = "prism"} or \code{type = "net"} or \code{type = "tree"}:
#'   Risk information is plotted in a network diagram of frequencies and probabilities (default).
#'   See \code{\link{plot_prism}} for further options.
#'
#'   \item \code{type = "tab"} or \code{type = "ftab"}:
#'   Risk information is plotted as a 2-by-2 frequency or contingency table.
#'   See \code{\link{plot_tab}} for further options.
#'
#'   \item \code{type = "area"} or \code{type = "mosaic"}:
#'   Risk information is plotted as a mosaic plot (scaled area).
#'   See \code{\link{plot_area}} for further options.
#'
#'   \item \code{type = "bar"} or \code{type = "fbar"}:
#'   Risk information is plotted as a bar chart.
#'   See \code{\link{plot_bar}} for further options.
#'
#'   \item \code{type = "icons"} or \code{type = "iconarray"}:
#'   The underlying population is plotted as an array of icons.
#'   See \code{\link{plot_icons}} for further options.
#'
#'   \item \code{type = "curve"} or \code{type = "curves"}:
#'   Draws curves of selected values (including \code{\link{PPV}}, \code{\link{NPV}})
#'   See \code{\link{plot_curve}} for further options.
#'
#'   \item \code{type = "plane"} or \code{type = "planes"}:
#'   Draws a 3D-plane of selected values (e.g., predictive values \code{\link{PPV}} or \code{\link{NPV}})
#'   See \code{\link{plot_plane}} for further options.
#'
#' }
#'
#' @param ... Additional parameters to be passed to the underlying plotting functions.
#'
#' @examples
#' # Select a scenario (from list of scenarios):
#' s1 <- scenarios$n1  # select scenario 1 from scenarios
#' plot(s1)  # default plot (type = "prism")
#'
#' # Plot types currently available:
#' plot(s1, type = "prism")                # prism/network diagram (default)
#' plot(s1, type = "tree", by = "cd")      # tree diagram (only 1 perspective)
#' plot(s1, type = "area")                 # area/mosaic plot
#' plot(s1, type = "tab")                  # 2x2 frequency/contingency table
#' plot(s1, type = "bar", dir = 2)         # bar plot
#' plot(s1, type = "icons")                # icon array
#' plot(s1, type = "curve", what = "all")  # curves as fn. of prev
#' plot(s1, type = "plane", what = "NPV")  # plane as function of sens & spec
#' plot(s1, type = "default")              # unknown type: use default plot
#'
#' @family visualization functions
#' @family riskyr scenario functions
#'
#' @seealso
#' \code{\link{riskyr}} initializes a \code{riskyr} scenario.
#'
#' @export

## plot.riskyr Definition: ------

plot.riskyr <- function(x = NULL,        # a riskyr scenario
                        type = "prism",  # default type
                        # by = "cddc",   # default perspective
                        main = NULL,     # main title
                        sub = NULL,      # subtitle
                        ...              # other parameters (passed to plot_xxx functions)
) {

  ## (1) Increase robustness: ------

  type <- tolower(type)  # ensure lowercase

  # Test type argument:
  if (all(type %in% c(#
    # plot_prism:
    "prism", "fprism", "tree", "ftree", "dtree", "double tree",
    # plot_fnet:
    "net", "fnet", "network", "fnetwork",
    # plot_area:
    "area", "farea", "mosaic", "unit square",
    # plot_tab:
    "tab", "table", "ftab", "ctab", "matrix",
    # plot_icons:
    "icon", "icons", "iconarray",
    # plot_bar:
    "bar", "bars", "barplot", "fbar",
    # plot_curve:
    "curve", "curves",
    # plot_plane:
    "plane", "planes", "cube")) == FALSE) {

    message("Unknown plot type (in plot.riskyr): Using type = 'prism'.")
    type <- "prism"

  }

  # # If type == "tree" and current by contains more than 1 perspective:
  # if ( (substr(type, 1, 4) == "tree") || (substr(type, 1, 5) == "ftree") ) {
  #   if ( !is.null(by) ) {
  #     if ( is.character(by) && (nchar(by) > 2) ) {
  #
  #       # print(paste0("1. by = ", by))  # debugging
  #       by <- substr(by, 1, 2)           # use only 1st perspective of by
  #       # print(paste0("2. by = ", by))  # debugging
  #
  #     }
  #   }
  # }


  ## (2) Use text info of scenario x for current txt information: ------

  x_txt <- init_txt(scen_lbl = x$scen_lbl,

                    popu_lbl = x$popu_lbl,
                    N_lbl = x$N_lbl,

                    cond_lbl = x$cond_lbl,
                    cond_true_lbl = x$cond_true_lbl,
                    cond_false_lbl = x$cond_false_lbl,

                    dec_lbl  = x$dec_lbl,
                    dec_pos_lbl = x$dec_pos_lbl,
                    dec_neg_lbl = x$dec_neg_lbl,

                    acc_lbl = x$acc_lbl,
                    dec_cor_lbl = x$dec_cor_lbl,
                    dec_err_lbl = x$dec_err_lbl,

                    sdt_lbl = x$sdt_lbl,
                    hi_lbl = x$hi_lbl,
                    mi_lbl = x$mi_lbl,
                    fa_lbl = x$fa_lbl,
                    cr_lbl = x$cr_lbl,

                    scen_txt = x$scen_txt,
                    scen_src = x$scen_src,
                    scen_apa = x$scen_apa,
                    scen_lng = x$scen_lng
  )

  # Set default main title:
  if (is.null(main) & (nchar(x$scen_lbl) > 0)){
    main <- x$scen_lbl
  }

  ## (3) Call plotting functions: ------

  # 1. Table/contingency/confusion/frequency table/tab plot: ----
  if ((substr(type, 1, 3) == "tab") || (type == "ftab") || (type == "ctab")) {

    plot_tab(prev = x$prev,
             sens = x$sens, mirt = NA,
             spec = x$spec, fart = NA,
             N = x$N,
             # Options:
             lbl_txt = x_txt,
             main = main,
             sub = sub,
             ...
    )

  } # if (type == "tab")


  # 2. Area / mosaic plot / unit square: ----
  if ((substr(type, 1, 4) == "area") || (type == "farea") ||
      (substr(type, 1, 6) == "mosaic") ||
      (substr(type, 1, 4) == "unit")) {  # "mosaic"

    plot_area(prev = x$prev,
              sens = x$sens, mirt = NA,
              spec = x$spec, fart = NA,
              N = x$N,
              # Options:
              lbl_txt = x_txt,
              main = main,
              sub = sub,
              ...
    )

  } # if (type == "area")


  # 3. Icon array: ----
  if (substr(type, 1, 4) == "icon") {

    plot_icons(prev = x$prev,             # probabilities
               sens = x$sens, mirt = NA,
               spec = x$spec, fart = NA,  # was: num$fart,
               N = x$N,    # ONLY freq used (so far)
               # Options:
               # lbl_txt = x_txt, # Do NOT pass for icons!
               main = main,
               sub = sub,
               type_lbls = x[c("hi_lbl", "mi_lbl", "fa_lbl", "cr_lbl")],
               ...
    )

  } #  if (type == "icon")


  # 4. Prism plot (tree/double tree): ----
  if ((substr(type, 1, 5) == "prism") || (substr(type, 1, 6) == "fprism") ||
      (substr(type, 1, 4) == "tree")  || (substr(type, 1, 5) == "ftree") ||
      (substr(type, 1, 6) == "double") || (substr(type, 1, 5) == "dtree")) {

    plot_prism(prev = x$prev,
               sens = x$sens, mirt = NA,
               spec = x$spec, fart = NA,
               N = x$N,
               # Options:
               lbl_txt = x_txt,
               main = main,
               sub = sub,
               ...
    )

  } # if (type == "prism")


  # 5. Frequency net plot (fnet): ----
  if ((type == "frequency net") ||
      (substr(type, 1, 3) == "net") || (substr(type, 1, 4) == "fnet")) {

    plot_fnet(prev = x$prev,
              sens = x$sens, mirt = NA,
              spec = x$spec, fart = NA,
              N = x$N,
              # Options:
              lbl_txt = x_txt,
              main = main,
              sub = sub,
              ...
    )

  } # if (type == "fnet")

  # 6. Bar plot / frequency bars: ----
  if ((substr(type, 1, 3) == "bar") || (substr(type, 1, 4) == "fbar")) {

    plot_bar(prev = x$prev,
             sens = x$sens, mirt = NA,
             spec = x$spec, fart = NA,
             N = x$N,
             # Options:
             lbl_txt = x_txt,
             main = main,
             sub = sub,
             ...
    )

  } # if (type == "bar")


  # 7. Curve of probabilities: ----
  if (substr(type, 1, 5) == "curve") {

    plot_curve(prev = x$prev,             # probabilities (3 essential, 2 optional)
               sens = x$sens, mirt = NA,
               spec = x$spec, fart = NA,
               # Options:
               lbl_txt = x_txt,
               main = main,
               sub = sub,
               ...
    )
  } # if (type == "curve")


  # 8. Plane/cube of probabilities: ----
  if ((substr(type, 1, 5) == "plane") || (substr(type, 1, 4) == "cube")) {

    plot_plane(prev = x$prev,             # probabilities (3 essential, 2 optional)
               sens = x$sens, mirt = NA,
               spec = x$spec, fart = NA,
               # Options:
               lbl_txt = x_txt,
               main = main,
               sub = sub,
               ...
    )
  } # if (type == "plane")


  ## Currently NO output. ------

} # plot.riskyr().


## Check: ------
## (A) with example scenarios (defined above):
# plot(scenario2, type = "icons")
# plot(scenario3, type = "tree")

## (B) with scenarios from scenarios (defined BELOW):
#
# s25 <- scenarios$n25  # select Scenario 25 from scenarios
#
# plot(s25)  # => default plot (prism/net)
# plot(s25, type = "ofnet")  # => old network diagram (old default)
# plot(s25, type = "tree", area = "vr") # => tree diagram (with vertical rectangles)
# plot(s25, type = "icons")
# plot(s25, type = "icons", type = "mosaic")  # passing on additional parameters.
# plot(s25, type = "curve", what = "all")
# plot(s25, type = "plane", what = "npv")
# plot(s25, type = "wetwork")
#
# New plots:
# plot(s25, type = "prism", by = "cddc", f_lbl = "num")
# plot(s25, type = "tree", by = "ac", f_lbl = "num")
# plot(s25, type = "area", by = "cddc")
# plot(s25, type = "tab", by = "cddc", f_lwd = 2)
# plot(s25, type = "bar", dir = 2)
#
# # Plot types currently available:
# plot(s25, type = "prism")                # prism/network diagram (default)
# plot(s25, type = "tree", by = "cd")      # tree diagram (only 1 perspective)
# plot(s25, type = "area")                 # area/mosaic plot
# plot(s25, type = "tab")                  # 2x2 frequency/contingency table
# plot(s25, type = "bar", dir = 2)         # bar plot
# plot(s25, type = "icons")                # icon array
# plot(s25, type = "curve", what = "all")  # curves as fn. of prev
# plot(s25, type = "plane", what = "NPV")  # plane as function of sens & spec
#
# # Older plot types (obsolete and RETIRED):
# plot(s25, type = "onet")     # plot_fnet (replaced by plot_prism)
# plot(s25, type = "otree")    # plot_tree (replaced by plot_prism)
# plot(s25, type = "omosaic")  # plot_mosaic (replaced by plot_area)



## Legacy code snippets: ------

## From documentation of older plot.riskyr:
# Older plot types (replaced by the above, retired since version 0.1.0.950):
#
# \item \code{type = "onet"} or \code{type = "ofnet"}:
# Risk information is plotted in a network diagram of frequencies and probabilities (default).
# Replaced by \code{\link{plot_prism}}, but see \code{\link{plot_fnet}} for further options.
#
# \item \code{type = "otree"} or \code{type = "oftree"}:
# Risk information is plotted in a frequency tree.
# Replaced by \code{\link{plot_prism}}, but see \code{\link{plot_tree}} for further options.
#
# \item \code{type = "omosaic"} or \code{type = "omosaicplot"}:
# Risk information is plotted as a mosaicplot.
# Replaced by \code{\link{plot_area}}, but see \code{\link{plot_mosaic}} for further options.




## 2a. summary.riskyr function: ------------------

## (a) Create a summary object:

## summary.riskyr Documentation: ------

#' Summarize a riskyr scenario.
#'
#' \code{summary.riskyr} provides a \code{summary} method for objects of class "riskyr".
#'
#' @format An object of class \code{summary.riskyr} with up to 9 entries.
#'
#' @return A summary list \code{obj.sum}
#' with up to 9 entries, dependent on which information is requested by \code{summarize}.
#'
#' Scenario name, relevant condition \code{}, and \code{N}
#' are summarized by default.
#'
#' @param object  An object of class "riskyr", usually a result of a call to \code{\link{riskyr}}.
#' Inbuilt \code{scenarios} are also of type "riskyr".
#'
#' @param summarize What is summarized as a vector consisting of \code{c("freq", "prob", "accu")}
#' for frequencies, probabilities, and accuracy respectively.
#' The default "all" is an alias to all three.
#'
#' @param ... Additional parameters (to be passed to summary functions).
#'
#' @examples
#' summary(scenarios$n4)
#'
#' @family summary functions
#' @family riskyr scenario functions
#'
#' @seealso
#' \code{\link{riskyr}} initializes a \code{riskyr} scenario.
#'
#' @export

## summary.riskyr Definition: ------

summary.riskyr <- function(object = NULL, summarize = "all", ...) {

  obj.sum <- list()  # initialize as list

  obj.sum$scen_lbl <- object$scen_lbl

  obj.sum$cond_lbl <- object$cond_lbl  # condition
  obj.sum$dec_lbl <- object$dec_lbl    # decision
  obj.sum$popu_lbl <- object$popu_lbl  # population
  obj.sum$N <- object$N                # N
  obj.sum$scen_src <- object$scen_src  # source (short)

  ## (0) If all should be summarized: ----------

  if (summarize == "all") summarize <- c("prob", "freq", "accu")


  ## (A) Probability information: ----------

  if (("prob" %in% summarize) || ("probs" %in% summarize) || ("probabilities" %in% summarize)) {

    # calculate all probabilities:
    probs <- comp_prob_prob(prev = object$prev, sens = object$sens, spec = object$spec)

    probs.ess <- unlist(probs[c("prev", "sens", "mirt", "spec", "fart")])  # essential probabilities.

    probs.ness <- unlist(probs[c("ppod", "PPV", "NPV", "FDR", "FOR", "acc")])  # non-essential probabilities.

    obj.sum$probs <- list(probs.ess = probs.ess, probs.ness = probs.ness)

  } # if "prob"


  ## (B) Frequency information: ----------

  if (("freq" %in% summarize) || ("freqs" %in% summarize) || ("frequencies" %in% summarize)) {

    # calculate frequencies:
    freqs <- comp_freq(prev = object$prev, sens = object$sens, spec = object$spec,
                       N = object$N)

    ## (a) Frequencies by condition:
    cond_freqs <- unlist(freqs[c("cond_true", "cond_false")])

    ## (b) Frequencies by decision:
    dec_freqs <- unlist(freqs[c("dec_pos", "dec_neg")])

    ## (c) Frequencies by accuracy (i.e., correspondence of decision to condition):
    acc.freqs <- unlist(freqs[c("dec_cor", "dec_err")])

    ## (d) SDT frequencies:
    sdt.freqs <- unlist(freqs[c("hi", "mi", "fa", "cr")])  # == "essential" frequencies.

    ## (+) Add to summary object:
    obj.sum$freqs <- list(cond_freqs = cond_freqs,
                          dec_freqs = dec_freqs,
                          acc.freqs = acc.freqs,
                          sdt.freqs = sdt.freqs)

  } # if "freq"


  ## (C) Accuracy information: ----------

  if (("acc" %in% summarize) || ("accu" %in% summarize) || ("accuracy" %in% summarize)) {

    ## Overall accuracy acc:
    obj.sum$acc <- comp_acc(prev = object$prev, sens = object$sens, spec = object$spec)

    ## ToDo: ALL accuracy metrics:
    # accu <- comp_accu_prob(prev = obj$prev, sens = obj$sens, spec = obj$spec)

  } # if "acc"



  ## Add class to summary object: ----------

  class(obj.sum) <- c("summary.riskyr")

  return(obj.sum)

} # summary.riskyr().


## 2b. print.summary.riskyr function: ------------------

## Create print function corresponding to summary object:

## print.summary.riskyr Documentation: ------

#' Print summary information of a riskyr scenario.
#'
#' \code{print.summary.riskyr} provides a \code{print} method for objects of class "summary.riskyr".
#'
#' @format Printed output of a "summary.riskyr" object.
#'
#' @param x An object of class "summary.riskyr", usually a result of a call to \code{summary.riskyr}.
#'
#' @param ... Additional parameters (to be passed to generic print function).
#'
#' @examples
#' summary(scenarios$n4)
#'
#' @family print functions
#'
#' @seealso
#' \code{\link{riskyr}} initializes a \code{riskyr} scenario.
#'
#' @export

## print.summary.riskyr Definition: ------

print.summary.riskyr <- function(x = NULL, ...) {

  ## 1. Always print header: ----------

  cat("Scenario: ",   x$scen_lbl, "\n\n")  # always show scenario name.

  cat("Condition: ",  x$cond_lbl, "\n")  # always show current condition.
  cat("Decision: ",   x$dec_lbl,  "\n")  # always show current decision.
  cat("Population: ", x$popu_lbl, "\n")  # always show current condition.
  cat("N = ", x$N, "\n")                 # always show population size N.
  cat("Source: ", x$scen_src, "\n")      # always show (short) source info

  ## 2. Print only on demand: ----------

  n <- names(x)  # save names.

  ## (A) Probabilities: ----------

  if ("probs" %in% n) {

    cat("\nProbabilities:\n\n")

    cat(" Essential probabilities:\n")
    # names(x$probs$probs.ess) <- c("Prevalence", "Sensitivity", "Miss rate", "Specificity", "False alarm rate")  # explicit
    names(x$probs$probs.ess) <- c("prev", "sens", "mirt", "spec", "fart")  # shorter
    print(x$probs$probs.ess)

    cat("\n Other probabilities:\n")
    print(round(x$probs$probs.ness, 3))  # no naming for non-essential probs.

  }

  ## (B) Frequencies: ----------

  if ("freqs" %in% n) {

    cat("\nFrequencies:\n")

    cat("\n by conditions:\n")
    # names(x$freqs$cond_freqs) <- c("True", "False")  # explicit
    names(x$freqs$cond_freqs) <- c("cond_true", "cond_false")  # more explicit
    print(x$freqs$cond_freqs)

    cat("\n by decision:\n")
    names(x$freqs$dec_freqs) <- c("Positive", "Negative")  # explicit
    names(x$freqs$dec_freqs) <- c("dec_pos", "dec_neg")  # more explicit
    print(x$freqs$dec_freqs)

    cat("\n by correspondence (of decision to condition):\n")
    # names(x$freqs$acc.freqs) <- c("Correct cases", "Incorrect cases")  # explicit
    names(x$freqs$acc.freqs) <- c("dec_cor", "dec_err")  # implicit
    print(x$freqs$acc.freqs)

    cat("\n 4 essential (SDT) frequencies:\n")
    # names(x$freqs$sdt.freqs) <- c("Hits", "Misses", "False alarms", "Correct rejections")  # explicit
    names(x$freqs$sdt.freqs) <- c("hi", "mi", "fa", "cr")  # implicit
    print(x$freqs$sdt.freqs)

  }

  ## (C) Accuracy: ----------

  if (("acc" %in% n) || ("accu" %in% n) || ("accuracy" %in% n)) {

    cat("\nAccuracy:\n\n")

    cat(" acc:\n")
    cat(x$acc)  # overall accuracy acc only!

    ## ToDo: Include ALL other accuracy metrics (accu).

  }

}


## Check: ------
# scenario2 <- df_scenarios[2, ]  # get scenario 2 of df_scenarios
# summary(scenario2)  # => all summaries
# summary(scenario2, summarize = "freq")
# summary(scenario2, summarize = "prob")
# summary(scenario2, summarize = "accu")


## (D) Demo: Typical user interaction / session: ----------

## 1. Defining and viewing a scenario: -----------
## see Vignette of riskyr primer

## 2. Exploring pre-defined scenarios: -----------

## Example 1: Bowel cancer (FOB screening): ------
## Source: https://en.wikipedia.org/wiki/Positive_and_negative_predictive_values#Worked_example

## Select scenario 1:
# s1 <- scenarios$n1   # select by name:   $ndd

# summary(s1)
# summary(s1, summarize = "freq")

# plot(s1, type = "tree", area = "hr") # => tree diagram (with vertical rectangles)
# plot(s1, type = "curve", what = "all", uc = .05)
# plot(s1, type = "icons")
# plot(s1, type = "icons", arr_type = "mosaic")  # passing on additional parameters.
# plot(s1, type = "mosaic")
# plot(s1, type = "plane", what = "NPV")
# plot(s1, type = "wetwork")


## Example 2: Mammography screening (standard example) -----
## Source: Hoffrage et al. (2015), p. 3

## (a) Select scenario 6:
# s6 <- scenarios[[6]] # select by number: [[dd]]
# s6 <- scenarios$n6   # select by name:   $ndd

## (b) Summary:
# summary(s6)

## (c) Visualization:
# plot(s6)  # => default plot (fnet)
# plot(s6, type = "icons")
# plot(s6, type = "curve")


## Example 3: PSA screening ------
## Source: Arkes & Gaissmaier (2012), p. 550

## A. Low prevalence version (Scenario 9: prev = 6.3%):

# summary(scenarios$n9)

# plot(scenarios$n9, type = "tree", area = "sq")
# plot(scenarios$n9, type = "icons")
# plot(scenarios$n9, type = "curves", what = "all")
# plot(scenarios$n9, type = "planes", what = "PPV")

## B. Contrast with high prevalence version (Scenario 10: prev = 50%):

# summary(scenarios$n10)

# plot(scenarios$n10, type = "tree", area = "sq")
# plot(scenarios$n10, type = "icons")
# plot(scenarios$n10, type = "curves", what = "all")
# plot(scenarios$n10, type = "planes", what = "PPV")



## (*) Done: ----------

## - etc.

## (+) ToDo: ----------

## - Flexible inputs: Allow riskyr() to take different kinds of inputs (prob, freq, ...),
##   until a full object is created.

## - Allow incomplete riskyr scenarios (e.g., 1 branch of tree)

## eof. ------------------------------------------

Try the riskyr package in your browser

Any scripts or data that you put into this service are public.

riskyr documentation built on Aug. 15, 2022, 9:09 a.m.