Small Area Estimation of Age-Specific and Total Fertility Rates in Bangladesh

Bangladesh has experienced a rapid national decline in fertility in recent decades, however, fertility rates vary considerably at the sub-national level (i.e., division). These variations are expected to be more pronounced at lower levels of geography (e.g., district level). However, routinely conducted demographic health surveys are designed for national estimates and do not have adequate samples to produce reliable estimate of fertility rates at lower levels of administrative units, particular when considering district level age-specific fertility rates. Data extracted from the Bangladesh Demographic Health Survey 2014 are used to derive direct estimates of age-specific fertility rates and associated smoothed standard errors. These are used as inputs for developing a small area model, which is expressed in a hierarchical Bayesian framework and fitted by Markov Chain Monte Carlo simulation. The model accounts for variation at different levels—women age-group, division, and district. The modeling results show large reductions in the estimated standards errors and provide consistent estimates of fertility at the detailed district age-specific level. There are significant differences in the fertility levels within and between districts and at the division level. Fertility rates are observed to be higher for Sylhet division and for women aged 20–24 years. We use geo-spatial maps of the fertility rates to visualize the variations over districts, and identify hot and cold-spots to have better targeted local level planning and policy decision making for further reductions in fertility rates in Bangladesh.


Introduction
Bangladesh has experienced a rapid national decline in fertility in past decades from a peak of more than six births per woman in mid 1970s to 2.3 births per woman in 2017-2018 (Bairagi & Datta, 2001;NIPORT et al., 2020). However, the declines in fertility rates do not adequately capture the geographical disparities in fertility and in fact, it is worthwhile to mention that fertility rates are unevenly distributed across the country. The 2017-2018 Bangladesh Demographic and Health Survey (BDHS) reveal that Sylhet (in the north-east) has the highest fertility rates (2.6 births per woman) following 2.5 births per woman in Mymensingh (north-central) and Chittagong (south-eastern coast) divisions (NIPORT et al., 2020). These disparities in divisional fertility rates indicate the inequality of fertility levels in Bangladesh. Moreover, these differences are expected to be more pronounced at lower levels of geography (e.g., district level). Thus, estimating age-specific fertility rates and total fertility rates at lower levels remain important in order to better understand the conditions, and attempts for strengthening family planning services and programmatic development to reach replacement-level fertility (2.1 births per woman) in Bangladesh. Earlier studies on the demographic transition (Bairagi & Datta, 2001;Caldwell et al., 1999;Cleland et al., 1994;Kabir et al., 2009;Khuda & Barkat, 2012) including recently published studies by Bongaarts and O'Neill (2018), and Bora et al. (2022) posited that the family planning program and the use of contraceptives were pivotal to the steep decline in fertility in Bangladesh.
Point estimates of fertility rates at the national and sub-national (i.e., division) levels are possible because the BDHS collects information on a wide range of indicators which is representative at these higher levels of geography. The BDHS, like most large nationally representative surveys, is designed to provide precise estimates at these (planned) areas. However, at lower levels of geography, for example the district, the survey data cannot be used to derive estimates of age-specific fertility as a result of small sample sizes. For policy planning purposes, it is often required to have data at much finer levels (e.g., by geographic area, ethnicity or other groups relevant to national context). These unplanned domains are referred to as 'small areas' because the sample sizes are too small to reliably estimate the parameter of interest (Das et al., 2019;Ghosh & Rao, 1994;Rao, 2003;Rao & Molina, 2015). To address this, survey methodologists derived small area estimation techniques to combine information from different sources using statistical modelling and survey sampling to increase the effective sample size, and thereby derive reliable estimates for these small (unplanned) areas (Ghosh & Rao, 1994;Pfeffermann, 2013). The basic concept of small area estimation is to borrow strength from related or similar small areas through statistical models that connect the small areas via supplementary data (e.g., census or administrative data). For this study, we use small area estimation to derive district level estimates of age-specific fertility and total fertility rates.
Small area estimation techniques are classified into two broad groups, based on whether the supplementary data used in the statistical modelling is available at the individual level (unit-level models) or aggregate level (area-level models). At an area-level, the Fay-Herriot model (Fay & Herriot, 1979) was developed to estimate per-capita income in municipalities (with populations less than 1000). In their model, they use auxiliary data such as municipal level per capita income, tax returns data and housing data to fit a linear model with random effects that improves the efficiency of the survey-based estimates. For this study, we adopted an area-level modelling approach since the auxiliary were available at the district level. There have been several authors, such as Prasad and Rao (1990), Rao and Yu (1994), Datta and Lahiri (2000), González-Manteiga et al. (2010), Benavent and Morales (2016), Pfeffermann (2002), without being exhaustive, have proposed extensions of the Fay-Herriot model in different applications of area-level small area estimation. In our paper, we fit the Fay-Heriot model in a hierarchical Bayesian framework using Markov Chain Monte Carlo (MCMC) simulations. The formulation of the Fay-Herriot model in the Bayesian framework has advantages in terms of its flexibility and computational efficiency. Furthermore, an R package "mcmcsae", developed by Boonstra (2020) is readily available to fit the small area models using MCMC simulations, and provide model assessment, selection and diagnostics. It has been used recently to estimate the mobility trends for small domains (i.e. cross-classification by sex, age purpose and mode) in the Netherlands (Boonstra et al., 2021). Das et al. (2022) also applied a similar methodology to estimate prevalence of child malnutrition for spatio-demographic domains (cross-classified by sex, age, place of residence and district) in Bangladesh. The advancement in SAE methodology in addition to freely available R package "mcmcsae", and the importance of disaggregate levels estimations on age-specific and total fertility rates are the motivation of this study. The results will help to identify the inequality and high fertility areas in Bangladesh.
As contextual background, there are three basic administrative layers of the country: seven divisions (Barisal, Chittagong, Dhaka, Khulna, Rajshahi, Rangpur, and Sylhet) as of 2014 BDHS, 64 districts and about 545 sub-districts. Subdistricts are further divided by rural-urban classification, such as wards (mohallas) in urban areas and union parishads (mouzas) in rural areas. Effective policy planning and service delivery on reproductive health and fertility requires information at the lowest level of administrative geography. Further, local authority administrators face challenges when implementing policy without relevant statistics. In this study, we derive robust model-based estimates of age-specific fertility rates and total fertility rates using data from round seven 2014 BDHS for the 64 administrative districts. This study builds on the existing demographic literature and contributes to a better understanding of the geographic variation in fertility in Bangladesh across districts.
The outline of the paper is as follows. Section 2 gives a brief description of the outcomes and data. Section 3 provides details about the small area estimation modelling, fixed effects, random effects and the estimation approach. Section 4 provides results, comprised of the parameter estimates (fixed and random effects) and estimations of TFRs and ASFRs, and maps examining the regional-level spatial patterns in fertility. Finally, the paper concludes with a discussion in Sect. 5. The background of this study, a number of tables and figures related to the described results, and a detailed description of the modelling approach and other related materials are provided as a separate Supplementary file.

