fit_MLiLee: A function to fit the stochastic mortality model MLiLee.

In BayesMoFo: Bayesian Mortality Forecasting

A function to fit the stochastic mortality model MLiLee.

Description

Carry out Bayesian estimation of the stochastic mortality model MLiLee (Li and Lee, 2005). Note that if the number of strata is one, results from this model are essentially the same as the Lee-Carter model, fit_LC().

Usage

fit_MLiLee(
  death,
  expo,
  n_iter = 10000,
  family = "nb",
  share_alpha = 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.
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="log", 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}k_{t,p}+B_xK_t ,

where 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}k_{t,p}+B_xK_t ,

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}k_{t,p}+B_xK_t ,

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} b_{x} = 1, \sum_{t,p} k_{t,p} = 0 .

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

a_{x}+b_{x,p}k_{t,p}+B_xK_t .

If forecast=TRUE, then a time series model (an AR(1) with linear drift) will be fitted on both k_{t,p} and K_t as follows:

k_{t,p} = \eta^k_1+\eta^k_2 t +\rho_k (k_{t-1,p}-(\eta^k_1+\eta^k_2 (t-1))) + \epsilon^k_{t,p} \text{ for }p=1,\ldots,P-1 \text{ and } t=1,\ldots,T,

and

K_{t} = \eta^K_1+\eta^K_2 t +\rho^K (K_{t-1}-(\eta^K_1+\eta^K_2 (t-1))) + \epsilon^K_{t} \text{ for }t=1,\ldots,T,

where \epsilon^k_{t,p}\sim N(0,\sigma_k^2), \epsilon^K_{t}\sim N(0,\sigma_K^2), while \eta^k_1,\eta^k_2,\rho_k,\sigma_k^2, \eta^K_1,\eta^K_2,\rho_K,\sigma_K^2 are additional parameters to be estimated. In principle, there are many other options for forecasting the mortality time trend. But currently, we assume that this serves as a general purpose forecasting model for simplicity.

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

Li N., & Lee R. (2005). Coherent mortality forecasts for a group of populations: an extension of the Lee-Carter method. Demography. 42(3):575-94. \Sexpr[results=rd]{tools:::Rd_expr_doi("https://doi.org/10.1353/dem.2005.0021")}

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. \Sexpr[results=rd]{tools:::Rd_expr_doi("https://doi.org/10.1016/j.insmatheco.2017.09.023")}

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_MLiLee_result<-fit_MLiLee(death=death,expo=expo,n_iter=50,n.adapt=50)
head(fit_MLiLee_result)


#if sharing the alphas (poisson family)
fit_MLiLee_result2<-fit_MLiLee(death=death,expo=expo,n_iter=1000,family="poisson",share_alpha=TRUE)
head(fit_MLiLee_result2)

#if forecast (poisson family)
fit_MLiLee_result3<-fit_MLiLee(death=death,expo=expo,n_iter=1000,family="poisson",forecast=TRUE)
plot_rates_fn(fit_MLiLee_result3)
plot_param_fn(fit_MLiLee_result3)

BayesMoFo documentation built on Aug. 11, 2025, 1:07 a.m.