R/DIC_fn.R

In BayesMoFo: Bayesian Mortality Forecasting

#' A function to calculate the DIC of stochastic mortality models using posterior samples
#'
#' Compute the Deviance Information Criterion (DIC) of stochastic mortality models using posterior samples stored in "fit_result" object.
#' @param result object of type either "fit_result" or "BayesMoFo".
#' @return A numeric value representing the DIC of the mortality model.
#' @keywords bayesian model selection
#' @concept DIC
#' @importFrom stats dbinom dpois quantile sd
#' @export
#' @examples
#' #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)
#' 
#' #a toy example
#' runBayesMoFo_result<-runBayesMoFo(death=death,expo=expo,models="M1A",n_iter=50,n.adapt=50)
#' 
#' #compute the DIC
#' DIC_fn(runBayesMoFo_result)

DIC_fn<-function(result){
  
  if ("BayesMoFo_obj" %in% names(result)){
    fit_result<-result$result$best
  } else {
    fit_result<-result
  }
  
  death<-fit_result$death
  expo<-fit_result$expo
  mcmc_object<-coda::as.mcmc(as.matrix(fit_result$post_sample))
  
  family<-fit_result$family
  
  p<-death$n_strat
  A<-death$n_ages
  T<-death$n_years  
  
  deaths_combined_vectorised<-as.vector(death$data)
  expo_combined_vectorised<-as.vector(expo$data)
  expo_initial_combined_vectorised<-as.vector(round(expo$data+0.5*death$data))
  
  forecast<-fit_result$forecast
  
  n<-dim(mcmc_object)[1] #number of posterior samples
  rates_mat<-matrix(0,nrow=n,ncol=A*T)
  
  deviance_vec<-vector(length=n)
  deviance_mean<-NULL
  
  rates_mat<-mcmc_object[,startsWith(colnames(mcmc_object),"q")]
  if (forecast){rates_mat<-rates_mat[,1:(A*T)]}
  rates_mean<-apply(rates_mat,2,mean)
  invisible(gc())
  
  for (s in 1:n){
    if (family=="binomial"){
      deviance_vec[s]<--2*sum(dbinom(x=deaths_combined_vectorised,size=expo_initial_combined_vectorised,prob=rates_mat[s,],log=TRUE))
    }
    if (family=="poisson"){
      deviance_vec[s]<--2*sum(dpois(x=deaths_combined_vectorised,lambda=(expo_combined_vectorised*rates_mat[s,]),log=TRUE))
    }
    if (family=="nb"){
      phi<-mcmc_object[,"phi"]
      deviance_vec[s]<--2*sum(dnbinom(x=deaths_combined_vectorised,size=phi[s],prob=(phi[s]/(phi[s]+expo_combined_vectorised*rates_mat[s,])),log=TRUE))
    }
  }
  
  if (family=="binomial"){
    deviance_mean<--2*sum(dbinom(x=deaths_combined_vectorised,size=expo_initial_combined_vectorised,prob=rates_mean,log=TRUE))
  }
  if (family=="poisson"){
    deviance_mean<--2*sum(dpois(x=deaths_combined_vectorised,lambda=(expo_combined_vectorised*rates_mean),log=TRUE))
  }
  if (family=="nb"){
    phi_mean<-mean(mcmc_object[,"phi"])
    deviance_mean<--2*sum(dnbinom(x=deaths_combined_vectorised,size=phi_mean,prob=(phi_mean/(phi_mean+expo_combined_vectorised*rates_mean)),log=TRUE))
  }
  
  invisible(gc())
  
  DIC<-2*mean(deviance_vec)-deviance_mean
  DIC
}