Package 'BayesMoFo'

BayesMoFo: Bayesian Mortality Forecasting

Description

The BayesMoFo package performs Bayesian inference in a wide variety of mortality models (see vignette for a full list of supported models). Inference is performed using JAGS via the rjags interface.

Details

Additional author and maintainer information is available in the DESCRIPTION file.

Author(s)

Maintainer: Jackie Siaw Tze Wong [email protected] (ORCID)

Other contributors:

• Alex Diana [email protected] [contributor]

• Aniketh Pittea [email protected] [contributor]

---

A function to check for overall convergence of fitted stochastic mortality models for monitoring purposes

Description

Compute several common MCMC convergence diagnostics of parameters fitted under stochastic mortality models using posterior samples stored in "fit_result" object.

Usage

converge_diag_fn(result, tol = 0.2, plot_gelman = FALSE, plot_geweke = FALSE)

Arguments

result	object of type either "fit_result" or "BayesMoFo".
tol	A numeric value between 0 and 1 specifying the tolerance percentage of the convergence diagnostics, i.e. if the proportion of convergence diagnostic checks/tests (either Gelman's or Heidel's) exceeds tol, warning messages will be triggered. Default is tol=0.20.
plot_gelman	A logical value indicating whether to produce a plot of Gelman's R statistic, PSRF ("potential scale reduction factor") for visualisation.
plot_geweke	A logical value indicating whether to produce a plot of Geweke's statistic for visualisation.

Value

A message of recommendations and (possibly) convergence diagnostics plots. Additionally, a list with components:

gelman_diag

A gelman.diag object which is a list containing the estimated R statistic (PSRF) along with the upper confidence intervals, and also the multivariate PSRF (Brooks and Gelman, 1998). See Gelman and Rubin (1992) for more details.

geweke_diag

A geweke.diag object which is a list containing the estimated Geweke's Z scores and the corresponding fractions used for equality of means test. See Geweke (1992) for more details.

heidel_diag

A heidel.diag object which is a matrix containing results from both Stationarity and Half-width tests. See Heidelberger and Welch (1981), Heidelberger and Welch (1983) for more details.

n

The total number of posterior samples (if more than one chain were generated, they are merged into one long chain).

References

Andrew Gelman, Donald B. Rubin. (1992). "Inference from Iterative Simulation Using Multiple Sequences," Statistical Science, Statist. Sci. 7(4), 457-472. doi: 10.1214/ss/1177011136

Brooks, SP. and Gelman, A. (1998). General methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics, 7, 434-455. doi:10.1080/10618600.1998.10474787

Geweke, J. (1992) Evaluating the accuracy of sampling-based approaches to calculating posterior moments. In Bayesian Statistics 4 (ed JM Bernado, JO Berger, AP Dawid and AFM Smith). Clarendon Press, Oxford, UK. doi:10.21034/sr.148

Heidelberger P and Welch PD. (1981). A spectral method for confidence interval generation and run length control in simulations. Comm. ACM. 24, 233-245. doi:10.1145/358598.358630

Heidelberger P and Welch PD. (1983) Simulation run length control in the presence of an initial transient. Opns Res., 31, 1109-44. doi:10.1287/opre.31.6.1109

Examples

#load and prepare data
data("dxt_array_product");data("Ext_array_product")
death<-preparedata_fn(dxt_array_product,strat_name = c("ACI","DB","SCI"),ages=35:65)
expo<-preparedata_fn(Ext_array_product,strat_name = c("ACI","DB","SCI"),ages=35:65)

#fit any mortality model
runBayesMoFo_result<-runBayesMoFo(death=death,expo=expo,models="M1A",n_iter=1000,
  family="poisson",n.chain=5,thin=1)

#compute convergence diagnostics (with plots)
converge_runBayesMoFo_result<-converge_diag_fn(runBayesMoFo_result,
  plot_gelman=TRUE,plot_geweke=TRUE)

#Gelman's R (PSRF)
converge_runBayesMoFo_result$gelman_diag

#Geweke's Z statistics
converge_runBayesMoFo_result$geweke_diag$z

#Heidel's Stationary and Halfwidth Tests
converge_runBayesMoFo_result$heidel_diag

---

A function to assess convergence of the posterior sampling of fitted parameters for monitoring purposes

Description

Produce several convergence diagnostic tools (e.g. trace/density/acf plots and effective sample sizes) from the posterior samples of fitted parameters.

Usage

converge_diag_param_fn(
  result,
  plot_params = NULL,
  trace = TRUE,
  density = TRUE,
  acf_plot = FALSE,
  ESS_all = FALSE
)

Arguments

result	object of type either "fit_result" or "BayesMoFo".
plot_params	A vector of character strings specifying which set of parameters to plot for visualisation. If not specified, a random selection of the parameters will be included in the plots (see fit_result$param or runBayesMoFo_result$result$best$param for a full list) will be chosen. Note that a specific combination of alpha, beta, kappa, and gamma, is to be plotted, then users need to specify the exact indices of them, e.g. plot_params="gamma[1,2]". Otherwise, only three randomly selected of them will be plotted. To see a complete list of parameters, i.e. colnames(fit_result$post_sample[[1]])[!startsWith(colnames(fit_result$post_sample[[1]]),"q[")] or colnames(runBayesMoFo_result$result$best$post_sample[[1]])[!startsWith(colnames(runBayesMoFo_result$result$best$post_sample[[1]]),"q[")].
trace	A logical value to indicate if trace plots of posterior samples of death rates should be shown (default) or suppressed (e.g. to aid visibility).
density	A logical value to indicate if density plots of posterior samples of death rates should be shown (default) or suppressed (e.g. to aid visibility).
acf_plot	A logical value to indicate if auto-correlation plots should be shown or suppressed (default).
ESS_all	A logical value indicating if effective sample sizes are to be computed for all parameters. The default is FALSE where only chosen parameters will be evaluated, if TRUE all parameters will be assessed.

Value

Some convergence-related plots of posterior samples of fitted parameters.

ESS

The effective sample sizes of the chosen parameters.

Examples