Outcomes and Sample Survey Data
The main outcomes of interest are to estimate district level point-in-time estimates of age-specific fertility rates (ASFRs) and total fertility rates (TFRs). The data are obtained from a national representative survey 2014 BDHS sample. Birth history and an array of health and demographic information at household and individual levels of reproductive aged ever-married women collected through face-to-face interview are documented in the 'individual record file'. The birth history information is used to calculate the number of births by seven women age-groups, namely, 15-19, 20-24, 25-29, 30-34, 35-39, 40-44 and 45-49 years. The BDHS 2014 survey was based on two-stage stratified sample of households. The survey used as a sampling frame the list of enumeration areas prepared for the 2011 Population and Housing Census. In the first stage, 600 EAs were selected with probability proportional to the estimation area size, with 207 clusters in urban areas and 393 in rural areas. In the second stage of sampling, a systematic sample of 30 households on average was selected per estimation area to provide statistically reliable estimates for the key demographic and health variables for the country as a whole, for urban and rural areas separately, and for each of the seven divisions (for details see NIPORT et al., 2016). This two-stage sampling design clearly demonstrates that districts (second level administrative units) are unplanned domains for estimating any demographic indicators and so district level estimates of ASFRs (even for TFRs) obtained directly from the survey data are expected to unreliable due to small number of observations. Therefore, it justifies our aim to fit small area estimation models for the district level ASFRs in Bayesian framework by performing Markov Chain Monte Carlo simulations (Boonstra et al., 2021). As input of the small area estimation model, district level direct estimates (DIR) of ASFRs along with their associated standard errors and weighted number of married women are calculated using the package "DHS.rates" (Elkasabi, 2021) in the R environment (R Core Team, 2019) utilizing the microdata of reproductive aged women categorized into seven agegroups, described above and birth history information available in the 'individual records' file of BDHS, 2014. The computation accounts for the survey design and sampling weight, to adjust for selection, non-response and non-coverage biases in the direct estimates. The TFRs for each district are then obtained by summing the ASFRs, at five yearly intervals, and multiplying by five. The number of ever married women for each age-specific domain was extracted from the 5% census data of Bangladesh Population and Housing Census 2011 and this information was merged with the input file for the small area estimation model estimation. Since the BDHS is a sample survey, census data is used to ensure that the number of reproductive aged women (used as a denominator) in the computation of the fertility rates was numerically consistent (in a benchmarking step as part of the modelling process). The benchmarking ensured that the model-based ASFRs and TFRs summed up at a district, division and national level. These benchmarking adjustments correct for any inconsistencies between the model-based estimates at the various levels of aggregation.
It is worthwhile to note here that the direct estimates of the TFRs were found to be zero for three hilly districts (Bandarban, Khagrachari, and Rangamati) due to small sample size (less than 30), which also lead to zero values for the ASFRs. Furthermore, the older age-groups in some districts also reveal zero ASFRs (Please see the sparse distributions of the district level ASFRs with their standard errors and the domain-specific number of women sampled in the survey data in Supplementary Table S.2.1). However, the small area estimation model using the Markov Chain Monte Carlo (MCMC) simulations will help to calibrate the missing values. This is a recent development and an advantage of the use of MCMC simulations within the small area estimation modelling framework (Boonstra, 2020). Since the district location is not disclosed in the recent 2017-18 BDHS, this study used 2014 BDHS data.

