Bayesian Population Forecasting: Extending the Lee-Carter Method

In this article, we develop a fully integrated and dynamic Bayesian approach to forecast populations by age and sex. The approach embeds the Lee-Carter type models for forecasting the age patterns, with associated measures of uncertainty, of fertility, mortality, immigration, and emigration within a cohort projection model. The methodology may be adapted to handle different data types and sources of information. To illustrate, we analyze time series data for the United Kingdom and forecast the components of population change to the year 2024. We also compare the results obtained from different forecast models for age-specific fertility, mortality, and migration. In doing so, we demonstrate the flexibility and advantages of adopting the Bayesian approach for population forecasting and highlight areas where this work could be extended.


Introduction
This work is guided by two aims. The first is to have a flexible platform for forecasting populations. Most statistical offices in developed countries utilize data obtained from different sources, including administrative registers, surveys, and censuses with varying levels of quality and measurement. Cohort component models have long been the standard apparatus for producing projections, but wide differences remain in the underlying assumptions, especially regarding the treatment of migration and the level of detail provided. The forecasting approach we have developed is one that can be adapted to include all demographic components of change by age and sex, and provides measures of the accuracy of the forecasts.
As we move away from deterministic population projections to those that provide measures of uncertainty, we believe it is important to integrate the various sources of uncertainty into the modeling framework. The rationale for considering a Bayesian approach is that it offers a natural probabilistic framework to predict future populations. Variability in the data and uncertainties in the parameters and model choice can be explicitly incorporated by using probability distributions, and the predictive distributions follow directly from the probabilistic model applied. The approach also allows the inclusion of expert judgments, together with their uncertainty, in the model framework.
The second aim of this article is to provide a flexible and consistent method for forecasting the age patterns of fertility, mortality, and migration that drive our forecast results. A vast literature focuses on modeling demographic events (e.g., Bongaarts and Bulatao 2000), but approaches for forecasting the age patterns are less abundant. Methods for forecasting migration, in particular, represent a persistent weakness in population forecast models (Bijak 2010). Our population forecasting model is developed with these two aims in mind. We focus on generalizing and extending the Lee-Carter model for forecasting mortality to age-specific fertility and migration (Lee and Carter 1992), and integrating these into the cohort projection mechanism. One of the contributions of the proposed approach is forecasting age-specific emigration rates and immigration volumes, following the suggestions of Rees (1986:148).
Because the age patterns of immigration and emigration are more regular than those observed for net migration, they are better amenable to modeling by using Lee-Carter models.

Forecasting the Age Patterns of Fertility, Mortality, and Migration
There is a long history of modeling the age patterns of fertility, mortality, and migration events (Booth 2006). This work has demonstrated the persistent and strong regularities in the age patterns over time and across space (Preston et al. 2001:191-210;Rogers 1986;Rogers and Little 1994) driven by biological and social life course mechanisms (Courgeau 1985). The age regularities exhibited in demographic patterns allow population forecasters to simplify their underlying assumptions and models. Indeed, some forecasts focus on indicator variables, such as the total fertility rate (TFR), life expectancy, or net migration rate, which are then converted into an assumed age distribution (see, e.g., Raftery et al. 2012;Wilson and Bell 2007).
The main approaches for modeling the age patterns of demographic components include the imposition of empirical tables obtained from other countries and historical settings (Coale and Demeney 1966), parametric model schedules (Coale and McNeil 1972;Coale and Trussell 1974;Heligman and Pollard 1980;Knudsen et al. 1993;Rogers 1986;Rogers and Castro 1981;Rogers et al. 1978), relational models (Brass 1974), functional models (De Beer 2011;Hyndman and Booth 2008;Hyndman and Ullah 2007;Lee and Carter 1992), and hierarchical Bayesian models (Czado et al. 2005;Girosi and King 2008). Of these, the most successful and widely used for forecasting future patterns and uncertainty have been the relational and functional modelsparticularly, the Lee-Carter model for mortality (Booth and Tickle 2008).
Relational and functional models include standards or time-invariant patterns, which are then perturbed based on a set of parameters. Because the shapes of age-specific fertility, mortality, and migration largely remain the same over time, these approaches provide a powerful tool for forecasting. In the Lee-Carter case, single-or five-year age groups are altered on the basis of time series equations. The main critique of the original Lee-Carter approach is that it can produce implausible forecasts for particular age groups (see Girosi and King 2008:38-42). As a consequence, several extensions have been developed to accommodate cohort effects, correlation between sexes, and smoothing (Booth and Tickle 2008;Lee 2000), as well as to forecast fertility rates (Lee 1993).