#load and prepare data
data("dxt_array_product");data("Ext_array_product")
death<-preparedata_fn(dxt_array_product,strat_name = c("ACI","DB","SCI"),ages=35:65)
expo<-preparedata_fn(Ext_array_product,strat_name = c("ACI","DB","SCI"),ages=35:65)

#fit any mortality model
runBayesMoFo_result<-runBayesMoFo(death=death,expo=expo,models="APCI")

#default plot
converge_runBayesMoFo_result<-converge_diag_param_fn(runBayesMoFo_result)

#ESS
converge_runBayesMoFo_result$ESS

#plot specific parameters 

runBayesMoFo_result$result$best$param  #run this line to check parameters of the model
converge_diag_param_fn(runBayesMoFo_result,plot_params=c("rho","sigma2_kappa","beta")) 
#note only three betas were plotted

colnames(runBayesMoFo_result$result$best$post_sample[[1]])[!startsWith(
  colnames(runBayesMoFo_result$result$best$post_sample[[1]]),"q[")]  
#run the above line to check full list of parameters of the model
converge_diag_param_fn(runBayesMoFo_result,plot_params=c("beta[1,2]","gamma[3,2]")) 

#ACF plot 
converge_diag_param_fn(runBayesMoFo_result,plot_params=c("beta[1,2]","gamma[3,2]"),
trace=FALSE,density=FALSE,acf_plot=TRUE)

---

A function to assess convergence of the posterior sampling of death rates for monitoring purposes

Description

Produce several convergence diagnostic tools (e.g. trace/density/acf plots and effective sample sizes) from the posterior samples of death rates.

Usage

converge_diag_rates_fn(
  result,
  plot_ages = NULL,
  plot_years = NULL,
  plot_strata = NULL,
  trace = TRUE,
  density = TRUE,
  acf_plot = FALSE,
  ESS_all = FALSE
)

Arguments

result	object of type either "fit_result" or "BayesMoFo".
plot_ages	A numeric vector specifying which range of ages to plot for visualisation. If not specified, a random selection of ages that were used to fit the model (i.e. fit_result$death$ages) will be chosen.
plot_years	A numeric vector specifying which range of years to plot for visualisation. If not specified, a random selection of years that were used to fit the model (i.e. fit_result$death$years) will be chosen.
plot_strata	A vector of character strings specifying which strata to plot for visualisation. If not specified, a random selection of the strata used to fit the model (i.e. fit_result$death$strat_name) will be chosen.
trace	A logical value to indicate if trace plots of posterior samples of death rates should be shown (default) or suppressed (e.g. to aid visibility).
density	A logical value to indicate if density plots of posterior samples of death rates should be shown (default) or suppressed (e.g. to aid visibility).
acf_plot	A logical value to indicate if auto-correlation plots should be shown or suppressed (default).
ESS_all	A logical value indicating if effective sample sizes are to be computed for all rates. The default is FALSE where only chosen rates will be evaluated, if TRUE all rates will be assessed.

Value

Some convergence-related plots of posterior samples of mortality rates.

ESS

The effective sample sizes of the chosen parameters.

Examples

#load and prepare data
data("dxt_array_product");data("Ext_array_product")
death<-preparedata_fn(dxt_array_product,strat_name = c("ACI","DB","SCI"),ages=35:65)
expo<-preparedata_fn(Ext_array_product,strat_name = c("ACI","DB","SCI"),ages=35:65)

#fit any mortality model
runBayesMoFo_result<-runBayesMoFo(death=death,expo=expo,models="APCI")

#default plot
converge_runBayesMoFo_result<-converge_diag_rates_fn(runBayesMoFo_result)

#ESS
converge_runBayesMoFo_result$ESS

#plot by age and changing pre-specified arguments 
converge_diag_rates_fn(runBayesMoFo_result,plot_ages=c(40,50,60),plot_years=c(2017,2020))

#ACF plot 
converge_diag_rates_fn(runBayesMoFo_result,plot_ages=c(40,50,60),plot_years=c(2017,2020),
trace=FALSE,density=FALSE,acf_plot=TRUE)

---

Sample mortality data stratified by insurance products

Description

This is a sample data set used for demonstration purposes.

Usage

data("data_summarised")

Format

A data frame with 1278 rows of observations and 9 variables:

Product

Character. The name of the insurance product associated with the observation. There are in total 4 types of products considered in the dataset:
"ACI": ;
"DB": ;
"SCI": ;
"Annuities": Note that this product contains a lot of missing values.

Age

Numeric. The claim age $x$ associated with the observation, ranging between 18-100.

Year

Numeric. The claim year $t$ associated with the observation, spanning years 2016-2020.

Exposure

Numeric. The central exposure to risk, $E_x^c$, associated with the observation.

Claim

Numeric. The number of claims ("deaths") associated with the observation.

ExpClaim

Numeric. The expected number of claims associated with the observation.

Qx

Numeric. The crude mortality rate associated with the observation. It can be computed as $\frac{\text{Claim}}{\text{Exposure}}$.

ExpQx

Numeric. The expected crude mortality rate associated with the observation. It can be computed as $\frac{\text{ExpClaim}}{\text{Exposure}}$.

StdQx

Numeric. The standard deviation of the crude mortality rate associated with the observation. It can be computed as $\sqrt{\frac{\text{Qx} (1-\text{Qx})}{\text{Exposure}}}$.

Examples

data("data_summarised")
str(data_summarised)
head(data_summarised)

#extracting a subset of the data (3 products)
data_summarised[data_summarised$Product==c("ACI","DB","SCI"),]

#extracting a subset of the data (ages 35-65)
data_summarised[(data_summarised$Age>=35 & data_summarised$Age<=65),]

---

A function to calculate the DIC of stochastic mortality models using posterior samples

Description

Compute the Deviance Information Criterion (DIC) of stochastic mortality models using posterior samples stored in "fit_result" object.

Usage

DIC_fn(result)

Arguments

result

object of type either "fit_result" or "BayesMoFo".

Value

A numeric value representing the DIC of the mortality model.

Examples