Method
The multilevel models considered are defined at the most detailed level, which is the full cross-classification of women age-groups by district. The number of domains (or small areas) is therefore given as the multiplication of the 64 (districts) × 7 (ever married women age-groups) = 448 small areas. At a larger division level, there are 7 × 7 = 49 small areas. There are various considerations that are required before producing the final model-based small area estimates. For instance, the transformation of the model inputs, smoothing of the standard errors, and the inclusion of fixed effects as well as random effects.
In addition to fixed effects (such as district, division, or women age-group), we consider four types of random effects. These random effects that account for the unobserved differences in the fertility behaviour by (i) women age-groups in districts, (ii) women age-groups, or (iii) women age-groups in divisions. Further details are explained in Sect. 4.1. Models that lack to include potential explanatory variables (known as omitting variables problem) lead to such unobserved heterogeneity in the data. We consider another type of random effects component that accounts for the spatial variation due to neighbouring districts exhibiting patterned fertility behaviour. These correlated random effects allow smoothing over space and other similar small areas and are easily fitted using MCMC simulation. Using the direct estimates of ASFRs and associated standard errors, as modeling inputs, we aim to provide small-area model-based estimates of ASFRs with their precisions for all small areas separately for districts and divisions, respectively. We also compute TFRs for districts and divisions by following the underlying demographic aggregation method. We obtain robust standard errors using the repeated draws under the MCMC simulation. Finally, the different candidate models are evaluated by comparing the Widely Applicable Information Criterion (WAIC) (Watanabe, 2013) and the Deviance Information Criterion (DIC) (Spiegelhalter et al., 2002) to consider different competing models and determine which one is the best fit for the data. The best model is the one that explains the greatest amount of variation using the fewest possible variables, with some penalty for model complexity. Further, graphical comparison of predictions of the ASFRs and TFRs at various aggregation levels such as at the district, subnational (division), and national levels will be contrasted with the available direct estimates to visualize and investigate any discrepancies.
The model development also has different estimation procedures depending on whether district level age-specific fertility data are (a) observed with non-zero estimates (314 district-by-age domains), (b) observed with zero estimates (107 domains), and (c) non-sampled (27 remaining domains = 448-314-107). For the observed district-byage domains, the BDHS sample has fertility data from 3 to 714 ever married women, allowing the computation of the direct ASFRs. For the domains with zero ASFRs, it is not possible to estimate their standard errors, and so a smoothing approach has been implemented to estimate their standard errors. Finally, the model predicts the estimates for the remaining non-sampled areas.

