Nothing
R/fit_APCI.R
In BayesMoFo: Bayesian Mortality Forecasting
Defines functions
fit_APCI_shareall
fit_APCI_sharebeta_sharegamma
fit_APCI_sharealpha_sharegamma
fit_APCI_sharealpha_sharebeta
fit_APCI_sharegamma
fit_APCI_sharebeta
fit_APCI_sharealpha
fit_APCI
Documented in
fit_APCI
#' A function to fit the APCI stochastic mortality model.
#'
#' Carry out Bayesian estimation of the stochastic mortality model, the \bold{Age-Period-Cohort-Improvement (APCI) model}.
#' Note that this model has been studied extensively by Richards et al. (2019) and Wong et al. (2023).
#'
#' 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}(t-\bar{t})+k_{t,p} + \gamma_{c,p} , }
#' where \eqn{c=t-x} is the cohort index, \eqn{\bar{t}} is the mean of \eqn{t},
#' \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}(t-\bar{t})+k_{t,p} + \gamma_{c,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}(t-\bar{t})+k_{t,p} + \gamma_{c,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_{t} k_{t,p} = \sum_{t} tk_{t,p} = 0, \sum_{c} \gamma_{c,p} = \sum_{c}c\gamma_{c,p} = \sum_{c}c^2\gamma_{c,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}(t-\bar{t})+k_{t,p}+ \gamma_{c,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}(t-\bar{t})+k_{t,p}+ \gamma_{c,p} .}
#' Similarly if \code{share_gamma=TRUE}, then the cohort parameter is the same across all strata \eqn{p}, i. e.
#' \deqn{a_{x,p}+b_{x,p}(t-\bar{t})+k_{t,p}+ \gamma_{c} .}
#' If \code{forecast=TRUE}, then the following time series models are fitted on \eqn{k_{t,p}} and \eqn{\gamma_{c,p}} as follows:
#' \deqn{k_{t,p} = \rho k_{t-1,p} + \epsilon_{t,p} \text{ for }p=1,\ldots,P \text{ and } t=2,\ldots,T,}
#' \deqn{k_{1,p}=\epsilon_{1,p}}
#' and
#' \deqn{\gamma_{c,p} = \gamma_{c-1,p}+ \rho_\gamma (\gamma_{c-1,p}-\gamma_{c-2,p}) + \epsilon^\gamma_{c,p} \text{ for }p=1,\ldots,P \text{ and } c=3,\ldots,C,}
#' \deqn{\gamma_2=\gamma_1+\frac{1}{\sqrt{1-\rho_\gamma^2}}\epsilon^\gamma_{2,p},}
#' \deqn{\gamma_1=100\epsilon^\gamma_{1,p}}
#' where \eqn{\epsilon_{t,p}\sim N(0,\sigma_k^2)} for \eqn{t=1,\ldots,T}, \eqn{\epsilon^\gamma_{c,p}\sim N(0,\sigma_\gamma^2)} for \eqn{c=1,\ldots,C}, while \eqn{\eta_1,\eta_2,\rho,\sigma_k^2,\rho_\gamma,\sigma_\gamma^2} are additional parameters to be estimated.
#' Note that the forecasting models are inspired by Wong et al. (2023).
#'
#' @references Richards S. J., Currie I. D., Kleinow T., and Ritchie G. P. (2019). A stochastic implementation of the APCI model for mortality projections. British Actuarial Journal, 24, E13. \doi{https://doi.org/10.1017/S1357321718000260}
#' @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 share_gamma a logical value indicating if the cohort parameter \eqn{\gamma_{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 APCI
#' @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_APCI_result<-fit_APCI(death=death,expo=expo,n_iter=50,n.adapt=50)
#' head(fit_APCI_result)
#'
#' \donttest{
#' #if sharing the cohorts only (poisson family)
#' fit_APCI_result2<-fit_APCI(death=death,expo=expo,n_iter=1000,family="poisson",share_gamma=TRUE)
#' head(fit_APCI_result2)
#'
#' #if sharing all alphas, betas, and cohorts (poisson family)
#' fit_APCI_result3<-fit_APCI(death=death,expo=expo,n_iter=1000,family="poisson",share_alpha=TRUE,
#' share_beta=TRUE,share_gamma=TRUE)
#' head(fit_APCI_result3)
#'
#' #if forecast
#' fit_APCI_result4<-fit_APCI(death=death,expo=expo,n_iter=1000,family="poisson",forecast=TRUE)
#' plot_rates_fn(fit_APCI_result4)
#' plot_param_fn(fit_APCI_result4)
#' }
fit_APCI<-function(death,expo,n_iter=10000,family="nb",share_alpha=FALSE,share_beta=FALSE,share_gamma=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
C<-A+T-1
cohorts<-1:C;cohorts_rev<-rev(cohorts)
t<-1:T;t_rev<-rev(t)
prior_mean_beta=rep(0,A)
sigma2_beta<-1
prior_prec_beta<-solve(diag(rep(1,A))*sigma2_beta)
prior_mean_kappa=rep(0,T-2)
sigma2_kappa<-1000
matrix_kappa_A<-diag(rep(1,T))
matrix_kappa_A<-rbind(rep(1/T,T),((1:T)-mean((1:T)))/sd((1:T)),matrix_kappa_A[(2:(T-1)),])
matrix_kappa_B<-matrix_kappa_A%*%t(matrix_kappa_A)
prior_prec_kappa=solve((matrix_kappa_B[3:T,3:T]-matrix_kappa_B[3:T,1:2]%*%solve(matrix_kappa_B[1:2,1:2])%*%matrix_kappa_B[1:2,3:T])*sigma2_kappa)
prior_mean_cohort=rep(0,C-3)
sigma2_cohort<-1
matrix_cohort_A<-diag(rep(1,C))
matrix_cohort_A<-rbind(rep(1/C,C),(cohorts-mean(cohorts))/sd(cohorts),(cohorts^2-mean(cohorts^2))/sd(cohorts^2),matrix_cohort_A[(2:(C-2)),])
matrix_cohort_B<-matrix_cohort_A%*%t(matrix_cohort_A)
prior_prec_cohort=solve((matrix_cohort_B[4:C,4:C]-matrix_cohort_B[4:C,1:3]%*%solve(matrix_cohort_B[1:3,1:3])%*%matrix_cohort_B[1:3,4:C])*sigma2_cohort)
if (!forecast){if (family=="binomial"){
expo_initial<-round(expo$data+0.5*death$data)
if (share_alpha){
if (share_beta){
if (share_gamma){
#8
data<-list(dxt=death$data,ext=expo_initial,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_shareall_jags<-rjags::jags.model(system.file("models/logit_APCI_shareall.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_shareall_jags,vars,n.iter=n_iter,thin=thin)} else {
#5
data<-list(dxt=death$data,ext=expo_initial,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharealpha_sharebeta_jags<-rjags::jags.model(system.file("models/logit_APCI_sharealpha_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sharealpha_sharebeta_jags,vars,n.iter=n_iter,thin=thin)
}
}else{
if (share_gamma){
#6
data<-list(dxt=death$data,ext=expo_initial,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharealpha_sharegamma_jags<-rjags::jags.model(system.file("models/logit_APCI_sharealpha_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sharealpha_sharegamma_jags,vars,n.iter=n_iter,thin=thin)
} else {
#2
data<-list(dxt=death$data,ext=expo_initial,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharealpha_jags<-rjags::jags.model(system.file("models/logit_APCI_sharealpha.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sharealpha_jags,vars,n.iter=n_iter,thin=thin)
}
}
} else {
if (share_beta){
if (share_gamma){
#7
data<-list(dxt=death$data,ext=expo_initial,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharebeta_sharegamma_jags<-rjags::jags.model(system.file("models/logit_APCI_sharebeta_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sharebeta_sharegamma_jags,vars,n.iter=n_iter,thin=thin)
} else {
#3
data<-list(dxt=death$data,ext=expo_initial,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharebeta_jags<-rjags::jags.model(system.file("models/logit_APCI_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sharebeta_jags,vars,n.iter=n_iter,thin=thin)
}
} else {
if (share_gamma){
#4
data<-list(dxt=death$data,ext=expo_initial,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharegamma_jags<-rjags::jags.model(system.file("models/logit_APCI_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sharegamma_jags,vars,n.iter=n_iter,thin=thin)
}else{
#1
data<-list(dxt=death$data,ext=expo_initial,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sep_jags<-rjags::jags.model(system.file("models/logit_APCI_sep.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sep_jags,vars,n.iter=n_iter,thin=thin)
}
}
}
}
if (family=="nb"){
if (share_alpha){
if (share_beta){
if (share_gamma){
#8
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_shareall_jags<-rjags::jags.model(system.file("models/nb_APCI_shareall.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_shareall_jags,vars,n.iter=n_iter,thin=thin)} else {
#5
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharealpha_sharebeta_jags<-rjags::jags.model(system.file("models/nb_APCI_sharealpha_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sharealpha_sharebeta_jags,vars,n.iter=n_iter,thin=thin)
}
}else{
if (share_gamma){
#6
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharealpha_sharegamma_jags<-rjags::jags.model(system.file("models/nb_APCI_sharealpha_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sharealpha_sharegamma_jags,vars,n.iter=n_iter,thin=thin)
} else {
#2
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharealpha_jags<-rjags::jags.model(system.file("models/nb_APCI_sharealpha.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sharealpha_jags,vars,n.iter=n_iter,thin=thin)
}
}
} else {
if (share_beta){
if (share_gamma){
#7
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharebeta_sharegamma_jags<-rjags::jags.model(system.file("models/nb_APCI_sharebeta_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sharebeta_sharegamma_jags,vars,n.iter=n_iter,thin=thin)
} else {
#3
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharebeta_jags<-rjags::jags.model(system.file("models/nb_APCI_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sharebeta_jags,vars,n.iter=n_iter,thin=thin)
}
} else {
if (share_gamma){
#4
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharegamma_jags<-rjags::jags.model(system.file("models/nb_APCI_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sharegamma_jags,vars,n.iter=n_iter,thin=thin)
}else{
#1
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sep_jags<-rjags::jags.model(system.file("models/nb_APCI_sep.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sep_jags,vars,n.iter=n_iter,thin=thin)
}
}
}
}
if (family=="poisson"){
if (share_alpha){
if (share_beta){
if (share_gamma){
#8
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_shareall_jags<-rjags::jags.model(system.file("models/log_APCI_shareall.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_shareall_jags,vars,n.iter=n_iter,thin=thin)} else {
#5
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharealpha_sharebeta_jags<-rjags::jags.model(system.file("models/log_APCI_sharealpha_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sharealpha_sharebeta_jags,vars,n.iter=n_iter,thin=thin)
}
}else{
if (share_gamma){
#6
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharealpha_sharegamma_jags<-rjags::jags.model(system.file("models/log_APCI_sharealpha_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sharealpha_sharegamma_jags,vars,n.iter=n_iter,thin=thin)
} else {
#2
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharealpha_jags<-rjags::jags.model(system.file("models/log_APCI_sharealpha.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sharealpha_jags,vars,n.iter=n_iter,thin=thin)
}
}
} else {
if (share_beta){
if (share_gamma){
#7
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharebeta_sharegamma_jags<-rjags::jags.model(system.file("models/log_APCI_sharebeta_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sharebeta_sharegamma_jags,vars,n.iter=n_iter,thin=thin)
} else {
#3
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharebeta_jags<-rjags::jags.model(system.file("models/log_APCI_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sharebeta_jags,vars,n.iter=n_iter,thin=thin)
}
} else {
if (share_gamma){
#4
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharegamma_jags<-rjags::jags.model(system.file("models/log_APCI_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sharegamma_jags,vars,n.iter=n_iter,thin=thin)
}else{
#1
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sep_jags<-rjags::jags.model(system.file("models/log_APCI_sep.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sep_jags,vars,n.iter=n_iter,thin=thin)
}
}
}
}
}
if (forecast){
t_forecast<-c(t,t[T]+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){
if (share_gamma){
#8
data<-list(dxt=death_forecast,ext=expo_forecast_initial,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_shareall_alt_forecast_jags<-rjags::jags.model(system.file("models/logit_APCI_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_APCI_shareall_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
} else {
#5
data<-list(dxt=death_forecast,ext=expo_forecast_initial,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharealpha_sharebeta_alt_forecast_jags<-rjags::jags.model(system.file("models/logit_APCI_alt_forecast_sharealpha_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sharealpha_sharebeta_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
}else{
if (share_gamma){
#6
data<-list(dxt=death_forecast,ext=expo_forecast_initial,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharealpha_sharegamma_alt_forecast_jags<-rjags::jags.model(system.file("models/logit_APCI_alt_forecast_sharealpha_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sharealpha_sharegamma_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
} else {
#2
data<-list(dxt=death_forecast,ext=expo_forecast_initial,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharealpha_alt_forecast_jags<-rjags::jags.model(system.file("models/logit_APCI_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_APCI_sharealpha_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
}
} else {
if (share_beta){
if (share_gamma){
#7
data<-list(dxt=death_forecast,ext=expo_forecast_initial,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharebeta_sharegamma_alt_forecast_jags<-rjags::jags.model(system.file("models/logit_APCI_alt_forecast_sharebeta_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sharebeta_sharegamma_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
} else {
#3
data<-list(dxt=death_forecast,ext=expo_forecast_initial,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharebeta_alt_forecast_jags<-rjags::jags.model(system.file("models/logit_APCI_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_APCI_sharebeta_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
} else {
if (share_gamma){
#4
data<-list(dxt=death_forecast,ext=expo_forecast_initial,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharegamma_alt_forecast_jags<-rjags::jags.model(system.file("models/logit_APCI_alt_forecast_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sharegamma_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}else{
#1
data<-list(dxt=death_forecast,ext=expo_forecast_initial,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sep_alt_forecast_jags<-rjags::jags.model(system.file("models/logit_APCI_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_APCI_sep_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
}
}
}
if (family=="nb"){
if (share_alpha){
if (share_beta){
if (share_gamma){
#8
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_shareall_alt_forecast_jags<-rjags::jags.model(system.file("models/nb_APCI_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_APCI_shareall_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
} else {
#5
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharealpha_sharebeta_alt_forecast_jags<-rjags::jags.model(system.file("models/nb_APCI_alt_forecast_sharealpha_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sharealpha_sharebeta_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
}else{
if (share_gamma){
#6
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharealpha_sharegamma_alt_forecast_jags<-rjags::jags.model(system.file("models/nb_APCI_alt_forecast_sharealpha_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sharealpha_sharegamma_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
} else {
#2
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharealpha_alt_forecast_jags<-rjags::jags.model(system.file("models/nb_APCI_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_APCI_sharealpha_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
}
} else {
if (share_beta){
if (share_gamma){
#7
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharebeta_sharegamma_alt_forecast_jags<-rjags::jags.model(system.file("models/nb_APCI_alt_forecast_sharebeta_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sharebeta_sharegamma_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
} else {
#3
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharebeta_alt_forecast_jags<-rjags::jags.model(system.file("models/nb_APCI_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_APCI_sharebeta_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
} else {
if (share_gamma){
#4
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharegamma_alt_forecast_jags<-rjags::jags.model(system.file("models/nb_APCI_alt_forecast_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sharegamma_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}else{
#1
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sep_alt_forecast_jags<-rjags::jags.model(system.file("models/nb_APCI_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_APCI_sep_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
}
}
}
if (family=="poisson"){
if (share_alpha){
if (share_beta){
if (share_gamma){
#8
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_shareall_alt_forecast_jags<-rjags::jags.model(system.file("models/log_APCI_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_APCI_shareall_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
} else {
#5
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharealpha_sharebeta_alt_forecast_jags<-rjags::jags.model(system.file("models/log_APCI_alt_forecast_sharealpha_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sharealpha_sharebeta_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
}else{
if (share_gamma){
#6
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharealpha_sharegamma_alt_forecast_jags<-rjags::jags.model(system.file("models/log_APCI_alt_forecast_sharealpha_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sharealpha_sharegamma_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
} else {
#2
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharealpha_alt_forecast_jags<-rjags::jags.model(system.file("models/log_APCI_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_APCI_sharealpha_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
}
} else {
if (share_beta){
if (share_gamma){
#7
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharebeta_sharegamma_alt_forecast_jags<-rjags::jags.model(system.file("models/log_APCI_alt_forecast_sharebeta_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sharebeta_sharegamma_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
} else {
#3
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharebeta_alt_forecast_jags<-rjags::jags.model(system.file("models/log_APCI_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_APCI_sharebeta_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
} else {
if (share_gamma){
#4
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharegamma_alt_forecast_jags<-rjags::jags.model(system.file("models/log_APCI_alt_forecast_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sharegamma_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}else{
#1
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sep_alt_forecast_jags<-rjags::jags.model(system.file("models/log_APCI_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_APCI_sep_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_APCI_sharealpha<-function(death,expo,n_iter=10000,family="nb",n.chain=1,thin=1,n.adapt=1000,forecast=FALSE,h=5,quiet=FALSE){
fit_APCI(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_APCI_sharebeta<-function(death,expo,n_iter=10000,family="nb",n.chain=1,thin=1,n.adapt=1000,forecast=FALSE,h=5,quiet=FALSE){
fit_APCI(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_APCI_sharegamma<-function(death,expo,n_iter=10000,family="nb",n.chain=1,thin=1,n.adapt=1000,forecast=FALSE,h=5,quiet=FALSE){
fit_APCI(death=death,expo=expo,share_gamma = TRUE,n_iter=n_iter,family=family,n.chain=n.chain,thin=thin,n.adapt=n.adapt,quiet=quiet,forecast=forecast,h=h)
}
fit_APCI_sharealpha_sharebeta<-function(death,expo,n_iter=10000,family="nb",n.chain=1,thin=1,n.adapt=1000,forecast=FALSE,h=5,quiet=FALSE){
fit_APCI(death=death,expo=expo,share_alpha = TRUE,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_APCI_sharealpha_sharegamma<-function(death,expo,n_iter=10000,family="nb",n.chain=1,thin=1,n.adapt=1000,forecast=FALSE,h=5,quiet=FALSE){
fit_APCI(death=death,expo=expo,share_alpha = TRUE,share_gamma = TRUE,n_iter=n_iter,family=family,n.chain=n.chain,thin=thin,n.adapt=n.adapt,quiet=quiet,forecast=forecast,h=h)
}
fit_APCI_sharebeta_sharegamma<-function(death,expo,n_iter=10000,family="nb",n.chain=1,thin=1,n.adapt=1000,forecast=FALSE,h=5,quiet=FALSE){
fit_APCI(death=death,expo=expo,share_beta = TRUE,share_gamma = TRUE,n_iter=n_iter,family=family,n.chain=n.chain,thin=thin,n.adapt=n.adapt,quiet=quiet,forecast=forecast,h=h)
}
fit_APCI_shareall<-function(death,expo,n_iter=10000,family="nb",n.chain=1,thin=1,n.adapt=1000,forecast=FALSE,h=5,quiet=FALSE){
fit_APCI(death=death,expo=expo,share_alpha = TRUE,share_beta = TRUE,share_gamma = TRUE,n_iter=n_iter,family=family,n.chain=n.chain,thin=thin,n.adapt=n.adapt,quiet=quiet,forecast=forecast,h=h)
}
Try the BayesMoFo package in your browser
Any scripts or data that you put into this service are public.
BayesMoFo documentation built on Aug. 11, 2025, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions fit_APCI_shareall fit_APCI_sharebeta_sharegamma fit_APCI_sharealpha_sharegamma fit_APCI_sharealpha_sharebeta fit_APCI_sharegamma fit_APCI_sharebeta fit_APCI_sharealpha fit_APCI
Documented in fit_APCI
#' A function to fit the APCI stochastic mortality model.
#'
#' Carry out Bayesian estimation of the stochastic mortality model, the \bold{Age-Period-Cohort-Improvement (APCI) model}.
#' Note that this model has been studied extensively by Richards et al. (2019) and Wong et al. (2023).
#'
#' 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}(t-\bar{t})+k_{t,p} + \gamma_{c,p} , }
#' where \eqn{c=t-x} is the cohort index, \eqn{\bar{t}} is the mean of \eqn{t},
#' \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}(t-\bar{t})+k_{t,p} + \gamma_{c,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}(t-\bar{t})+k_{t,p} + \gamma_{c,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_{t} k_{t,p} = \sum_{t} tk_{t,p} = 0, \sum_{c} \gamma_{c,p} = \sum_{c}c\gamma_{c,p} = \sum_{c}c^2\gamma_{c,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}(t-\bar{t})+k_{t,p}+ \gamma_{c,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}(t-\bar{t})+k_{t,p}+ \gamma_{c,p} .}
#' Similarly if \code{share_gamma=TRUE}, then the cohort parameter is the same across all strata \eqn{p}, i. e.
#' \deqn{a_{x,p}+b_{x,p}(t-\bar{t})+k_{t,p}+ \gamma_{c} .}
#' If \code{forecast=TRUE}, then the following time series models are fitted on \eqn{k_{t,p}} and \eqn{\gamma_{c,p}} as follows:
#' \deqn{k_{t,p} = \rho k_{t-1,p} + \epsilon_{t,p} \text{ for }p=1,\ldots,P \text{ and } t=2,\ldots,T,}
#' \deqn{k_{1,p}=\epsilon_{1,p}}
#' and
#' \deqn{\gamma_{c,p} = \gamma_{c-1,p}+ \rho_\gamma (\gamma_{c-1,p}-\gamma_{c-2,p}) + \epsilon^\gamma_{c,p} \text{ for }p=1,\ldots,P \text{ and } c=3,\ldots,C,}
#' \deqn{\gamma_2=\gamma_1+\frac{1}{\sqrt{1-\rho_\gamma^2}}\epsilon^\gamma_{2,p},}
#' \deqn{\gamma_1=100\epsilon^\gamma_{1,p}}
#' where \eqn{\epsilon_{t,p}\sim N(0,\sigma_k^2)} for \eqn{t=1,\ldots,T}, \eqn{\epsilon^\gamma_{c,p}\sim N(0,\sigma_\gamma^2)} for \eqn{c=1,\ldots,C}, while \eqn{\eta_1,\eta_2,\rho,\sigma_k^2,\rho_\gamma,\sigma_\gamma^2} are additional parameters to be estimated.
#' Note that the forecasting models are inspired by Wong et al. (2023).
#'
#' @references Richards S. J., Currie I. D., Kleinow T., and Ritchie G. P. (2019). A stochastic implementation of the APCI model for mortality projections. British Actuarial Journal, 24, E13. \doi{https://doi.org/10.1017/S1357321718000260}
#' @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 share_gamma a logical value indicating if the cohort parameter \eqn{\gamma_{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 APCI
#' @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_APCI_result<-fit_APCI(death=death,expo=expo,n_iter=50,n.adapt=50)
#' head(fit_APCI_result)
#'
#' \donttest{
#' #if sharing the cohorts only (poisson family)
#' fit_APCI_result2<-fit_APCI(death=death,expo=expo,n_iter=1000,family="poisson",share_gamma=TRUE)
#' head(fit_APCI_result2)
#'
#' #if sharing all alphas, betas, and cohorts (poisson family)
#' fit_APCI_result3<-fit_APCI(death=death,expo=expo,n_iter=1000,family="poisson",share_alpha=TRUE,
#' share_beta=TRUE,share_gamma=TRUE)
#' head(fit_APCI_result3)
#'
#' #if forecast
#' fit_APCI_result4<-fit_APCI(death=death,expo=expo,n_iter=1000,family="poisson",forecast=TRUE)
#' plot_rates_fn(fit_APCI_result4)
#' plot_param_fn(fit_APCI_result4)
#' }
fit_APCI<-function(death,expo,n_iter=10000,family="nb",share_alpha=FALSE,share_beta=FALSE,share_gamma=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
C<-A+T-1
cohorts<-1:C;cohorts_rev<-rev(cohorts)
t<-1:T;t_rev<-rev(t)
prior_mean_beta=rep(0,A)
sigma2_beta<-1
prior_prec_beta<-solve(diag(rep(1,A))*sigma2_beta)
prior_mean_kappa=rep(0,T-2)
sigma2_kappa<-1000
matrix_kappa_A<-diag(rep(1,T))
matrix_kappa_A<-rbind(rep(1/T,T),((1:T)-mean((1:T)))/sd((1:T)),matrix_kappa_A[(2:(T-1)),])
matrix_kappa_B<-matrix_kappa_A%*%t(matrix_kappa_A)
prior_prec_kappa=solve((matrix_kappa_B[3:T,3:T]-matrix_kappa_B[3:T,1:2]%*%solve(matrix_kappa_B[1:2,1:2])%*%matrix_kappa_B[1:2,3:T])*sigma2_kappa)
prior_mean_cohort=rep(0,C-3)
sigma2_cohort<-1
matrix_cohort_A<-diag(rep(1,C))
matrix_cohort_A<-rbind(rep(1/C,C),(cohorts-mean(cohorts))/sd(cohorts),(cohorts^2-mean(cohorts^2))/sd(cohorts^2),matrix_cohort_A[(2:(C-2)),])
matrix_cohort_B<-matrix_cohort_A%*%t(matrix_cohort_A)
prior_prec_cohort=solve((matrix_cohort_B[4:C,4:C]-matrix_cohort_B[4:C,1:3]%*%solve(matrix_cohort_B[1:3,1:3])%*%matrix_cohort_B[1:3,4:C])*sigma2_cohort)
if (!forecast){if (family=="binomial"){
expo_initial<-round(expo$data+0.5*death$data)
if (share_alpha){
if (share_beta){
if (share_gamma){
#8
data<-list(dxt=death$data,ext=expo_initial,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_shareall_jags<-rjags::jags.model(system.file("models/logit_APCI_shareall.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_shareall_jags,vars,n.iter=n_iter,thin=thin)} else {
#5
data<-list(dxt=death$data,ext=expo_initial,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharealpha_sharebeta_jags<-rjags::jags.model(system.file("models/logit_APCI_sharealpha_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sharealpha_sharebeta_jags,vars,n.iter=n_iter,thin=thin)
}
}else{
if (share_gamma){
#6
data<-list(dxt=death$data,ext=expo_initial,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharealpha_sharegamma_jags<-rjags::jags.model(system.file("models/logit_APCI_sharealpha_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sharealpha_sharegamma_jags,vars,n.iter=n_iter,thin=thin)
} else {
#2
data<-list(dxt=death$data,ext=expo_initial,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharealpha_jags<-rjags::jags.model(system.file("models/logit_APCI_sharealpha.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sharealpha_jags,vars,n.iter=n_iter,thin=thin)
}
}
} else {
if (share_beta){
if (share_gamma){
#7
data<-list(dxt=death$data,ext=expo_initial,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharebeta_sharegamma_jags<-rjags::jags.model(system.file("models/logit_APCI_sharebeta_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sharebeta_sharegamma_jags,vars,n.iter=n_iter,thin=thin)
} else {
#3
data<-list(dxt=death$data,ext=expo_initial,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharebeta_jags<-rjags::jags.model(system.file("models/logit_APCI_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sharebeta_jags,vars,n.iter=n_iter,thin=thin)
}
} else {
if (share_gamma){
#4
data<-list(dxt=death$data,ext=expo_initial,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharegamma_jags<-rjags::jags.model(system.file("models/logit_APCI_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sharegamma_jags,vars,n.iter=n_iter,thin=thin)
}else{
#1
data<-list(dxt=death$data,ext=expo_initial,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sep_jags<-rjags::jags.model(system.file("models/logit_APCI_sep.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sep_jags,vars,n.iter=n_iter,thin=thin)
}
}
}
}
if (family=="nb"){
if (share_alpha){
if (share_beta){
if (share_gamma){
#8
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_shareall_jags<-rjags::jags.model(system.file("models/nb_APCI_shareall.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_shareall_jags,vars,n.iter=n_iter,thin=thin)} else {
#5
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharealpha_sharebeta_jags<-rjags::jags.model(system.file("models/nb_APCI_sharealpha_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sharealpha_sharebeta_jags,vars,n.iter=n_iter,thin=thin)
}
}else{
if (share_gamma){
#6
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharealpha_sharegamma_jags<-rjags::jags.model(system.file("models/nb_APCI_sharealpha_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sharealpha_sharegamma_jags,vars,n.iter=n_iter,thin=thin)
} else {
#2
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharealpha_jags<-rjags::jags.model(system.file("models/nb_APCI_sharealpha.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sharealpha_jags,vars,n.iter=n_iter,thin=thin)
}
}
} else {
if (share_beta){
if (share_gamma){
#7
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharebeta_sharegamma_jags<-rjags::jags.model(system.file("models/nb_APCI_sharebeta_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sharebeta_sharegamma_jags,vars,n.iter=n_iter,thin=thin)
} else {
#3
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharebeta_jags<-rjags::jags.model(system.file("models/nb_APCI_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sharebeta_jags,vars,n.iter=n_iter,thin=thin)
}
} else {
if (share_gamma){
#4
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharegamma_jags<-rjags::jags.model(system.file("models/nb_APCI_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sharegamma_jags,vars,n.iter=n_iter,thin=thin)
}else{
#1
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sep_jags<-rjags::jags.model(system.file("models/nb_APCI_sep.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sep_jags,vars,n.iter=n_iter,thin=thin)
}
}
}
}
if (family=="poisson"){
if (share_alpha){
if (share_beta){
if (share_gamma){
#8
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_shareall_jags<-rjags::jags.model(system.file("models/log_APCI_shareall.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_shareall_jags,vars,n.iter=n_iter,thin=thin)} else {
#5
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharealpha_sharebeta_jags<-rjags::jags.model(system.file("models/log_APCI_sharealpha_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sharealpha_sharebeta_jags,vars,n.iter=n_iter,thin=thin)
}
}else{
if (share_gamma){
#6
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharealpha_sharegamma_jags<-rjags::jags.model(system.file("models/log_APCI_sharealpha_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sharealpha_sharegamma_jags,vars,n.iter=n_iter,thin=thin)
} else {
#2
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharealpha_jags<-rjags::jags.model(system.file("models/log_APCI_sharealpha.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sharealpha_jags,vars,n.iter=n_iter,thin=thin)
}
}
} else {
if (share_beta){
if (share_gamma){
#7
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharebeta_sharegamma_jags<-rjags::jags.model(system.file("models/log_APCI_sharebeta_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sharebeta_sharegamma_jags,vars,n.iter=n_iter,thin=thin)
} else {
#3
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharebeta_jags<-rjags::jags.model(system.file("models/log_APCI_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sharebeta_jags,vars,n.iter=n_iter,thin=thin)
}
} else {
if (share_gamma){
#4
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharegamma_jags<-rjags::jags.model(system.file("models/log_APCI_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sharegamma_jags,vars,n.iter=n_iter,thin=thin)
}else{
#1
data<-list(dxt=death$data,ext=expo$data,A=A,T=T,C=C,p=p,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=T-2),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sep_jags<-rjags::jags.model(system.file("models/log_APCI_sep.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt=n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sep_jags,vars,n.iter=n_iter,thin=thin)
}
}
}
}
}
if (forecast){
t_forecast<-c(t,t[T]+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){
if (share_gamma){
#8
data<-list(dxt=death_forecast,ext=expo_forecast_initial,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_shareall_alt_forecast_jags<-rjags::jags.model(system.file("models/logit_APCI_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_APCI_shareall_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
} else {
#5
data<-list(dxt=death_forecast,ext=expo_forecast_initial,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharealpha_sharebeta_alt_forecast_jags<-rjags::jags.model(system.file("models/logit_APCI_alt_forecast_sharealpha_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sharealpha_sharebeta_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
}else{
if (share_gamma){
#6
data<-list(dxt=death_forecast,ext=expo_forecast_initial,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharealpha_sharegamma_alt_forecast_jags<-rjags::jags.model(system.file("models/logit_APCI_alt_forecast_sharealpha_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sharealpha_sharegamma_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
} else {
#2
data<-list(dxt=death_forecast,ext=expo_forecast_initial,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharealpha_alt_forecast_jags<-rjags::jags.model(system.file("models/logit_APCI_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_APCI_sharealpha_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
}
} else {
if (share_beta){
if (share_gamma){
#7
data<-list(dxt=death_forecast,ext=expo_forecast_initial,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharebeta_sharegamma_alt_forecast_jags<-rjags::jags.model(system.file("models/logit_APCI_alt_forecast_sharebeta_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sharebeta_sharegamma_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
} else {
#3
data<-list(dxt=death_forecast,ext=expo_forecast_initial,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharebeta_alt_forecast_jags<-rjags::jags.model(system.file("models/logit_APCI_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_APCI_sharebeta_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
} else {
if (share_gamma){
#4
data<-list(dxt=death_forecast,ext=expo_forecast_initial,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sharegamma_alt_forecast_jags<-rjags::jags.model(system.file("models/logit_APCI_alt_forecast_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(logit_APCI_sharegamma_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}else{
#1
data<-list(dxt=death_forecast,ext=expo_forecast_initial,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
logit_APCI_sep_alt_forecast_jags<-rjags::jags.model(system.file("models/logit_APCI_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_APCI_sep_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
}
}
}
if (family=="nb"){
if (share_alpha){
if (share_beta){
if (share_gamma){
#8
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_shareall_alt_forecast_jags<-rjags::jags.model(system.file("models/nb_APCI_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_APCI_shareall_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
} else {
#5
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharealpha_sharebeta_alt_forecast_jags<-rjags::jags.model(system.file("models/nb_APCI_alt_forecast_sharealpha_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sharealpha_sharebeta_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
}else{
if (share_gamma){
#6
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharealpha_sharegamma_alt_forecast_jags<-rjags::jags.model(system.file("models/nb_APCI_alt_forecast_sharealpha_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sharealpha_sharegamma_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
} else {
#2
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharealpha_alt_forecast_jags<-rjags::jags.model(system.file("models/nb_APCI_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_APCI_sharealpha_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
}
} else {
if (share_beta){
if (share_gamma){
#7
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharebeta_sharegamma_alt_forecast_jags<-rjags::jags.model(system.file("models/nb_APCI_alt_forecast_sharebeta_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sharebeta_sharegamma_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
} else {
#3
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharebeta_alt_forecast_jags<-rjags::jags.model(system.file("models/nb_APCI_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_APCI_sharebeta_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
} else {
if (share_gamma){
#4
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sharegamma_alt_forecast_jags<-rjags::jags.model(system.file("models/nb_APCI_alt_forecast_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(nb_APCI_sharegamma_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}else{
#1
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1,phi=100))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma","phi")
nb_APCI_sep_alt_forecast_jags<-rjags::jags.model(system.file("models/nb_APCI_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_APCI_sep_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
}
}
}
if (family=="poisson"){
if (share_alpha){
if (share_beta){
if (share_gamma){
#8
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_shareall_alt_forecast_jags<-rjags::jags.model(system.file("models/log_APCI_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_APCI_shareall_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
} else {
#5
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharealpha_sharebeta_alt_forecast_jags<-rjags::jags.model(system.file("models/log_APCI_alt_forecast_sharealpha_sharebeta.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sharealpha_sharebeta_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
}else{
if (share_gamma){
#6
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharealpha_sharegamma_alt_forecast_jags<-rjags::jags.model(system.file("models/log_APCI_alt_forecast_sharealpha_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sharealpha_sharegamma_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
} else {
#2
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=rep(0,A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharealpha_alt_forecast_jags<-rjags::jags.model(system.file("models/log_APCI_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_APCI_sharealpha_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
}
} else {
if (share_beta){
if (share_gamma){
#7
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharebeta_sharegamma_alt_forecast_jags<-rjags::jags.model(system.file("models/log_APCI_alt_forecast_sharebeta_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sharebeta_sharegamma_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
} else {
#3
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=rep(0,A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharebeta_alt_forecast_jags<-rjags::jags.model(system.file("models/log_APCI_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_APCI_sharebeta_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}
} else {
if (share_gamma){
#4
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest=rep(0,C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sharegamma_alt_forecast_jags<-rjags::jags.model(system.file("models/log_APCI_alt_forecast_sharegamma.jags",package="BayesMoFo"),data=data,inits=inits,n.chain=n.chain,n.adapt = n.adapt,quiet=quiet)
result_jags<-rjags::coda.samples(log_APCI_sharegamma_alt_forecast_jags,vars,n.iter=n_iter,thin=thin)
}else{
#1
data<-list(dxt=death_forecast,ext=expo_forecast,A=A,T=T,C=C,p=p,h=h,prior_mean_beta=prior_mean_beta,prior_prec_beta=prior_prec_beta,t=t,tbar=mean(t),t_forecast=t_forecast,t_rev=t_rev,cohorts=cohorts,cohorts_rev=cohorts_rev)
inits<-function() (list(alpha=matrix(0,nrow=p,ncol=A),beta=matrix(0,nrow=p,ncol=A),kappa_rest_mat=matrix(0,nrow=p,ncol=(T-2)),gamma_rest_mat=matrix(0,nrow=p,ncol=C-3),rho=0.5,i_sigma2_kappa=0.1,rho_gamma=0.1,sigma_gamma=0.1))
vars<-c("q","alpha","beta","kappa","gamma","rho","sigma2_kappa","rho_gamma","sigma2_gamma")
log_APCI_sep_alt_forecast_jags<-rjags::jags.model(system.file("models/log_APCI_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_APCI_sep_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_APCI_sharealpha<-function(death,expo,n_iter=10000,family="nb",n.chain=1,thin=1,n.adapt=1000,forecast=FALSE,h=5,quiet=FALSE){
fit_APCI(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_APCI_sharebeta<-function(death,expo,n_iter=10000,family="nb",n.chain=1,thin=1,n.adapt=1000,forecast=FALSE,h=5,quiet=FALSE){
fit_APCI(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_APCI_sharegamma<-function(death,expo,n_iter=10000,family="nb",n.chain=1,thin=1,n.adapt=1000,forecast=FALSE,h=5,quiet=FALSE){
fit_APCI(death=death,expo=expo,share_gamma = TRUE,n_iter=n_iter,family=family,n.chain=n.chain,thin=thin,n.adapt=n.adapt,quiet=quiet,forecast=forecast,h=h)
}
fit_APCI_sharealpha_sharebeta<-function(death,expo,n_iter=10000,family="nb",n.chain=1,thin=1,n.adapt=1000,forecast=FALSE,h=5,quiet=FALSE){
fit_APCI(death=death,expo=expo,share_alpha = TRUE,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_APCI_sharealpha_sharegamma<-function(death,expo,n_iter=10000,family="nb",n.chain=1,thin=1,n.adapt=1000,forecast=FALSE,h=5,quiet=FALSE){
fit_APCI(death=death,expo=expo,share_alpha = TRUE,share_gamma = TRUE,n_iter=n_iter,family=family,n.chain=n.chain,thin=thin,n.adapt=n.adapt,quiet=quiet,forecast=forecast,h=h)
}
fit_APCI_sharebeta_sharegamma<-function(death,expo,n_iter=10000,family="nb",n.chain=1,thin=1,n.adapt=1000,forecast=FALSE,h=5,quiet=FALSE){
fit_APCI(death=death,expo=expo,share_beta = TRUE,share_gamma = TRUE,n_iter=n_iter,family=family,n.chain=n.chain,thin=thin,n.adapt=n.adapt,quiet=quiet,forecast=forecast,h=h)
}
fit_APCI_shareall<-function(death,expo,n_iter=10000,family="nb",n.chain=1,thin=1,n.adapt=1000,forecast=FALSE,h=5,quiet=FALSE){
fit_APCI(death=death,expo=expo,share_alpha = TRUE,share_beta = TRUE,share_gamma = TRUE,n_iter=n_iter,family=family,n.chain=n.chain,thin=thin,n.adapt=n.adapt,quiet=quiet,forecast=forecast,h=h)
}
Try the BayesMoFo package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.