Bayesian Population Forecasting
The need to incorporate probabilistic uncertainty into population estimates and forecasts is well known. Probabilistic forecasts have the advantage over variant style projections in that they specify the chances or probability that a particular future population value will be within any given range (Ahlburg and Land 1992;Alho andSpencer 1985, 2005;Bongaarts and Bulatao 2000;Keilman 1990;Lee and Tuljapurkar 1994;Lutz 1996). With variant projections, on the other hand, the user has no idea how likely future population values are, but only that they are plausible scenarios representing the "most likely" and the "extreme" high and low possibilities.
However, despite the known advantages of probabilistic forecasts, they have yet to be widely adopted by statistical agencies (Lutz and Goldstein 2004). The reason is that there are many types of uncertainties to consider, and including them in projections is not always straightforward.
According to Alho (1999), probabilistic population forecasting within the Bayesian framework has a tradition dating back over 60 years to the seminal works of Leo Törnqvist and colleagues (Hyppölä et al. 1949). However, it was not until the 1980s when probabilistic methods began entering mainstream demography. These included time series extrapolations, expert-based forecasting, and past error propagation (for a detailed overview of different approaches, see Bijak 2010). Examples of Bayesian models for population forecasting were practically non-existent until the past few years, with the notable exception of Daponte et al. (1997). Recent advances in fast computation and numerical methods have enabled a more widespread use of the Bayesian approach in many fields of application, including population forecasting. At a global level, Raftery et al. (2012) proposed a generic model for all countries of the world, based on aggregate indicators (TFRs and life expectancies) and model life tables. At the opposite end of the data spectrum, Bryant and Graham (2013) suggested a comprehensive Bayesian approach to reconciling different data sources for New Zealand, a country with very good availability of population statistics. Bryant and Graham's (2013) approach uses an accounting framework for estimating New Zealand's current population disaggregated by regions, age, sex and time. It combines various data sources, including vital events registers, censuses, and school and electoral rolls. The model constrains the true values of the demographic components by the population accounting equation.
Also, age, time, sex and regional patterns are specified by main effect and two-way interaction terms within a Poisson-gamma model (similar to non-Bayesian approaches of Smith et al. (2010) and Raymer et al. (2011b)).

