Compare 559d018 ... +8 ... 7d49199

#

#   This file contains the mcmcRocPrc() S3 generic, which constructs objects

#   of class "mcmcRocPrc". For methods for this class, see mcmcRocPrc-methods.R

#   S3 methods for the mcmcRocPrc() generic handle different types of input

#   e.g. "rjags" input produced by R2jags.

#

#' probit regression using [rstan::stan()], [rstanarm::stan_glm()], 

#' [R2jags::jags()], [R2WinBUGS::bugs()], [MCMCpack::MCMClogit()], or other 

#' functions that provide samples from a posterior density. 

#' @param object A fitted binary choice model, e.g. "rjags" object 

#'   (see [R2jags::jags()]), or a `[N, iter]` matrix of predicted probabilites.

#' @param yvec A `numeric(N)` vector of observed outcomes. 

#' @param yname (`character(1)`)\cr

#'   The name of the dependent variable, should match the variable name in the 

#'   JAGS data object.

#' @param xnames ([base::character()])\cr

#'   A character vector of the independent variable names, should match the 

#'   corresponding names in the JAGS data object.

#' @param posterior_samples a "mcmc" object with the posterior samples

#' @details If only the average AUC-ROC and PR are of interest, setting 

#'   `curves = FALSE` and `fullsims = FALSE` can greatly speed up calculation 

#'   time. The curve data (`curves = TRUE`) is needed for plotting. The plot

#'   method will always plot both the ROC and PR curves, but the underlying

#'   data can easily be extracted from the output for your own plotting; 

#'   see the documentation of the value returned below. 

#'   

#'   The default method works with a matrix of predicted probabilities and the 

#'   vector of observed incomes as input. Other methods accommodate some of the 

#'   common Bayesian modeling packages like rstan (which returns class "stanfit"),

#'   rstanarm ("stanreg"), R2jags ("jags"), R2WinBUGS ("bugs"), and 

#'   MCMCpack ("mcmc"). Even if a package-specific method is not implemented, 

#'   the default method can always be used as a fallback by manually calculating

#'   the matrix of predicted probabilities for each posterior sample. 

#'   

#'   Note that MCMCpack returns generic "mcmc" output that is annotated with 

#'   some additional information as attributes, including the original function

#'   call. There is no inherent way to distinguish any other kind of "mcmc" 

#'   object from one generated by a proper MCMCpack modeling function, but as a

#'   basic precaution, `mcmcRocPrc()` will check the saved call and return an 

#'   error if the function called was not `MCMClogit()` or `MCMCprobit()`. 

#'   This behavior can be suppressed by setting `force = TRUE`. 

#' 

#'   - "area_under_roc": `numeric()`; either length 1 if `fullsims = FALSE`, or 

#'     one value for each posterior sample otherwise

#'   - "area_under_prc": `numeric()`; either length 1 if `fullsims = FALSE`, or 

#'     one value for each posterior sample otherwise

#'   - "prc_dat": only if `curves = TRUE`; a list with length 1 if 

#'     `fullsims = FALSE`, longer otherwise

#'   - "roc_dat": only if `curves = TRUE`; a list with length 1 if 

#'     `fullsims = FALSE`, longer otherwise

#' 

#' # Equivalently, we can calculate the matrix of predicted probabilities 

#' # ourselves; using the example from ?jags_logit:

#' library(R2jags)

#' 

#' data("sim_data")

#' yvec <- sim_data$Y

#' xmat <- sim_data[, c("X1", "X2")]

#' 

#' # add intercept to the X data

#' xmat <- as.matrix(cbind(Intercept = 1L, xmat))

#' 

#' beta <- as.matrix(as.mcmc(jags_logit))[, c("b[1]", "b[2]", "b[3]")]

#' pred_mat <- plogis(xmat %*% t(beta)) 

#' 

#' # the matrix of predictions has rows matching the number of rows in the data;

#' # the column are the predictions for each of the 2,000 posterior samples

#' nrow(sim_data)

#' dim(pred_mat)

#' 

#' # now we can call mcmcRocPrc; the default method works with the matrix

#' # of predictions and vector of outcomes as input

#' mcmcRocPrc(object = pred_mat, curves = TRUE, fullsims = FALSE, yvec = yvec)

#' 

mcmcRocPrc <- function(object, curves = FALSE, fullsims = FALSE, ...) {

}

#' Constructor for mcmcRocPrc objects

#' 

#' This function actually does the heavy lifting once we have a matrix of 

#' predicted probabilities from a model, plus the vector of observed outcomes.

#' The reason to have it here in a single function is that we don't replicate 

#' it in each function that accomodates a JAGS, BUGS, RStan, etc. object.

#' 

#' @param pred_prob a `\[N, iter\]` matrix of predicted probabilities 

#' @param yvec a `numeric(N)` vector of observed outcomes

#' @param curves include curve data in output?

#' @param fullsims collapse posterior samples into single summary?

#' 

#' @md

#' @keywords internal

new_mcmcRocPrc <- function(pred_prob, yvec, curves, fullsims) {

  )

#' @rdname mcmcRocPrc

#' 

#' @md

#' @export

mcmcRocPrc.default <- function(object, curves, fullsims, yvec, ...) {

  )

}

# Under the hood ROC/PRC calculations -------------------------------------

# auc_roc and auc_pr are not really used, but keep around just in case

auc_roc <- function(obs, pred) {

}

auc_pr <- function(obs, pred) {

}

# JAGS-like input (rjags, R2jags, runjags) --------------------------------

