R/plot_deaths_fn.R

In BayesMoFo: Bayesian Mortality Forecasting

#' A function to plot the fitted and forecast number of deaths, accompanied by credible intervals,
#' from posterior samples generated for stochastic mortality models
#'
#' Plot the median fitted and forecast number of deaths, accompanied by credible intervals (user-specified level), using posterior samples stored in "fit_result" object.
#' @param result object of type either "fit_result" or "BayesMoFo".
#' @param expo_forecast An optional 3-dimensional array (of dimensions \eqn{p \times A \times h}) containing exposure data for the forecast period. If not provided, the exposure data from the most recent year will be used for forecasting.
#' @param pred_int A numeric value (between 0 and 1) specifying the credible level of uncertainty bands. Default is \code{pred_int=0.95} (\eqn{95\%} intervals). 
#' @param plot_type A character string (\code{c("age","time")}) to indicate whether to plot by age (default) or by time/year.
#' @param plot_ages A numeric vector specifying which range of ages to plot for visualisation. If not specified, use whatever ages that were used to fit the model (i.e. \code{fit_result$death$ages}). One panel will be constructed per age when \code{plot_type="time"}, with a maximum of nine panels. If exceeded, only the first nine ages will be plotted.
#' @param plot_years A numeric vector specifying which range of years to plot for visualisation. If not specified, use whatever years that were used to fit the model (i.e. \code{fit_result$death$years}). One panel will be constructed per year when \code{plot_type="age"}, with a maximum of nine panels. If exceeded, only the first nine years will be plotted.
#' @param legends A logical value to indicate if legends of the plots should be shown (default) or suppressed (e.g. to aid visibility).
#' @return A plot illustrating the median fitted and forecast number of deaths, accompanied by credible intervals.
#' @keywords graphics visualization plots
#' @concept fitted deaths
#' @concept forecast deaths
#' @concept credible intervals 
#' @importFrom graphics legend lines par points
#' @importFrom stats dbinom dpois quantile sd rbinom rnbinom rpois
#' @export
#' @examples
#' \donttest{
#' #load and prepare 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)
#' 
#' #fit any mortality model
#' runBayesMoFo_result<-runBayesMoFo(death=death,expo=expo,models="APCI",n_iter=1000,forecast=TRUE)
#' 
#' #default plot
#' plot_deaths_fn(runBayesMoFo_result)
#' 
#' #plot by age and changing pre-specified arguments 
#' plot_deaths_fn(runBayesMoFo_result,pred_int=0.8,plot_ages=40:60,plot_years=c(2017,2020))
#' 
#' #plot by time/year
#' plot_deaths_fn(runBayesMoFo_result,plot_type="time",plot_ages=c(40,50,60))
#' }

plot_deaths_fn<-function(result,expo_forecast=NULL,pred_int=0.95,plot_type="age",
                         plot_ages=NULL,plot_years=NULL,legends=TRUE){
  
  deaths_qtl <- predict_deaths_fn(result, expo_forecast, pred_int)
  
  if ("BayesMoFo_obj" %in% names(result)){
    fit_result<-result$result$best
  } else {
    fit_result<-result
  }
  
  death<-fit_result$death
  p<-death$n_strat
  A<-death$n_ages
  T<-death$n_years  
  h<-fit_result$h
  forecast<-fit_result$forecast
  
  p_names<-dimnames(death$data)[[1]]
  ages_names<-dimnames(death$data)[[2]]
  years_names<-dimnames(death$data)[[3]]
  
  if (is.null(plot_ages)) {
    plot_ages = death$ages
  }
  
  if(!forecast) h <- 0
  
  if (is.null(plot_years)) {
    plot_years = c(death$years, death$years[T] + seq_len(h))
  }
  
  oldpar <- par(no.readonly = TRUE) 
  on.exit(par(oldpar)) 
  
  #plot by age
  if (plot_type=="age"){
    
    length_years <- length(plot_years)
    if (length_years <= 3) {
      par(mfrow = c(1, length_years))
    } else if (length_years > 3 & length_years <= 6) {
      par(mfrow = c(2, ceiling(length_years / 2)))
    } else if (length_years > 6 & length_years <= 9) {
      par(mfrow = c(3, 3))
    } else{
      par(mfrow = c(3, 3))
      warning("Too many years selected, only printing the first 9 years.")
      plot_years <- plot_years[1:9]
    }
    
    for (i in 1:length(plot_years)){
      
      yrange_plot <- range(deaths_qtl[,,as.character(plot_ages),as.character(plot_years[i])])
      
      plot(NULL,xlim=range(plot_ages),ylim=yrange_plot,main=plot_years[i],xlab="age",ylab="Predicted deaths")
      for (j in 1:p){
        
        lines(plot_ages,deaths_qtl[2,j,as.character(plot_ages),as.character(plot_years[i])],type="l",col=(j+1));
        lines(plot_ages,deaths_qtl[1,j,as.character(plot_ages),as.character(plot_years[i])],lty=2,col=(j+1));
        lines(plot_ages,deaths_qtl[3,j,as.character(plot_ages),as.character(plot_years[i])],lty=2,col=(j+1))
        
        if (legends){legend("bottomright",p_names,lty=1,col=((1:p)+1))}
        
      }
    }
  }
  
  #plot by time/year
  if (plot_type=="time"){
    
    length_ages<-length(plot_ages)
    if (length_ages<=3) {
      par(mfrow = c(1, length_ages))
    } else if (length_ages > 3 & length_ages <= 6) {
      par(mfrow = c(2, ceiling(length_ages / 2)))
    } else if (length_ages > 6 & length_ages <= 9) {
      par(mfrow = c(3, 3))
    } else{
      par(mfrow = c(3, 3))
      warning("Too many ages selected, only printing the first 9 ages.")
      plot_ages <- plot_ages[1:9]
    }
    
    for (i in 1:length(plot_ages)){
      yrange_plot<-range(deaths_qtl[,,as.character(plot_ages[i]),])
      plot(NULL,xlim=range(plot_years),ylim=yrange_plot,main=paste0("Age=",plot_ages[i]),xlab="year",ylab="Predicted deaths")
      for (j in 1:p){
        lines(plot_years,deaths_qtl[2,j,as.character(plot_ages[i]),as.character(plot_years)],type="l",col=(j+1));
        lines(plot_years,deaths_qtl[1,j,as.character(plot_ages[i]),as.character(plot_years)],lty=2,col=(j+1));
        lines(plot_years,deaths_qtl[3,j,as.character(plot_ages[i]),as.character(plot_years)],lty=2,col=(j+1))
        if (legends){legend("bottomright",p_names,lty=1,col=((1:p)+1))}
      }
    }
  }
  
  invisible(gc())
  
}