#load and prepare data
data("dxt_array_product");data("Ext_array_product")
death<-preparedata_fn(dxt_array_product,strat_name = c("ACI","DB","SCI"),ages=35:65)
expo<-preparedata_fn(Ext_array_product,strat_name = c("ACI","DB","SCI"),ages=35:65)

#a toy example
runBayesMoFo_result<-runBayesMoFo(death=death,expo=expo,models="M1A",n_iter=50,n.adapt=50)

#compute the DIC
DIC_fn(runBayesMoFo_result)

---

Sample mortality data stratified by country

Description

This is a sample data set used for demonstration purposes. They consist of two 3-dimensional arrays, one for number of deaths (dxt_array_country) and another for the the corresponding central exposures to risk (Ext_array_country).

Usage

data("Ext_array_country")

data("dxt_array_country")

Format

An object of class array of dimension 5 x 96 x 50.

A 3-dimensional data array (dim 1: strata, dim 2: ages, dim 3: years) containing 5 countries, 96 ages, and 50 years:

Strata

Character. Indicate the country (in total 5) origin of the data, labelled as:
"AUS": Australia;
"ITALY": Italy;
"JAPAN": Japan;
"UK": United Kingdom ;
"US": United States.

Ages

Numeric. Ages at deaths, ranging between 0-95.

Years

Numeric. Years at deaths, spanning years 1951-2000.

Examples

##Load exposure data
data("Ext_array_country")
str(Ext_array_country)
head(Ext_array_country)

#extracting a subset of the data (2 countries, ages 0-90, years 1961-2000)
Ext_array_country[c("AUS","JAPAN"),as.character(0:90),as.character(1961:2000)]

##Load death data
data("dxt_array_country")
str(dxt_array_country)
head(dxt_array_country)

#extracting a subset of the data (2 countries, ages 0-90, years 1961-2000)
dxt_array_country[c("AUS","JAPAN"),as.character(0:90),as.character(1961:2000)]

---

Sample mortality data stratified by insurance products

Description

This is a sample data set used for demonstration purposes. They consist of two 3-dimensional arrays, one for number of deaths (dxt_array_product) and another for the the corresponding central exposures to risk (Ext_array_product).

Usage

data("Ext_array_product")

data("dxt_array_product")

Format

An object of class array of dimension 4 x 83 x 5.

A 3-dimensional data array (dim 1: strata, dim 2: ages, dim 3: years) with 4 strata (products), 83 ages, and 5 years:

Strata/Products

Character. Names of the insurance products considered in the dataset. There are in total 4 products:
"ACI": ;
"DB": ;
"SCI": ;
"Annuities": Note that this product contains a lot of missing values.

Ages

Numeric. Ages at claims, ranging between 18-100.

Years

Numeric. Years at claims, spanning years 2016-2020.

Examples

##Load exposure data
data("Ext_array_product")
str(Ext_array_product)
head(Ext_array_product)

#extracting a subset of the data (3 products, ages 35-65, years 2016-2020)
Ext_array_product[c("ACI","DB","SCI"),as.character(35:65),as.character(2016:2020)]

##Load death data
data("dxt_array_product")
str(dxt_array_product)
head(dxt_array_product)

#extracting a subset of the data (3 products, ages 35-65, years 2016-2020)
dxt_array_product[c("ACI","DB","SCI"),as.character(35:65),as.character(2016:2020)]

---

Sample mortality data stratified by gender/sex in the UK

Description

This is a sample data set used for demonstration purposes. They consist of two 3-dimensional arrays, one for number of deaths (dxt_array_sex) and another for the the corresponding central exposures to risk (Ext_array_sex).

Usage

data("Ext_array_sex")

data("dxt_array_sex")

Format

An object of class array of dimension 2 x 111 x 181.

A 3-dimensional data array (dim 1: strata, dim 2: ages, dim 3: years) with 2 strata (males and females), 111 ages, and 181 years:

Strata

Character. Indicate the gender/sex, labelled as:
"Male": ;
"Female": ..

Ages

Numeric. Ages at deaths, ranging between 0-109, and 110+.

Years

Numeric. Years at deaths, spanning years 1841-2021.

Examples

##Load exposure data
data("Ext_array_sex")
str(Ext_array_sex)
head(Ext_array_sex)

#extracting a subset of the data (2 genders, ages 0-100, years 1961-2000)
Ext_array_sex[c("Male","Female"),as.character(0:100),as.character(1961:2000)]

##Load death data
data("dxt_array_sex")
str(dxt_array_sex)
head(dxt_array_sex)

#extracting a subset of the data (2 genders, ages 0-100, years 1961-2000)
dxt_array_sex[c("Male","Female"),as.character(0:100),as.character(1961:2000)]

---

A function to fit the APCI stochastic mortality model.

Description

Carry out Bayesian estimation of the stochastic mortality model, the Age-Period-Cohort-Improvement (APCI) model. Note that this model has been studied extensively by Richards et al. (2019) and Wong et al. (2023).

Usage

fit_APCI(
  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
)

Arguments

death	death data that has been formatted through the function preparedata_fn.
expo	expo data that has been formatted through the function preparedata_fn.
n_iter	number of iterations to run. Default is n_iter=10000.
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.
share_alpha	a logical value indicating if $a_{x,p}$ should be shared across all strata (see details below). Default is FALSE.
share_beta	a logical value indicating if $b_{x,p}$ should be shared across all strata (see details below). Default is FALSE.
share_gamma	a logical value indicating if the cohort parameter $\gamma_{x,p}$ should be shared across all strata (see details below). Default is FALSE.
n.chain	number of parallel chains for the model.
thin	thinning interval for monitoring purposes.
n.adapt	the number of iterations for adaptation. See ?rjags::adapt for details.
forecast	a logical value indicating if forecast is to be performed (default is FALSE). See below for details.
h	a numeric value giving the number of years to forecast. Default is h=5.
quiet	if TRUE then messages generated during compilation will be suppressed, as well as the progress bar during adaptation.

Details

The model can be described mathematically as follows: If family="poisson", then

$d_{x,t,p} \sim \text{Poisson}(E^c_{x,t,p} m_{x,t,p}) ,$

