Probabilistic Estimation and Projection of the Annual Total Fertility Rate Accounting for Past Uncertainty: A Major Update of the bayesTFR R Package

The bayesTFR package for R provides a set of functions to produce probabilistic projections of the total fertility rates (TFR) for all countries, and is widely used, including as part of the basis for the UN's official population projections for all countries. Liu and Raftery (2020) extended the theoretical model by adding a layer that accounts for the past TFR estimation uncertainty. A major update of bayesTFR implements the new extension. Moreover, a new feature of producing annual TFR estimation and projections extends the existing functionality of estimating and projecting for five-year time periods. An additional autoregressive component has been developed in order to account for the larger autocorrelation in the annual version of the model. This article summarizes the updated model, describes the basic steps to generate probabilistic estimation and projections under different settings, compares performance, and provides instructions on how to summarize, visualize and diagnose the model results.


Introduction
In 2015 for the first time, the United Nations (UN) adopted the Bayesian method described by Alkema, Raftery, Gerland, Clark, Pelletier, Buettner, and Heilig (2011) for their official population projections for all countries, the World Population Prospects (WPP) 2015 (United Nations 2015). This method is probabilistic and based on a principled statistical footing, replacing the previous deterministic method. One of the major components is the projection of the total fertility rate (TFR) which is implemented in the bayesTFR R package (Ševčíková, The article is organized as follows. Section 2 summarizes the theoretical models developed by Alkema et al. (2011); Raftery, Alkema, and Gerland (2014); Liu and Raftery (2020), and the autoregressive model in the fertility transition phase. Section 4 describes how to use the package, using a step-by-step approach with different model settings. Section 5 presents experiments on the performance of the models and the selection of the various settings. The article concludes with a discussion in Section 6.

Annual TFR model with uncertainty about the past
Here, we first summarize the original TFR model developed for five-year time periods . We then review the new methodology for probabilistic estimation and projection of TFR for all countries of the world accounting for uncertainty about the past, as proposed by Liu and Raftery (2020). Finally, we describe the changes in the methodology to work for annual estimation and projections.
TFR can be defined as the number of children a woman would have if she were subject to the prevailing fertility rates at all ages from a single given year, and survived throughout her childbearing years. Alkema et al. (2011) defined a three-phase model for the evolution of TFR over time in a country: 1. Phase I: pre-transition phase with fluctuations at high fertility level.
2. Phase II: transition from high to low fertility, where decrements are modeled by a random walk with drift given by a double logistic function.
3. Phase III: post-transition phase where fertility fluctuates around the replacement level (a level close to 2.1), modeled by an autoregressive AR(1) process.
We will use the same notation asŠevčíková et al. (2011). Specifically, f c,t denotes the TFR in country c and time period t, τ c denotes the start period of Phase II for country c, λ c is the start period of Phase III for country c, while g(θ c , f c,t ) and θ c denote the parametric decline function and the corresponding country-specific parameters, respectively.

Existing model with five-year estimates
The pre-transition phase (Phase I) is not modeled, as all countries have already entered Phase II. Thus, for the purpose of projecting into the future it is not needed.
The fertility transition phase (Phase II) is modeled by a random walk with drift. This is specified by The decrement d c,t in (1) is modeled as the sum of a function of the level of the TFR and the noise, as follows: where g(θ c , f c,t ) are the double logistic decrements, which are determined by the countryspecific parameter vector θ c = (∆ c1 , ∆ c2 , ∆ c3 , ∆ c4 , d c ) and given by . ( The random distortions ε c,t in each period have normal distributions as follows: The quantity σ(f c,t ) is the standard deviation of the distortions during the later periods with σ(f c,t ) = c 1975 (t) σ 0 + (f c,t − S)(−aI [S,∞) (f c,t ) + bI [0,S] (f c,t )) .
The constant c 1975 (t) is added to model the higher error variance of the distortions before 1975. For further details about the model and its priors, seeŠevčíková et al. (2011). For the purpose of this article, we only point to the definition of two parameters, namely the country-specific maximum decrement d c , and the hyperparameter for the maximum standard deviation of the distortions σ 0 . The d c parameter is defined as The prior distribution of σ 0 is σ 0 ∼ U [0.01, 0.6].
The TFR in the post-transition phase (Phase III) is modeled by a first order autoregressive time series model (Raftery et al. 2014) as where µ c is the country-specific long-term mean fertility rate, and ρ c is the autoregressive parameter with ρ c ∈ (0, 1). In bayesTFR these parameters can be estimated via the Markov chain Monte Carlo (MCMC) method. Alternatively, country-independent values can be predefined or estimated by maximum likelihood.
The start period of Phase II, τ c , is defined as first estimation year, otherwise, where M c is the maximum observed TFR outcome in country c, and L c,t denote local maxima.
The start period of Phase III for country c, λ c , is defined as the period where two consecutive increases of TFR below 2 have been observed. More formally, λ c = min{t : f c,t > f c,t−1 , f c,t+1 > f c,t and f c,p < 2 for p = t − 1, t, t + 1} .

Probabilistic TFR estimation with uncertainty
The method described above uses observed TFR values as input to estimate the model parameters. In the UN context, these input values are taken from the latest revision of the WPP. Such TFR data are in fact estimates of the observed values, often derived from multiple data sources and involve varying amounts of uncertainty. The TFR model from the previous section however, treats these estimates as true values. Liu and Raftery (2020) developed a method that assesses the uncertainty around past estimates of the observed TFR values and propagates it into the projections. The medians of the resulting posterior distributions can be used as point estimates of the past TFR, reducing the need for manual analysis and assessments by individual UN analysts. In addition, TFR projections resulting from the method of Liu and Raftery (2020) show better out-of-sample validation, especially better coverage of the prediction intervals, than the existing method.
The method uses the World Fertility Data database (United Nations 2019a) for past raw TFR estimates from surveys, reports and vital registrations for most regions in the world. We denote these data points by y c,t,s , i.e., the raw TFR estimate for country c, time t and source s. The source s may refer to a census, a survey, vital registration statistics or other sources. For each observation, there are features x c,s that describe the sources, estimating methods, recall lags and other aspects of data collection and estimation, often measures of the quality of the data. The observed y c,t,s are modeled as: The δ c,s and ρ c,s are country-specific parameters which are estimated by maximum likelihood. In Liu and Raftery (2020), the features used are the sources of the data and the corresponding estimation methods, but the model allows for any user-specified features. This part is combined with the existing Bayesian hierarchical model implemented in bayesTFR. Here, the past TFRs are considered as unknown, and are part of the parameters to estimate. The complete model is described in the Appendix.
If we are using TFR for five-year intervals, as for example in the tfr dataset in the wpp2019 package (United Nations Population Division 2020), the true TFR at any time stamp is considered to be the linear interpolation of two consecutive five-year TFRs, namely If we are estimating from annual TFR, for each observation of the raw data, we take the floor of t. For example, if an observation in the raw data is recorded at 1975.5, we consider this observation as an estimate of the calendar year 1975.
Since we are now also modeling the past, not just the future as in the extant method, we need to model the pre-transition phase (Phase I), which is not necessary for projecting the future. The TFR in this phase will be modeled by a random walk, specified by where the distributions of the random distortions in each period are given by Here, we simplify the model by setting the variance to be the same as the variance in the first period of the TFR decrements. This is a reasonable assumption because the starting period of Phase II is linked to Phase I, and the expected decline of TFR at the starting period of Phase II is small. Thus, the distortions of TFR share similar behavior.
The estimation of all country-specific parameters and hyperparameters conditional on the TFRs, other than the TFRs themselves, in the Phase II model remains the same as described byŠevčíková et al. (2011). To estimate past TFR, the model for Phase III is estimated together with the Phase II model via an MCMC algorithm (Gelfand and Smith 1990). This algorithm is a combination of Gibbs sampling, Metropolis-Hastings (for ∆ ci in (3) and TFR), and slice sampling steps (Neal 2003).
The estimation yields a set of TFR trajectories about the past. To project into the future, we apply the existing projection method as described inŠevčíková et al. (2011) starting with the last estimation period of each trajectory. This is in contrast with the previous version, where the projection for each country starts from a single data point, namely the last observed TFR.

Annual version of bayesTFR
The original model described in Section 2.1 was designed to work with five-year data. Several modifications needed to be made in order to estimate and project annual data well.
Most importantly, we found that unlike in the five-year version, the residuals of the Phase II model are highly autocorrelated when using annual data. We found that the lag 1 autocorrelation coefficients are about 0.7 for residuals of Phase II model for some major countries. Thus, we modified the Phase II model defined in Equations (1)-(2) by adding an additional first-order autoregressive component. The random walk with drift model is then specified as The prior distribution of φ is set to be Uniform(0, 1) and is not country-specific. For the random distortions ε c,t , the distribution is considered to be the same as in Equations (4)- (5).
The same prior distributions as in the five-year version is used for most parameters. One exception is σ 0 where the lower bound was decreased by a factor of the square root of five, i.e., σ 0 ∼ U [0.0045, 0.6]. The upper bound was kept the same to allow for the possibility of additional correlation.
The definition of the maximum decrement defined in Equation (6) was changed to be one-fifth of that for the five-year model: No changes have been made to the model of the post-transition phase of TFR, Phase III. It is modeled by a first-order autoregressive time series model as defined in Equation (7).
The rule for determining the start period of Phase II, τ c , as defined in Equation (8), was unchanged. However, since the local maxima are calculated using annual TFR data, the results can differ from those obtained from a five-year dataset.
To determine the start periods of Phase III, λ c , as defined in Equation (9), we first obtain five-year averages of TFR and apply the same rule as in the five-year version, namely that Phase III starts when two consecutive increases of TFR below 2 are observed.

Changes in TFR Projections
There are three main differences in the TFR projections between the new implementation and the one described byŠevčíková et al. (2011).
The first difference (which we alluded to at the end of Section 2.2), relates to the fact that by accounting for the past TFR uncertainty (Equation 10), instead of observed point estimates, the model results in a set of TFR trajectories about the past which changes the starting values of the forecast. To project f c,T +1 where T is the last period of the estimation, the i-th sample from the MCMC output is given by f c,T . Thus, the uncertainty in the first forecast period will be wider than if we use a model without accounting for past uncertainty, in which case f (i) c,T = f c,T is the same for all trajectories. The second difference relates to the annual model described in Section 2.3. When the additional autocorrelation of Phase II is taken into account (Equation 11), the past noise is carried over to the next time period. Specifically, to project f c,t+1 for a country c that is in Phase II at time t, the i-th sample is given by f c,t )). For the first forecast, i.e., at the time period T + 1, the distortion of the last estimation period T is calculated before starting the projections.
Finally, the last difference regards the updated Phase III model as described in Raftery et al. (2014), where country-specific long-term means µ c and autocorrelation coefficients ρ c were incorporated into the model (Equation 7) and estimated by MCMC. However, this change has been available in bayesTFR since version 3.0-0 was published in 2013. Using this approach, to project f c,t+1 for a country c that is in Phase III at time t, the i-th MCMC sample is drawn from a Normal distribution N µ

