R/predict_deaths_fn.R

In BayesMoFo: Bayesian Mortality Forecasting

#' A function to compute the fitted and forecast number of deaths, accompanied by credible intervals, 
#' from posterior samples generated for stochastic mortality models
#'
#' Return 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). 
#' @return An array containing the lower, median, and upper quantiles of the number of deaths for both the fitted and forecast periods.
#' @keywords models forecasting
#' @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)
#' 
#' #default
#' predict_deaths_fn(runBayesMoFo_result)
#' 
#' #changing pre-specified arguments 
#' predict_deaths_fn(runBayesMoFo_result,pred_int=0.8)
#' }

predict_deaths_fn<-function(result,expo_forecast=NULL,pred_int=0.95){
  
  if ("BayesMoFo_obj" %in% names(result)){
    fit_result<-result$result$best
  } else {
    fit_result<-result
  }
  
  death<-fit_result$death
  expo_data<-fit_result$expo$data
  if (fit_result$family=="binomial"){expo_data<-round(expo_data+0.5*death$data)}
  
  mcmc_object<-coda::as.mcmc(as.matrix(fit_result$post_sample))
  
  probs<-c(0.5*(1-pred_int),0.5,1-0.5*(1-pred_int))
  
  p<-death$n_strat
  A<-death$n_ages
  T<-death$n_years  
  h<-fit_result$h
  
  forecast<-fit_result$forecast
  
  n<-dim(mcmc_object)[1]
  
  p_names<-dimnames(death$data)[[1]]
  ages_names<-dimnames(death$data)[[2]]
  years_names<-dimnames(death$data)[[3]]
  
  crude_rates<-death$data/expo_data
  crude_rates<-provideDimnames(crude_rates,base=list(p_names,ages_names,years_names))
  
  if(!forecast){
    h <- 0
  }
  
  rates_mat<-mcmc_object[,startsWith(colnames(mcmc_object),"q")]
  
  familyResults <- fit_result$family
  
  niter <- dim(rates_mat)[1]
  deaths_pred <- array(dim=c(p,A,T+h,niter))
  deaths_pred<-provideDimnames(deaths_pred,base=list(p_names,ages_names,c(years_names,as.character(death$years[T]+seq_len(h)))))
  
  if (forecast){
    if (is.null(expo_forecast)){
      expo_forecast<-array(dim=c(p,A,h))
      for (m in 1:h){
        expo_forecast[,,m]<-expo_data[,,T]
      }
      expo_data_temp<-array(dim=c(p,A,T+h))
      expo_data_temp[,,1:T]<-expo_data
      expo_data_temp[,,(T+1):(T+h)]<-expo_forecast
      expo_data<-expo_data_temp
      rm(expo_data_temp)
      invisible(gc())
    } else {
      if (!all(dim(expo_forecast)==c(p,A,h))){stop("Please provide exposure data with the correct dimensions.")}
      expo_data_temp<-array(dim=c(p,A,T+h))
      expo_data_temp[,,1:T]<-expo_data
      expo_data_temp[,,(T+1):(T+h)]<-expo_forecast
      expo_data<-expo_data_temp
      rm(expo_data_temp)
      invisible(gc())
    }
  }
  
  for (i in 1:p){
    for(j in 1:A){
      for(k in 1:(T+h)){
        
        index<-paste0("q[",i,",",j,",",k,"]")
        
        if(familyResults == "poisson"){
          
          deaths_pred[i,j,k,] <- sapply(1:niter, function(iter){
            rpois(1, rates_mat[iter,index] * expo_data[i,j,k])
          })
          
        } else if(familyResults == "binomial"){
          
          deaths_pred[i,j,k,] <- sapply(1:niter, function(iter){
            rbinom(1, expo_data[i,j,k], rates_mat[iter,index])
          })
          
        } else if(familyResults == "nb"){
          
          phi_mat <- mcmc_object[,startsWith(colnames(mcmc_object),"phi")]
          
          deaths_pred[i,j,k,] <- sapply(1:niter, function(iter){
            phi_current <- phi_mat[iter]
            prob_current <-
              phi_current / (phi_current + rates_mat[iter,index] * expo_data[i,j,k])
            rnbinom(1, size = phi_current, 
                    prob =  prob_current)
          })
          
        } 
      }
    }
  }
  
  deaths_qtl<-apply(deaths_pred,c(1,2,3),quantile_fn,probs=probs)

  return(deaths_qtl)
  
}