R/sf_output_mcmc_summaries.R

In AurMad/STOCfree: STOC free model: prediction of probability of freedom from longitudinal data

#' Summary quantiles for the STOC free model parameters
#'
#' @param x a data frame with MCMC samples for the different model parameters generated by the STOC free model
#' @param parameters model parameters to summarise
#' @param probs quantiles to compute for each parameter. The default values are the median as well as the 2.5 and 97.5 percentiles.
#' @param digits number of digits used to round the results
#'
#' @return this function returns the quantiles supplied with the probs argument for the different parameters of the STOC free model.
#' @export
param_summary <- function(x = data.frame(),
                          parameters = c("Se", "Sp", "tau", "pi_within", "theta"),
                          probs = c(.5, .025, .975),
                          digits = 4){

  ## checks if the object contains MCMC samples
  x_type_check <- unlist(sapply(c(".chain", ".iteration", ".draw"), grep, colnames(x)))
  if(length(x_type_check) < 3) stop("Wrong type of input")

  ## column in which each parameter is stored
  param_column <- unlist(sapply(parameters, grep, colnames(x)))

  ## creation of the output object
  param_summary <- data.frame(
    parameter = colnames(x)[param_column]
  )

  param_summary <- cbind(param_summary,
                         as.data.frame(matrix(rep(NA, nrow(param_summary) * length(probs)), ncol = length(probs)))
  )

  colnames(param_summary)[2:length(param_summary)] <- paste(probs * 100, "%")

  ## computing the summary values for each parameter
  for(i in 1:nrow(param_summary)){

    param_i <- as.character(param_summary$parameter[i])

    val <- round(quantile(unlist(x[, param_i]), probs), digits)

    param_summary[i, 2:length(param_summary)] <- val

  }

  param_summary

}


#' Summary of herd level probabilities of infection generated by the STOC free model
#'
#' @param x a data frame with MCMC samples for the herd level probabilities of being latent status positive
#' @param quantile use to summarise each herd's MCMC samples
#' @param cut_off cut-off value used to categorise herds as negative or positive
#'
#' @return for each herd, the quantile supplied is computed from the MCMC samples for the predicted probabilities of infection. So, if a value of 0.5 is supplied for quantile, the median probability of infection is computed. The cut_off argument is then used to categorise each herd as status positive or status negative.
#' @export
#'
herd_summary <- function(x = data.frame(),
                         quantile = .5,
                         cut_off = .5){

  ## checks if the object contains MCMC samples with herd-level predicted proba
  x_type_check <- unlist(sapply(c(".chain", ".iteration", ".draw", "predicted_proba"), grep, colnames(x)))
  if(length(x_type_check) < 4) stop("Wrong type of input")

  herd_col <- colnames(x)[1]
  colnames(x)[1] <- "herd_id"
  
  ## computing summary value for each herd
  hrd_smmry <- dplyr::group_by(x, herd_id)
  hrd_smmry <- dplyr::summarise(hrd_smmry,
                                proba = quantile(predicted_proba, quantile), 
                                .groups = 'drop')

  colnames(hrd_smmry)[match("herd_id", colnames(hrd_smmry))] <- herd_col

  hrd_smmry$predicted_status <- ifelse(hrd_smmry$proba > cut_off, 1, 0)

  hrd_smmry

}