fit_LC: A function to fit the stochastic mortality model by Lee and...
In BayesMoFo: Bayesian Mortality Forecasting
fit_LC R Documentation
A function to fit the stochastic mortality model by Lee and Carter (1992).
Description
Carry out Bayesian estimation of the stochastic mortality model, the Lee-Carter model (Lee and Carter, 1992).
Note that this model is equivalent to "M1" under the formulation of Cairns et al. (2009).
Usage
fit_LC(
death,
expo,
n_iter = 10000,
family = "nb",
share_alpha = FALSE,
share_beta = FALSE,
n.chain = 1,
thin = 1,
n.adapt = 1000,
forecast = FALSE,
h = 5,
quiet = FALSE
)
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.
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}k_{t,p} ,
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} ,
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} ,
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,p} = 1, \sum_{t} k_{t,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}k_{t,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}k_{t,p} .
If forecast=TRUE, then a time series model (an AR(1) with linear drift) will be fitted on k_{t,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,
where \epsilon_{t,p}\sim N(0,\sigma_k^2), \eta_1,\eta_2,\rho,\sigma_k^2 are additional parameters to be estimated.
In principle, there are many other options for forecasting the mortality time trend. But currently, we assume that this serves as a general purpose forecasting model for simplicity.
Value
A list with components:
post_sampleAn mcmc.list object containing the posterior samples generated.
paramA vector of character strings describing the names of model parameters.
deathThe death data that was used.
expoThe expo data that was used.
familyThe family function used.
forecastA logical value indicating if forecast has been performed.
hThe 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. \Sexpr[results=rd]{tools:::Rd_expr_doi("https://doi.org/10.1080/10920277.2009.10597538")}
Lee, R. D., and Carter, L. R. (1992). Modeling and Forecasting U.S. Mortality. Journal of the American Statistical Association, 87(419), 659–671. \Sexpr[results=rd]{tools:::Rd_expr_doi("https://doi.org/10.1080/01621459.1992.10475265")}
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")}
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. \Sexpr[results=rd]{tools:::Rd_expr_doi("https://doi.org/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_LC_result<-fit_LC(death=death,expo=expo,n_iter=50,n.adapt=50)
head(fit_LC_result)
#fit the model (poisson family)
fit_LC_result<-fit_LC(death=death,expo=expo,n_iter=1000,family="poisson")
head(fit_LC_result)
#if sharing the betas
fit_LC_result2<-fit_LC(death=death,expo=expo,n_iter=1000,family="poisson",share_beta=TRUE)
head(fit_LC_result2)
#if sharing both alphas and betas
fit_LC_result3<-fit_LC(death=death,expo=expo,n_iter=1000,family="poisson",
share_alpha=TRUE,share_beta=TRUE)
head(fit_LC_result3)
#if forecast
fit_LC_result4<-fit_LC(death=death,expo=expo,n_iter=1000,family="poisson",forecast=TRUE)
plot_rates_fn(fit_LC_result4)
plot_param_fn(fit_LC_result4)
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.
| fit_LC | R Documentation |
A function to fit the stochastic mortality model by Lee and Carter (1992).
Description
Carry out Bayesian estimation of the stochastic mortality model, the Lee-Carter model (Lee and Carter, 1992). Note that this model is equivalent to "M1" under the formulation of Cairns et al. (2009).
Usage
fit_LC(
death,
expo,
n_iter = 10000,
family = "nb",
share_alpha = FALSE,
share_beta = FALSE,
n.chain = 1,
thin = 1,
n.adapt = 1000,
forecast = FALSE,
h = 5,
quiet = FALSE
)
Arguments
death |
death data that has been formatted through the function |
expo |
expo data that has been formatted through the function |
n_iter |
number of iterations to run. Default is |
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 |
share_beta |
a logical value indicating if |
n.chain |
number of parallel chains for the model. |
thin |
thinning interval for monitoring purposes. |
n.adapt |
the number of iterations for adaptation. See |
forecast |
a logical value indicating if forecast is to be performed (default is |
h |
a numeric value giving the number of years to forecast. Default is |
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}k_{t,p} ,
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} ,
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} ,
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,p} = 1, \sum_{t} k_{t,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}k_{t,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}k_{t,p} .
If forecast=TRUE, then a time series model (an AR(1) with linear drift) will be fitted on k_{t,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,
where \epsilon_{t,p}\sim N(0,\sigma_k^2), \eta_1,\eta_2,\rho,\sigma_k^2 are additional parameters to be estimated.
In principle, there are many other options for forecasting the mortality time trend. But currently, we assume that this serves as a general purpose forecasting model for simplicity.
Value
A list with components:
post_sampleAn
mcmc.listobject containing the posterior samples generated.paramA vector of character strings describing the names of model parameters.
deathThe death data that was used.
expoThe expo data that was used.
familyThe family function used.
forecastA logical value indicating if forecast has been performed.
hThe 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. \Sexpr[results=rd]{tools:::Rd_expr_doi("https://doi.org/10.1080/10920277.2009.10597538")}
Lee, R. D., and Carter, L. R. (1992). Modeling and Forecasting U.S. Mortality. Journal of the American Statistical Association, 87(419), 659–671. \Sexpr[results=rd]{tools:::Rd_expr_doi("https://doi.org/10.1080/01621459.1992.10475265")}
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")}
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. \Sexpr[results=rd]{tools:::Rd_expr_doi("https://doi.org/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_LC_result<-fit_LC(death=death,expo=expo,n_iter=50,n.adapt=50)
head(fit_LC_result)
#fit the model (poisson family)
fit_LC_result<-fit_LC(death=death,expo=expo,n_iter=1000,family="poisson")
head(fit_LC_result)
#if sharing the betas
fit_LC_result2<-fit_LC(death=death,expo=expo,n_iter=1000,family="poisson",share_beta=TRUE)
head(fit_LC_result2)
#if sharing both alphas and betas
fit_LC_result3<-fit_LC(death=death,expo=expo,n_iter=1000,family="poisson",
share_alpha=TRUE,share_beta=TRUE)
head(fit_LC_result3)
#if forecast
fit_LC_result4<-fit_LC(death=death,expo=expo,n_iter=1000,family="poisson",forecast=TRUE)
plot_rates_fn(fit_LC_result4)
plot_param_fn(fit_LC_result4)
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.