Methodology
In this section, we first introduce the forecasting model proposed by Lee and Carter (1992) and then describe how it can be extended and applied within a Bayesian framework. The Lee-Carter model was originally designed to forecast age-specific mortality rates with the following specification: (1) where the logarithm of the age and time-specific mortality rate m x,t is decomposed into an overall age profile, a x , averaged over the entire period under consideration, and age-specific changes in mortality b x . The subscripts x and t denote age and time, respectively. The b x parameter describes how fast the rates decline over time in response to changes in the time-specific effect k t . The error term x x, t is assumed to be normally distributed with a mean of 0 and a constant variance. To forecast mortality rates into the future, a simple random walk with drift model for k t was proposed: (2) To ensure identifiability of the model parameters, b x are constrained to sum to unity and k t are constrained to sum to 0 (Lee and Carter 1992:661). Lee (1993) subsequently proposed a similar model to forecast age-specific fertility for the United States. In that model, several constraints were introduced to represent the prior information on fertility.
In this article, we extend and adapt the Lee-Carter model to create a general framework for forecasting mortality, fertility, emigration and immigration. In this framework, we first assume that count data on a given population component Y x,t g,k follow a Poisson distribution: where denotes a component of population change: D represents deaths (mortality), B is births (fertility), E is emigration, I is immigration, sex k = F denotes females, and k = M stands for males. Parameter m x,t g, k denotes a ratio of counts of demographic event scaled to the size of the population exposed to risk of these events, R x,t g,k . As in the original Lee-Carter model, t and x denote time and age, respectively. For mortality and emigration, the population at risk is the same; for fertility, it consists of women of reproductive age. For immigration, we forecast the counts rather than rates; hence, we assume that . Czado et al. (2005) developed the extension of the Lee-Carter model to incorporate Poisson variability of death counts within the Bayesian framework.
Second, we assume that the logarithm of the rate follows a normal distribution: where a x g,k , b x g,k , and k t g,k represent the same parameters as in the Lee-Carter model, and denotes a cohort effect. The cohort effect, introduced by Renshaw and Haberman (2006) for mortality, is incorporated in our framework for the sake of generality, but it may be omitted if not required. Throughout this article, denotes a normal distribution with a mean of m and precision (inverse variance) t . The normal distribution assumed for rates is an extension of the Czado et al. (2005) model. It allows capturing the overdispersion that is not explained by the variability resulting from the Poisson sampling of count data.
Third, for the time-specific effects, k t g,k , we require a time series model, which facilitates the forecasting. This model can be univariate, such as random walk with drift in the original Lee-Carter specification, for each component and each sex; alternatively, the model can be multivariate (e.g., vector autoregression (VAR)), which allows exploring correlations between sexes and components. A time series model is also utilized for a cohort effect, . In our application, we use a univariate model. However, more general multivariate frameworks can be used.
To ensure identification of the parameters a x g,k , b x g,k , k t g,k , and , the following constraints are imposed: where z denotes the oldest age group. These constraints suffice to identify the bilinear model in Eq. (4) as long as there is a clear differentiation in the b x -that is, as long as they are not all equal to 1 / z, in which case the model reduces to the linear age-period-cohort (APC) model. The problem of identification of the period and cohort effects in the APC model has long been discussed in the literature (see Luo 2013, with a comment by Fienberg 2013).
To learn about the model parameters, and specifically about the forecasts of the population components, we use Bayesian inference. Bayes theorem states that the posterior distribution of the model parameters (e.g., forecasts of the age-specific mortality rates) is proportional to the product of the likelihood for the data and the prior distribution. In our approach, Bayesian inference integrates uncertainty from the demographic event (count) data expressed by the Poisson likelihood with weakly informative prior distributions that give preference to the historical data. Subsequently, all four components of population change are combined within the cohort component projection model.
In the following subsections, we present specific adaptations of the preceding framework to forecast mortality, fertility, and migration. We adopt a convention of proposing a very simple model for the data (such as the original Lee-Carter model) with a more general one. Because our extensions of the Lee-Carter model lead to a relatively complex specification of the probabilistic model, the closed forms of the posterior distributions are difficult to obtain analytically. Hence, we sample from the posterior distributions by using the Markov chain Monte Carlo (MCMC) algorithms implemented in the OpenBUGS software (Lunn et al. 2009). The example code used for the simulations for fertility is available in Online Resource 1 (section A.3). The other codes are available from the corresponding author upon request.

Forecasting Mortality
To forecast the mortality of males and females, we consider the Bayesian version of the original Lee-Carter model, denoted as M1, and a general extension of this model, denoted by M2. The Lee-Carter model M1 is specified as in Eqs. (1) and (2), and the age-specific mortality rates are calculated as . Model M2 for death counts Y x,t D,k follows Eqs. (3) and (4).
In M1, the time-specific parameters k t k (we drop the superscript D for each parameter for clarity of notation) follow univariate random walk with drift models, as in the original Lee-Carter specification. In M2, k t k for both sexes follow a bivariate vector autoregressive VAR (1) process with drift: where MVN 2 denotes a two-dimensional multivariate normal distribution, with the precision matrix T k ; f ij , i, j = 1, 2 are the parameters of the VAR(1) model, with f 0 j , j = 1, 2 being drift terms. The instantaneous and lagged correlations assumed for males and females reflect the assumption of parallel improvements in health conditions over time 1 . Finally, the cohort effect for each sex k follows a univariate autoregressive process AR(1): (7)