Overview of the Package Updates
The updated package bayesTFR retains all the functionalities of the previous version, which implements the model of Alkema, Raftery, Gerland, Clark, and Pelletier (2012). Its new features allow the user to conduct estimation of past TFR while accounting for uncertainty as described in Liu and Raftery (2020), as well as to forecast TFR for both five-year and one-year time intervals, as requested by the UN.
These new functionalities are implemented in the form of either new functions or additional arguments to existing functions. For convenience, especially for users who are familiar with the previous version of the package, this section summarizes the various changes. For users who are new to the package we recommend skipping to Section 4 where the usage is explained in more detail.
The following are established bayesTFR functions that were updated.
• run.tfr.mcmc. This is the core function for MCMC estimation of fertility transition model parameters. The following optional arguments were added: 1. annual. Logical argument determining whether the model is trained based on annual TFR data (TRUE) or on the five-year data (FALSE). The default is FALSE. 2. ar.phase2. Logical argument. If TRUE, model (11) will be used in the estimation, and the parameter φ will be estimated through the MCMC process. This is relevant only if annual = TRUE. 3. uncertainty. Logical argument determining whether the model described in Liu and Raftery (2020)  • tfr.predict. This is the core function for TFR prediction. There was one optional argument added: uncertainty. Logical argument. If TRUE and the corresponding estimation was produced via run.tfr.mcmc(..., uncertainty = TRUE), then each prediction trajectory starts from a trajectory representing the past. Otherwise all prediction trajectories start from the same point, namely the last observed TFR.
• run.tfr.mcmc.extra. Originally, this function has been implemented in order to estimate the TFR transition model for very small countries or countries with unusual historical patterns. These countries were excluded from run.tfr.mcmc in order not to bias the world parameters, and were estimated separately via this function. In this update, the function has been extended to recompute past TFR estimates on a countryspecific basis. This allows users to analyze the impact of changes on the raw TFR of individual countries without needing to run a new simulation for the whole world. Added arguments to this function have the same meaning as for run.tfr.mcmc: uncertainty, my.tfr.raw.file, covariates and cont_covariates, iso.unbiased • tfr.trajectories.table, tfr.trajectories.plot. These functions give projection quantiles in tabular and graphical formats, respectively. They have been extended to include uncertainty information about the past, if such information exists.
The following are new functions added to the package. They are described in Section 4.4 in more detail.
• get.tfr.estimation. Allows exploring the estimated trajectories as well as any quantiles of the past TFR estimates.
• tfr.estimation.plot. Function for plotting estimates of past TFR for individual countries.
• tfr.bias.sd. Allows exploring the bias and standard deviation estimated for the raw TFR estimates.

