---
title: "BayesMoFo-vignette"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{BayesMoFo-vignette}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: BayesMoFo.bib  
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

This vignette is intended to show how an analysis of mortality data would work
using the **BayesMoFo** R package. We start by installing and loading the package.


```{r}
# install package
# Recommended installation
# install.packages("BayesMoFo")

# Development version (use only if needed)
# install.packages("devtools")
# devtools::install_github("jstw1g09/Rpackage-BayesMoFo")

#load package
library(BayesMoFo)
```


The package fits various mortality models to two types of data: **age-period (AP)** data and
**age-period-product (APP)** data. AP data refers to data structured by individuals' age and the calendar period (or year) in which events (such as deaths or exposures) are observed. Each observation corresponds to a specific age and period combination, summarising the number of events and the population at risk for that group. For example, an AP dataset might record the number of deaths and the exposure (population at risk) for individuals aged 50 in the year 2010, aged 51 in 2011, and so on.

APP data extends the AP framework by including an additional stratifying variable - the "product" - which in general can be any other stratifying variable (such as
cause of death, country, deprivation level, gender/sex, geographical location/region, insurance product, marital status, socioeconomic group, smoking behaviour, etc.). Each observation in an APP dataset corresponds to a specific combination of age, period, and product, capturing the number of events and the population at risk for that stratum. For example, an APP dataset might record deaths and exposures for individuals aged 50 in 2010, by cause of death or by insurance product type.

We will separately consider each analysis. 

# Age-Period (AP) data

## Data preparation

The package accepts several types of data formats. For AP data, users can supply a data-frame, a 3-dimensional (3D) data array, or a data matrix.

## Data supplied as a data-frame

The data needs to be formatted as a `data.frame` with columns name *Age*, *Year*, *Deaths*
and *Exposures*. An example is provided below.

You can access a dataset in this format for comparison by running

```{r}
data(uk_mortalitydata)
```

The first few lines of the data look like the following:
```{r}
head(uk_mortalitydata, n = 20)
```

Next, we need to format the data in the format necessary for the function **runBayesMoFo** to work properly. 
To do this, we pass the dataset (separately for death and exposure) to the function **preparedata_fn**, which takes
the arguments *ages*, *years*, and *data*. In the case of deaths, we pass the data 
using the column  *Age*, *Year*, and *Claim*; and in the case of exposures, we use *Exposure* in place of *Claim*.

```{r, echo = FALSE, eval = FALSE}
death <- preparedata_fn(data_summarised[, c("Age", "Year", "Claim")], 
                                  ages = 35:65, years = 2016:2020)
expo <- preparedata_fn(data_summarised[, c("Age", "Year", "Exposure")], 
                                ages = 35:65, years = 2016:2020)
```

```{r}
death <- preparedata_fn(uk_mortalitydata[, c("Age", "Year", "Deaths")], 
                                  ages = 30:60, years = 2000:2020)
expo <- preparedata_fn(uk_mortalitydata[, c("Age", "Year", "Exposures")], 
                                ages = 30:60, years = 2000:2020)
```

## Data supplied as a 3D data array

Alternatively, users can supply the data as a 3-dimensional (3D) data array. `dxt_array_product` is a 3D array containing mortality data stratified by insurance product (see `?dxt_array_product` for details), where dim one: 4 insurance products, dim two: 83 ages, dim three: 5 years.

```{r}
data("dxt_array_product")
data("Ext_array_product")

# preview of death data the 1st insurance product called "ACI"
str(dxt_array_product["ACI",,,drop = FALSE])
```

Similarly, users can prepare the data by either by inputting the data as a 3-way array, or by specifying the name of the stratum to load using the argument `strat_name`:

