R/summary_fn.R

In BayesMoFo: Bayesian Mortality Forecasting

#' A function to summarise the fitted mortality rates and parameters of stochastic mortality models 
#'
#' Provide summaries (means, standard errors, percentiles) of the fitted mortality rates and parameters, derived using posterior samples stored in "fit_result" object.
#' @param result object of type either "fit_result"  or "BayesMoFo".
#' @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 A list with components:
#' rates_summary=list(mean=rates_mean,std=rates_std),rates_pn=list(lower=rates_lower,median=rates_median,upper=rates_upper),param_summary=list(mean=param_mean,std=param_std),param_pn=param_pn
#' \describe{
#'   \item{\code{rates_summary}}{A list containing 2 components, respectively called "mean" (\code{$rates_summary$mean}) and "std" (\code{$rates_summary$std}). Both return a 3-dimensional data array (dim 1: strata, dim 2: ages, dim 3: years), with the former giving posterior means of fitted mortality rates while the latter giving standard errors.}
#'   \item{\code{rates_pn}}{A list containing 3 components, respectively called "lower" (\code{$rates_pn$lower}), "median" (\code{$rates_pn$median}), and "upper" (\code{$rates_pn$upper}). All return a 3-dimensional data array (dim 1: strata, dim 2: ages, dim 3: years), representing the respective percentiles for the fitted mortality rates.}
#'   \item{\code{param_summary}}{A list containing 2 components, respectively called "mean" (\code{$param_summary$mean}) and "std" (\code{$param_summary$std}). Both return a 3-dimensional data array (dim 1: strata, dim 2: ages, dim 3: years), with the former giving posterior means of fitted parameters while the latter giving standard errors.}
#'   \item{\code{param_pn}}{A 2-dimensional matrix containing percentiles of fitted parameters.}
#' }
#' @keywords statistics
#' @concept death rates
#' @concept parameters
#' @concept credible intervals
#' @concept summary 
#' @importFrom graphics legend lines par points
#' @importFrom stats dbinom dpois quantile sd
#' @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,n_iter=1000,models="APCI")
#' 
#' #default summary
#' summary_runBayesMoFo<-summary_fn(runBayesMoFo_result)
#' 
#' #mean of fitted mortality rates 
#' summary_runBayesMoFo$rates_summary$mean
#' 
#' #standard errors of fitted mortality rates 
#' summary_runBayesMoFo$rates_summary$std
#' 
#' #97.5th percentile of fitted mortality rates 
#' summary_runBayesMoFo$rates_pn$upper
#' 
#' #mean of fitted parameters 
#' summary_runBayesMoFo$param_summary$mean
#' 
#' #standard errors of fitted parameters 
#' summary_runBayesMoFo$param_summary$std
#' 
#' #97.5th percentile of fitted parameters 
#' summary_runBayesMoFo$param_pn[,"upper"]
#' }