mcmcRocPrc.jags <- function(object, curves = FALSE, fullsims = FALSE, yname, 

                            xnames, posterior_samples, ...) {

  )

  # add intercept by default, maybe revisit this

  # this is not very robust, assumes pars are 'b[x]'

  # for both this and the intercept addition above, maybe a more robust solution

  # down the road would be to dig into the object$model$model() string

  } 

#' @export

mcmcRocPrc.rjags <- function(object, curves = FALSE, fullsims = FALSE, yname, 

                             xnames, ...) {

    stop("Package \"R2jags\" is needed for this function to work. Please install it.", call. = FALSE)  # nocov

  }

  # pass it on to the "jags" method

}

#' @rdname mcmcRocPrc

mcmcRocPrc.runjags <- function(object, curves = FALSE, fullsims = FALSE, yname, 

                               xnames, ...) {

  # as.mcmc.runjags will issue a warning when converting multiple chains

  # because it combines them

  # pass it on to the "jags" method

}

# STAN-like input (rstan, rstanarm, brms) ---------------------------------

#' @rdname mcmcRocPrc

#' 

#' @param data the data that was used in the `stan(data = ?, ...)` call

#' 

#' @export

mcmcRocPrc.stanfit <- function(object, curves = FALSE, fullsims = FALSE, data, 

                               xnames, yname, ...) {

    stop("Package \"rstan\" is needed for this function to work. Please install it.", call. = FALSE)  # nocov

  }

  }

  }

  # add intercept by default, maybe revisit this

  # this is not very robust, assumes pars are 'b[x]'

  } 

}

#' Try to identify if a stanfit model is a binary choice model

#' 

#' @param obj stanfit object

#' 

#' @keywords internal

is_binary_model <- function(obj) {

}

#' Try to identify link function 

#' 

#' @param obj stanfit object

#' 

#' @return Either "logit" or "probit"; if neither can be identified the function

#' will return `NA_character_`. 

#' 

#' @keywords internal

identify_link_function <- function(obj) {

  NA_character_

}

#' @rdname mcmcRocPrc

#' 

#' @export

mcmcRocPrc.stanreg <- function(object, curves = FALSE, fullsims = FALSE, ...) {

    stop("Package \"rstanarm\" is needed for this function to work. Please install it.", call. = FALSE)  # nocov

  }

  # posterior_linepred returns a matrix in which data cases are columns, and 

  # MCMC samples are row; we need to transpose this so that columns are samples

#' @export

mcmcRocPrc.brmsfit <- function(object, curves = FALSE, fullsims = FALSE, ...) {

    stop("Package \"brms\" is needed for this function to work. Please install it.", call. = FALSE)  # nocov

  }

  }

  # posterior_epred returns a matrix in which data cases are columns, and 

  # MCMC samples are row; we need to transpose this so that columns are samples

}

# Other input types (MCMCpack, ...) ---------------------------------------

#' #' @rdname mcmcRocPrc

#' #'

#' #' @export

#' mcmcRocPrc.bugs <- function(object, curves = FALSE, fullsims = FALSE, ...) {

#'   stop("not implemented yet")

#' }

#' @rdname mcmcRocPrc

#' @param type "logit" or "probit"

#' @param force for MCMCpack models, suppress warning if the model does not 

#'   appear to be a binary choice model?

#'

mcmcRocPrc.mcmc <- function(object, curves = FALSE, fullsims = FALSE, data, 

                            xnames, yname, type = c("logit", "probit"), 

                            force = FALSE, ...) {

    } else {

      }

  }

  )

  # add intercept by default, maybe revisit this

#

#   Methods for class "mcmcRocPrc", generated by mcmcRocPrc()

#

#' @rdname mcmcRocPrc

#' 

#' @export

print.mcmcRocPrc <- function(x, ...) {

  } else {

  }

}

#' @rdname mcmcRocPrc

#' 

#' @param n plot method: if `fullsims = TRUE`, how many sample curves to draw?

#' @param alpha plot method: alpha value for plotting sampled curves; between 0 and 1

#' 

#' @export

plot.mcmcRocPrc <- function(x, n = 40, alpha = .5, ...) {

  )

    # use first non-NaN y-value for y[1]

  } else {

      )

    }

    # PRC

    # The elements of prc_dat have different lengths, unlike roc_dat, so we

    # have to do the central curve differently.

    # Instead of combining the list of curve coordinates from each sample into

    # two x and y matrices, we can first make a point cloud with all curve 

    # points from all samples, and then average the y values at all distinct

    # x coordinates. The x-axis plots recall (TPR), which will only have as 

    # many distinct values as there are positives in the data, so this does 

    # not lose any information about the x coordinates. 

    )

                                    # factor implicitly encodes distinct values only,

                                    # since they will get the same labels

      )

    }

  }

}

#' @rdname mcmcRocPrc

#' 

#' @param row.names see [base::as.data.frame()] 

#' @param optional see [base::as.data.frame()]

#' @param what which information to extract and convert to a data frame?

#' 

#' @export

as.data.frame.mcmcRocPrc <- function(x, row.names = NULL, optional = FALSE,

                                     what = c("auc", "roc", "prc"), ...) {

    # all 4 output types have AUC, so this should work across the board

    # if curves was FALSE, there will be no curve data...

    }

    # Otherwise, there will be either one set of coordinates if mcmcmRegPrc()

    # was called with fullsims = FALSE, or else N_sims curve data sets.

    # If the latter, we can return a long data frame with an identifying 

    # "sim" column to delineate the sim sets. To ensure consistency in output,

    # also add this column when fullsims = FALSE.

    # averaged, single coordinate set

    }

    # full sims

    # add a unique ID to each coordinate set

    # combine into long data frame

  } 

  stop("Developer error (I should not be here): please file an issue on GitHub")  # nocov

}