R/runBayesMoFo.R

In BayesMoFo: Bayesian Mortality Forecasting

#' A function to fit mortality models.
#'
#' Carry out Bayesian estimation a selection of stochastic mortality models considered in the paper. \cr
#' DIC (Spiegelhalter et al., 2002) is used to perform model selection in determining the best/worst model for the specified (stratified) mortality data.
#'
#' The standard set of models used (25 in total) is as follows: \cr 
#' \code{M1A}, \code{M1M}, \code{M1U}, \code{M2A1}, \code{M2A2}, \code{M2Y1}, \code{M2Y2}, \cr
#' \code{MLiLee}, \code{MLiLee_sharealpha}, \cr
#' \code{CBD_M1}, \cr
#' \code{CBD_M2}, \code{CBD_M2_sharegamma}, \cr 
#' \code{CBD_M3}, \code{CBD_M3_sharegamma}, \cr 
#' \code{CBD_M5}, \cr
#' \code{CBD_M6}, \code{CBD_M6_sharegamma}, \cr
#' \code{CBD_M7}, \code{CBD_M7_sharegamma}, \cr
#' \code{CBD_M8}, \code{CBD_M8_sharegamma}, \cr
#' \code{APCI}, \code{APCI_sharegamma}, \cr
#' \code{PLAT}, \code{PLAT_sharegamma}.
#' 
#' The full list of mortality models fitted (44 in total) is as follows: \cr
#' \code{M1A}, \code{M1M}, \code{M1U}, \code{M2A1}, \code{M2A2}, \code{M2Y1}, \code{M2Y2}, \cr
#' \code{MLiLee}, \code{MLiLee_sharealpha}, \cr
#' \code{CBD_M1}, \code{CBD_M1_sharealpha}, \code{CBD_M1_sharebeta}, \code{CBD_M1_shareall}, \cr
#' \code{CBD_M2}, \code{CBD_M2_sharealpha}, \code{CBD_M2_sharebeta}, \code{CBD_M2_sharegamma}, \code{CBD_M2_sharealpha_sharebeta}, \code{CBD_M2_sharealpha_sharegamma}, \code{CBD_M2_sharebeta_sharegamma}, \code{CBD_M2_shareall}, \cr
#' \code{CBD_M3}, \code{CBD_M3_sharealpha}, \code{CBD_M3_sharegamma}, \code{CBD_M3_shareall}, \cr
#' \code{CBD_M5}, \cr
#' \code{CBD_M6}, \code{CBD_M6_sharegamma}, \cr
#' \code{CBD_M7}, \code{CBD_M7_sharegamma}, \cr
#' \code{CBD_M8}, \code{CBD_M8_sharegamma}, \cr
#' \code{APCI}, \code{APCI_sharealpha}, \code{APCI_sharebeta}, \code{APCI_sharegamma}, \code{APCI_sharealpha_sharebeta}, \code{APCI_sharealpha_sharegamma}, \code{APCI_sharebeta_sharegamma}, \code{APCI_shareall}, \cr
#' \code{PLAT}, \code{PLAT_sharealpha}, \code{PLAT_sharegamma}, \code{PLAT_shareall}.
#'
#' 
#' @references Spiegelhalter, David J., Best, Nicola G., Carlin, Bradley P., and van der Linde, Angelika. (2002). "Bayesian measures of model complexity and fit (with discussion)". Journal of the Royal Statistical Society, Series B. 64 (4): 583–639.\doi{https://doi.org/10.1111/1467-9868.00353}
#' 
#' @param death death data that has been formatted through the function \code{preparedata_fn}.
#' @param expo expo data that has been formatted through the function \code{preparedata_fn}.
#' @param family a string of characters that defines the family function associated with the mortality model. "poisson" would assume that deaths follow a Poisson distribution and use a log family; "binomial" would assume that deaths follow a Binomial distribution and a logit family; "nb" (default) would assume that deaths follow a Negative-Binomial distribution and a log family.
#' @param forecast a logical value indicating if forecast is to be performed (default is \code{FALSE}). See below for details.
#' @param h a numeric value giving the number of years to forecast. Default is \code{h=5}.
#' @param models a vector of strings specifying the models to run. If not specified, a standard set of models is run. If we set \code{models="all"}, all the possible models considered will be run.
#' @param n_iter number of iterations to run. Default is \code{n_iter=1000}.
#' @param n.chain number of parallel chains for the model. Default is \code{n.chain=1}.
#' @param thin thinning interval for monitoring purposes.
#' @param n.adapt the number of iterations for adaptation. See \code{?rjags::adapt} for details.
#' @param quiet if TRUE then messages generated during compilation will be suppressed, as well as the progress bar during adaptation.
#' @return A list with components:
#' \describe{
#'   \item{\code{result}}{A list containing 2 lists, respectively called "best" (\code{$result$best}) and "worst" (\code{$result$best}). Both return the similar output as \code{fit_result}, with the former giving those of the best model while the latter giving those of the worst model.}
#'   \item{\code{DIC}}{A table containing the numeric values of the DIC of all mortality models fitted.}
#'   \item{\code{best_model}}{A character string indicating the best model (lowest DIC).}
#'   \item{\code{worst_model}}{A character string indicating the worst model (highest DIC). If only one model was specified, then this is the same as `best_model`.}
#'   \item{\code{BayesMoFo_obj}}{A logical value indicating whether the result has been generated using the function`runBayesMoFo` (default=TRUE).}
#' }
#' @rdname runBayesMoFo
#' @keywords bayesian estimation model selection
#' @concept DIC
#' @concept stochastic mortality models
#' @importFrom stats dnbinom dbinom dpois quantile sd
#' @import tidyverse
#' @export
#' @examples
#' \donttest{
#' #load and prepare mortality data
#' data("dxt_array_product");data("Ext_array_product")
#' death<-preparedata_fn(dxt_array_product,strat_name = c("ACI","DB","SCI"),ages=35:65)
#' expo<-preparedata_fn(Ext_array_product,strat_name = c("ACI","DB","SCI"),ages=35:65)
#' 
#' ##automatically find the best model from a standard set
#' runBayesMoFo_result<-runBayesMoFo(death=death,expo=expo,n_iter=50,n.adapt=50)
#' head(runBayesMoFo_result$result$best)
#' runBayesMoFo_result$DIC
#' runBayesMoFo_result$best_model
#' 
#' ##if fit all the models
#' runBayesMoFo_result<-runBayesMoFo(death=death,expo=expo,models="all",
#'   n_iter=50,n.adapt=50)
#'
#' # fit a subset of the models and forecast for 10 years 
#' runBayesMoFo_result<-runBayesMoFo(death=death,expo=expo,models=c("APCI","LC","PLAT"),
#'   n_iter=1000,n.adapt=1000,n.chain=2,forecast=TRUE,h=10)
#'
#' ##plot the best model
#' plot_rates_fn(runBayesMoFo_result)
#' plot_rates_fn(runBayesMoFo_result,plot_type="time",plot_ages=c(40,50,60))
#' plot_param_fn(runBayesMoFo_result)
#' 
#' ##convergence diagnostics plots
#' 
#' #trace and density plots of death rates
#' converge_diag_rates_fn(runBayesMoFo_result)
#' 
#' #trace and density plots of parameters
#' converge_diag_param_fn(runBayesMoFo_result)
#' 
#' #ACF plots of death rates
#' converge_diag_rates_fn(runBayesMoFo_result, trace = FALSE, density = FALSE, acf_plot = TRUE)
#' 
#' #ACF plots of parameters
#' converge_diag_param_fn(runBayesMoFo_result, trace = FALSE, density = FALSE, acf_plot = TRUE)
#' 
#' #Some MCMC diagnostics (Gelman, Geweke, Heidel diagnostics etc.)
#' converge_diag_fn(runBayesMoFo_result)
#' }

