R/ppc.R

In BeeGUTS: General Unified Threshold Model of Survival for Bees using Bayesian Inference

#' Generates an object to be used in posterior predictive check for \code{beeSurvFit}, \code{beeSurvPred}
#'
#' @param x  an object used to select a method \code{ppc}
#'
#' @return a \code{data.frame} of class \code{ppc}
#' @export
#'
ppc <- function(x){
  UseMethod("ppc")
}



#' Posterior predictive check method for \code{beeSurvFit} objects
#'
#' @param x an object of class \code{beeSurvFit}
#'
#'
#' @return a \code{data.frame} of class \code{ppc}
#'
#' @examples
#' data(fitBetacyfluthrin_Chronic)
#' out <- ppc(fitBetacyfluthrin_Chronic)
#'
#' @export
#'
ppc.beeSurvFit <- function(x){
  NsurvPred_all<- as.data.frame(x$stanFit, pars = "Nsurv_ppc")
  NsurvPred_quantiles<- NsurvPred_all%>%
    tidyr::pivot_longer(cols = tidyr::starts_with('Nsurv'),
                        names_to = "ppc",
                        values_to = "value")%>%
    dplyr::group_by(ppc)%>%
    dplyr::summarise(median = stats::quantile(value,  0.5, na.rm = TRUE),
                     q_0.025=stats::quantile(value,  0.025, na.rm = TRUE),
                     q_0.975=stats::quantile(value,  0.975, na.rm = TRUE))

  NsurvData_all<- data.frame(value=x$dataFit$Nsurv, id=seq(1,x$dataFit$nData_Nsurv, 1))%>%
    dplyr::mutate(ppc=paste0("Nsurv_ppc[",id, "]"))

  Nsurv_ppc<- dplyr::full_join( NsurvPred_quantiles,  NsurvData_all, by="ppc")%>%
    dplyr::mutate(col=ifelse(value<q_0.025|value>q_0.975, "red", "green")) %>%
    dplyr::arrange(id)

  Nsurv_ppc$data<-"Survival"

  class(Nsurv_ppc) <- c("ppc", class(Nsurv_ppc))
  return(Nsurv_ppc)
}


#' Posterior predictive check method for \code{beeSurvValidation} objects
#'
#' @param x an object of class \code{beeSurvValidation}
#'
#'
#' @return a \code{data.frame} of class \code{ppc}
#'
#' @examples
#' data(fitBetacyfluthrin_Chronic)
#' data(betacyfluthrinChronic)
#' valid <- validate(fitBetacyfluthrin_Chronic,betacyfluthrinChronic)
#' out <- ppc(valid)
#'
#' @export
#'
ppc.beeSurvValidation <- function(x){
  NSurv_ppc = data.frame(median = x$sim$Nsurv_q50_valid,
                         q_0.025= x$sim$Nsurv_qinf95_valid,
                         q_0.975= x$sim$Nsurv_qsup95_valid,
                         value  = x$sim$Nsurv,
                         data   = rep("Survival",length(x$sim$Nsurv_q50_valid)))

  NSurv_ppc = NSurv_ppc %>% dplyr::mutate(col = ifelse(value<q_0.025|value>q_0.975, "red", "green"))
  class(NSurv_ppc) <- c("ppc", class(NSurv_ppc))
  return(NSurv_ppc)
}