```{r, eval = FALSE}

# inputting the data as a 3-way array
death <- preparedata_fn(dxt_array_product["ACI",,,drop = FALSE], ages = 35:65, years = 2016:2020)
expo <- preparedata_fn(Ext_array_product["ACI",,,drop = FALSE], ages = 35:65, years = 2016:2020)

# specifying the name of the stratum to load using `strat_name`
death <- preparedata_fn(dxt_array_product,strat_name="ACI", ages = 35:65, years = 2016:2020)
expo <- preparedata_fn(Ext_array_product,strat_name="ACI", ages = 35:65, years = 2016:2020)
```

Data array types  are less conventional but can be useful if data has been stored as it is. Preserving this data structure is useful for JAGS implementation later (for package development purposes).

## Data supplied as a data matrix

Suppose the data is provided in a 2-dimensional matrix format by age and year, commonly used in the literature. For example, the following is an illustration:

```{r}
# preview of death data the 1st insurance product called "ACI"
str(dxt_array_product["ACI",,,drop = TRUE])
```

To prepare data of type matrix, users need to specify the argument `data_matrix=TRUE`.

```{r, eval = FALSE}
death<-preparedata_fn(dxt_array_product["ACI",,,drop = TRUE],data_matrix=TRUE,ages=35:65)
expo<-preparedata_fn(Ext_array_product["ACI",,,drop = TRUE],data_matrix=TRUE,ages=35:65)
```

## Running the model

For illustrative purposes, we chose the UK mortality data.

Once the data have been prepared, they can be passed to **runBayesMoFo**, which is the core function in the package for estimating mortality models. 

As this package is built on top of the `rjags` package, it is capable of handling missing values in the death data, provided they are coded as `NA`. However, if missing values are present in the exposures data, these will be 
automatically replaced with a default value of 100, and predictions will be performed using that value.
If the user wishes to perform prediction for a specific exposure value, they can manually
set the desired value of the exposure (leaving the value of death count as `NA`, of course).

Users also have the option to perform model selection, depending on their needs. If more than one model is provided in the `models` argument, model selection is performed by default using deviance information criterion (DIC). The argument `models` can be set equal to:

* **"all"**, in which case, all models are fitted; 
* a vector with the name of the models, in which case the models in the list are
fitted;
* If the argument is left unspecified, a default set of models is fitted.

For example, the code below fit the `LC`, `CMB_M3`, and `APCI` models to the data.

```{r, eval = FALSE}
fitmodel <- runBayesMoFo(death, expo,
                         models = c("LC",
                                    "CBD_M3",
                                    "APCI")
                         )
```

Note that one can also run the individual functions rather than using the function **runBayesMoFo**. For example,

```{r, eval = FALSE}
fitmodel <- fit_LC(death, expo)
```

All other functions for analysing the output (see later) would work equally. That being said, users are highly recommend to use the function **runBayesMoFo** even when only one model is needed. That is,

```{r, eval = FALSE}
fitmodel <- runBayesMoFo(death, expo, models = "LC")
```

The full list of models, with the specific names, is available by checking `?runBayesMoFo`. Alternatively, one can query the model details through the documentation within the package, i.e. `?fit_LC`.

The argument `family` defines the specification for the distribution of death. A summary is as below (note that for AP data, just suppress the subscript $p$):

* If `family="poisson"`, then as proposed by @poissonlca,
$$d_{x,t,p} \sim \text{Poisson}(E^c_{x,t,p} m_{x,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$. The specification is used in conjunction with the log link function, i.e. 
$$\log(m_{x,t,p}) = \eta_{x,t,p} $$
where $\eta_{x,t,p}$ is the predictor that depends on the functional form of the mortality model, a full list of which is provided in Appendices A and B.

* Similarly, if `family="nb"` (default), then a negative binomial distribution is fitted with the log link function, 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}) = \eta_{x,t,p} , $$
where $\phi$ is the overdispersion parameter. This specification of the death distribution is standard practice for incorporating overdispersion, a phenomenon commonly observed in mortality modelling. See @wong_2023.