summary_fn<-function(result,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  
  C<-A+T-1
  
  h<-fit_result$h
  forecast<-fit_result$forecast
  
  n<-dim(mcmc_object)[1]
  
  if (!forecast){
  ages<-death$ages
  years<-death$years
  cohorts<-(1:C)+years[1]-ages[A]-1
  
  p_names<-dimnames(death$data)[[1]]
  ages_names<-dimnames(death$data)[[2]]
  years_names<-dimnames(death$data)[[3]]
  cohorts_names<-as.character(cohorts)
  
  param<-fit_result$param
  k<-length(param)
  
  crude_rates<-death$data/expo_data
  crude_rates<-provideDimnames(crude_rates,base=list(p_names,ages_names,years_names))
  
  rates_lower<-array(dim=c(p,A,T))
  rates_median<-array(dim=c(p,A,T))
  rates_upper<-array(dim=c(p,A,T))
  rates_mean<-array(dim=c(p,A,T))
  rates_std<-array(dim=c(p,A,T))
  rates_lower<-provideDimnames(rates_lower,base=list(p_names,ages_names,years_names))
  rates_median<-provideDimnames(rates_median,base=list(p_names,ages_names,years_names))
  rates_upper<-provideDimnames(rates_upper,base=list(p_names,ages_names,years_names))
  rates_mean<-provideDimnames(rates_mean,base=list(p_names,ages_names,years_names))
  rates_std<-provideDimnames(rates_std,base=list(p_names,ages_names,years_names))
  interval_name<-c("lower","median","upper")
  
  rates_mat<-matrix(0,nrow=n,ncol=A*T)
  
  rates_mat<-mcmc_object[,startsWith(colnames(mcmc_object),"q")]
  
  rates_pn<-apply(rates_mat,2,quantile_fn,probs=probs)
  rates_mean_vec<-apply(rates_mat,2,mean)
  rates_std_vec<-apply(rates_mat,2,sd)
  
  for (i in 1:p){
    for(j in 1:A){
      for(k in 1:T){
        index<-paste0("q[",i,",",j,",",k,"]")
        rates_lower[i,j,k]<-rates_pn[1,index]
        rates_median[i,j,k]<-rates_pn[2,index]
        rates_upper[i,j,k]<-rates_pn[3,index]
        rates_mean[i,j,k]<-rates_mean_vec[index]
        rates_std[i,j,k]<-rates_std_vec[index]
      }
    }
  }
  } else {
    
    ages<-death$ages
    years<-c(death$years,death$years[T]+(1:h))
    cohorts<-(1:C)+years[1]-ages[A]-1
    
    p_names<-dimnames(death$data)[[1]]
    ages_names<-dimnames(death$data)[[2]]
    years_names<-c(dimnames(death$data)[[3]],as.character(death$years[T]+(1:h)))
    cohorts_names<-as.character(cohorts)
    
    param<-fit_result$param
    k<-length(param)
    
    crude_rates<-death$data/expo_data
    crude_rates<-provideDimnames(crude_rates,base=list(p_names,ages_names,years_names))
    
    rates_lower<-array(dim=c(p,A,T+h))
    rates_median<-array(dim=c(p,A,T+h))
    rates_upper<-array(dim=c(p,A,T+h))
    rates_mean<-array(dim=c(p,A,T+h))
    rates_std<-array(dim=c(p,A,T+h))
    rates_lower<-provideDimnames(rates_lower,base=list(p_names,ages_names,years_names))
    rates_median<-provideDimnames(rates_median,base=list(p_names,ages_names,years_names))
    rates_upper<-provideDimnames(rates_upper,base=list(p_names,ages_names,years_names))
    rates_mean<-provideDimnames(rates_mean,base=list(p_names,ages_names,years_names))
    rates_std<-provideDimnames(rates_std,base=list(p_names,ages_names,years_names))
    interval_name<-c("lower","median","upper")
    
    rates_mat<-matrix(0,nrow=n,ncol=A*(T+h))
    
    rates_mat<-mcmc_object[,startsWith(colnames(mcmc_object),"q")]
    
    rates_pn<-apply(rates_mat,2,quantile_fn,probs=probs)
    rates_mean_vec<-apply(rates_mat,2,mean)
    rates_std_vec<-apply(rates_mat,2,sd)
    
    for (i in 1:p){
      for(j in 1:A){
        for(k in 1:(T+h)){
          index<-paste0("q[",i,",",j,",",k,"]")
          rates_lower[i,j,k]<-rates_pn[1,index]
          rates_median[i,j,k]<-rates_pn[2,index]
          rates_upper[i,j,k]<-rates_pn[3,index]
          rates_mean[i,j,k]<-rates_mean_vec[index]
          rates_std[i,j,k]<-rates_std_vec[index]
        }
      }
    }
  }
  
  param_samples<-mcmc_object[,!startsWith(colnames(mcmc_object),"q")]
  param_pn<-apply(param_samples,2,quantile_fn,probs=probs)
  rownames(param_pn)<-c("lower","median","upper")
  param_pn<-t(param_pn)
  param_mean<-apply(param_samples,2,mean)
  param_std<-apply(param_samples,2,sd)
  
  list(rates_summary=list(mean=rates_mean,std=rates_std),rates_pn=list(lower=rates_lower,median=rates_median,upper=rates_upper),param_summary=list(mean=param_mean,std=param_std),param_pn=param_pn)
  
}