runBayesMoFo <- function(death,
                         expo,
                         models = NULL,
                         family = "nb",
                         forecast = FALSE,
                         h = 5,
                         n_iter = 1000,
                         n.chain = 1,
                         n.adapt = 1000,
                         thin = 1,
                         quiet = FALSE) {
  
  # convert NA exposures
  expo$data[is.na(expo$data)] <- 100
  
  p <- death$n_strat
  
  if (is.null(models)) {
    if (p == 1) {
      models <- c("LC",
                  "RH",
                  "CBD_M3",
                  "CBD_M5",
                  "CBD_M6",
                  "CBD_M7",
                  "CBD_M8",
                  "APCI",
                  "MLiLee",
                  "PLAT")
    } else {
      models <- c(
        "M1A",
        "M1U",
        "M1M",
        "M2A1",
        "M2A2",
        "M2Y1",
        "M2Y2",
        "MLiLee",
        "MLiLee_sharealpha",
        "LC",
        "RH_sharegamma",
        "RH",
        "CBD_M3_sharegamma",
        "CBD_M3",
        "CBD_M5",
        "CBD_M6_sharegamma",
        "CBD_M6",
        "CBD_M7_sharegamma",
        "CBD_M7",
        "CBD_M8_sharegamma",
        "CBD_M8",
        "APCI_sharegamma",
        "APCI",
        "PLAT_sharegamma",
        "PLAT"
      )
    }
    
  } else if(length(models) == 1 && models == "all"){
    if (p == 1) {
      models <- c("LC",
                  "RH",
                  "CBD_M3",
                  "CBD_M5",
                  "CBD_M6",
                  "CBD_M7",
                  "CBD_M8",
                  "APCI",
                  "MLiLee",
                  "PLAT")
    } else {
      models <- c(
        "M1A",
        "M1U",
        "M1M",
        "M2A1",
        "M2A2",
        "M2Y1",
        "M2Y2",
        "MLiLee",
        "MLiLee_sharealpha",
        "LC",
        "LC_sharealpha",
        "LC_sharebeta",
        "LC_shareall",
        "RH",
        "RH_sharealpha",
        "RH_sharebeta",
        "RH_sharegamma",
        "RH_sharealpha_sharebeta",
        "RH_sharealpha_sharegamma",
        "RH_sharebeta_sharegamma",
        "RH_shareall",
        "CBD_M3",
        "CBD_M3_sharealpha",
        "CBD_M3_sharegamma",
        "CBD_M3_shareall",
        "CBD_M5",
        "CBD_M6",
        "CBD_M6_sharegamma",
        "CBD_M7",
        "CBD_M7_sharegamma",
        "CBD_M8",
        "CBD_M8_sharegamma",
        "APCI",
        "APCI_sharealpha",
        "APCI_sharebeta",
        "APCI_sharegamma",
        "APCI_sharealpha_sharebeta",
        "APCI_sharealpha_sharegamma",
        "APCI_sharebeta_sharegamma",
        "APCI_shareall",
        "PLAT",
        "PLAT_sharealpha",
        "PLAT_sharegamma",
        "PLAT_shareall"
      )
    }
  }
  
  # if (p == 1) {
  
  M <- length(models)
  DIC_table <- matrix(0, nrow = 1, ncol = M)
  colnames(DIC_table) <- models
  fit_results_list <- list()
  
  args <- list(
    "death" = death,
    "expo" = expo,
    "n_iter" = n_iter,
    "family" = family,
    "n.chain" = n.chain,
    "thin" = thin,
    "n.adapt" = n.adapt,
    "forecast" = forecast,
    "h" = h,
    "quiet" = quiet
  )
  
  iter = 1
  for (idx_model in seq_along(models)) {
    
    modelName <- models[idx_model]
    
    # currentFun <- paste0("fit_", modelName)
    if(p == 1){
      currentFun <- modelsName_ap[modelsName_ap[,1] == modelName,2]  
    } else {
      currentFun <- modelsName_app[modelsName_app[,1] == modelName,2]  
    }
    
    fit_fun <- get(currentFun, envir = asNamespace("BayesMoFo"))
    arg_names <- names(formals(fit_fun))
    args_filtered <- args[names(args) %in% arg_names]
    fit_result <- do.call(fit_fun, args=args_filtered)
    
    fit_results_list[[modelName]] <- fit_result
    
    DIC_table[1, modelName] <- DIC_fn(fit_result)
    insight::print_colour(paste0("Completed: ", modelName, " (", iter, "/", M, ")", "\n"),
                          "red")
    iter <- iter + 1
  }
  
  # gc()
  
  best_model <- colnames(DIC_table)[which.min(DIC_table)]
  worst_model <- colnames(DIC_table)[which.max(DIC_table)]
  
  result <- list(best = fit_results_list[[best_model]], 
                 worst = fit_results_list[[worst_model]])
  
  list(
    result = result,
    DIC = DIC_table,
    best_model = best_model,
    worst_model = worst_model,
    BayesMoFo_obj=TRUE
  )
  
}