Forecasting Fertility
For forecasting age-specific fertility rates, we apply a simple version of the Lee (1993) model F1.
The extended model, denoted by F2, includes a cohort effect (see also Cheng and Lin 2010;Lee 1993Lee , 2000. The population at risk represents all women of reproductive age. The time component in F1 follows an ARMA(1,1) process (again we suppress the superscripts B and g because the model relates only to females): In Model F2, we use a simple univariate autoregressive process AR(1) for the time component k t and for : (9)

Forecasting Immigration Counts and Emigration Rates
To forecast immigration counts and emigration rates, we introduce two models: (1) a univariate model, denoted by IE1, which assumes no correlation between emigration and immigration of males and females; and (2) a multivariate model, denoted by IE2, in which we assume correlation between the time parameters k t for both sexes and both directions of migration. In both models, we incorporate smoothing, built in into the prior distributions for the age-specific model parameters a x and b x . Unlike mortality and fertility, there is no clear rationale for including a cohort effect parameter in forecasting immigration and emigration. Also, because immigration does not have an easily defined population at risk, we model the counts, which is a common practice in population projections (McDonald and Kippen 2002;Rees 1986).
In IE1, we assume a random walk without drift for emigration rates and immigration counts for both sexes: This simple model leads to forecasts with a constant expectation (as the last observation) and increasing uncertainty.
For IE2, we assume instantaneous correlation in the time parameters for emigration and immigration for both sexes: where , ,

Prior Distributions
In the absence of any prior information, we suggest using weakly informative distributions that will allow the data to drive the estimation. For the general forecasting model, we propose a set of priors, the hyperparameters of which may differ in particular applications. For both sexes k, the specifications of the prior distributions for the model parameters are as follows: where G(a,b) denotes the gamma distribution with a mean of a / b and variance a / b 2 ; U denotes uniform distribution; and Y b k is a precision matrix for a conditional distribution of b x , given that they sum to 1 (details are presented in Online Resource 1, section A.2). The preceding prior distributions imply weak information a priori about the model parameters. For the age-specific parameters a x and b x , we avoid specifications based on the data (i.e., the empirical Bayes approach) suggested by Czado et al. (2005). For the precision parameters t k and t g k , we follow Gelman's (2006) suggestions of avoiding conjugate gamma distributions, and we use uniform priors for standard deviations over a large range. The Wishart distribution is a standard prior distribution for precision matrix in the multivariate time series model and is easy to implement in the OpenBUGS software. The subscript l denotes the number of series in the model: males and females for mortality (i.e., l = 2) or male and female migration in both directions (i.e., l = 4). In the univariate model (e.g., for fertility or for a single sex), we replace the Wishart prior with a uniform prior for the standard deviation; that is, .
For low-quality data, smoothing of the age profile may be required. In this case, we propose a smoothing technique that is embedded in the specification of the prior distributions for parameters a x and b x . Smoothing prevents the artificial age patterns resulting from the sample data from being propagated in the forecasts. Our smoothing method is based on the spatial autoregressive processes (see, e.g., Besag 1986).
Prior densities for smoothing the age-specific parameters a x are constructed in the following way, with the sex and population component superscripts omitted for clarity. For the youngest and oldest age groups, we assume that the mean of the prior distribution depends on the second-youngest and second-oldest group, respectively: For the remaining age groups, we assume that their means depend on the average of the two neighboring age groups x -1 and x + 1: The preceding construction of the conditional precisions for each age group ensures that the unconditional precision is constant for all age groups. Finally, we assume a priori that , which is based on our assessment of results obtained with various values for this parameter and from visual inspections of the model fits. It implies a moderate degree of smoothing for a x . For b x , we assume the same pattern of smoothing as for a x , but we derive a distribution conditional on . The resulting multivariate normal distribution is This prior distribution is similar to the one in Eq. (12). However, the matrix Y b here is derived analogously to Eqs. (13) and (14) by assuming that all elements of b x follow an autoregressive process that in the limit tends to a random walk, but also conditional on (for more details see Online Resource 1, section A.2). The smoothing parameter t b can be sex-specific and direction-of-flow-specific, or one parameter can be used for all four flows, which allows borrowing of strength. The gamma distribution assumed for this parameter is characterized by a very heavy tail; thus, it allows this parameter to explore regions of large values that lead to "smoother" age profiles. Because the smoothing parameter has a vague prior distribution, the whole smoothing procedure is driven by the data rather than by subjective judgment. However, the results exhibit sensitivity to the specification of the prior for the smoothing parameter. Other prior distributions, such as truncated t or Cauchy, can be used. The degree of smoothing may also be controlled by fixing t b at some value, which can be found by a grid search, for example.

Population Projection Model
The results of forecasting the four components of population change-that is, samples from the posterior distributions of mortality, fertility, and emigration rates, as well as immigration counts-are subsequently combined into a cohort component projection model (see Preston et al. 2001:117-137;Rogers 1995). The projection model is specified as (16) where is vector of midyear population sizes by age and sex k and They include the standard transformation of the mortality and emigration rates, both of which result in the decrease of the population, into the survival rates. For the last year of age x = z, we assume the standard formula following Rogers (1995:107). Finally, the birth rates are constructed by using the fertility rates obtained from our model for forecasting fertility and survivorship of infants from the mortality forecasting model:

Model Validation and Selection
The models for the population components that underlie the population forecast are selected from the models proposed in the previous sections. The selection process is based on (1)  The DIC is a tool for assessing the goodness-of-fit of a model to the data, which enables selecting the best performing model. It is often considered a generalization of the Akaike information criterion (AIC) for comparing complex Bayesian hierarchical models. It utilizes a deviance of the likelihood evaluated at the mean of the posterior distribution of the likelihood as the goodness-of-fit measure, corrected with the "effective number of parameters" in a model (for definitions and formal derivation, see Spiegelhalter et al. 2002). The requirement for using DIC is that the posterior distribution is approximately multivariate normal. Selecting the best performing model is similar for the AIC-namely, the lower the value of the criterion, the better fit of the model.

Data
In this section, we illustrate our forecasting method with the data for the United Kingdom (UK).
These data represent a case in which the counts of all population components by single year of age and sex are available but their quality is varying. Because the data on vital events are recorded by the registers, they are considered to be of relatively good quality and precise.
Immigration and emigration counts are, however, produced by using the International Passenger Survey (IPS) and include both sampling and nonsampling errors, especially in the age profiles by

Results
In this section, we present the results of forecasting the population components with the models described in the previous section. For each component, we discuss the model's goodness-of-fit to the data and forecasts of the future patterns, and select the underlying model to be used for the population forecast.

Forecasts of Mortality
In the first row of Fig. 3 2010). An alternative argument for the inclusion of the cohort effect in the model is lifelong processes that might affect mortality, such as smoking (e.g., Doll et al. 2004).
To support our rationale for selecting model M2, we analyze the in-sample forecasts from