$\log(m_{x,t,p})=a_{x,p}+b_{x,p}(t-\bar{t})+k_{t,p} + \gamma_{c,p} ,$

where $c=t-x$ is the cohort index, $\bar{t}$ is the mean of $t$, $d_{x,t,p}$ represents the number of deaths at age $x$ in year $t$ of stratum $p$, while $E^c_{x,t,p}$ and $m_{x,t,p}$ represents respectively the corresponding central exposed to risk and central mortality rate at age $x$ in year $t$ of stratum $p$. Similarly, if family="nb", then a negative binomial distribution is fitted, i.e.

$d_{x,t,p} \sim \text{Negative-Binomial}(\phi,\frac{\phi}{\phi+E^c_{x,t,p} m_{x,t,p}}) ,$

$\log(m_{x,t,p})=a_{x,p}+b_{x,p}(t-\bar{t})+k_{t,p} + \gamma_{c,p} ,$

where $\phi$ is the overdispersion parameter. See Wong et al. (2018). But if family="binomial", then

$d_{x,t,p} \sim \text{Binomial}(E^0_{x,t,p} , q_{x,t,p}) ,$

$\text{logit}(q_{x,t,p})=a_{x,p}+b_{x,p}(t-\bar{t})+k_{t,p} + \gamma_{c,p} ,$

where $q_{x,t,p}$ represents the initial mortality rate at age $x$ in year $t$ of stratum $p$, while $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:

$\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 share_alpha=TRUE, then the additive age-specific parameter is the same across all strata $p$, i. e.

$a_{x}+b_{x,p}(t-\bar{t})+k_{t,p}+ \gamma_{c,p} .$

Similarly if share_beta=TRUE, then the multiplicative age-specific parameter is the same across all strata $p$, i. e.

$a_{x,p}+b_{x}(t-\bar{t})+k_{t,p}+ \gamma_{c,p} .$

Similarly if share_gamma=TRUE, then the cohort parameter is the same across all strata $p$, i. e.

$a_{x,p}+b_{x,p}(t-\bar{t})+k_{t,p}+ \gamma_{c} .$

If forecast=TRUE, then the following time series models are fitted on $k_{t,p}$ and $\gamma_{c,p}$ as follows:

$k_{t,p} = \rho k_{t-1,p} + \epsilon_{t,p} \text{ for }p=1,\ldots,P \text{ and } t=2,\ldots,T,$

$k_{1,p}=\epsilon_{1,p}$

and

$\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,$

$\gamma_2=\gamma_1+\frac{1}{\sqrt{1-\rho_\gamma^2}}\epsilon^\gamma_{2,p},$

$\gamma_1=100\epsilon^\gamma_{1,p}$

where $\epsilon_{t,p}\sim N(0,\sigma_k^2)$ for $t=1,\ldots,T$, $\epsilon^\gamma_{c,p}\sim N(0,\sigma_\gamma^2)$ for $c=1,\ldots,C$, while $\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).

Value

A list with components:

post_sample

An mcmc.list object containing the posterior samples generated.

param

A vector of character strings describing the names of model parameters.

death

The death data that was used.

expo

The expo data that was used.

family

The family function used.

forecast

A logical value indicating if forecast has been performed.

h

The forecast horizon used.

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:10.1017/S1357321718000260

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:10.1016/j.insmatheco.2017.09.023

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:10.1093/jrsssc/qlad021

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)


#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)

---

A function to fit the Age-Period-Cohort (APC) stochastic mortality model.

Description

Carry out Bayesian estimation of the stochastic mortality model, the Age-Period-Cohort (APC) model (Jacobsen et al., 2002). Note that this model is equivalent to "M3" under the formulation of Cairns et al. (2009).

Usage

fit_CBD_M3(
  death,
  expo,
  n_iter = 10000,
  family = "nb",
  share_alpha = FALSE,
  share_gamma = FALSE,
  n.chain = 1,
  thin = 1,
  n.adapt = 1000,
  forecast = FALSE,
  h = 5,
  quiet = FALSE
)

Arguments

death	death data that has been formatted through the function preparedata_fn.
expo	expo data that has been formatted through the function preparedata_fn.
n_iter	number of iterations to run. Default is n_iter=10000.
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.
share_alpha	a logical value indicating if $a_{x,p}$ should be shared across all strata (see details below). Default is FALSE.
share_gamma	a logical value indicating if the cohort parameter $\gamma_{x,p}$ should be shared across all strata (see details below). Default is FALSE.
n.chain	number of parallel chains for the model.
thin	thinning interval for monitoring purposes.
n.adapt	the number of iterations for adaptation. See ?rjags::adapt for details.
forecast	a logical value indicating if forecast is to be performed (default is FALSE). See below for details.
h	a numeric value giving the number of years to forecast. Default is h=5.
quiet	if TRUE then messages generated during compilation will be suppressed, as well as the progress bar during adaptation.

Details

The model can be described mathematically as follows: If family="poisson", then

$d_{x,t,p} \sim \text{Poisson}(E^c_{x,t,p} m_{x,t,p}) ,$

$\log(m_{x,t,p})=a_{x,p}+k_{t,p} + \gamma_{c,p} ,$

where $c=t-x$ is the cohort index, $d_{x,t,p}$ represents the number of deaths at age $x$ in year $t$ of stratum $p$, while $E^c_{x,t,p}$ and $m_{x,t,p}$ represents respectively the corresponding central exposed to risk and central mortality rate at age $x$ in year $t$ of stratum $p$. Similarly, if family="nb", then a negative binomial distribution is fitted, i.e.

$d_{x,t,p} \sim \text{Negative-Binomial}(\phi,\frac{\phi}{\phi+E^c_{x,t,p} m_{x,t,p}}) ,$

$\log(m_{x,t,p})=a_{x,p}+k_{t,p} + \gamma_{c,p} ,$

where $\phi$ is the overdispersion parameter. See Wong et al. (2018). But if family="binomial", then

$d_{x,t,p} \sim \text{Binomial}(E^0_{x,t,p} , q_{x,t,p}) ,$