modelsName_ap <- rbind(c("LC","fit_LC"),
                       c("RH","fit_RH"),
                       c("CBD_M3","fit_CBD_M3"),
                       c("CBD_M5","fit_CBD_M5"),
                       c("CBD_M6","fit_CBD_M6"),
                       c("CBD_M7","fit_CBD_M7"),
                       c("CBD_M8","fit_CBD_M8_try"),
                       c("APCI","fit_APCI"),
                       c("MLiLee","fit_MLiLee"),
                       c("PLAT","fit_PLAT")
)

modelsName_app <- rbind(
  c("M1A","fit_M1A"),
  c("M1U","fit_M1U"),
  c("M1M","fit_M1M"),
  c("M2A1","fit_M2A1"),
  c("M2A2","fit_M2A2"),
  c("M2Y1","fit_M2Y1"),
  c("M2Y2","fit_M2Y2"),
  c("MLiLee","fit_MLiLee"),
  c("MLiLee_sharealpha","fit_MLiLee_sharealpha"),
  c("LC","fit_LC"),
  c("LC_sharealpha","fit_LC_sharealpha"),
  c("LC_sharebeta","fit_LC_sharebeta"),
  c("LC_shareall","fit_LC_shareall"),
  c("RH","fit_RH"),
  c("RH_sharealpha","fit_RH_sharealpha"),
  c("RH_sharebeta","fit_RH_sharebeta"),
  c("RH_sharegamma","fit_RH_sharegamma"),
  c("RH_sharealpha_sharebeta","fit_RH_sharealpha_sharebeta"),
  c("RH_sharealpha_sharegamma","fit_RH_sharealpha_sharegamma"),
  c("RH_sharebeta_sharegamma","fit_RH_sharebeta_sharegamma"),
  c("RH_shareall","fit_RH_shareall"),
  
  c("CBD_M3","fit_CBD_M3"),
  c("CBD_M3_sharealpha","fit_CBD_M3_sharealpha"),
  c("CBD_M3_sharegamma","fit_CBD_M3_sharegamma"),
  c("CBD_M3_shareall","fit_CBD_M3_shareall"),
  
  c("CBD_M5","fit_CBD_M5"),
  
  c("CBD_M6","fit_CBD_M6"),
  c("CBD_M6_sharegamma","fit_CBD_M6_sharegamma"),
  
  c("CBD_M7","fit_CBD_M7"),
  c("CBD_M7_sharegamma","fit_CBD_M7_sharegamma"),
  
  c("CBD_M8","fit_CBD_M8_try"),
  c("CBD_M8_sharegamma","fit_CBD_M8_sharegamma"),
  
  c("APCI","fit_APCI"),
  c("APCI_sharealpha","fit_APCI_sharealpha"),
  c("APCI_sharebeta","fit_APCI_sharebeta"),
  c("APCI_sharegamma","fit_APCI_sharegamma"),
  c("APCI_sharealpha_sharebeta","fit_APCI_sharealpha_sharebeta"),
  c("APCI_sharealpha_sharegamma","fit_APCI_sharealpha_sharegamma"),
  c("APCI_sharebeta_sharegamma","fit_APCI_sharebeta_sharegamma"),
  c("APCI_shareall","fit_APCI_shareall"),
  
  c("PLAT","fit_PLAT"),
  c("PLAT_sharealpha","fit_PLAT_sharealpha"),
  c("PLAT_sharegamma","fit_PLAT_sharegamma"),
  c("PLAT_shareall","fit_PLAT_shareall")
  
)