Forecasts of Fertility
The age-specific forecasts for fertility are presented in Fig. 4. In the first row, we observe the fit of the models F1 and F2 to the 2009 data. The model with the cohort effect (F2) provides a better fit with lower uncertainty. Also, the 2024 forecast (second row) appears more plausible than the forecast based on the F1 model, which produces an unrealistic median fertility rate of 0.3 for females aged 33-35 years.
The resulting TFRs are presented in the third row of Fig. 4. It is clear from the plots that F2 fits the data better than F1. Moreover, the projected TFR from F1 shows an explosive pattern that we consider unrealistic, with an explosive predicted TFR. Hence, we believe that the pattern of gradual diminishing of the recently increasing TFR produced by F2 better reflects our expectations about future fertility in the UK.
[ Fig. 4 about here.] The in-sample forecasts of the fertility rates confirm our rationale for choosing F2 as the foundation of the population forecast. Again, F2 appears to fit the data better (see the fourth row of Fig. 4). The resulting forecasts of TFR under F2 seem to be more uncertain than those of F1.  (Bongaarts and Feeney 1998). In particular, we refer to the recent postponement and subsequent recuperation of fertility in many developed countries, where the cohort effects are the most profound (see, e.g., Sobotka et al. 2011). In our results, slightly declining but still uncertain fertility rates may indicate a possibility of yet another period of postponement.

Forecasts of Emigration and Immigration
To forecast emigration rates and immigration counts for the UK we need to reflect on the quality of the survey data on migration and their implications for modeling. To overcome the irregularities described in the earlier Data section, we modified the models described in the section Forecasting Immigration Counts and Emigration Rates to account for rounding to the nearest thousand and included smoothing in the model. Detailed equations are presented in Online Resource 1 (section A.1).
The results of forecasting emigration rates and immigration counts for females are summarized in Figs. 5 and 6, respectively (the similar patterns for males are not shown because of space constraints). In the first row, we present the IE1 and IE2 forecasts for 2009. We observe that both models fit the data reasonably well. As expected, the univariate random walk model in IE1 leads to stable forecasts over time, with increasing uncertainty for both emigration rates and immigration counts. The drift term and log-linear trend incorporated in the IE2 lead to everincreasing immigration and emigration.