Multilevel Small Area Estimation Model
The model is defined in Eq. (1) as a multilevel small area model, can also be an extension of the well-known Fay and Herriot area-level small area model (Fay & Herriot, 1979). The model is multilevel or hierarchical in the sense that the outcome values (ASFRs and associated standard errors) are nested at different levels for example, women age-groups and administrative regions (divisions or districts).
To define the model specifically for districts, let us to consider ŷ ad denote the ASFRs for age-groups a(1, 2, .., 7) over districts d(1, 2, .., 64) . The district level ASFRs are combined in one large vector as Ŷ = (ŷ 11 ,ŷ 12 , ..,ŷ 1d , … ,ŷ a1 ,ŷ a2 , … ,ŷ ad ) of size M = ad = 7 × 64 = 448 . To define the multilevel model for district level ASFRs, the vector Ŷ can be expressed as a general linear additive form of where X is a design matrix of size M × p , is a p-vector of fixed effects, Z is a design matrix of size M × q , is a q -vector of random effects, and is a M-vector of sampling errors. The sampling errors are assumed to follow normal distribution (with a mean zero vector and the standard deviation given by the standard errors of the estimated direct estimates). Also, a Student-t distribution instead of the normal distribution can be examined to give smaller weight to more outlying observations (West, 1984). The model (1) indicates that several independent vectors of random effects ( ) can be added to the linear predictor to obtain a better fitted multilevel model. Though vectors of random effects are assumed independent, the components within a vector can be correlated to accommodate cross-sectional and spatial correlations. The advantage of this multilevel small area estimation approach is the flexibility to include random effects and the distributional assumptions of the distribution of over the standard Fay-Herriot model. Each vector of random effects (say, ) is assumed to follow a normal distribution with a covariance matrix which is a Kronecker product of two covariance matrices A and D (say) of size a × a and d × d . The vector consists of ad random effects which can be thought of as women age-group level (a = 7) random effects varying over d = 64 districts. The covariance matrix A can be either unstructured or diagonal with equal or unequal variances. The covariance matrix D is assumed to be known and its precision matrix is used for computational efficiency (Rue & Held, 2005). Details of these random effects vectors and their covariance structure are referred to Boonstra and van den Brakel (2022); and Boonstra et al. (2021). As the model regression parameters (fixed effects), we have investigated the effects of women's age categories (Age) and administrative regions (Division) separately, and their interaction effects, denoted by where Age × Division indicates interaction effects of the two factor variables Age and Division . District level contextual variables (such as educational attainment, labour force status, or religion) calculated from the 5% sample of census data were not included to avoid the associated measurement errors (Ybarra & Lohr, 2008).
The random effects components are: (i) The first random effects component assumes that all the detailed level domains (district-by-age) have random intercepts which follow normal distribution with a common variance denoted by 2 1 . The model with these random effects component in addition to the fixed effects component is a standard Fay-Herriot model (Fay & Herriot, 1979). These random effects component accounts for unstructured variation at the most detailed level (i.e., due to unobserved differences in the fertility behavior of women). (ii) The second component assumes that random effects for the women age-groups vary over 64 districts with diagonal covariance structure with unequal variances. The variances for each agegroups are denoted by 2 2a , (a = 1, 2, .., 7 ). These random effects component accounts for cross-sectional correlation between the women age-groups and districts. (iii) We also include normally distributed random effects for domains cross-classified by women age-groups and higher levels of administrative region 'division' (divisionby-age) with a common variance ( 2 3 ). (iv) Finally, spatial strength has been brought by including district level spatial random effects, which are assumed to vary with equal variance denoted by 2 4 . The spatial random effects are assumed to follow an intrinsic conditional autoregressive (ICAR) model (Besag & Kooperberg, 1995). We estimate the model specified in Eq. (1) using the sqrt transformed smoothed standard errors obtained from generalized variance function (GVF) model (Wolter, 2007). This gives better results than applying the GVF models to the untransformed sampling error variances of the point estimates (Boonstra et al., 2021). For detail about transformation and GVF model we refer to Supplementary file section S.2. First, we run the model based on the assumption that the systematic idiosyncratic errors (sampling errors) are distributed normally. However, as the distribution of the ASFRs are observed to be skewed, we also run the model with the assumption that sampling errors follow the student-t distribution to provide less weight to the outlying values. The results of parameter estimates are given in Table 1 (see Sect. 4.1).

Simulations
Using the distributions of estimations obtain from the multilevel Bayesian models, we perform simulations of 10,000 MCMC draws. However, 2000 simulations 2 Page 8 of 24 of these 10,000 were randomly drawn to estimate our targeted means and standard errors of ASFR and TFR at different levels of aggregation. The simulation vectors of the linear predictions based on the fitted model are formed as where superscript (r) indexes the retained MCMC draws of estimated domain specific (over districts and age-groups) ASFRs. The means and standard deviations of the corresponding MCMC draws ̂ (r) are used to estimate ASFRs and associated standard errors, respectively at the district level. Further, the MCMC draws are aggregated to obtain the targeted district level TFRs and associated standard errors. Aggregation of MCMC draws of ASFRs and TFRs also provides national and division level ASFRs and TFRs, respectively. Since the linear predictions are on squareroot scale of ASFRs, the target estimates are obtained by back transformation of MCMC draws with bias correction using the smoothed standard errors following Boonstra et al. (2021) and Das et al. (2022). More detail regarding the model specification and estimation is provided in Section S.3 of the Supplementary file.

Results
The parameter estimates obtained from the final model (fixed and random effects) are explained in Sect. 4.1, the diagnostics of the model performance is discussed in Sect. 4.2. Finally, the model-based estimates of fertility rates are presented in Sect. 4.3, along with spatial maps identifying district and division level differences in fertility rates.

