R/fit_LC.R

In BayesMoFo: Bayesian Mortality Forecasting

#' A function to fit the stochastic mortality model by Lee and Carter (1992).
#'
#' Carry out Bayesian estimation of the stochastic mortality model, the \bold{Lee-Carter model} (Lee and Carter, 1992).
#' Note that this model is equivalent to "M1" under the formulation of Cairns et al. (2009).
#' 
#' The model can be described mathematically as follows:
#' If \code{family="poisson"}, then
#' \deqn{d_{x,t,p} \sim \text{Poisson}(E^c_{x,t,p} m_{x,t,p}) , }
#' \deqn{\log(m_{x,t,p})=a_{x,p}+b_{x,p}k_{t,p} , }
#' where \eqn{d_{x,t,p}} represents the number of deaths at age \eqn{x} in year \eqn{t} of stratum \eqn{p},
#' while \eqn{E^c_{x,t,p}} and \eqn{m_{x,t,p}} represents respectively the corresponding central exposed to risk and central mortality rate at age \eqn{x} in year \eqn{t} of stratum \eqn{p}.
#' Similarly, if \code{family="nb"}, then a negative binomial distribution is fitted, i.e.
#' \deqn{d_{x,t,p} \sim \text{Negative-Binomial}(\phi,\frac{\phi}{\phi+E^c_{x,t,p} m_{x,t,p}}) , }
#' \deqn{\log(m_{x,t,p})=a_{x,p}+b_{x,p}k_{t,p} , } 
#' where \eqn{\phi} is the overdispersion parameter. See Wong et al. (2018).
#' But if \code{family="binomial"}, then 
#' \deqn{d_{x,t,p} \sim \text{Binomial}(E^0_{x,t,p} , q_{x,t,p}) , }
#' \deqn{\text{logit}(q_{x,t,p})=a_{x,p}+b_{x,p}k_{t,p} , }
#' where \eqn{q_{x,t,p}} represents the initial mortality rate at age \eqn{x} in year \eqn{t} of stratum \eqn{p},
#' while \eqn{E^0_{x,t,p}\approx E^c_{x,t,p}+\frac{1}{2}d_{x,t,p}} is the corresponding initial exposed to risk.
#' Constraints used are:
#' \deqn{\sum_{x} b_{x,p} = 1, \sum_{t} k_{t,p} = 0 \text{\ \ \ for each stratum }p.}
#' If \code{share_alpha=TRUE}, then the additive age-specific parameter is the same across all strata \eqn{p}, i. e.  
#' \deqn{a_{x}+b_{x,p}k_{t,p} .}
#' Similarly if \code{share_beta=TRUE}, then the multiplicative age-specific parameter is the same across all strata \eqn{p}, i. e.  
#' \deqn{a_{x,p}+b_{x}k_{t,p} .}
#' If \code{forecast=TRUE}, then a time series model (an AR(1) with linear drift) will be fitted on \eqn{k_{t,p}} as follows:
#' \deqn{k_{t,p} = \eta_1+\eta_2 t +\rho (k_{t-1,p}-(\eta_1+\eta_2 (t-1))) + \epsilon_{t,p} \text{ for }p=1,\ldots,P \text{ and } t=1,\ldots,T,}
#' where \eqn{\epsilon_{t,p}\sim N(0,\sigma_k^2)}, \eqn{\eta_1,\eta_2,\rho,\sigma_k^2} are additional parameters to be estimated.
#' In principle, there are many other options for forecasting the mortality time trend. But currently, we assume that this serves as a general purpose forecasting model for simplicity.
#' @references Cairns, A. J. G., Blake, D., Dowd, K., Coughlan, G. D., Epstein, D., Ong, A., and Balevich, I. (2009). A Quantitative Comparison of Stochastic Mortality Models Using Data From England and Wales and the United States. North American Actuarial Journal, 13(1), 1–35. \doi{https://doi.org/10.1080/10920277.2009.10597538}
#' @references Lee, R. D., and Carter, L. R. (1992). Modeling and Forecasting U.S. Mortality. Journal of the American Statistical Association, 87(419), 659–671. \doi{https://doi.org/10.1080/01621459.1992.10475265}
#' @references Jackie S. T. Wong, Jonathan J. Forster, and Peter W. F. Smith. (2018). Bayesian mortality forecasting with overdispersion, Insurance: Mathematics and Economics, Volume 2018, Issue 3, 206-221. \doi{https://doi.org/10.1016/j.insmatheco.2017.09.023}
#' @references Jackie S. T. Wong, Jonathan J. Forster, and Peter W. F. Smith. (2023). Bayesian model comparison for mortality forecasting, Journal of the Royal Statistical Society Series C: Applied Statistics, Volume 72, Issue 3, 566–586. \doi{https://doi.org/10.1093/jrsssc/qlad021}
#' 
#' @param death death data that has been formatted through the function \code{preparedata_fn}.
#' @param expo expo data that has been formatted through the function \code{preparedata_fn}.
#' @param n_iter number of iterations to run. Default is \code{n_iter=10000}. 
#' @param family a string of characters that defines the family function associated with the mortality model. "poisson" would assume that deaths follow a Poisson distribution and use a log link; "binomial" would assume that deaths follow a Binomial distribution and a logit link; "nb" (default) would assume that deaths follow a Negative-Binomial distribution and a log link.
#' @param share_alpha a logical value indicating if \eqn{a_{x,p}} should be shared across all strata (see details below). Default is \code{FALSE}.
#' @param share_beta a logical value indicating if \eqn{b_{x,p}} should be shared across all strata (see details below). Default is \code{FALSE}.
#' @param n.chain number of parallel chains for the model.
#' @param thin thinning interval for monitoring purposes.
#' @param n.adapt the number of iterations for adaptation. See \code{?rjags::adapt} for details.
#' @param forecast a logical value indicating if forecast is to be performed (default is \code{FALSE}). See below for details.
#' @param h a numeric value giving the number of years to forecast. Default is \code{h=5}.
#' @param quiet if TRUE then messages generated during compilation will be suppressed, as well as the progress bar during adaptation.
#' @return A list with components:
#' \describe{
#'   \item{\code{post_sample}}{An \code{mcmc.list} object containing the posterior samples generated.}
#'   \item{\code{param}}{A vector of character strings describing the names of model parameters.}
#'   \item{\code{death}}{The death data that was used.}
#'   \item{\code{expo}}{The expo data that was used.}
#'   \item{\code{family}}{The family function used.}
#'   \item{\code{forecast}}{A logical value indicating if forecast has been performed.}
#'   \item{\code{h}}{The forecast horizon used.}
#' }
#' @keywords bayesian estimation models
#' @concept stochastic mortality models
#' @concept parameter estimation
#' @concept Lee-Carter
#' @importFrom stats dnbinom dbinom dpois quantile sd
#' @export
#' @examples
#' #load and prepare mortality 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 the model (negative-binomial family)
#' #NOTE: This is a toy example, please run it longer in practice.
#' fit_LC_result<-fit_LC(death=death,expo=expo,n_iter=50,n.adapt=50)
#' head(fit_LC_result)
#' 
#' \donttest{
#' 
#' #fit the model (poisson family)
#' fit_LC_result<-fit_LC(death=death,expo=expo,n_iter=1000,family="poisson")
#' head(fit_LC_result)
#' 
#' #if sharing the betas
#' fit_LC_result2<-fit_LC(death=death,expo=expo,n_iter=1000,family="poisson",share_beta=TRUE)
#' head(fit_LC_result2)
#' 
#' #if sharing both alphas and betas
#' fit_LC_result3<-fit_LC(death=death,expo=expo,n_iter=1000,family="poisson",
#' share_alpha=TRUE,share_beta=TRUE)
#' head(fit_LC_result3)
#' 
#' #if forecast
#' fit_LC_result4<-fit_LC(death=death,expo=expo,n_iter=1000,family="poisson",forecast=TRUE)
#' plot_rates_fn(fit_LC_result4)
#' plot_param_fn(fit_LC_result4)
#' }

