R/plot_param_fn.R

In BayesMoFo: Bayesian Mortality Forecasting

#' A function to plot the fitted parameters of stochastic mortality models, accompanied by credible intervals 
#'
#' Plot the fitted parameters, 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 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 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 parameters, accompanied by credible intervals.
#' @keywords graphics visualization plots
#' @concept fitted parameters
#' @concept forecast parameters
#' @concept credible intervals
#' @importFrom graphics legend lines par points boxplot
#' @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",
#'   family="poisson",forecast=TRUE)
#' 
#' #default plot
#' plot_param_fn(runBayesMoFo_result)
#' 
#' #with 80% credible intervals 
#' plot_param_fn(runBayesMoFo_result,pred_int=0.8)
#' }

plot_param_fn<-function(result,pred_int=0.95,legends=TRUE){
  
  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))
  param<-fit_result$param
  
  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]
  
  ages<-death$ages
  if (forecast){years<-c(death$years,death$years[T]+(1:h))}else{
  years<-death$years}
  if (forecast){
    cohorts<-(1:(C+h))+years[1]-ages[A]-1
  }else{
    cohorts<-(1:C)+years[1]-ages[A]-1
  }
  
  p_names<-dimnames(death$data)[[1]]
  ages_names<-dimnames(death$data)[[2]]
  if (forecast){
    years_names<-c(dimnames(death$data)[[3]],as.character(death$years[T]+(1:h)))
  } else {
  years_names<-dimnames(death$data)[[3]]}
  cohorts_names<-as.character(cohorts)
  
  k<-length(param)
  
  oldpar <- par(no.readonly = TRUE) 
  on.exit(par(oldpar)) 
  
  if (k<=3){
    par(mfrow=c(1,k))}else if(k>3 & k<=6){
      par(mfrow=c(2,ceiling(k/2)))
    }else if(k>6 & k<=9){
      par(mfrow=c(3,3))
    }else{
      par(mfrow=c(3,3))
    }
  
  plot_dummy<-0
  
  for (l in 1:k){
    
    if (((plot_dummy %% 9)==0) & plot_dummy!=0 && interactive()){readline(prompt = "Press [Enter] to see the next plot...")}
    
    param_l_names<-colnames(mcmc_object)[startsWith(colnames(mcmc_object),param[l])]
    
    param_l_names_comma<-grep(",",param_l_names)
    
    if (length(param_l_names_comma)==0){
      temp<-mcmc_object[,(startsWith(colnames(mcmc_object),param[l]))]
      if (!is.null(dim(temp))){
      assign(paste0(param[l],"_pn"),apply(temp,2,quantile_fn,probs=probs))}
      
      if (startsWith(param[l],"alpha") | startsWith(param[l],"beta") | startsWith(param[l],"Beta")){
        xlab<-"ages"
        x<-ages
      } else if(startsWith(param[l],"gamma")){
        xlab<-"cohorts"
        x<-cohorts
      } else if(startsWith(param[l],"kappa") | startsWith(param[l],"Kappa") ){
        xlab<-"years"
        x<-years
      } else{
        if (is.null(dim(temp))){
          xlab<-""
          x<-1
        } else {
          xlab<-""
          x<-1:(dim(temp)[2])}
      }
      
      if (startsWith(param[l],"alpha") | startsWith(param[l],"beta") | startsWith(param[l],"Beta") | startsWith(param[l],"gamma") | startsWith(param[l],"kappa") | startsWith(param[l],"Kappa")){
        plot(x,get(paste0(param[l],"_pn"))[2,],type="l",xlab=xlab,cex=0.5,main=param[l],ylim=range(temp),ylab="")
        lines(x,get(paste0(param[l],"_pn"))[1,],xlab=xlab,lty=2)
        lines(x,get(paste0(param[l],"_pn"))[3,],xlab=xlab,lty=2)
        if (forecast & (startsWith(param[l],"kappa") | startsWith(param[l],"Kappa"))){abline(v=death$years[T],lty=3)}
        if (forecast & startsWith(param[l],"gamma")){abline(v=cohorts[C],lty=3)}
        if ((startsWith(param[l],"kappa") | startsWith(param[l],"Kappa"))){
          if (legends){legend("bottomleft",c("shared"),lty=1,col=(1))}
        } else {
          if (legends){legend("bottomright",c("shared"),lty=1,col=(1))}}}
      else {
        boxplot(as.matrix(temp),xlab=xlab,main=param[l])
      }
    }
    
    if (length(param_l_names_comma)>0){
      temp<-mcmc_object[,(startsWith(colnames(mcmc_object),param[l]))]
      
      if (startsWith(param[l],"alpha") | startsWith(param[l],"beta") | startsWith(param[l],"Beta")){
        pp<-dim(temp)[2]/A
        xlab<-"ages"
        x<-ages
      } else if(startsWith(param[l],"gamma")){
        if (forecast){
          pp<-dim(temp)[2]/(C+h)
        } else {
          pp<-dim(temp)[2]/C
        }
        xlab<-"cohorts"
        x<-cohorts
      } else if(startsWith(param[l],"kappa") | startsWith(param[l],"Kappa") ){
        if (forecast){
          pp<-dim(temp)[2]/(T+h)
        } else {
          pp<-dim(temp)[2]/T
        }
        xlab<-"years"
        x<-years
      } else{
        if (is.null(dim(temp))){
          xlab=""
          x<-1
        } else {
          xlab<-""
          x<-1:(dim(temp)[2])}
      }
      
      for (i in 1:pp){
        assign(paste0(param[l],i,"_pn"),apply(temp[,seq(i,dim(temp)[2],by=pp)],2,quantile_fn,probs=probs))
      }
      
      color<-(1:p)+1
      if (pp==(p-1)){color<-color[-p]}
      
      plot(x,get(paste0(param[l],1,"_pn"))[2,],xlab=xlab,type="l",cex=0.5,main=param[l],ylim=range(temp),ylab="",col=color[1])
      lines(x,get(paste0(param[l],1,"_pn"))[1,],xlab=xlab,lty=2,col=color[1])
      lines(x,get(paste0(param[l],1,"_pn"))[3,],xlab=xlab,lty=2,col=color[1])
      
      if (forecast & (startsWith(param[l],"kappa") | startsWith(param[l],"Kappa"))){abline(v=death$years[T],lty=3)}
      if (forecast & startsWith(param[l],"gamma")){abline(v=cohorts[C],lty=3)}
      if ((startsWith(param[l],"kappa") | startsWith(param[l],"Kappa"))){
        if (legends){legend("bottomleft",c(p_names),lty=1,col=((1:p)+1))}
      } else {
        if (legends){legend("bottomright",c(p_names),lty=1,col=((1:p)+1))}
      }
      
      if (pp>1){for (j in 2:pp){
        lines(x,get(paste0(param[l],j,"_pn"))[2,],xlab=xlab,col=color[j])
        lines(x,get(paste0(param[l],j,"_pn"))[1,],xlab=xlab,lty=2,col=color[j])
        lines(x,get(paste0(param[l],j,"_pn"))[3,],xlab=xlab,lty=2,col=color[j])
      }}
      
    }
    
    plot_dummy<-plot_dummy+1
    
  }
  
  invisible(gc())
}