Parameter Estimates
The estimated regression coefficients and variance parameters of the considered four random effects components along with the DIC and WAIC values are given in Table 1. The number of degrees of freedom of the student-t distribution is estimated as approximately 7.5. This model is found to be better fitting in terms of WAIC, and DIC when compared to the model with normally distributed sampling errors. Thus, we select the model with student-t distributed sampling errors as our benchmark model. We examined how the inclusion of the last three random effects components improved the model performance based on the DIC and WAIC information criteria (please see Table S.3.1 in the Supplementary file). However, the parameter estimation of the alternative models that account for normally distributed sampling errors did not improve the model fit (Table S.3.2 in the Supplementary file). The estimated regression coefficients by women age-groups and divisions reveal higher levels of fertility in the age-group 20-24 and in the Sylhet division. Further, the estimated variance parameters indicate that fertility levels significantly vary over districts, but this is more apparent for the women in age-group 35-39 2 25 . However, as number of births of this age-group is relatively small, controlling for such large effects may not contribute much in improvement of the model estimation. The variance components corresponding to division-by-age level random effects ( 2 3 ) and district level spatial random effects ( 2 4 ) are found to be large and thus have significantly higher contribution in the improvement of model estimation. However, the random effects parameter by district and women age-groups ( 2 1 ) was found to be small.

Model Diagnostics and Performance
For assessing the validity and the reliability of the model-based small area estimates, we used a set of bias diagnostics. Model-based small area estimates should be (a) consistent with unbiased direct estimates and (b) more precise than direct estimates. As a bias diagnostic, the multilevel model-based estimates of division level ASFRs, district level ASFRs and TFRs (denoted by SAE) are plotted against the extreme than the DIR estimates. The fitted (OLS) regression (blue) line between the direct estimates and the model-based estimates indicates that they are highly correlated and linearly related. The middle-left panel plot also shows how the model-based estimator provides unbiased estimates for district level ASFRs. Additionally, it is also clear that the SAE estimator provides reasonable estimates for the domains with zero district level ASFRs. In particular, the SAE estimator provides smaller ASFRs for the older age-groups and higher ASFRs for the younger age-groups taking into account of the cross-sectional and spatial correlations. The estimated standard errors (SEs) of the SAE and DIR estimates for ASFRs are plotted in the right panel, which clearly shows that the SAE estimates are consistently estimated compared to the DIR even for the division level ASFRs. Similarly, the bias diagnostic plot (bottom-left panel) for district level TFRs indicates that the model-based estimates are approximately unbiased compared to the direct estimates in addition to reasonable TFR estimates for three districts with zero TFR. The bottom-right panel for the SEs also confirms that the model-based district level SAE estimates of TFRs are more consistent than the direct estimates. These diagnostic plots show that the model-based small area estimates derived from the fitted model should be consistent with the unbiased direct survey estimates. Furthermore, the model-based estimates should provide an approximation to the direct estimates that is consistent with these values being 'close' to the expected values of the direct estimates (shown by the proximity of estimates to the Y = X line). Also, the model-based small area estimates should have significantly lower variances than the corresponding direct survey estimates.
We also examined the performance of the modelling approach by considering whether the small area models were sensitive to the data inputs. We accomplished this by developing a simple Fay-Herriot model fitted to the district level estimates of TFRs, which used 64 small area domains (based on the 64 districts). Next, we developed the model using direct estimates of the district level ASFRs using 448 (= 64 × 7) small area domains. In both cases, the developed small area models, and derived SAE estimators, performed better than the direct estimates. While the model with the TFRs is less variable (due to having larger district level sample sizes), we see that there are bigger improvements at the district-by-age level when we examine the SAE and DIR estimates. We provided further detail on the TFR based model in the Supplementary section S.4. Figure 2 presents the error bars using the credible intervals (CIs) across women age-groups and divisions. This is to indicate the uncertainty in the estimates. We observe that for older age-groups, there is greater precision (depicted by narrower credible intervals) for the SAE estimators across all divisions except in Sylhet when compared to the DIR estimators. This is the distinct advantage of the use of advanced modeling through small area estimation. However, for the younger agegroups (i.e., women aged 15-19 and 20-24 years old), the SAE and DIR estimate shown to be close, with similar credible intervals. This is another property of the small area modeling that shows that where the sampled data are adequate then the direct and model-based estimates will be consistent. It is a confirmation of the law of large samples that an observed sample average from a large sample will be close to the true population average. Sylhet showed the highest levels of fertility across all Fig. 2 Age-specific fertility rates (ASFRs) at division level with 95% credible interval estimated by direct (DIR, black error-bar) and model-based (SAE, green error-bar) estimators age-groups and wide credible intervals; however, modeling results suggest precise estimates compared to the DIR estimates.
In Table 2 we investigated other measures of model diagnostics for example, relative bias (RB) (which measures the difference between the estimates SAE and DIR relative to DIR), absolute relative bias (ARB) (i.e., the absolute difference between SAE and DIR relative to DIR), percentage of relative reduction standard error (RRSE) (which measures the percentage reduction in the standard error of SAE compared to the standard error of DIR) and the ratio of coefficient of variation (CV). These model diagnostics are examined to show how the model-based SAE estimates outperform the DIR estimates (see Boonstra et al., 2021 for the mathematical definitions). These measures are helpful to establish the validity and reliability of the model estimators in terms of precision. Mainly, the smaller the values of the measures RB and ARB are, the better the SAE estimates are. However, larger the values of RRSE indicate the SAE estimates are more accurate than the direct estimates (Boonstra et al., 2021). Coverage of 95% credible interval of the SAE estimates is also examined by checking whether the direct estimates line within the estimated Table 2 Summary statistics of relative bias (RB, in %), absolute relative bias (ARB, in %), relative reduction standard errors (RRSE, in %), coefficient of variation ration (CV ratio), and coverage rate (CR, in %) of the model-based estimates based on multilevel model assuming sampling errors follow student-t distribution The performances are measured for age-specific fertility rates (ASFRs) and total fertility rates (TFRs) at district and division levels credible interval, see Fig. 2. Further, we investigated the aggregation properties of the model-based district level estimations at sub-national (division) and national levels (see Table 3). The summary statistics of the considered discrepancy measures indicate that model-based ASFRs and TFRs reveal consistent estimates in terms of lower mean values of RB and ARB at the division level and higher mean values of RRSE at the district level. The summary of CV ratios also shows that the model-based estimator was more efficient at the disaggregated level. The coverage rates of the small area model-based estimates are higher in the case of the TFRs compared to the ASFRs as the DIR estimates of TFRs are estimated consistently due to the larger sample size. By contrast, the estimates of the ASFRs are more variable as observations of some age-groups were very small or even zero. However, overall, the summary statistics indicate that the small area modeling estimates are found to be consistent and unbiased.