* If `\code{family="binomial"}`, then 
$$d_{x,t,p} \sim \text{Binomial}(E^0_{x,t,p} , q_{x,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. 
The binomial specification is used in conjunction with the logit link function, i.e.
$$\text{logit}(q_{x,t,p})= \log(\frac{q_{x,t,p}}{1-q_{x,t,p}}) = \eta_{x,t,p}. $$

For example, the following fit the same set of models using the Poisson distribution for modelling number of death.

```{r, eval = FALSE}
fitmodel <- runBayesMoFo(death, expo,
                         models = c("LC",
                                    "CBD_M3",
                                    "APCI"),
                         family="poisson"
                         )
```

There are also other arguments for customising the MCMC sampling of the posterior distributions, as below:

* `n.adapt` specifies the number of iterations for the adaptation phase. See `?rjags::adapt` for more details.
* `n.chain` specifies the number of parallel chains for the posterior  sampling under each model.
* `n_iter` specifies the number of iterations in the posterior sampling.
* `thin` specifies the thinning interval for the posterior sampling

## Forecast

As part of the function **runBayesMoFo**, forecasting can be performed by setting the argument `forecast=TRUE`, with the parameter `h` specifying the forecast horizon. For example, the code below fit the `LC`, `CMB_M3`, and `APCI` models to the data, and forecast the models for `h=6` time points into the future.

```{r, message = FALSE, warning = FALSE}
fitmodel_forecast <- runBayesMoFo(death, expo,
                         models = c("LC",
                                    "CBD_M3",
                                    "APCI"),
                         forecast = TRUE,
                         h = 6, 
                         n.chain = 2
                         )
```


## Analyzing the output

After running the model, users can then query the best and worst models (in terms of DIC) among the competing models.

```{r}
fitmodel_forecast$best_model
fitmodel_forecast$worst_model
```

A table showing the DIC of all models fitted can be returned too.

```{r}
fitmodel_forecast$DIC
```

One can retrieve the fitted results for the best and worst performing models, both of which are of type `fit_result`.

```{r, eval = FALSE}
fitmodel_forecast$result$best
fitmodel_forecast$result$worst
```

The function `plot_param_fn` plots all the fitted parameters of the model specified, using posterior samples generated. If more than one models were specified previously when running **runBayesMoFo**, then only the best model will be illustrated.

```{r, fig.width = 10, fig.height = 10}
plot_param_fn(fitmodel_forecast)
```

As evident in the plot, all fitted and forecasted parameters will be included, with solid lines indicating the medians and dashed lines representing the credible intervals generated from the posterior samples. By default, the intervals are constructed based on 95\% credibility, but can be changed using the argument `pred_int`. For instance, for 80\% credible intervals,

```{r, eval = FALSE}
plot_param_fn(fitmodel_forecast, pred_int = 0.80)
```

The argument `legends` argument can be used to suppress the legends for better visualisation (e.g. if visibility is blocked by the legend boxes). 

```{r, fig.width = 10, fig.height = 10}
plot_param_fn(fitmodel_forecast, pred_int = 0.80, legends = FALSE)
```

The function `plot_rates_fn` plots the fitted death rates of the model specified for specific ages and years, using posterior samples generated. Again, if more than one models were specified previously when running **runBayesMoFo**, then only the best model will be illustrated. By default, the (log) death rates will be plotted against age for the first nine years.

As before, both fitted and forecasted death rates will be included, with solid lines indicating the medians and dashed lines representing the credible intervals (95\% by default but can be changed using the argument `pred_int`) generated from the posterior samples. Also, observed crude death rates will also be included as coloured dots.

```{r, fig.width = 10, fig.height = 10}
plot_rates_fn(fitmodel_forecast)
```

For better visualisation, one may customise the argument `plot_years` to plot only selected years. Note that if more than nine years have been specified, then only the first nine years will be plotted.

```{r, fig.width = 10, fig.height = 5}
plot_rates_fn(fitmodel_forecast, plot_years = c(2016,2020,2024))
```

The argument `plot_type` allows users to plot death rates against year instead to better visualise temporal variations in death rates. The argument `plot_ages` can be used accordingly to specify which ages to plot.

```{r, fig.width = 10, fig.height = 5}
plot_rates_fn(fitmodel_forecast, plot_type = "time", plot_ages = c(35,45,55))
```

The function `summary_fn` produces a summary of the model results, including posterior means, standard deviations, medians, lower and upper quantiles based on the credibility specified using `pred_int`.

```{r}
summary_fitmodel<-summary_fn(fitmodel_forecast)
```

To obtain the posterior means and standard deviations of all death rates,
```{r, eval = FALSE}
#posterior means
summary_fitmodel$rates_summary$mean

#posterior standard deviations
summary_fitmodel$rates_summary$std
```

To obtain the posterior medians and lower/upper quantiles of all death rates,
```{r, eval = FALSE}
#posterior medians
summary_fitmodel$rates_pn$median

#lower quantiles
summary_fitmodel$rates_pn$lower

#upper quantiles
summary_fitmodel$rates_pn$upper
```

Correspondingly, for model parameters,

```{r, eval = FALSE}
#posterior means
summary_fitmodel$param_summary$mean

#posterior standard deviations
summary_fitmodel$param_summary$std

#posterior medians
summary_fitmodel$param_pn$median

#lower quantiles
summary_fitmodel$param_pn$lower

#upper quantiles
summary_fitmodel$param_pn$upper
```

## Convergence diagnostics

Users can assess if convergence has been attained by the MCMC posterior sampling procedure. The functions `diag_rates_fn` produces (by default) trace plots, density plots, as well as effective sample sizes of the posterior samples of death rates under the best model. For example, 

```{r, fig.width = 10, fig.height = 10}
diagnostics_rates_result<-converge_diag_rates_fn(fitmodel_forecast)
```

The effective sample sizes can be viewed as:

```{r}
diagnostics_rates_result$ESS
```

The arguments `plot_strata`, `plot_ages`, `plot_years` can be used to specify which death rates to examine, as follows.

```{r, fig.width = 10, fig.height = 10}
converge_diag_rates_fn(fitmodel_forecast, plot_ages = c(35,45,55), plot_years = c(2016,2020,2024))
```

The arguments `trace` and `density` can be used to specify only plotting one of them or none according to personal preferences. 

```{r, eval = FALSE}
#for only trace plots
converge_diag_rates_fn(fitmodel_forecast, plot_ages = c(35,45,55), plot_years = c(2016,2020,2024), trace = TRUE, density = FALSE)
```

Auto-correlation plots can also displayed to check if posterior samples are too correlated. 

```{r, fig.width = 10, fig.height = 10}
#for only acf plots
converge_diag_rates_fn(fitmodel_forecast, plot_ages = c(35,45,55), plot_years = c(2016,2020,2024), trace = FALSE, density = FALSE, acf_plot = TRUE)
```

Similarly, convergence can be assessed for fitted parameters using the function `converge_diag_rates_fn`.

```{r, fig.width = 10, fig.height = 10}
diagnostics_param_result<-converge_diag_param_fn(fitmodel_forecast)
```

The effective sample sizes can be viewed as:

```{r}
diagnostics_param_result$ESS
```

By default, the function examines a selection of the parameters. But the arguments `plot_params` can be used to specify which set of parameters to examine, as follows.

```{r, fig.width = 10, fig.height = 10}
converge_diag_param_fn(fitmodel_forecast, plot_params = c("kappa","gamma","rho","phi","sigma2_kappa"))
```

To check the names of parameters available for examining:

```{r}
fitmodel_forecast$result$best$param
```

For the rate-related parameters such as alpha, beta, kappa, gamma etc., only three of the randomly selected subset will be examined when specified. If a particular parameter is to be assessed, users need to specify clearly the indices of the parameters to be examined. For instance, the following will assess the beta parameters for the first stratum and the third age, as well as kappa for the first stratum and the fourth year.

```{r, fig.width = 9, fig.height = 5}
converge_diag_param_fn(fitmodel_forecast, plot_params = c("kappa[1,4]","gamma[1,2]"))
```

To check the full list of parameters available for examining:

```{r}
colnames(fitmodel_forecast$result$best$post_sample[[1]])[!startsWith(colnames(fitmodel_forecast$result$best$post_sample[[1]]),"q[")]
```

The arguments `trace` and `density` can be used to specify only plotting one of them or none according to personal preferences. 

```{r, eval = FALSE}
#for only trace plots
converge_diag_param_fn(fitmodel_forecast, plot_params = c("kappa[1,4]","gamma[1,2]"), trace = TRUE, density = FALSE)
```

Auto-correlation plots can also displayed to check if posterior samples are correlated. 

```{r, fig.width = 5, fig.height = 5}
#for only acf plots
converge_diag_param_fn(fitmodel_forecast, plot_params = c("kappa[1,4]","gamma[1,2]"), trace = FALSE, density = FALSE, acf_plot = TRUE)
```

Several other commonly used MCMC convergence diagnostics, such as Gelman's R statistic (@rubin), Geweke's Z scores (@geweke1991evaluating), Heidel's Stationarity and Half-width tests (see @heidelberger1981spectral and @heidelberger1983 for more details), can be computed and illustrated using the function `converge_diag_fn`.

```{r, fig.width = 9, fig.height = 5}
converge_diag_result<-converge_diag_fn(fitmodel_forecast, plot_gelman = TRUE, plot_geweke = TRUE)
```

```{r}
#Gelman's R
head(converge_diag_result$gelman_diag$psrf)

#Geweke's Z
head(converge_diag_result$geweke_diag$z)

#Heidel's Stationarity and Half-width tests
head(converge_diag_result$heidel_diag)
```

# Age-period-product (APP) data

## Data preparation

Similarly to before, the data needs to be formatted in a data.frame with columns name *Age*, *Year*, *Deaths*,*Exposures* and *Cause*. Some examples are provided below.

## Data supplied as a data-frame 

```{r}
data(uk_deathscausedata)
head(uk_deathscausedata, n = 10)
```

```{r}
death <- preparedata_fn(uk_deathscausedata[,c("Age","Year","Deaths","Cause")])
expo <- preparedata_fn(uk_deathscausedata[,c("Age","Year","Exposures","Cause")])

#or if require a subset of the data
death <- preparedata_fn(uk_deathscausedata[,c("Age","Year","Deaths","Cause")], 
                                 ages = seq(45,90,by=5), years = 2001:2020)
expo <- preparedata_fn(uk_deathscausedata[,c("Age","Year","Exposures","Cause")], 
                                ages = seq(45,90,by=5), years = 2001:2020)

str(death)
str(expo)
```

As shown above, the cause of death data consists of numbers of death and exposures for five causes of death, spanning years 2001-2020 and 5-year age groups between 45-90. 

## Data supplied as a 3D data array

Alternatively, users may wish to supply data which is already sorted as a 3D array (dim 1: strata, dim 2: ages, dim 3: years). 

```{r, eval = FALSE}
data("dxt_array_product");data("Ext_array_product")
str(dxt_array_product) # 3D data array

death<-preparedata_fn(dxt_array_product,ages=35:65)
expo<-preparedata_fn(Ext_array_product,ages=35:65)
```

# Running the model

The syntax is similar to the case of the age-period data. The function automatically recognises the structure of the data after being processed by `preparedata_fn`. For illustration, we fit the model by @li2005coherent on the cause of death data.

```{r, message = FALSE, warning = FALSE}
fitmodel_forecast <- runBayesMoFo(death, expo,
                         models = "MLiLee",
                         forecast = TRUE,
                         h = 5,
                         quiet = TRUE,
                         n.chain = 2
                         )
```

The argument `quiet=TRUE` was used to suppress messages generated during model compilation stage.

```{r, eval = FALSE, echo = FALSE}
fitmodel_forecast$best_model
fitmodel_forecast$worst_model
fitmodel_forecast$DIC
```

## Plot the output

```{r, fig.width = 10, fig.height = 10}
plot_param_fn(fitmodel_forecast)
```

```{r, fig.width = 10, fig.height = 5}
plot_rates_fn(fitmodel_forecast, plot_years = c(2005,2020,2025))
```

The argument `plot_type` allows users to plot death rates against year instead to better visualise temporal variations in death rates. The argument `plot_ages` can be used accordingly to specify which ages to plot.

```{r, fig.width = 10, fig.height = 5}
plot_rates_fn(fitmodel_forecast, plot_type = "time", plot_ages = c(45,65,85))
```

## Convergence diagnostics

```{r, fig.width = 10, fig.height = 10}
diagnostics_param_result<-converge_diag_rates_fn(fitmodel_forecast, plot_ages = c(45,65,85), plot_years = c(2005,2020,2025))
```


```{r, fig.width = 10, fig.height = 10}
#for only acf plots
converge_diag_rates_fn(fitmodel_forecast, plot_ages = c(45,65,85), plot_years = c(2005,2020,2025), trace = FALSE, density = FALSE, acf_plot = TRUE)
```

The effective sample sizes can be viewed as:

```{r}
diagnostics_param_result$ESS
```

By default, the function examines a selection of the parameters. But the arguments `plot_params` can be used to specify which set of parameters to examine, as follows.

```{r, fig.width = 10, fig.height = 10}
converge_diag_param_fn(fitmodel_forecast, plot_params = c("beta","kappa","rho","phi","sigma2_kappa"))
```

To check the names of parameters available for examining:

```{r}
fitmodel_forecast$result$best$param
```

For the predictor-related parameters such as alpha, beta, kappa, gamma etc., only three of the randomly selected subset will be examined when specified. If a particular parameter is to be assessed, users need to specify clearly the indices of the parameters to be examined. For instance, the following will assess the beta parameters for the first stratum and the third age, as well as kappa for the second stratum and the fourth year.

```{r, fig.width = 9, fig.height = 5}
converge_diag_param_fn(fitmodel_forecast, plot_params = c("beta[1,3]","kappa[2,4]"))
```

To check the full list of parameters available for examining:

```{r}
colnames(fitmodel_forecast$result$best$post_sample[[1]])[!startsWith(colnames(fitmodel_forecast$result$best$post_sample[[1]]),"q[")]
```

The arguments `trace` and `density` can be used to specify only plotting one of them or none according to personal preferences. 

```{r, eval= FALSE}
#for only trace plots
converge_diag_param_fn(fitmodel_forecast, plot_params = c("beta[1,3]","kappa[2,4]"), trace = TRUE, density = FALSE)
```

Auto-correlation plots can also displayed to check if posterior samples are correlated. 

```{r, fig.width = 5, fig.height = 5}
#for only acf plots
converge_diag_param_fn(fitmodel_forecast, plot_params = c("beta[1,3]","kappa[2,4]"), trace = FALSE, density = FALSE, acf_plot = TRUE)
```

Convergence diagnostics can be applied as usual, but are not run in the example below.

```{r, eval= FALSE}
converge_diag_result<-converge_diag_fn(fitmodel_forecast, plot_gelman = TRUE, plot_geweke = TRUE)
```

```{r, eval= FALSE}
#Gelman's R
head(converge_diag_result$gelman_diag$psrf)

#Geweke's Z
head(converge_diag_result$geweke_diag$z)

#Heidel's Stationarity and Half-width tests
head(converge_diag_result$heidel_diag)
```

Interestingly, there is an article discussing the use of common (shared) cohort effects for modelling cause of death data as described by @cairns2023common. Thus, we can fit some of the models that incorporate shared cohort effects as below (NOT RUN).

```{r, eval = FALSE}
fitmodel_forecast <- runBayesMoFo(death, expo,
                         models = c("APCI_sharegamma",
                                    "RH_sharegamma"),
                         forecast = TRUE,
                         h = 5,
                         quiet = TRUE
                         )
```


```{r, eval = FALSE}
fitmodel_forecast$DIC

fitmodel_forecast$best_model

plot_param_fn(fitmodel_forecast)
```

```{r, eval = FALSE}

plot_rates_fn(fitmodel_forecast, plot_years = c(2005,2020,2025))

plot_rates_fn(fitmodel_forecast, plot_type = "time", plot_ages = c(45,65,85))

```

```{r, eval = FALSE}
converge_diag_rates_fn(fitmodel_forecast, plot_ages = c(45,65,85), plot_years = c(2005,2020,2025))

converge_diag_param_fn(fitmodel_forecast, plot_params = c("kappa","rho","phi","sigma2_kappa"))
```

# Appendix A: full list of age-period models considered

| Model    | Predictor, $\eta_{x,t}$ |
|----------|-------------------------|
| APCI     | $a_{x}+b_{x}(t-\bar{t})+k_{t} + \gamma_{c}$ |
| CBD_M3   | $a_{x}+k_{t} + \gamma_{c}$  |
| CBD_M5   | $k^1_{t} + k^2_{t}(x-\bar{x})$  |
| CBD_M6   | $k^1_{t} + k^2_{t}(x-\bar{x}) +\gamma_{c}$  |
| CBD_M7   | $k^1_{t} + k^2_{t}(x-\bar{x}) + k^3_{t}((x-\bar{x})^2-\hat{\sigma}_x^2) +\gamma_{c}$  |
| CBD_M8   | $k^1_{t} + k^2_{t}(x-\bar{x}) +\gamma_{c}(constant-x)$  |
| LC       | $a_{x}+b_{x}k_{t}$  |
| MLiLee   | $a_{x}+B_xK_t$  |
| PLAT     | $a_{x}+k^1_{t} + k^2_{t}(\bar{x}-x) + k^3_{t}(\bar{x}-x)^+ +\gamma_{c}$  |
| RH       | $a_{x}+b_{x}k_{t} + \gamma_{c}$  |


# Appendix B: full list of age-period-product models considered

| Model                        | Predictor, $\eta_{x,t,p}$ |
|------------------------------|---------------------------|
| APCI                         | $a_{x,p}+b_{x,p}(t-\bar{t})+k_{t,p} + \gamma_{c,p}$ |
| APCI_sharealpha              | $a_{x}+b_{x,p}(t-\bar{t})+k_{t,p} + \gamma_{c,p}$ |
| APCI_sharebeta               | $a_{x,p}+b_{x}(t-\bar{t})+k_{t,p} + \gamma_{c,p}$ |
| APCI_sharegamma              | $a_{x,p}+b_{x,p}(t-\bar{t})+k_{t,p} + \gamma_{c}$ |
| APCI_sharealpha_sharebeta    | $a_{x}+b_{x}(t-\bar{t})+k_{t,p} + \gamma_{c,p}$ |
| APCI_sharealpha_sharegamma   | $a_{x}+b_{x,p}(t-\bar{t})+k_{t,p} + \gamma_{c}$ |
| APCI_sharebeta_sharegamma    | $a_{x,p}+b_{x}(t-\bar{t})+k_{t,p} + \gamma_{c}$ |
| APCI_shareall                | $a_{x}+b_{x}(t-\bar{t})+k_{t,p} + \gamma_{c}$ |
| CBD_M3                       | $a_{x,p}+k_{t,p} + \gamma_{c,p}$  |
| CBD_M3_sharealpha            | $a_{x}+k_{t,p} + \gamma_{c,p}$  |
| CBD_M3_sharegamma            | $a_{x,p}+k_{t,p} + \gamma_{c}$  |
| CBD_M3_shareall              | $a_{x}+k_{t,p} + \gamma_{c}$  |
| CBD_M5                       | $k^1_{t,p} + k^2_{t,p}(x-\bar{x})$  |
| CBD_M6                       | $k^1_{t,p} + k^2_{t,p}(x-\bar{x}) +\gamma_{c,p}$  |
| CBD_M6_sharegamma            | $k^1_{t,p} + k^2_{t,p}(x-\bar{x}) +\gamma_{c}$  |
| CBD_M7                       | $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}$  |
| CBD_M7_sharegamma            | $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}$  |
| CBD_M8                       | $k^1_{t,p} + k^2_{t,p}(x-\bar{x}) +\gamma_{c,p}(c_p-x)$  |
| CBD_M8_sharegamma            | $k^1_{t,p} + k^2_{t,p}(x-\bar{x}) +\gamma_{c}(c_p-x)$  |
| LC                           | $a_{x,p}+b_{x,p}k_{t,p}$  |
| LC_sharealpha                | $a_{x}+b_{x,p}k_{t,p}$  |
| LC_sharebeta                 | $a_{x,p}+b_{x}k_{t,p}$  |
| LC_shareall                  | $a_{x}+b_{x}k_{t,p}$  |
| M1A                          | $a_{x}+c_p+b_xk_t$  |
| M1U                          | $a_{x,p}+b_xk_t$  |
| M1M                          | $a_{x}c_p+b_xk_t$  |
| M2A1                         | $a_{x}+(c_p+b_x)k_t$  |
| M2A2                         | $a_{x}+b_{x,p}k_t$  |
| M2Y1                         | $a_{x}+b_x(k_t+c_p)$  |
| M2Y2                         | $a_{x}+b_{x}k_{t,p}$  |
| MLiLee                       | $a_{x,p}+b_{x,p}k_{t,p}+B_xK_t$  |
| MLiLee_sharealpha            | $a_{x}+b_{x,p}k_{t,p}+B_xK_t$  |
| PLAT                         | $a_{x,p}+k^1_{t,p} + k^2_{t,p}(\bar{x}-x) + k^3_{t,p}(\bar{x}-x)^+ +\gamma_{c,p}$  |
| PLAT_sharealpha              | $a_{x}+k^1_{t,p} + k^2_{t,p}(\bar{x}-x) + k^3_{t,p}(\bar{x}-x)^+ +\gamma_{c,p}$  |
| PLAT_sharegamma              | $a_{x,p}+k^1_{t,p} + k^2_{t,p}(\bar{x}-x) + k^3_{t,p}(\bar{x}-x)^+ +\gamma_{c}$  |
| PLAT_shareall                | $a_{x}+k^1_{t,p} + k^2_{t,p}(\bar{x}-x) + k^3_{t,p}(\bar{x}-x)^+ +\gamma_{c}$  |
| RH                           | $a_{x,p}+b_{x,p}k_{t,p} + \gamma_{c,p}$  |
| RH_sharealpha                | $a_{x}+b_{x,p}k_{t,p} + \gamma_{c,p}$  |
| RH_sharebeta                 | $a_{x,p}+b_{x}k_{t,p} + \gamma_{c,p}$  |
| RH_sharegamma                | $a_{x,p}+b_{x,p}k_{t,p} + \gamma_{c}$  |
| RH_sharealpha_sharebeta      | $a_{x}+b_{x}k_{t,p} + \gamma_{c,p}$  |
| RH_sharealpha_sharegamma     | $a_{x}+b_{x,p}k_{t,p} + \gamma_{c}$  |
| RH_sharebeta_sharegamma      | $a_{x,p}+b_{x}k_{t,p} + \gamma_{c}$  |
| RH_shareall                  | $a_{x}+b_{x}k_{t,p} + \gamma_{c}$  |

# References