Population Forecasts
A final step in forecasting population is combining the population components within the cohort component projection model. As an illustration, we select models M2 for mortality, F2 for fertility, and IE2 for migration.
The age composition of the predicted UK population in 2024 is presented in the first row of Fig. 7. Forecasts of the total population for females and males are presented in the second row.
We observe that the age profile of the 2024 population is shaped mostly by future migration and, to a lesser extent, fertility. The largest uncertainty concerns the youngest population, as well as the population aged 20-45, for both males and females. These findings are in line with Keyfitz's (1981) observation on the plausible limits of population forecasting, which were set to about 20 years ahead. On the other hand, the forecasts suggest that the number of the elderly persons will be larger in 2024 but the population aged 20-45 will be most numerous.
[ Fig. 7 about here.] Our forecast can be compared with the official deterministic projections for 2024 prepared by the Office for National Statistics (2011). The dashed line in Fig. 7 represents the principal projection, whereas the dotted lines are low and high population scenarios. The main drivers of differences between these projections are assumptions about fertility and migration.
We observe that the population aged 20-35 is substantially smaller in all scenarios of the Office for National Statistics projections compared with our forecast. This results from the assumption that net migration stays constant at the levels observed in recent years: at the level of 200,000, 140,000, and 260,000 persons annually in the principal, low, and high scenarios, respectively (Office for National Statistics 2011). Moreover, our results suggest that the uncertainty about the age profile of future population is considerably larger than it is reflected in the projections based on the deterministic scenario.
As far as the total population size is concerned, the forecast indicates that there will be only 49,000 fewer females than males in 2024, whereas the difference was more than 1 million in 2009 and was 1.5 million in favor of females in 1975. This is most likely due to the larger proportion of male migration and a gradual closing of the life expectancy gap between the sexes.
The median size of the 2024 population is 70.8 million, which is around 9 million larger than the population size observed in 2009. In the principal projection for 2024 prepared by the Office for National Statistics (2011), the predicted total population is 69.0 million (i.e., 1.8 million lower), whereas the low and high scenarios are 66.4 and 71.1 million, respectively. However, only high scenario falls into our 80 % predictive interval for the total population. Also, the population predicted by the Office for National Statistics seems to be growing more slowly than in our forecast. This results from the aforementioned more-conservative prediction of international net migration.
As a further validation step in population forecasting, we present an in-sample population forecast in Fig. 8

Discussion
In this article, we have demonstrated that extension of the Lee-Carter model can serve as a general platform for estimating age schedules of the four demographic components of population.
We then combined these components into a single forecast by means of a cohort-component projection model. We also explored the correlation of each of the components in time, as well as between sexes and components (for emigration and immigration), which is embedded in our extension of the Lee-Carter method. For emigration and immigration, we provided a tool for smoothing irregularities in the data. This tool, however, can be easily extended to fertility or mortality. Finally, we have illustrated the use of the forecasting model on the UK's data.
This research makes two contributions to the literature. The first is the development of a new approach for integrating demographic components to provide stochastic population forecasts by age and sex. The Bayesian approach that we adopted accounts for the uncertainties embedded in births, deaths, and emigration and immigration, as well as across age and sexes. We show that the same general framework of the Lee-Carter approach for modeling age and sex patterns of mortality and fertility can be coherently applied to model corresponding patterns of migration.
Irregularities in the data, such as those observed for the UK, can also be accounted for within the model.
The second contribution is the application of the approach to a situation of relatively good yet imperfect data availability. In this way, we position our work between the generic global approach with far fewer data requirements, which has been proposed by Raftery et al. (2012), and a specific multiple-data situation discussed by Bryant and Graham (2013). Where possible, population forecasting should follow a bottom-up approach, in which the age-specific rates of the demographic components are utilized. The rates describe the underlying processes more comprehensively than summary aggregates, such as TFRs or life expectancies, in the topdown approach.
Further research should explore other models for forecasting age patterns of demographic components, such as the functional models developed by Hyndman and Booth (2008).
Analogously, various specifications for the time component models (such as ARIMA or VAR models of higher order) should be investigated. Next, the underlying models of components for the population forecast can be selected by using various techniques, of which Bayesian model averaging (Raftery et al. 1997) seems to be most appealing. In this way, the model uncertainty would be accounted for coherently. Finally, the uncertainty of the baseline population size used for projections could be incorporated into the projection. We believe this work provides a strong foundation for such extensions.