$\text{logit}(q_{x,t,p})=a_{x,p}+k_{t,p} + \gamma_{c,p} ,$

where $q_{x,t,p}$ represents the initial mortality rate at age $x$ in year $t$ of stratum $p$, while $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:

$\sum_{t} k_{t,p} = 0, \gamma_{1,p}=0, \sum_{x,t} \gamma_{c,p} = 0 \text{\ \ \ for each stratum }p.$

If share_alpha=TRUE, then the additive age-specific parameter is the same across all strata $p$, i. e.

$a_{x}+k_{t,p}+ \gamma_{c,p} .$

Similarly if share_gamma=TRUE, then the cohort parameter is the same across all strata $p$, i. e.

$a_{x,p}+k_{t,p}+ \gamma_{c} .$

If forecast=TRUE, then the following time series models are fitted on $k_{t,p}$ and $\gamma_{c,p}$ as follows:

$k_{t,p} = \eta_1+\eta_2 t +\rho (k_{t-1,p}-(\eta_1+\eta_2 (t-1))) + \epsilon_{t,p} \text{ for }p=1,\ldots,P \text{ and } t=1,\ldots,T,$

and

$\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,$

$\gamma_2=\gamma_1+\frac{1}{\sqrt{1-\rho_\gamma^2}}\epsilon^\gamma_{2,p},$

$\gamma_1=100\epsilon^\gamma_{1,p}$

where $\epsilon_{t,p}\sim N(0,\sigma_k^2)$ for $t=1,\ldots,T$, $\epsilon^\gamma_{c,p}\sim N(0,\sigma_\gamma^2)$ for $c=1,\ldots,C$, while $\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).

Value

A list with components:

post_sample

An mcmc.list object containing the posterior samples generated.

param

A vector of character strings describing the names of model parameters.

death

The death data that was used.

expo

The expo data that was used.

family

The family function used.

forecast

A logical value indicating if forecast has been performed.

h

The forecast horizon used.

References

Cairns, A. J. G., Blake, D., Dowd, K., Coughlan, G. D., Epstein, D., Ong, A., and Balevich, I. (2009). A Quantitative Comparison of Stochastic Mortality Models Using Data From England and Wales and the United States. North American Actuarial Journal, 13(1), 1–35. doi:10.1080/10920277.2009.10597538

Jacobsen R, Keiding N, and Lynge E. (2002). Long term mortality trends behind low life expectancy of Danish women. J Epidemiol Community Health. 56(3): 205-8. doi: 10.1136/jech.56.3.205

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:10.1016/j.insmatheco.2017.09.023

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:10.1093/jrsssc/qlad021

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_CBD_M3_result<-fit_CBD_M3(death=death,expo=expo,n_iter=50,n.adapt=50)
head(fit_CBD_M3_result)


#if sharing the cohorts only (poisson family)
fit_CBD_M3_result2<-fit_CBD_M3(death=death,expo=expo,n_iter=1000,family="poisson",share_gamma=TRUE)
head(fit_CBD_M3_result2)

#if sharing both alphas and cohorts (poisson family)
fit_CBD_M3_result3<-fit_CBD_M3(death=death,expo=expo,n_iter=1000,family="poisson",
share_alpha=TRUE,share_gamma=TRUE)
head(fit_CBD_M3_result3)

---

A function to fit the stochastic mortality model "M5".

Description

Carry out Bayesian estimation of the stochastic mortality model referred to as "M5" under the formulation of Cairns et al. (2009).

Usage

fit_CBD_M5(
  death,
  expo,
  n_iter = 10000,
  family = "nb",
  n.chain = 1,
  thin = 1,
  n.adapt = 1000,
  forecast = FALSE,
  h = 5,
  quiet = FALSE
)

Arguments

death	death data that has been formatted through the function preparedata_fn.
expo	expo data that has been formatted through the function preparedata_fn.
n_iter	number of iterations to run. Default is n_iter=10000.
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.
n.chain	number of parallel chains for the model.
thin	thinning interval for monitoring purposes.
n.adapt	the number of iterations for adaptation. See ?rjags::adapt for details.
forecast	a logical value indicating if forecast is to be performed (default is FALSE). See below for details.
h	a numeric value giving the number of years to forecast. Default is h=5.
quiet	if TRUE then messages generated during compilation will be suppressed, as well as the progress bar during adaptation.

Details

The model can be described mathematically as follows: If family="poisson", then

$d_{x,t,p} \sim \text{Poisson}(E^c_{x,t,p} m_{x,t,p}) ,$

$\log(m_{x,t,p})=k^1_{t,p} + k^2_{t,p}(x-\bar{x}) ,$

where $\bar{x}$ is the mean of $x$, $d_{x,t,p}$ represents the number of deaths at age $x$ in year $t$ of stratum $p$, while $E^c_{x,t,p}$ and $m_{x,t,p}$ represents respectively the corresponding central exposed to risk and central mortality rate at age $x$ in year $t$ of stratum $p$. Similarly, if family="nb", then a negative binomial distribution is fitted, i.e.

$d_{x,t,p} \sim \text{Negative-Binomial}(\phi,\frac{\phi}{\phi+E^c_{x,t,p} m_{x,t,p}}) ,$

$\log(m_{x,t,p})=k^1_{t,p} + k^2_{t,p}(x-\bar{x}) ,$

where $\phi$ is the overdispersion parameter. See Wong et al. (2018). But if family="binomial", then

$d_{x,t,p} \sim \text{Binomial}(E^0_{x,t,p} , q_{x,t,p}) ,$

$\text{logit}(q_{x,t,p})=k^1_{t,p} + k^2_{t,p}(x-\bar{x}) ,$

where $q_{x,t,p}$ represents the initial mortality rate at age $x$ in year $t$ of stratum $p$, while $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. No constraints are needed. If forecast=TRUE, then the following time series models are fitted on $k^1_{t,p}$ as follows:

$k^1_{t,p} = \eta^1_1+\eta^1_2 t +\rho (k^1_{t-1,p}-(\eta^1_1+\eta^1_2 (t-1))) + \epsilon^1_{t,p} \text{ for }p=1,\ldots,P \text{ and } t=1,\ldots,T,$