Divisional and National Level ASFRs and TFRs
The national and division level estimates shown in Table 3 are calculated from the aggregation of district level estimates of ASFRs and TFRs. The SAE estimates of both TFRs and ASFRs at division and national levels are observed consistent and are very close to the DIR estimates. The national level TFR is 2.3, and this tells us about the average number of children born to a woman in Bangladesh (this is the same for both the SAE and DIR, as a result of the consistency property at aggregate national level). The lowest TFR is estimated for Rangpur division (DIR: 1.81 vs SAE: 1.90) following Khulna (DIR: 1.87 vs SAE: 1.92), Rajshahi (DIR: 2.02 vs SAE: 2.05), Barisal (DIR: 2.13 vs SAE: 2.12), Dhaka (DIR: 2.31 vs SAE: 2.23), Table 3 Age-specific fertility rates (ASFRs) and total fertility rates (TFRs) at division and national levels estimated by direct (DIR) and model-based (SAE) estimators Estimates and standard errors are denoted by "est" and "se" respectively  Fig. 1) particularly for the older age-groups (See Supplementary Table S.5.1). Detailed results on the division level ASFRs are given in Supplementary Section S.5.

District Level ASFRs and TFRs
The spatial distribution of district level TFRs estimated from SAE estimator is mapped in Fig. 3 Fig. 4 as map as well as given in supplementary Tables S.6.1-S.6.7, and these show remarkable variability in fertility for different women ages. Highest fertility rates are estimated for women in age-group 20-24 in the majority of the districts, with a few exceptions where the highest rates are estimated for the age-group 15-19 (in districts Dinajpur, Gazipur, Munshiganj, Narayanganj, Panchagarh, and Shariatpur) (see Fig. 4). District level ASFRs are displayed by division in tabular format as well as in map for each age-group in Supplementary Section S.6.