Using the Updated bayesTFR
Previous versions of the bayesTFR package which implemented the model described in Section 2.1 have been used by UN analysts and others to train TFR projection models based on past five-year estimates. New UN requirements added the need to update the package so that analysts can conduct estimation of past TFR accounting for uncertainty, and make corresponding forecasts for both five-year and annual time periods.
To make probabilistic TFR projections accounting for past TFR uncertainty, the user will follow four steps in the following order: 1. Data assembly (optional): (a) Prepare a dataset of raw TFR values. By default, the World Fertility Data 2019 (United Nations 2019a) is used.
(b) Assemble a dataset of reference (initial) TFR values for all countries and time periods. By default, the the World Population Prospects (United Nations 2019b) in the wpp2019 package is used.

Model estimation:
(a) Train linear models to estimate systematic bias and standard deviation for each observation from the raw TFR dataset.
(b) Given the reference TFR, calculate the start period of Phase II and the start period of Phase III for each country (τ c and λ c ).
(c) Run the MCMC process to obtain posterior samples of the Phase II and Phase III model parameters, and posterior samples of the past TFR for all countries.
3. Generate future TFR trajectories as discussed in Section 2.4.
4. Analyze the outputs using a set of functions that summarize, plot, diagnose and export the results of the three steps above.
As described byŠevčíková et al. (2011), steps 2 and 3 are relatively time-consuming. Adding the estimation uncertainty feature as well as working with annual estimates adds even more run time. Even though we optimized the code wherever possible, it takes several hours to complete these steps in a production-like setting.
The following sections describe the steps above in more detail, especially the parts that are different fromŠevčíková et al. (2011). We will elaborate on how to use the new features, as well as how to use the package in the original way. We will demonstrate the package on a realistic example with a relatively large number of MCMC iterations, which might take several hours to process. Therefore, users who wish to explore the functionality quickly should reduce the number of iterations to the order of magnitude of 10. However note that since the Metropolis-Hastings step for the TFR updates will have an acceptance rate of around 30%, a small number of iterations will result in estimation plots that are less smooth than what we will present in this article.