where $\epsilon^1_{t,p}\sim N(0,(\sigma^1)_k^2)$ for $t=1,\ldots,T$, while $\eta^1_1,\eta^1_2,\rho^1,(\sigma^1)_k^2$ are additional parameters to be estimated. Similarly for $k^2_{t,p}$. Note that the forecasting models are inspired by Wong et al. (2023).

Value

A list with components:

post_sample

An mcmc.list object containing the posterior samples generated.

param

A vector of character strings describing the names of model parameters.

death

The death data that was used.

expo

The expo data that was used.

family

The family function used.

forecast

A logical value indicating if forecast has been performed.

h

The forecast horizon used.

References

Cairns, A. J. G., Blake, D., Dowd, K., Coughlan, G. D., Epstein, D., Ong, A., and Balevich, I. (2009). A Quantitative Comparison of Stochastic Mortality Models Using Data From England and Wales and the United States. North American Actuarial Journal, 13(1), 1–35. doi:10.1080/10920277.2009.10597538

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:10.1016/j.insmatheco.2017.09.023

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:10.1093/jrsssc/qlad021

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 (poisson family)
#NOTE: This is a toy example, please run it longer in practice.
fit_CBD_M5_result<-fit_CBD_M5(death=death,expo=expo,n_iter=50,n.adapt=50,family="poisson")
head(fit_CBD_M5_result)

---

A function to fit the stochastic mortality model "M6".

Description

Carry out Bayesian estimation of the stochastic mortality model referred to as "M6" under the formulation of Cairns et al. (2009).

Usage

fit_CBD_M6(
  death,
  expo,
  share_gamma = FALSE,
  n_iter = 10000,
  family = "nb",
  n.chain = 1,
  thin = 1,
  n.adapt = 1000,
  forecast = FALSE,
  h = 5,
  quiet = FALSE
)

Arguments

death	death data that has been formatted through the function preparedata_fn.
expo	expo data that has been formatted through the function preparedata_fn.
share_gamma	a logical value indicating if the cohort parameter $\gamma_{c,p}$ should be shared across all strata (see details below). Default is FALSE.
n_iter	number of iterations to run. Default is n_iter=10000.
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.
n.chain	number of parallel chains for the model.
thin	thinning interval for monitoring purposes.
n.adapt	the number of iterations for adaptation. See ?rjags::adapt for details.
forecast	a logical value indicating if forecast is to be performed (default is FALSE). See below for details.
h	a numeric value giving the number of years to forecast. Default is h=5.
quiet	if TRUE then messages generated during compilation will be suppressed, as well as the progress bar during adaptation.

Details

The model can be described mathematically as follows: If family="poisson", then

$d_{x,t,p} \sim \text{Poisson}(E^c_{x,t,p} m_{x,t,p}) ,$

$\log(m_{x,t,p})=k^1_{t,p} + k^2_{t,p}(x-\bar{x}) +\gamma_{c,p} ,$

where $c=t-x$ is the cohort index, $\bar{x}$ is the mean of $x$, $d_{x,t,p}$ represents the number of deaths at age $x$ in year $t$ of stratum $p$, while $E^c_{x,t,p}$ and $m_{x,t,p}$ represents respectively the corresponding central exposed to risk and central mortality rate at age $x$ in year $t$ of stratum $p$. Similarly, if family="nb", then a negative binomial distribution is fitted, i.e.

$d_{x,t,p} \sim \text{Negative-Binomial}(\phi,\frac{\phi}{\phi+E^c_{x,t,p} m_{x,t,p}}) ,$

$\log(m_{x,t,p})=k^1_{t,p} + k^2_{t,p}(x-\bar{x}) +\gamma_{c,p} ,$

where $\phi$ is the overdispersion parameter. See Wong et al. (2018). But if family="binomial", then

$d_{x,t,p} \sim \text{Binomial}(E^0_{x,t,p} , q_{x,t,p}) ,$

$\text{logit}(q_{x,t,p})=k^1_{t,p} + k^2_{t,p}(x-\bar{x}) +\gamma_{c,p} ,$

where $q_{x,t,p}$ represents the initial mortality rate at age $x$ in year $t$ of stratum $p$, while $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:

$\sum_{c} \gamma_{c,p} = \sum_{c} c\gamma_{c,p} = 0 \text{\ \ \ for each stratum }p.$

If share_gamma=TRUE, then the cohort parameter is the same across all strata $p$, i. e.

$k^1_{t,p} + k^2_{t,p}(x-\bar{x}) +\gamma_{c} .$

If forecast=TRUE, then the following time series models are fitted on $k^1_{t,p}$ as follows:

$k^1_{t,p} = \eta^1_1+\eta^1_2 t +\rho (k^1_{t-1,p}-(\eta^1_1+\eta^1_2 (t-1))) + \epsilon^1_{t,p} \text{ for }p=1,\ldots,P \text{ and } t=1,\ldots,T,$

where $\epsilon^1_{t,p}\sim N(0,(\sigma^1)_k^2)$ for $t=1,\ldots,T$, while $\eta^1_1,\eta^1_2,\rho^1,(\sigma^1)_k^2$ are additional parameters to be estimated. Similarly for $k^2_{t,p}$. Also,

$\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,$

$\gamma_2=\gamma_1+\frac{1}{\sqrt{1-\rho_\gamma^2}}\epsilon^\gamma_{2,p},$

$\gamma_1=100\epsilon^\gamma_{1,p}$

where $\epsilon^\gamma_{c,p}\sim N(0,\sigma_\gamma^2)$ for $c=1,\ldots,C$, while $\rho_\gamma,\sigma_\gamma^2$ are additional parameters to be estimated. Note that the forecasting models are inspired by Wong et al. (2023).

Value

A list with components:

post_sample

An mcmc.list object containing the posterior samples generated.

param

A vector of character strings describing the names of model parameters.

death

The death data that was used.

expo

The expo data that was used.

family

The family function used.

forecast

A logical value indicating if forecast has been performed.

h

The forecast horizon used.

References