Discussion and Conclusion
Deriving detailed fertility levels such as age-specific fertility rates and total fertility rates are particularly important for Bangladesh. This is due to the fact that the national government lacks the resources or the civil registration system to create a complete census or collect data of such detailed local outcomes. While there are nationally representative surveys like the DHS that are designed to provide reliable estimates at national and sub-national (division) levels (NIPORT et al., 2016), due to small sample size, the estimates are not reliable at disaggregated (e.g., district) levels. Moreover, the observed samples in some districts are found to be zero. Furthermore, when considering age-specific fertility there are a substantial number of districts with no women in particular age-groups, making it difficult to properly estimate the ASFRs using only information from the BDHS sample. This study has attempted to generate district level estimates of ASFRs and TFRs by applying Markov Chain Monte Carlo simulations within the hierarchical Bayesian framework. The data used in this study was sourced from 2014 BDHS. The model that has been developed helps to simulate the values (i.e., ASFRs and standard errors) based on the direct estimates of ASFRs and associated standard errors as input of the model. This small area model achieves higher accuracy by accounting for the unobserved heterogeneity and spatial heterogeneity through random effects. This process is done by means of cross-classification of women age-groups in 64 districts and 7 divisions. This analysis of fertility behaviour at District level map of age-specific fertility rates (ASFRs) in Bangladesh district level allows for better inference-making and policy planning, which in turn illustrates the model's applicability and its benefits. Reliable estimates of fertility at the district level can highlight the sub-regional differences that are often overlooked at higher levels of geographical aggregation, such as at national or divisional levels. For assessing validity and reliability of the model-based small area estimations, the small area estimation literature suggest two types of assessments (for example, see Brown et al., 2001;Chandra et al., 2011;Baffour et al., 2019). Firstly, the model-based estimates should be consistent with the unbiased direct survey estimates; they should provide an approximation to the direct estimates by giving values which are close to the expected values of the direct estimates. Secondly, the estimates will be more precise than the direct estimates as the aim of the model-based estimates is to reduce the mean squared errors.
We found that there was significant improvement in the model-based district level estimates especially in the domains with sparse samples. Furthermore, the developed model captured the variation at women age-groups, district and division levels. The lack of significance in the variance estimates by women age-group and fertility levels at smaller areas (i.e., districts) ( 2 1 ) as can be observed in Table 1, is indicative of where possible improvements in the modelling can be made for future SAE analysis that can better account for variations at smaller areas. Investigating Fig. 2, it is observed that the variation in ASFRs was highest for Sylhet division, among older age-groups. However, the model-based estimates were found to be better with smaller standard errors, particularly for the older age-groups when contrasted with the direct estimates, indicative of the outperformance of the direct estimation by the small area modeling. We also investigated the aggregation properties of the model-based district level estimates at higher (division and national) level and found that the model-based estimates were very close to the direct estimates of the divisions for both fertility outcomes, i.e., ASFRs and TFRs. By examining the CVs (Table 2), it was also observed that the model-based estimator of ASFRs and TFRs outperformed the DIR estimator both at the district and division level. Again, the diagnostic measures explained in Fig. 1 confirm that the model-based estimates are robust enough to provide reliable estimates of the district level total and age-specific fertility rates. However, notably it is observed that model estimations of ASFRs better performed compared to the model estimation of TFRs in terms of precision. This implies that ASFRs model captured the differences in fertility levels by women age-groups, which is important. For further detail, we refer Supplementary section S.4.
We use geo-spatial maps of the fertility rates to visualize the variations over districts. The areas with high fertility (hot-spots), or low fertility (cold-spots) are indicative of localized behavioural and cultural factors, such as poor contraceptive use or early child marriage (Caldwell et al., 1999). This can aid policy makers and program managers in formulating appropriate policies at the local level. Our results given in Table 1 suggest that division level estimates of TFRs in Sylhet and Chittagong divisions have the highest TFR followed by other division like Dhaka, Khulna, Rajshahi. These fertility rates are analyzed by women age-groups (ASFRs). The results show higher fertility rates among the younger women age-groups (15-19, 20-24 years), and this pattern is similar for all the administrative hierarchies (Fig. 2). However, a considerable amount of differences in the ASFRs are observed between the eastern and western parts of the country. The spatial distribution of TFRs and ASFRs (Figs. 3, 4) at district level confirms this difference between the eastern and western districts. Higher level of TFRs and ASFRs are observed in the north-eastern districts, while comparatively lower levels in the north-western districts. The district level maps of ASFRs and TFRs further confirm this geo-spatial correlation in fertility. The north-western districts of Bangladesh share a border with the Bengali speaking districts of India. Earlier analysis highlights the importance of common language which helps to spread cross-cultural knowledge concerning development factors like contraception and family planning. This exchange of knowledge may explain the lower fertility rates in the north-western region of Bangladesh compared to the north-eastern region (Amin et al., 2002;Basu & Amin, 2000). Basu and Amin (2000) have stressed that cultural boundaries and social effects need to be taken into account when analyzing the causes for Bangladesh's fertility decline.
Research suggests that the high fertility regions (such as Syhlet and Chittagong in the east) have the higher levels of unmet need and method discontinuation, highest desired levels of fertility, lesser husband-wife agreement about the number of desired children compared to the lowest TFR regions Khulna and Rajshahi located in the western part of Bangladesh (Islam et al., 2003). The same study also indicates that availability and utilization of the family planning outreach services are poor in the highest TFR regions. This is evident from the coverage and visits to satellite clinics and a percentage of visits by field workers (Islam et al., 2003). However, there are observed inconsistencies between total fertility rates and contraceptive prevalence rates in Bangladesh (Saha & Bairagi, 2007). Such differences were explained by other factors like reductions in breastfeeding, son preference etc. Therefore, policies need to target multifactor interventions including implementing comprehensive and high-quality family planning programs targeted at the low performing areas, especially in the eastern regions of Bangladesh. Further, Saha et al. (2020) highlights the role of the husband to augment their knowledge in family planning and fertility reductions in Sylhet division. Adolescent girls often experience early childbearing and have poor utilization of reproductive, maternal, and neonatal health services. Also, there is an expansive amount of evidence that has showed that early marriage increases fertility rates (Nahar et al., 2013;Yaya et al., 2019;Kamal & Ulas, 2021). The recommended minimum of four skilled antenatal care visits was often not completed by adolescent girls who experienced unintended pregnancies (Khan et al., 2020). Unintended pregnancies occur among the women who want to avoid the pregnancy but do not use effective contraceptive method. These results can aid to address the unmet need among adolescents and youth by operationalizing the national adolescent health strategy.
On a final note, human geographers and demographers have studied the effects of climate change on health, migration, mortality and early child marriage in Bangladesh (Ahmed et al., 2019;Alston et al., 2014;Babalola et al., 2018;Carrico & Donato, 2019;Carrico et al., 2020;Chen et al., 2021). For example, Ahmed et al. (2019) revealed that two districts, Brahmanbaria and Sunamganj of Bangladesh, are more vulnerable to extreme weather events in terms of cyclone and flooding. This finding is in line with our findings of higher fertility estimates in Brahmanbaria (2.79) and Sunamganj (3.47). However, Carrico et al. (2020) revealed that the relationship between heat waves and marriage is strongest for women aged 18-23. Our estimates demonstrate high fertility rates among agegroup 20-24 years for most of the districts with an exception of six districts where fertility estimates were found high in the age-group 15-19 years. Again, Ahmed et al. (2019) suggests that early marriage occurs as a coping strategy in the aftermath of adverse weather events, which often lead to sexual violence against girls, subsequently degrading the family's reputation and the girl's future marriage prospects. Therefore, our results of high fertility regions highlight the importance of population policies to adopt strategies to cope with adverse climate weather that may help to reduce early child marriage and teenaged fertility. Accordingly, amongst other factors like increasing girls' education, the strategy should include strengthening social security to prevent sexual violence against girls and reduce early marriage in Bangladesh. Yet another set of research in Bangladesh suggests the pathway of impact of climate change on reproductive behavior and fertility (Chen et al., 2021) and importance to control for climate change parameters in the fertility analysis, in addition to other factors.
In conclusion, we see that our resulted model-based estimates can be used by the director of general health and family planning services (DGH&FPS) of the Bangladesh government to better understand the regional differences (at the lower levels of geography i.e., districts) in the fertility levels and teenaged fertility. These can in turn be used to support the Sustainable Development Goal (SDG) target 17.18 that emphasizes the need for disaggregated data for the local level planning and may help to guide family planning policy in Bangladesh to subsequently improve maternal and child health to meet the SDG goals 3 & 5. Nevertheless, the modeling technique aids the implication of small area estimation in the area of demographic research in a resource constrained country like Bangladesh. Further, we find evidence of the importance of accounting for variances by small domains like age-groups, districts, and divisions etc., which are common and important differences in many health outcomes. However, we notice that there are some limitations. In the main, our small area model fails to include district level covariates such as religious affiliation of the households, socio-economic development, family planning services, contraception use, educational attainment and autonomy of female, which are factors frequently associated with the fertility transition (Basu & Amin, 2000;Bongaarts & O'Neill, 2018;Gupta & Narayana, 1997). Johnson et al. (2012) estimated district level contraceptive use and unmet need for contraception using covariates in their small area model (using data from the census and DHS in Ghana). Their results revealed increasing rates of contraceptive use among women who had access to health care index. However, our model incorporated several random effects components to capture such unobserved differences in the data. As a result, the model fit was improved and obtained better precision in the fertility estimates for small areas compared to the direct estimates. However, again it remains an area for future research to include important area level covariates, for example, climate change and global warming variables in the fertility analysis. This will allow researchers to draw covariate-specific inference related to population policies and interventions on fertility decline.