fit_LC<-function(death,expo,n_iter=10000,family="nb",share_alpha=FALSE,share_beta=FALSE,n.chain=1,thin=1,n.adapt=1000,forecast=FALSE,h=5,quiet=FALSE){
  
  p<-death$n_strat
  A<-death$n_ages
  T<-death$n_years
  
  prior_mean_beta=rep(1/A,A-1)
  sigma2_beta<-0.1
  prior_prec_beta=solve(sigma2_beta*(diag(rep(1,A-1))-1/A*(matrix(1,nrow=A-1,ncol=A-1))))
  
  prior_mean_kappa=rep(0,T-1)
  sigma2_kappa<-1000
  matrix_kappa_A<-diag(rep(1,T));matrix_kappa_A[1,]<-1
  matrix_kappa_B<-matrix_kappa_A%*%t(matrix_kappa_A)
  prior_prec_kappa=solve((matrix_kappa_B[-1,-1]-1/matrix_kappa_B[1,1]*matrix_kappa_B[-1,1]%*%t(matrix_kappa_B[1,-1]))*sigma2_kappa)
  
  if (!forecast){
    
    t<-(1:T)-mean(1:T)
    matrix_kappa_X<-matrix(c(rep(1,T),t),byrow=F,ncol=2)
    
    prior_prec_eta<-solve(matrix(c(400,0,0,2),nrow=2));prior_mean_eta<-c(0,0)
    
  if (family=="binomial"){
    expo_initial<-round(expo$data+0.5*death$data)
    if (share_alpha){
      if (share_beta){
        #4.
        data<-list(dxt=death$data,ext=expo_initial,A=A,T=T,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
        inits<-function() (list(alpha=rep(0,A),beta_rest=rep(1/A,(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1))
        vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa")
        logit_LC_CBD_M1_shareall_jags<-rjags::jags.model(system.file("models/logit_LC_CBD_M1_shareall.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
        result_jags<-rjags::coda.samples(logit_LC_CBD_M1_shareall_jags,vars,n.iter=n_iter,thin=thin)
      }else{
        #2.
        data<-list(dxt=death$data,ext=expo_initial,A=A,T=T,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
        inits<-function() (list(alpha=rep(0,A),beta_rest_mat=matrix(1/A,nrow=p,ncol=(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1))
        vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa")
        logit_LC_CBD_M1_sharealpha_jags<-rjags::jags.model(system.file("models/logit_LC_CBD_M1_sharealpha.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
        result_jags<-rjags::coda.samples(logit_LC_CBD_M1_sharealpha_jags,vars,n.iter=n_iter,thin=thin)
      }
    } else {
      if (share_beta){
        #3.
        data<-list(dxt=death$data,ext=expo_initial,A=A,T=T,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
        inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta_rest=rep(1/A,(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1))
        vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa")
        logit_LC_CBD_M1_sharebeta_jags<-rjags::jags.model(system.file("models/logit_LC_CBD_M1_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
        result_jags<-rjags::coda.samples(logit_LC_CBD_M1_sharebeta_jags,vars,n.iter=n_iter,thin=thin)
      } else {
        #1.
        data<-list(dxt=death$data,ext=expo_initial,A=A,T=T,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
        inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta_rest_mat=matrix(1/A,nrow=p,ncol=(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1))
        vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa")
        logit_LC_all_separately_jags<-rjags::jags.model(system.file("models/logit_LC_CBD_M1_sep.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
        result_jags<-rjags::coda.samples(logit_LC_all_separately_jags,vars,n.iter=n_iter,thin=thin)
      }
    }
    
  }
  if (family=="poisson"){
    if (share_alpha){
      if (share_beta){
        #4.
        data<-list(dxt=death$data,ext=expo$data,A=A,T=T,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
        inits<-function() (list(alpha=rep(0,A),beta_rest=rep(1/A,(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1))
        vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa")
        log_LC_CBD_M1_shareall_jags<-rjags::jags.model(system.file("models/log_LC_CBD_M1_shareall.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
        result_jags<-rjags::coda.samples(log_LC_CBD_M1_shareall_jags,vars,n.iter=n_iter,thin=thin)
      }else{
        #2.
        data<-list(dxt=death$data,ext=expo$data,A=A,T=T,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
        inits<-function() (list(alpha=rep(0,A),beta_rest_mat=matrix(1/A,nrow=p,ncol=(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1))
        vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa")
        log_LC_CBD_M1_sharealpha_jags<-rjags::jags.model(system.file("models/log_LC_CBD_M1_sharealpha.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
        result_jags<-rjags::coda.samples(log_LC_CBD_M1_sharealpha_jags,vars,n.iter=n_iter,thin=thin)
      }
    } else {
      if (share_beta){
        #3.
        data<-list(dxt=death$data,ext=expo$data,A=A,T=T,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
        inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta_rest=rep(1/A,(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1))
        vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa")
        log_LC_CBD_M1_sharebeta_jags<-rjags::jags.model(system.file("models/log_LC_CBD_M1_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
        result_jags<-rjags::coda.samples(log_LC_CBD_M1_sharebeta_jags,vars,n.iter=n_iter,thin=thin)
      } else {
        #1.
        data<-list(dxt=death$data,ext=expo$data,A=A,T=T,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
        inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta_rest_mat=matrix(1/A,nrow=p,ncol=(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1))
        vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa")
        log_LC_all_separately_jags<-rjags::jags.model(system.file("models/log_LC_CBD_M1_sep.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
        result_jags<-rjags::coda.samples(log_LC_all_separately_jags,vars,n.iter=n_iter,thin=thin)
      }
    }
    
  }
    
    if (family=="nb"){
      if (share_alpha){
        if (share_beta){
          #4.
          data<-list(dxt=death$data,ext=expo$data,A=A,T=T,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
          inits<-function() (list(alpha=rep(0,A),beta_rest=rep(1/A,(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1,phi=100))
          vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa","phi")
          nb_LC_CBD_M1_shareall_jags<-rjags::jags.model(system.file("models/nb_LC_CBD_M1_shareall.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
          result_jags<-rjags::coda.samples(nb_LC_CBD_M1_shareall_jags,vars,n.iter=n_iter,thin=thin)
        }else{
          #2.
          data<-list(dxt=death$data,ext=expo$data,A=A,T=T,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
          inits<-function() (list(alpha=rep(0,A),beta_rest_mat=matrix(1/A,nrow=p,ncol=(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1,phi=100))
          vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa","phi")
          nb_LC_CBD_M1_sharealpha_jags<-rjags::jags.model(system.file("models/nb_LC_CBD_M1_sharealpha.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
          result_jags<-rjags::coda.samples(nb_LC_CBD_M1_sharealpha_jags,vars,n.iter=n_iter,thin=thin)
        }
      } else {
        if (share_beta){
          #3.
          data<-list(dxt=death$data,ext=expo$data,A=A,T=T,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
          inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta_rest=rep(1/A,(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1,phi=100))
          vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa","phi")
          nb_LC_CBD_M1_sharebeta_jags<-rjags::jags.model(system.file("models/nb_LC_CBD_M1_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
          result_jags<-rjags::coda.samples(nb_LC_CBD_M1_sharebeta_jags,vars,n.iter=n_iter,thin=thin)
        } else {
          #1.
          data<-list(dxt=death$data,ext=expo$data,A=A,T=T,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
          inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta_rest_mat=matrix(1/A,nrow=p,ncol=(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1,phi=100))
          vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa","phi")
          nb_LC_all_separately_jags<-rjags::jags.model(system.file("models/nb_LC_CBD_M1_sep.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
          result_jags<-rjags::coda.samples(nb_LC_all_separately_jags,vars,n.iter=n_iter,thin=thin)
        }
      }
      
    }
  }
  
  if (forecast){
    
    t<-(1:T)-mean(1:T)
    matrix_kappa_X<-matrix(c(rep(1,T+h),c(t,t[T]+1:h)),byrow=F,ncol=2)
    
    prior_prec_eta<-solve(matrix(c(400,0,0,2),nrow=2));prior_mean_eta<-c(0,0)
    
    #d_names<-dimnames(death)
    #d_names[[3]]<-c(d_names[[3]],paste0("F_year_",1:h))
    
    death_forecast<-array(dim=c(p,A,T+h));expo_forecast<-array(dim=c(p,A,T+h))
    death_forecast[,,1:T]<-death$data
    death_forecast[,,(T+1):(T+h)]<-NA
    expo_forecast[,,1:T]<-expo$data
    expo_forecast[,,(T+1):(T+h)]<-expo$data[,,T]
    
    if (family=="binomial"){
      expo_forecast_initial<-expo_forecast
      expo_forecast_initial[,,1:T]<-round(expo_forecast[,,1:T,drop=FALSE]+0.5*death$data)
      expo_forecast_initial[,,(T+1):(T+h)]<-round(expo$data[,,T]+0.5*death$data[,,T])
      
      if (share_alpha){
        if (share_beta){
          data<-list(dxt=death_forecast,ext=expo_forecast_initial,A=A,T=T,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
          inits<-function() (list(alpha=rep(0,A),beta_rest=rep(1/A,(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1))
          vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa")
          logit_LC_shareall_AR1_alt_forecast_jags<-rjags::jags.model(system.file("models/logit_LC_CBD_M1_AR1_alt_forecast_shareall.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
          result_jags<-rjags::coda.samples(logit_LC_shareall_AR1_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
        } else {
          data<-list(dxt=death_forecast,ext=expo_forecast_initial,A=A,T=T,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
          inits<-function() (list(alpha=rep(0,A),beta_rest_mat=matrix(1/A,nrow=p,ncol=(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1))
          vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa")
          logit_LC_sharealpha_AR1_alt_forecast_jags<-rjags::jags.model(system.file("models/logit_LC_CBD_M1_AR1_alt_forecast_sharealpha.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
          result_jags<-rjags::coda.samples(logit_LC_sharealpha_AR1_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
        }
      } else {
        if (share_beta){
          data<-list(dxt=death_forecast,ext=expo_forecast_initial,A=A,T=T,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
          inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta_rest=rep(1/A,(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1))
          vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa")
          logit_LC_sharebeta_AR1_alt_forecast_jags<-rjags::jags.model(system.file("models/logit_LC_CBD_M1_AR1_alt_forecast_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
          result_jags<-rjags::coda.samples(logit_LC_sharebeta_AR1_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
        } else {
          data<-list(dxt=death_forecast,ext=expo_forecast_initial,A=A,T=T,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
          inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta_rest_mat=matrix(1/A,nrow=p,ncol=(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1))
          vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa")
          logit_LC_all_separately_AR1_alt_forecast_jags<-rjags::jags.model(system.file("models/logit_LC_CBD_M1_AR1_alt_forecast_sep.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
          result_jags<-rjags::coda.samples(logit_LC_all_separately_AR1_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
        }
      }
      
    }
    if (family=="poisson"){
      
      if (share_alpha){
        if (share_beta){
          data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
          inits<-function() (list(alpha=rep(0,A),beta_rest=rep(1/A,(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1))
          vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa")
          log_LC_shareall_AR1_alt_forecast_jags<-rjags::jags.model(system.file("models/log_LC_CBD_M1_AR1_alt_forecast_shareall.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
          result_jags<-rjags::coda.samples(log_LC_shareall_AR1_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
        } else {
          data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
          inits<-function() (list(alpha=rep(0,A),beta_rest_mat=matrix(1/A,nrow=p,ncol=(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1))
          vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa")
          log_LC_sharealpha_AR1_alt_forecast_jags<-rjags::jags.model(system.file("models/log_LC_CBD_M1_AR1_alt_forecast_sharealpha.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
          result_jags<-rjags::coda.samples(log_LC_sharealpha_AR1_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
        }
      } else {
        if (share_beta){
          data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
          inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta_rest=rep(1/A,(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1))
          vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa")
          log_LC_sharebeta_AR1_alt_forecast_jags<-rjags::jags.model(system.file("models/log_LC_CBD_M1_AR1_alt_forecast_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
          result_jags<-rjags::coda.samples(log_LC_sharebeta_AR1_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
        } else {
          data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
          inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta_rest_mat=matrix(1/A,nrow=p,ncol=(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1))
          vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa")
          log_LC_all_separately_AR1_alt_forecast_jags<-rjags::jags.model(system.file("models/log_LC_CBD_M1_AR1_alt_forecast_sep.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
          result_jags<-rjags::coda.samples(log_LC_all_separately_AR1_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
        }
      }
      
    }
    
    if (family=="nb"){
      
      if (share_alpha){
        if (share_beta){
          data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
          inits<-function() (list(alpha=rep(0,A),beta_rest=rep(1/A,(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1,phi=100))
          vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa","phi")
          nb_LC_shareall_AR1_alt_forecast_jags<-rjags::jags.model(system.file("models/nb_LC_CBD_M1_AR1_alt_forecast_shareall.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
          result_jags<-rjags::coda.samples(nb_LC_shareall_AR1_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
        } else {
          data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
          inits<-function() (list(alpha=rep(0,A),beta_rest_mat=matrix(1/A,nrow=p,ncol=(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1,phi=100))
          vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa","phi")
          nb_LC_sharealpha_AR1_alt_forecast_jags<-rjags::jags.model(system.file("models/nb_LC_CBD_M1_AR1_alt_forecast_sharealpha.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
          result_jags<-rjags::coda.samples(nb_LC_sharealpha_AR1_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
        }
      } else {
        if (share_beta){
          data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
          inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta_rest=rep(1/A,(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1,phi=100))
          vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa","phi")
          nb_LC_sharebeta_AR1_alt_forecast_jags<-rjags::jags.model(system.file("models/nb_LC_CBD_M1_AR1_alt_forecast_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
          result_jags<-rjags::coda.samples(nb_LC_sharebeta_AR1_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
        } else {
          data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,matrix_kappa_X=matrix_kappa_X,prior_mean_eta=prior_mean_eta,prior_prec_eta=prior_prec_eta)
          inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta_rest_mat=matrix(1/A,nrow=p,ncol=(A-1)),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-1)),rho=0.5,eta=c(0,0),i_sigma2_kappa=0.1,phi=100))
          vars<-c("q","alpha","beta","kappa","eta","rho","sigma2_kappa","phi")
          nb_LC_all_separately_AR1_alt_forecast_jags<-rjags::jags.model(system.file("models/nb_LC_CBD_M1_AR1_alt_forecast_sep.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
          result_jags<-rjags::coda.samples(nb_LC_all_separately_AR1_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
        }
      }
      
    }
      
  }
  
  invisible(gc())
  list(post_sample=result_jags,param=vars[-1],death=death,expo=expo,family=family,forecast=forecast,h=h)
  
}

fit_LC_sharealpha<-function(death,expo,n_iter=10000,family="nb",n.chain=1,thin=1,n.adapt=1000,forecast=FALSE,h=5,quiet=FALSE){
  
  fit_LC(death=death,expo=expo,share_alpha = TRUE,n_iter=n_iter,family=family,n.chain=n.chain,thin=thin,n.adapt=n.adapt,quiet=quiet,forecast=forecast,h=h)
  
}

fit_LC_sharebeta<-function(death,expo,n_iter=10000,family="nb",n.chain=1,thin=1,n.adapt=1000,forecast=FALSE,h=5,quiet=FALSE){
  
  fit_LC(death=death,expo=expo,share_beta = TRUE,n_iter=n_iter,family=family,n.chain=n.chain,thin=thin,n.adapt=n.adapt,quiet=quiet,forecast=forecast,h=h)
  
}

fit_LC_shareall<-function(death,expo,n_iter=10000,family="nb",n.chain=1,thin=1,n.adapt=1000,forecast=FALSE,h=5,quiet=FALSE){
  
  fit_LC(death=death,expo=expo,share_beta = TRUE,share_alpha = TRUE,n_iter=n_iter,family=family,n.chain=n.chain,thin=thin,n.adapt=n.adapt,quiet=quiet,forecast=forecast,h=h)
  
}