Data assembly and estimation settings
The datasets assembled in this step will be passed to the main estimation function, run.tfr.mcmc, which now has additional arguments for this purpose. It can be specified what raw TFR data to use, whether to estimate and predict annually (logical argument annual), and whether to use the AR(1) model in Phase II as defined in Equation 11 (logical argument ar.phase2).
The argument uncertainty = TRUE specifies that uncertainty about the past is incorporated into the estimation. In this case, a raw TFR dataset can be provided. By default, the World Fertility Data 2019 (United Nations 2019a) is used. This dataset contains 12,079 records for 201 countries, each of which includes the corresponding estimation method and data source. These are then used by the model as data quality covariates in Equation (10 The default covariates are c("source", "method"). Users can provide their own file and covariates of their choice. Required columns are "country code", "year" and "tfr". The name of this file is passed to the argument my.tfr.raw.file, names of categorical variables to the argument covariates, and names of continuous variables to the argument cont_covariates.
An additional option allows an analyst to specify that vital registration data from selected countries are unbiased, if there is a belief that these data are not systematically biased in a particular direction. (Note that this is not the same as saying that these data are perfect.) The UN country codes of those countries can be passed into the argument iso.unbiased.
In this case, the bias and standard deviation of the records of those countries for which the source column specifies "VR" (as Vital Registration) are forced to be equal to 0 and to be near 0, respectively (in practice the standard deviation is set to 0.0161). This option targets finetuning of the estimation of developed countries, especially because the annual TFR estimates are often not openly accessible.
The second dataset to assemble is a dataset on a reference, or initial, TFR. Its file name is passed into the argument my.tfr.file. If uncertainty = FALSE, this dataset is considered in the estimation as the true observed TFR. Otherwise, it is used as the starting points of TFR in the MCMC process. By default, if my.tfr.file is not given, the tfr dataset from the wpp2019 package is used for this purpose, which is a five-year dataset. Thus, if annual = TRUE, a linear interpolation of the default dataset is computed.

Fitting the TFR model
Most arguments of run.tfr.mcmc remain the same as described inŠevčíková et al. (2011). Importantly, start.year and present.year set the first and the last year of the time series included in the computation, respectively. The arguments nr.chains, iters and output.dir determine the number of chains, the number of iterations and the directory to store the MCMC simulated values, respectively.
In the prior version of bayesTFR, the function run.tfr.mcmc was designed to obtain a posterior sample of Phase II model parameters. The estimation of Phase III parameters (as defined in Equation 7) is implemented in the function run.tfr3.mcmc. When building a full probabilistic model as described in Liu and Raftery (2020), the MCMC steps for updating TFR will affect both phases. Thus, if uncertainty = TRUE, the new run.tfr.mcmc function combines the estimation of both phases, and there is no need to invoke the run.tfr3.mcmc function explicitly. We call this method a "one-step estimation". However, this is not the case if uncertainty about the past is not taken into account. In this case, the workflow of estimating Phase II and Phase III separately remains and is referred to as a a "two-step estimation".
The various combinations of the possible settings of the arguments annual and uncertainty are summarized in Table 1. We have marked each cell with a letter which will be referred to in the remainder of this section.
annual uncertainty TRUE FALSE TRUE A B one-step estimation; two-step estimation; Phase II -AR(1) allowed Phase II -AR(1) allowed FALSE C D one-step estimation two-step estimation

Starting a new simulation with two-step estimation
The two-step estimation should be performed if uncertainty about the past is not taken into account (cells B and D in Table 1). The main differences between cells B and D are the setting of prior distributions as described in Section 2.3, and whether the autoregressive component can be included in the model. Here we give an example of a simulation with annual data (cell B) without the autoregressive component. However, we will not use this example further in the text, as the main focus of the article is on cell A which will be discussed in the next section.
Our example simulation consists of three MCMC chains, each of which is 5000 iterations long where thinning is disabled. (As noted earlier, the user is advised to decrease the number of iterations to the order of ten for faster processing.) We will save the simulation results to a directory called "annual".
R> annual <-TRUE R> nr.chains <-3 R> total.iter <-5000 R> thin <-1 R> simu.dir <-file.path(getwd(), "annual") The first step is to launch an estimation of Phase II: R> m2 <-run.tfr.mcmc(output.dir = simu.dir, nr.chains = nr.chains, + iter = total.iter, thin = thin, annual = annual) The second step is to start an estimation of Phase III: R> m3 <-run.tfr3.mcmc(sim.dir = simu.dir, nr.chains = nr.chains, + iter = total.iter, thin = thin) Here, we are using the same number of chains and iterations for Phase II and Phase III. However, this is not a requirement, but rather depends on the MCMC convergence. Even the 3 × 5000 iterations might be not enough to reach convergence, but will usually give realistic outputs. Setting total.iter = 62000 or total.iter = "auto" will most likely result in full convergence.

Starting a new simulation with one-step estimation
We now show an example of a simulation with uncertainty which is performed with one step only (cells A and C in Table 1). In particular, here we set annual to TRUE (cell A), but the same function would be used if annual is FALSE (cell C). We will also include the Phase II -AR(1) model (ar.phase2) which would not have any effect in cell C. The results will be saved in the directory "annual unc". We will use this simulation throughout the article.
R> annual <-TRUE R> ar.phase2 <-TRUE R> nr.chains <-3 R> total.iter <-5000 R> thin <-1 R> simu.dir.unc <-file.path(getwd(), "annual_unc") As in the previous case, this setting may not be enough to yield fully converged MCMC simulations, but will still give realistic outputs. The processing time is within a range of a couple of hours. For faster processing, set total.iter = 50 for a toy simulation. In addition, the parallel argument can be set to TRUE, in which case the three chains will be processed in parallel. In Section 5.2, we will give recommendations for settings that yield fully converged MCMC simulations. When appropriate, we will use results from such converged simulations to present outputs.
As mentioned in Section 4.1, additional arguments of run.tfr.mcmc allow one to pass userspecific raw TFR data (my.tfr.raw.file), names of categorical covariates (covariates), names of continuous covariates (cont_covariates), or to specify countries with unbiased vital registration data (iso.unbiased). If the iso.unbiased option is used, the covariates argument should include the variable "source", or more specifically, the variable defined by the argument source.col.name which is "source" be default. In such a case, the function reduces the bias and standard deviation of records where the "source" column specifies "VR". In our example we will specify that the vital registration data of Canada and the USA (codes 124 and 840) are unbiased. In comparison to the two-step model, the training process here has an extra Metropolis-Hastings step per iteration for generating posterior TFR samples.

Continuing an existing simulation
If an existing simulation needs to be extended by more iterations, one would proceed as in the previous version of bayesTFR: • Use continue.tfr.mcmc if the MCMCs were originally created via run.tfr.mcmc, regardless of whether one is in the one-step or the two-step estimation mode.
Now suppose we want to extend the simulation in the previous section by 100 iterations. Then we could do R> m <-continue.tfr.mcmc(output.dir = simu.dir.unc, iter = 100) (Set the iter argument to 10 if working with a toy simulation.) This will continue running both TFR phases in a one-step estimation while inheriting all settings from the original simulation, including uncertainty, annual, ar.phase2 or settings about the raw data. At the end of the processing, each chain will be 5,100 iterations long.

Generating projections
The main function for generating projections is called tfr.predict. The new version of the package adds the argument uncertainty. If it is TRUE and the model was estimated taking uncertainty about the past into consideration, then that past uncertainty will be carried over to the projections. In this case, each future trajectory starts from a trajectory representing the past.
Suppose we want to generate projections represented by 1,000 posterior trajectories until the year 2100 based on the simulation stored in the directory given by simu.dir.unc, with burn-in of the first 2100 iterations for each chain. This can be done using the following command: R> pred <-tfr.predict(sim.dir = simu.dir.unc, end.year = 2100, + burnin = 2100, nr.traj = 1000, uncertainty = TRUE) The function takes the existing 5,100 iterations in each chain, removes the first 2,100 values and generates 1,000 TFR trajectories based on 1,000 equally spaced parameter values and past TFR, out of the remaining 3, 000×3 = 9, 000 iterations. For a toy simulation, use burnin = 20 and nr.traj = 50.
If uncertainty is set to FALSE, all future trajectories start from the last observed data point. If the estimation process accounted for uncertainty, but the projection does not, the starting value of the projections is the initial TFR value for the last observed time period. This is however not recommended but may serve the purpose of apples-to-apples comparisons.

Analyzing Results
If results are to be explored at a later time point, one can load the estimation object from disk using the command

Summary functions
For the summary statistics of the estimation object in this section, we will use the following thinning and burn in settings:
To view a summary of country-independent parameters, one can use R> summary(m, thin = thin, burnin = burnin) Since the object m was got by one-step estimation, the output includes estimation summaries for both phases. In comparison to the previous version of the package, Phase II contains one additional parameter, namely "rho phase2" which represents φ in model (11). As with any other parameter, the name, or multiple parameter names, can be passed to the function to view summary statistics for those selected parameters.
R> summary(m, par.names = c("rho_phase2", "sigma0"), + thin = thin, burnin = burnin) [1] "alpha" "alphat" "delta" "Triangle4" "delta4" "psi" "chi" [8] "a_sd" "b_sd" "const_sd" "S_sd" "sigma0" "mean_eps_tau" "sd_eps_tau" [15] "rho_phase2" Passing the meta argument is needed to identify that the simulation contains a Phase II -AR(1) estimation, and thus it contains the "rho phase2" parameter. Phase III parameter names are not dependent on the simulation, thus no meta argument is needed: R> tfr3.parameter.names() [1] "mu" "rho" "sigma.mu" "sigma.rho" "sigma.eps" Specifying a country in the summary function will show results for the country-specific parameters of that country. This is done via the country argument which accepts either the name or the numerical code of the country, as well as an ISO-2 or ISO-3 character code. This is the case for any function in the package that accepts the the country argument, as is shown throughout the paper.
As for the parameters to summarize, functions tfr.parameter.names.cs() and tfr3.parameter.names.cs() list the allowed parameter names for Phase II and Phase III, respectively. For a simulation that took into account uncertainty about the past, there is an additional country-specific parameter, called "tfr," capturing that uncertainty. It is not listed explicitly via the above functions, but it can be explored like any other parameter. For the summary function it means passing it to the par.names.cs argument. For example, to view summary statistics of TFR estimation for Nigeria, we can do R> summary(m, country = "Nigeria", par.names.cs = "tfr", + thin = thin, burnin = burnin) The tabular sections of the output contain one row per past observed period each (by default 71, i.e., from 1950 to 2020). To select a subset we can specify which time periods we are interested in as tfr_time. For example, to view results for time periods 1, 30 and 70 (corresponding to 1950, 1979 and 2019) we do R> summary(m, country = "Nigeria", + par.names.cs = c("tfr_1", "tfr_30", "tfr_70"), + thin = thin, burnin = burnin) 2.5% 25% 50% 75% 97.5% tfr_1_c566 5.765 6.135 6.290 6.438 6.731 tfr_30_c566 6.577 6.655 6.704 6.760 6.879 tfr_70_c566 4.827 5.113 5.765 5.916 6.320
Several arguments in this function need to be clarified: 1. sim.dir: Users can specify the location of the simulation set, or use the mcmc.list argument to pass the m object directly. For example tfr.estimation.plot(mcmc.list = m, ...).
2. pis: Specifies which probability intervals will be plotted. It is a vector of at most two elements.
3. plot.raw: Logical argument which determines whether the raw data used for estimating past uncertainty are plotted. If TRUE and the estimation process was not based on the default data, it is recommended to provide the argument grouping, which should be one of the categorical covariates in the raw data set. This covariate defines the color and shape of the points, as seen in Figure 1 where the default grouping is the "source" column of the rawTFR dataset.
4. save.image: (not used in the call above) If TRUE, the plot will be saved as a pdf file in the directory specified in the argument plot.dir, named "tfr countrycode.pdf".

Exploring bias and standard deviation of observations
Information about the bias and standard deviation of observations will give users an indication of the quality of the observations and whether these quantities were poorly estimated. Now suppose we are interested in the bias and standard deviation estimates of the observations for Nigeria. Then we could use R> bias_sd <-tfr.bias.sd(sim.dir = simu.dir.unc, country = 566) The function will return a list with elements model_bias, model_sd and table. The model_bias and model_sd objects are of class lm and contain the linear models used to estimate the bias and standard deviation, respectively, while the table object includes the observed data points, data quality covariates, and the actual estimates for the specified country, here for Nigeria.

R> summary(bias_sd$model_bias) R> head(bias_sd$table)
The results are shown in Tables 2 and 3   To generate the estimates in the table object, the predict S3 method is applied to the model_* objects. Then the following steps are performed: 1. For some countries, the number of data points is very small for several groups. This could lead to a large bias, but a very small variance. As a result, the estimation will be unreasonably concentrated on the bias-adjusted data points. In this case, the standard deviation estimates were adjusted as max(0.1, |bias|/2). The reason for such an adjustment is that it is unlikely that one observation is very biased but with a very small relative standard deviation. It is also unlikely that there is a source of data that is very precise (with standard deviation less than 0.1), but is only collected once.
2. For countries included in iso.unbiased, the model estimates are overwritten with zero or close to zero values as explained in Section 4.1.
3. Duplicates are dropped so that the combinations of data quality covariates are unique.
The output can help to detect problematic estimates on certain data points so that adjustments can be made by the analyst if necessary. In the example above, the estimated bias and standard deviation for the Indirect method and nationwide DHS surveys were 0.06 and 0.09, respectively. These estimates were derived based on only three data points in this category, and all of them were very close to the UN estimates (three of the brown dots in Figure 1 in year 1969, 1974 and 1979). Since the number of data points from nationwide DHS estimates is small (3 data points), the estimated bias (0.06) and standard deviation (0.09) may be too small.

Exploring TFR prediction
Plotting Here, the parameter uncertainty is used to specify whether the uncertainty about the past TFR should be plotted together with the prediction. If uncertainty is TRUE, the optional parameters thin, burnin, col_unc can be used to define the burn-in, thinning and the color for the past uncertainty plot.
The code above applied to a converged simulation results in the plots shown in Figure 2. If the user selects uncertainty = FALSE for a simulation where past uncertainty was taken into account (similarly to the right panel of Figure 2), the past TFR used for the initialization of the model is shown as the observed TFR. In this case, there could be a discontinuity between the last observed and the first projected time period.
These new arguments are also accepted by the tfr.trajectories.plot.all function which generates projection plots for all countries at once, as described byŠevčíková et al. (2011).
TFR predictions in a tabular format can be explored using the tfr.trajectories.table and summary functions which work the same way as in the previous versions of the package, except that in the former case, the output now contains uncertainty information about the past. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q median PI 95 PI 80 Phase II data Phase I data

Exploring double logistic function
The double logistic function defined in (3)  Results can be seen in the left panel of Figure 3, while the right panel shows the result of the same call with country = "Thailand".
If a simulation contains information about past uncertainty, then the Phase II and I data (black dots and squares) represent decrements of the estimated TFR median. In case of an annual simulation, these are annual decrements, otherwise they would correspond to five-year decrements. If the projections were produced without taking past uncertainty into account, then the data points represent the observed decrements.
This also applies to the DLcurve.plot.all function which plots the double logistic curves for all countries at once.

MCMC traces, density and diagnosis
To explore traces of the MCMC parameters, the existing functions tfr.partraces.plot (for country-independent parameters) and tfr.partraces.cs.plot (for country-specific parameters) can be used. Similarly, for density plots, tfr.pardensity.plot and tfr.pardensity.cs.plot are available.
As mentioned previously, there are two additional parameters in this version of the package, namely "rho phase2", which is country-independent and defined in model (11), and "tfr" which is a country-specific parameter. These two parameters can be used within the aforementioned functions, like any other parameters . For example, the trace plots and the density plots of φ and Nigeria's TFR estimate in year 1985 (as shown in Figures 4 and 5 ) can be visualized via R> tfr.partraces.plot(m, par.names = "rho_phase2", nr.points = 200) R> tfr.partraces.cs.plot(m, country = "Nigeria", + par.names = "tfr_36", nr.points = 200) R> tfr.pardensity.plot(m, par.names = "rho_phase2", burnin = burnin) R> tfr.pardensity.cs.plot(m, country = "Nigeria", + par.names = "tfr_36", burnin = burnin) To check if the MCMC algorithm has converged and adequately explored the parameter space, the tfr.diagnose function can be used; seeŠevčíková et al. (2011) for more details.
In the case of one-step estimation, the function checks parameters from Phase II as well as Phase III. In the case of two-step estimation, one would use tfr.diagnose for assessing the convergence of Phase II parameters, and tfr3.diagnose for assessing the convergence of Phase III parameters. Both functions accept a logical argument express which can disable or reduce the checking of country-specific parameters in order to speed up the process.
If the estimation includes uncertainty about the past, the assessment of country-specific parameters include the "tfr" parameter for each country and time period, in our case more than 14200 "tfr" parameters. In practice, it is often impossible to achieve convergence for all of them. Thus, we introduced the rule of accepting the "tfr" parameters as having converged if 95% of them have converged.
To apply the convergence diagnostics to our simulation, one could do R> tfr.diagnose(simu.dir.unc, thin = thin, + burnin = burnin, express = TRUE) As mentioned earlier, in our illustrative code examples the MCMC algorithm has not been run for long enough to achieve full convergence. See Section 5.2 for alternative settings. Note that the toy simulation we proposed earlier cannot be checked for convergence, as there is a requirement of a minimum number of iterations per chain, which the toy simulation does not satisfy.

Estimating a small set of countries
The Bayesian framework we have shown so far is designed to estimate all countries of the world at once, where the historical experience of an individual country influences the distribution of its own parameters as well as of the world parameters, while using the same settings for all countries. However, this is not always practical for several reasons: 1. Analysts might want to experiment with settings for individual countries without waiting several hours for a simulation of the whole world to finish.
2. Different sets of covariates might be needed to estimate different countries.
3. Countries with unusual historical patterns or very small countries might be excluded from the simulation in order not to bias the world parameters.
It was the last reason, as well as the need to include aggregations in the estimation, that motivated us to implement the run.tfr.mcmc.extra function in the original version of the package. The idea is that, while run.tfr.mcmc updates all parameters, the run.tfr.mcmc.extra function updates only the country-specific parameters of the specified countries, while re-using the existing distribution of the global parameters.
Since the function was designed for special cases of countries or aggregations, the original implementation allowed the user to process only the locations that had not been included in the world simulation. With the two additional use cases above, we have now relaxed that restriction and made it possible to rerun and overwrite existing estimations of country-specific parameters and past TFR estimates for individual countries, while allowing the user to change various estimation settings. However, several global settings are not subject to change, such as switching between annual and five-year estimation, or changing the ar.phase2 argument.
Suppose that after running the simulation with the default data from the World Fertility Data, the user wishes to experiment with their own data that exclude Nigeria's questionable data points, such as the Indirect DHS-NS data points identified in Table 3 as having unreasonably low standard deviations and biases. Unlike in the main simulation, the experiment will not force the VR data of the United States to have zero bias and variance. For that purpose, we will extract data for Nigeria (code 566) and the USA (code 840) from the default raw dataset discussed in Section 4.1, remove the Indirect DHS-NS points for Nigeria and store them into a file called "raw tfr user.csv": R> countries <-c(566, 840) R> myrawTFR <-subset(rawTFR, country_code %in% countries) R> myrawTFR <-subset(myrawTFR, !(country_code == 566 + & method == "Indirect" & source == "DHS-NS")) R> write.csv(myrawTFR, file = "raw_tfr_user.csv", row.names = FALSE) For experimentation with the run.tfr.mcmc.extra function, we recommend copying the main simulation into a different directory and applying the function to the copy. This is because the processing overwrites the existing estimation results, and thus there is no way back to the original results in case the experiments do not yield satisfactory outputs. Here we will append " extra" to the directory name stored in simu.dir.unc and copy the content from simu.dir.unc into it. This step is equivalent to the command "cp -r annual_unc annual_unc_extra" on unix-based systems: R> simu.dir.extra <-paste0(simu.dir.unc, "_extra") R> dir.create(simu.dir.extra) R> file.copy(list.files(simu.dir.unc, full.names = TRUE), + simu.dir.extra, recursive = TRUE) To run the new estimation for the two selected countries, we can do R> run.tfr.mcmc.extra(sim.dir = simu.dir.extra, countries = countries, + iter = total.iter, burnin = burnin, uncertainty = TRUE, + my.tfr.raw.file = "raw_tfr_user.csv", + covariates = c("source", "method")) We recommend using the same values of total.iter and burnin as in the main simulation.
To compare the new estimation results to those shown in Figure 1 we again use the tfr.estimation.plot function, now passing simu.dir.extra into the sim.dir argument. It can be seen in the left panel of Figure 6 that excluding the Indirect DHS-NS data points for Nigeria changed the estimates, especially for 1979. The uncertainty increased for the United States (right panel of Figure 6), since it was removed from the iso.unbiased set.
Finally, the option uncertainty = TRUE can be used even in two-step estimation where uncertainty about the past was not taken into account. This is possible because we do not expect the global parameters to be significantly different in the two situations (i.e., with and without uncertainty).

Structure of the output directory
Having a look at the simulation directory, here "annual unc", one should see a structure similar to the following: annual_unc bayesTFR.mcmc.meta.rda diagnostics mc1 mc2 mc3 phaseIII bayesTFR.mcmc.meta.rda mc1 mc2 mc3 predictions thinned_mcmc_9_2100 bayesTFR.mcmc.meta.rda mc1 The directories "mc1", "mc2" and "mc3" on the first level are generated by the run.tfr.mcmc function and contain results from the three chains of the Phase II estimation. Each of the directories contains one text file per parameter. The names of the hyperparameters and their corresponding notation are the same as described in Table 1 inŠevčíková et al. (2011). In addition, the parameter "rho phase2" representing φ from Equation 11 is also stored as a hyperparameter if the Phase II-AR(1) is considered. The names of the files storing countryindependent parameters consist of the parameter name and the suffix ".txt", while in the case of the files storing country-specific parameters the parameter name is followed by the suffix " countrycode.txt".
If uncertainty is taken into account, the MCMC algorithm also generates estimates for the past TFR data. These samples are considered as country-specific parameters, called "tfr", and thus stored in files "tfr countrycode.txt". They contain matrices of size the number of (thinned) iteration times the number of time periods. In the example above, the default starting year is 1950, and the present year is 2020, i.e., 71 years. Therefore, each file contains TFR estimates in 5100 rows and 71 columns.
The file "bayesTFR.mcmc.meta.rda" on the first level stores meta information about the Phase II estimation, which is contained in the m$meta object. If uncertainty is taken into account, the raw data used to obtain the estimates of TFR are stored as an additional element, called raw_data.original. A logical element ar.phase2 indicates whether the autoregressive component of Phase II is considered in the estimation. In order to allow users to work with different subsets of countries with the same base of global estimates, information indicating whether the countries were processed separately has been also stored in the meta object. It is accessible via the extra element, created only if the run.tfr.mcmc.extra function has been invoked and if uncertainty is TRUE. Here, extra_iter and extra_thin are used to retrieve the settings for specific countries. The raw data in this case are stored in a list called raw_data_extra. It is overwritten every time run.tfr.mcmc.extra is called for the same country.
The results of Phase III are stored in the directory "phaseIII". It has the same structure as described above. It is generated either by the run.tfr.mcmc function in case of a one-step estimation, or by the run.tfr3.mcmc function, in case of a two-step estimation. The meta file contains meta information related to the Phase III estimation. In the "mcx " directories, the names of the hyperparameters and their notations for Phase III are listed in Table 4. Similarly, the country-specific parameters and their notations are listed in Table 5. All files in this case contain one value per (thinned) iteration. Note that the country-specific parameters for Phase III are only estimated for countries which are already in Phase III, which in our case is 41 countries.μρ σ µ σ ρ σ ε mu rho sigma.mu sigma.rho sigma.eps Table 4: Country-independent parameters for Phase III in model 7, with their corresponding names in the code. They can be obtained using tfr3.parameter.names(). µ c ρ c mu.c rho.c Table 5: Country-specific parameters for Phase III in model 7, with their corresponding names in the code. They can be obtained using tfr3.parameter.names.cs().
The "predictions" directory is created by the pop.predict function and it holds binary files, one per country, each containing the predicted TFR trajectories for that country.
Other convenience directories might have been created for speeding up processing. For example, the "thinned mcmc 9 2100" directory was created by pop.predict to hold the final chain for each parameter derived by applying the burnin, thinning and collapsing the three chains into one, in order to generate the predictions. Since we asked to generate 1,000 posterior TFR trajectories with burnin of 2,100 iterations, a thinning of 9 was applied to retrieve those trajectories: 3 · (5, 100 − 2, 100)/9 = 1, 000. Thus, the parameter files in the "mc1" subdirectory here all contain 1,000 rows. Note that these values will differ when working with a toy simulation.
If functions for convergence diagnostics have been used, the simulation directory contains a folder "diagnostics" which holds results from these runs, one file per unique combination of thin and burnin.

Experiments
We have shown how the updated bayesTFR package can handle different versions of the TFR projection model. In this section, we will present results of experiments under different settings and discuss the implications of these settings. Based on those experiments we will give recommendations for a reasonable configuration of the model. Finally, we will discuss future directions in the development of the package.

Experiments with Settings
The new version of bayesTFR allows the user to handle different types of modeling needs, summarized in Table 1. An analyst can choose between a five-year and an annual model, as well as between accounting for past uncertainty or not. Flexibility is added by allowing the user to treat vital registration (VR) records for selected countries as unbiased, as well as using the autoregressive component in Phase II.
However, a question of consistency of results between the various settings may arise. For example, a forecast should not change dramatically when switching from five-year to annual data. Currently, there are no annual observations collected for all countries, and only a few countries (such as New Zealand) have good annual vital registration data, the only available annual observations. Thus, if past uncertainty is not taken into account the model would be estimated on some version of interpolated data for most countries.

Countries in Phase III with good records
The first major difference can be seen for countries in Phase III, especially for countries with high quality VR records. We take Switzerland as an example. The left panel of Figure 7 shows TFR projections for a five-year model without accounting for past uncertainty (cell D in Table 1), while the right panel shows results from an annual model with uncertainty about the past (cell A in Table 1). It can be seen that the results on the right yield wider probability intervals. For countries like Switzerland, the bias and uncertainty of past estimation is very low. Since the estimating process takes the linear interpolated TFR as the reference, the process can add extra bias to these data. Even though this is not large, the uncertainty propagated from the beginning of the forecast period could lead to a large difference. It can be seen that when compared to results from a five-year model (left panel), the differences between the two sets of projections are negligible. It is important especially for countries with nearly perfect historical data, such as Switzerland, that similar results be obtained whether annual or five-year data are used.

Countries in Phase II
The second major difference relates to countries in Phase II, such as Nigeria. Figure 9 shows the difference between a projection resulting from a five-year model without accounting for past uncertainty (left panel) and from an annual model with uncertainty about the past without applying the Phase II-AR(1) component.
It can be seen that if we account for uncertainty and use annual data, the prediction shows a faster decline. Without performing an out-of-sample validation, it is impossible to say which of these projections is better. Nevertheless, a more detailed analysis revealed that the posterior median of the residuals ε c,t for all countries in model (2) is highly autocorrelated. Figure 10 summarizes the estimates.
This suggests including the autocorrelation process in the modeling as defined in Equation (11). Figure 11 summarizes the differences. The decline has become slower, which is more in line with the five-year projections.
It can be seen however, that the starting point of the projections (year 2020) is now lower, and in fact it is significantly lower than the current UN estimates. The standard deviation of ε in model (11) is less than 0.02, if the autoregressive component is included. This could be problematic, given that for developed countries with low TFR and relatively stable societies, the standard deviation of annual TFR changes is about twice as much as 0.02. This is likely a result of a possible smoothing of the data. To remedy that, we introduce a new lower bound on the σ 0 parameter (argument sigma0.min in run.tfr.mcmc) of 0.04, which becomes the new default. Figure 12 shows the relevant differences.   If the lower bound on σ 0 is applied, the prediction yields wider probability intervals as well as a higher median (top right panel), which better matches the five-year forecast. The estimation in this case (bottom right panel) also shows a better match with the raw data as well as with the UN estimates, which is another argument for using the new default for sigma0.min.

Recommendations
We have shown the flexibility of the new version of bayesTFR which can incorporate different variations of the TFR model as well as being compatible with the extant version of the model. As one of the key components in population projections currently adopted by the United Nations, this is a key step for migrating population projections from a five-year basis to an annual one. The package is designed to support UN analysts in this process, as well as to give other researchers and practitioners a tool to generate their own projections.
In addition to incorporating past uncertainty of TFR in the forecast, and performing annualbased projections, the package has introduced two other important components, namely the ability to specify vital registration data as unbiased, and the autoregressive component in Phase II. In Section 5.1, we have described the reasoning behind these two new options, as well as for setting a lower bound on the standard deviation.
You should expect a full simulation with these settings to run for several days. Thus, we recommend processing it by a batch script in the background, so that it can be left unattended.

Discussion
In this article, we have described the latest major update of the R package bayesTFR. This update significantly enriches the modeling framework in the previous version of the package, and gives analysts the flexibility to account for past TFR uncertainty, use annual data, and allow for an autoregressive model in Phase II. Moreover, by making use of the vectorization nature of R, the computational cost has been kept at a reasonable level while making the model more sophisticated. New functions for visualizing estimation results, as well as updated analysis tools will further support analysts in exploring the package outputs.
On the package development side, there are at least two major areas for future improvements. The first is modeling age-specific fertility rates with past uncertainty which is of interest to demographers. The second would be further vectorizing the MCMC process. If past uncertainty is included in the model, updating the estimates of TFR is the most time-consuming part of the process. Since we consider each past TFR per country and time period as a parameter, it adds over 14,000 parameters in the annual case. Thus, the speed of the Metropolis-Hastings step for updating TFR plays a big role in determining the overall speed of the method. If past uncertainty is not included, updates of country-specific parameters dominate the computing time, and thus are subject to further optimization.
On the modeling side, there are also two obvious directions for improvement. First, instead of modeling the bias and standard deviation based on linear regression for each country separately, these could be folded into the process, giving a fully united probabilistic model. A pooled version could yield more robust estimates, especially given the small amount of data in some surveys. Another direction is related to the completeness of the VR data. The completeness of VR coverage is the most important factor for how precise the VR records are, and this is an important consideration for VR but not for other surveys. Due to the low bias of high quality vital registration systems, more research could be done on how to incorporate