Cairns, A. J. G., Blake, D., Dowd, K., Coughlan, G. D., Epstein, D., Ong, A., and Balevich, I. (2009). A Quantitative Comparison of Stochastic Mortality Models Using Data From England and Wales and the United States. North American Actuarial Journal, 13(1), 1–35. doi:10.1080/10920277.2009.10597538

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:10.1016/j.insmatheco.2017.09.023

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:10.1093/jrsssc/qlad021

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_CBD_M6_result<-fit_CBD_M6(death=death,expo=expo,n_iter=50,n.adapt=50)
head(fit_CBD_M6_result)


#if sharing the cohorts only (poisson family)
fit_CBD_M6_result2<-fit_CBD_M6(death=death,expo=expo,n_iter=1000,family="poisson",share_gamma=TRUE)
head(fit_CBD_M6_result2)

---

A function to fit the stochastic mortality model "M7".

Description

Carry out Bayesian estimation of the stochastic mortality model referred to as "M7" under the formulation of Cairns et al. (2009).

Usage

fit_CBD_M7(
  death,
  expo,
  share_gamma = FALSE,
  n_iter = 10000,
  family = "nb",
  n.chain = 1,
  thin = 1,
  n.adapt = 1000,
  forecast = FALSE,
  h = 5,
  quiet = FALSE
)

Arguments

death	death data that has been formatted through the function preparedata_fn.
expo	expo data that has been formatted through the function preparedata_fn.
share_gamma	a logical value indicating if the cohort parameter $\gamma_{c,p}$ should be shared across all strata (see details below). Default is FALSE.
n_iter	number of iterations to run. Default is n_iter=10000.
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.
n.chain	number of parallel chains for the model.
thin	thinning interval for monitoring purposes.
n.adapt	the number of iterations for adaptation. See ?rjags::adapt for details.
forecast	a logical value indicating if forecast is to be performed (default is FALSE). See below for details.
h	a numeric value giving the number of years to forecast. Default is h=5.
quiet	if TRUE then messages generated during compilation will be suppressed, as well as the progress bar during adaptation.

Details

The model can be described mathematically as follows: If family="poisson", then

$d_{x,t,p} \sim \text{Poisson}(E^c_{x,t,p} m_{x,t,p}) ,$

$\log(m_{x,t,p})=k^1_{t,p} + k^2_{t,p}(x-\bar{x}) + k^3_{t,p}((x-\bar{x})^2-\hat{\sigma}_x^2) +\gamma_{c,p} ,$

where $c=t-x$ is the cohort index, $\bar{x}$ is the mean of $x$, $\hat{\sigma}_x^2$ is the mean of $(x-\bar{x})^2$, $d_{x,t,p}$ represents the number of deaths at age $x$ in year $t$ of stratum $p$, while $E^c_{x,t,p}$ and $m_{x,t,p}$ represents respectively the corresponding central exposed to risk and central mortality rate at age $x$ in year $t$ of stratum $p$. Similarly, if family="nb", then a negative binomial distribution is fitted, i.e.

$d_{x,t,p} \sim \text{Negative-Binomial}(\phi,\frac{\phi}{\phi+E^c_{x,t,p} m_{x,t,p}}) ,$

$\log(m_{x,t,p})=k^1_{t,p} + k^2_{t,p}(x-\bar{x}) + k^3_{t,p}((x-\bar{x})^2-\hat{\sigma}_x^2) +\gamma_{c,p} ,$

where $\phi$ is the overdispersion parameter. See Wong et al. (2018). But if family="binomial", then

$d_{x,t,p} \sim \text{Binomial}(E^0_{x,t,p} , q_{x,t,p}) ,$

$\text{logit}(q_{x,t,p})=k^1_{t,p} + k^2_{t,p}(x-\bar{x}) + k^3_{t,p}((x-\bar{x})^2-\hat{\sigma}_x^2) +\gamma_{c,p} ,$

where $q_{x,t,p}$ represents the initial mortality rate at age $x$ in year $t$ of stratum $p$, while $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:

$\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 share_gamma=TRUE, then the cohort parameter is the same across all strata $p$, i. e.

$k^1_{t,p} + k^2_{t,p}(x-\bar{x}) + k^3_{t,p}((x-\bar{x})^2-\hat{\sigma}_x^2) +\gamma_{c} .$

If forecast=TRUE, then the following time series models are fitted on $k^1_{t,p}$ as follows:

$k^1_{t,p} = \eta^1_1+\eta^1_2 t +\rho (k^1_{t-1,p}-(\eta^1_1+\eta^1_2 (t-1))) + \epsilon^1_{t,p} \text{ for }p=1,\ldots,P \text{ and } t=1,\ldots,T,$

where $\epsilon^1_{t,p}\sim N(0,(\sigma^1)_k^2)$ for $t=1,\ldots,T$, while $\eta^1_1,\eta^1_2,\rho^1,(\sigma^1)_k^2$ are additional parameters to be estimated. Similarly for $k^2_{t,p}$ and $k^3_{t,p}$. Also,

$\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,$

$\gamma_2=\gamma_1+\frac{1}{\sqrt{1-\rho_\gamma^2}}\epsilon^\gamma_{2,p},$

$\gamma_1=100\epsilon^\gamma_{1,p}$

where $\epsilon^\gamma_{c,p}\sim N(0,\sigma_\gamma^2)$ for $c=1,\ldots,C$, while $\rho_\gamma,\sigma_\gamma^2$ are additional parameters to be estimated. Note that the forecasting models are inspired by Wong et al. (2023).

Value

A list with components:

post_sample

An mcmc.list object containing the posterior samples generated.

param

A vector of character strings describing the names of model parameters.

death

The death data that was used.

expo

The expo data that was used.

family

The family function used.

forecast

A logical value indicating if forecast has been performed.

h

The forecast horizon used.

References

Cairns, A. J. G., Blake, D., Dowd, K., Coughlan, G. D., Epstein, D., Ong, A., and Balevich, I. (2009). A Quantitative Comparison of Stochastic Mortality Models Using Data From England and Wales and the United States. North American Actuarial Journal, 13(1), 1–35. doi:10.1080/10920277.2009.10597538

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:10.1016/j.insmatheco.2017.09.023

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:10.1093/jrsssc/qlad021

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_CBD_M7_result<-fit_CBD_M7(death=death,expo=expo,n_iter=50,n.adapt=50)
head(fit_CBD_M7_result)


#if sharing the cohorts only (poisson family)
fit_CBD_M7_result2<-fit_CBD_M7(death=death,expo=expo,n_iter=1000,family="poisson",share_gamma=TRUE)
head(fit_CBD_M7_result2)

---

A function to fit the stochastic mortality model "M8".

Description

Carry out Bayesian estimation of the stochastic mortality model referred to as "M8" under the formulation of Cairns et al. (2009).

Usage

fit_CBD_M8(
  death,
  expo,
  share_gamma = FALSE,
  n_iter = 10000,
  family = "nb",
  n.chain = 1,
  thin = 1,
  n.adapt = 1000,
  forecast = FALSE,
  h = 5,
  quiet = FALSE
)

Arguments

death	death data that has been formatted through the function preparedata_fn.
expo	expo data that has been formatted through the function preparedata_fn.
share_gamma	a logical value indicating if the cohort parameter $\gamma_{c,p}$ should be shared across all strata (see details below). Default is FALSE.
n_iter	number of iterations to run. Default is n_iter=10000.
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.
n.chain	number of parallel chains for the model.
thin	thinning interval for monitoring purposes.
n.adapt	the number of iterations for adaptation. See ?rjags::adapt for details.
forecast	a logical value indicating if forecast is to be performed (default is FALSE). See below for details.
h	a numeric value giving the number of years to forecast. Default is h=5.
quiet	if TRUE then messages generated during compilation will be suppressed, as well as the progress bar during adaptation.

Details

The model can be described mathematically as follows: If family="poisson", then

$d_{x,t,p} \sim \text{Poisson}(E^c_{x,t,p} m_{x,t,p}) ,$

$\log(m_{x,t,p})=k^1_{t,p} + k^2_{t,p}(x-\bar{x}) +\gamma_{c,p}(c_p-x) ,$

where $c=t-x$ is the cohort index, $\bar{x}$ is the mean of $x$, $c_p$ are stratum-specific parameters, $d_{x,t,p}$ represents the number of deaths at age $x$ in year $t$ of stratum $p$, while $E^c_{x,t,p}$ and $m_{x,t,p}$ represents respectively the corresponding central exposed to risk and central mortality rate at age $x$ in year $t$ of stratum $p$. Similarly, if family="nb", then a negative binomial distribution is fitted, i.e.

$d_{x,t,p} \sim \text{Negative-Binomial}(\phi,\frac{\phi}{\phi+E^c_{x,t,p} m_{x,t,p}}) ,$

$\log(m_{x,t,p})=k^1_{t,p} + k^2_{t,p}(x-\bar{x}) +\gamma_{c,p}(c_p-x) ,$

where $\phi$ is the overdispersion parameter. See Wong et al. (2018). But if family="binomial", then

$d_{x,t,p} \sim \text{Binomial}(E^0_{x,t,p} , q_{x,t,p}) ,$

$\text{logit}(q_{x,t,p})=k^1_{t,p} + k^2_{t,p}(x-\bar{x}) +\gamma_{c,p}(c_p-x) ,$

where $q_{x,t,p}$ represents the initial mortality rate at age $x$ in year $t$ of stratum $p$, while $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:

$\sum_{x,t} \gamma_{c,p} = 0 \text{\ \ \ for each stratum }p.$

If share_gamma=TRUE, then the cohort parameter is the same across all strata $p$, i. e.

$k^1_{t,p} + k^2_{t,p}(x-\bar{x}) +\gamma_{c}(c_p-x) .$

If forecast=TRUE, then the following time series models are fitted on $k^1_{t,p}$ as follows:

$k^1_{t,p} = \eta^1_1+\eta^1_2 t +\rho (k^1_{t-1,p}-(\eta^1_1+\eta^1_2 (t-1))) + \epsilon^1_{t,p} \text{ for }p=1,\ldots,P \text{ and } t=1,\ldots,T,$

where $\epsilon^1_{t,p}\sim N(0,(\sigma^1)_k^2)$ for $t=1,\ldots,T$, while $\eta^1_1,\eta^1_2,\rho^1,(\sigma^1)_k^2$ are additional parameters to be estimated. Similarly for $k^2_{t,p}$. Also,

$\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,$

$\gamma_2=\gamma_1+\frac{1}{\sqrt{1-\rho_\gamma^2}}\epsilon^\gamma_{2,p},$

$\gamma_1=100\epsilon^\gamma_{1,p}$

where $\epsilon^\gamma_{c,p}\sim N(0,\sigma_\gamma^2)$ for $c=1,\ldots,C$, while $\rho_\gamma,\sigma_\gamma^2$ are additional parameters to be estimated. Note that the forecasting models are inspired by Wong et al. (2023).

Value

A list with components:

post_sample

An mcmc.list object containing the posterior samples generated.

param

A vector of character strings describing the names of model parameters.

death

The death data that was used.

expo

The expo data that was used.

family

The family function used.

forecast

A logical value indicating if forecast has been performed.

h

The forecast horizon used.

References

Cairns, A. J. G., Blake, D., Dowd, K., Coughlan, G. D., Epstein, D., Ong, A., and Balevich, I. (2009). A Quantitative Comparison of Stochastic Mortality Models Using Data From England and Wales and the United States. North American Actuarial Journal, 13(1), 1–35. doi:10.1080/10920277.2009.10597538

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:10.1016/j.insmatheco.2017.09.023

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:10.1093/jrsssc/qlad021

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 (poisson family)
#Note that this model sometimes fails numerically.
try(fit_CBD_M8_result<-fit_CBD_M8(death=death,expo=expo,n_iter=1000,family="poisson"))

#if sharing the cohorts only (negative-binomial family)
try(fit_CBD_M8_result2<-fit_CBD_M8(death=death,expo=expo,n_iter=1000,share_gamma=TRUE))

---