Multilevel time series modelling of antenatal care coverage in Bangladesh at disaggregated administrative levels
Section 8. Discussion
In this study, multilevel time-series (MTS) models have been developed for the percentage of women receiving no antenatal consult (ANC0) and the percentage of women receiving at least 4 antenatal consults (ANC4) in Bangladesh, using only seven editions of the Bangladesh Demographic and Health Survey (BDHS) over the period of 1994-2014. Time series models are defined at an annual frequency where years without a survey edition are treated as missing. In this way, the model accounts for the varying time gaps between the subsequent editions of the BDHS and produce predictions in the years without sample surveys. Trends are produced at three regional levels, namely the national level, a break down in 7 divisions and a breakdown in 64 districts.
In the first model (MTS-I) year-domain-specific direct estimates and their standard errors are used as input in the MTS model. Trends obtained under this model, are rather volatile since the trend estimates tend to follow the direct estimates. Another drawback of the MTS-I model is that it hampers the use of auxiliary information from two available censuses, since values for the auxiliary information from a particular census does not change in two or three subsequent editions of the survey. To use this census information, it is proposed to develop cross-sectional Fay-Herriot (FH) models for each survey year separately. In a second MTS model, MTS-II, these FH estimates and their standard errors are used as input series. In a third model, called MTS-III, the FH estimates with their full covariance matrices, are used as input series. This MTS model properly accounts for the cross-sectional correlations between the FH estimates. The overall model for MTS-II and MTS-III is then a two-step non-iterated process, for which the first stage is producing the FH estimates.
The models are developed at the most detailed regional level of districts. Division and national level trends are estimated by aggregating predictions of the district level trends. In this way, figures at different aggregation levels are numerically consistent by definition.
Compared to other time series small area estimation models proposed in the literature, our models contain more structure, since dynamic trend models are specified at different aggregation levels. This is necessary to obtain accurate aggregated predictions for the divisions and the national level and is a more parsimonious way of modelling cross-sectional correlations. Further model regularization was considered by specifying global-local priors. This, however, did not further improve the model fits.
In small area estimation, domain estimates are often benchmarked to the direct estimates at the national level for numerical consistency and as an attempt to reduce the bias in the model based domain predictions. In this application the trend estimates at the national level under the MTS models are already very close to the direct estimates. Therefore we do not consider an additional benchmark step.
All three time series models provide estimates with improved accuracy. Because MTS-II ignores the predominantly positive correlations between the cross-sectional FH input series, the standard errors of the trends at aggregated levels are actually too small. Since MTS-III accounts for these correlations, the standard errors for national and division trends are larger but also more realistic. The MTS-II model, however, seems to provide most plausible trends for both response variables, particularly at the district level, by compromising volatility in the trends under the MTS-I model and flatness in the trends under the MTS-III model. This choice is supported by the fact that these variables are likely to be relatively smooth over time. Fitting these models to the series of ANC0 and ANC4 is therefore certainly suitable in concept. This also justifies the interpolation of the trends for the years without sample surveys. Our approach can be useful also for many developing countries with repeated DHS surveys, since these are typically observed with varying time lags and mainly depend on census information that is not updated within two or three subsequent editions of the survey.
Using predictions of the cross-sectional FH models as input series for the MTS models, is proposed as a practical solution to make better use of the available census information. The additional advantage of this approach is that it stabilizes the input series for the MTS models by removing large sampling errors from the direct estimates. This requires, however, a careful model selection and evaluation process for the cross-sectional FH models, since model miss-specification of the cross-sectional FH models can result in biased input series with estimated standard errors that underestimate the real uncertainty.
One limitation of this study is related to the bias correction for the square root transformation that is applied to ANC4. The bias correction can only be applied to the trend estimates in the survey years. This results in awkward increases of point estimates if the sampling error is not smoothed enough, particularly for the domains with small sample size. This hampers estimation of period-to-period changes between survey years and non-survey years. Therefore the bias correction is only applied to the cross-sectional FH models and not to the MTS models.
The prevalence of ANC0 and ANC4 visits are negatively correlated, so a multivariate model may be an interesting alternative to the univariate models used here. The two series could be combined with the series of the remainder category in a single multivariate model. This, however, requires a multinomial model that has the advantages that it may further improve the precision of the estimates and guarantees that the predictions take values in their admissible range, and that the predictions over the categories add up to hundred percent. The multinomial model is, however, not easy to implement. Particularly in this study the variance-covariance matrix can be difficult to estimate for the districts with small number of observation. Furthermore, Datta et al. (2002) shows that univariate models may provide as good results as multivariate models proposed in Ghosh, Nangia and Kim (1996), while being simpler to implement. The extension of our univariate models to a multinomial model is therefore left for further research.
For ANC0, the national level shows a downward trend. The decline in the trend temporarily stopped during 2004-2011. The trend of ANC4 shows steady increase over the considered study period. Division level trends for ANC0 show a steady decline for all the divisions except Dhaka, Chittagong and Sylhet divisions. The trends for these three divisions remained stable during the period of 2004-2011 which mainly causes the flat trend at the national level of ANC0. On the other hand, at the division level ANC4 shows almost linear upward trends for most of the divisions except Dhaka and Chittagong. The greatest improvement is observed for Khulna and Rangpur divisions where the trends of ANC4 reach to more than 40% in 2014. District-level trends help to identify highly vulnerable districts in terms of the two considered response variables. Though the national level trend of ANC0 declines to about 21% in 2014, a few districts get below 10% (Dhaka, Jhenaidaha, and Meherpur) while a considerable number of districts still have ANC0 higher than 35% (Bhola, cox’s Bazar, Kishoregonj, Noakhali, Sunamganj, Sirajgonj, and three Chittagong hill tract districts). For ANC4, a few districts have estimates above 50% (Dhaka, Nilphamari, and Panchagarh) and most of the districts with high ANC0 have ANC4 estimates less than 20%. These district level trends might help policy makers to focus on vulnerable hotspots where both ANC0 and ANC4 indicators are still poor. Obviously, detailed level trends might help policy makers to take actions for reducing disaggregated level inequalities in the race to achieve SDGs.
Acknowledgements
We wish to thank Measure Evaluation and National Institute of Population Research and Training (NIPORT) for making the BDHS data publicly available. In addition, IPUMS deserves thanks for providing the access to the sample data of Bangaldesh Census 1991, Census 2001, and Census 2011. The views expressed in this paper are those of the authors and do not necessarily reflect the policy of Statistics Netherlands. The authors are grateful to two anonymous reviewers and the Associate Editor for providing useful comments on a former draft of the paper.
Appendix
| Variable | Definition |
|---|---|
| Division | Barishal, Chittagong, Dhaka, Khulna, Rajshahi, Rangpur, Sylhet. |
| Region | (1) Densely populated Dhaka, Chittagong and Gazipur districts, |
| (2) 9 regional districts with big cities, | |
| (3) 3 hilly districts (Bandarban, Khagrachhari and Rangamati), | |
| (4) 49 other districts (less urbanized areas). | |
| Chittagong | Chittagong Division? |
| Dhaka | Dhaka Division? |
| Khulna | Khulna Division? |
| Rangpur | Rangpur Division? |
| Rajshahi | Rajshahi Division? |
| Proportion of Under-5 children. | |
| Proportion of women aged 15-49 years. | |
| Proportion of married women aged 15-49 years. | |
| Proportion of married women aged 15-49 years having primary education. | |
| Proportion of married women aged 15-49 years having at least secondary education. | |
| Proportion of household (HH) with illiterate women aged 15-49 years. | |
| Proportion of household (HH) with primary educated women aged 15-49 years. | |
| Proportion of household (HH) with higher educated women aged 15-49 years. | |
| Proportion of HH with at least secondary educated HH head. | |
| Proportion of rural HH of size 4 and more. | |
| Proportion of rural HH with electricity. | |
| Proportion of rural HH with single mother. | |
| Proportion of HH having under-5 children and women aged -49 years having at least secondary education. | |
| Proportion of HH with 2 or more under-5 children. | |
| Proportion of rural HH with under-5 children. | |
| Proportion of rural HH with 2 or more under-5 children. | |
| Proportion of urban HH with 2 or more under-5 children. |
| Survey Year | Transformation | Fixed Effects | Random Effect | Census Data |
|---|---|---|---|---|
| 1994 | No | RI: District level Random Intercept | 1991 | |
| 1997 | No | RI | 1991 | |
| 2000 | No | RI | 1991 | |
| 2004 | No | RI | 2001 | |
| 2007 | No | RI | 2001 | |
| 2011 | SQRT | RI | 2011 | |
| 2014 | SQRT | RI | 2011 |
| Survey Year | Transformation | Fixed Effects | Random Effect | Census Data |
|---|---|---|---|---|
| 1994 | SQRT | RI: District level Random Intercept | 1991 | |
| 1997 | SQRT | RI | 1991 | |
| 2000 | SQRT | RI | 1991 | |
| 2004 | SQRT | RI | 2001 | |
| 2007 | SQRT | RI | 2001 | |
| 2011 | SQRT | RI | 2011 | |
| 2014 | SQRT | RI | 2011 |
References
BBS (2020). SDG Tracker: Bangladesh Development Mirror. https://www.sdg.gov.bd/page/thirty_nine_ plus_one_indicator/5#1, [Online; accessed 18-December-2020].
BBS and UNICEF (2019). Progotir Pathey Bangladesh: Bangladesh Multiple Indicator Cluster Survey 2019. Technical report, Bangladesh Bureau of Statistics (BBS).
Besag, J., and Kooperberg, C. (1995). On conditional and intrinsic autoregression. Biometrika, 82(4), 733-746.
Bollineni-Balabay, O., van den Brakel, J., Palm, F. and Boonstra, H.J. (2017). Multilevel hierarchical Bayesian versus state space approach in time series small area estimation: The Dutch Travel Survey. Journal of the Royal Statistical Society: Series A (Statistics in Society), 180(4), 1281-1308.
Boonstra, H.J. (2021). mcmcsae: Markov Chain Monte Carlo Small Area Estimation. R package version 0.6.0.
Boonstra, H.J., and van den Brakel, J. (2019). Estimation of level and change for unemployment using structural time series models. Survey Methodology, 45, 3, 395-425. Paper available at https://www150.statcan.gc.ca/n1/en/pub/12-001-x/2019003/article/00005-eng.pdf.
Boonstra, H.J., and van den Brakel, J. (2022). Multilevel time series models for small area estimation at different frequencies and domain levels. Annals of Applied Statistics, 16, 4, 2314-2338, DOI: 10.1214/21-AOAS1592.
Boonstra, H.J., van den Brakel, J. and Das, S. (2021). Multilevel time series modelling of mobility trends in the Netherlands for small domains. Journal of the Royal Statistical Society: Series A (Statistics in Society), 184(3), 985-1007.
Carvalho, C.M., Polson, N.G. and Scott, J.G. (2010). The horseshoe estimator for sparse signals. Biometrika, 97(2), 465-480.
Das, S., Kumar, B. and Kawsar, L.A. (2020). Disaggregated level child morbidity in Bangladesh: An application of small area estimation method. PLoS ONE, 15(5), e0220164.
Das, S., van den Brakel, J., Boonstra, H.J. and Haslett, S. (2021). Multilevel time series modelling of antenatal care coverage in Bangladesh at disaggregated administrative levels. Discussion paper, Statistics Netherlands.
Datta, G.S., Lahiri, P. and Maiti, T. (2002). Empirical Bayes estimation of median income of four-person families by state using time series and cross-sectional data. Journal of Statistical Planning and Inference, 102(1), 83-97.
Datta, G.S., Lahiri, P., Maiti, T. and Lu, K.L. (1999). Hierarchical Bayes estimation of unemployment rates for the States of the U.S. Journal of the American Statistical Association, 94(448), 1074-1082.
DHS (2021). The DHS Program: Demographic and Health Surveys. https://dhsprogram.com/, accessed 2021-11-01.
Fay, R.E., and Herriot, R.A. (1979). Estimates of income for small places: An application of James-Stein procedures to Census data. Journal of the American Statistical Association, 74(366), 269-277.
Gelfand, A.E., and Smith, A.F.M. (1990). Sampling based approaches to calculating marginal densities. Journal of the American Statistical Association, 85, 398-409.
Gelman, A. (2006). Prior distributions for variance parameters in hierarchical models. Bayesian Analysis, 1(3), 515-533.
Gelman, A., and Hill, J. (2007). Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press.
Gelman, A., and Rubin, D.B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7(4), 457-472.
Geman, S., and Geman, D. (1984). Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images. IEEE Trans. Pattn Anal. Mach. Intell., 6, 721-741.
Ghosh, M., Nangia, N. and Kim, D.H. (1996). Estimation of median income of four-person families: A Bayesian time series approach. Journal of the American Statistical Association, 91(436), 1423-1431.
Haslett, S., and Jones, G. (2004). Local Estimation of Poverty and Malnutrition in Bangladesh. Technical report, Bangladesh Bureau of Statistics and United Nations World Food Programme.
Haslett, S., Jones, G. and Isidro, M. (2014). Local Estimation of Poverty and Malnutrition in Bangladesh. Technical report, United Nations World Food Programme and International Fund for Agricultural Development.
Hossain, M.J., Das, S. and Chandra, H. (2020). Disaggregate level estimates and spatial mapping of food insecurity in Bangladesh by linking survey and census data. PLoS ONE, 15(4), e0230906.
Marhuenda, Y., Molina, I. and Morales, D. (2013). Small area estimation with spatio-temporal FayHerriot models. Computational Statistics & Data Analysis, 58, 308-325.
Mrisho, M., Obrist, B., Schellenberg, J.A., Haws, R.A., Mushi, A.K., Mshinda, H., Tanner, M. and Schellenberg, D. (2009). The use of antenatal and postnatal care: Perspectives and experiences of women and health care providers in rural southern Tanzania. BMC Pregnancy and Childbirth, 9(1), 10.
NIPORT (2015a). Bangladesh Demographic and Health Survey 2014: Key Indicators. Technical report, National Institute of Population Research and Training.
NIPORT (2015b). Bangladesh Demographic and Health Survey 2014: Policy Briefs. Technical report, National Institute of Population Research and Training.
NIPORT (2016). Bangladesh Demographic and Health Survey 2014. Technical report, National Institute of Population Research and Training.
O’Malley, A.J., and Zaslavsky, A.M. (2008). Domain-level covariance analysis for multilevel survey data with structured nonresponse. Journal of the American Statistical Association, 103(484), 1405-1418.
Park, T., and Casella, G. (2008). The Bayesian Lasso. Journal of the American Statistical Association, 103(482), 681-686.
Pfeffermann, D., and Burck, L. (1990). Robust small area estimation combining time series and cross-sectional data. Survey Methodology, 16, 2, 217-237. Paper available at https://www150.statcan.gc.ca/n1/en/pub/12-001-x/1990002/article/14534-eng.pdf.
Pfeffermann, D., and Tiller, R. (2006). Small area estimation with state-space models subject to benchmark constraints. Journal of the American Statistical Association, 101, 1387-1397.
Pratesi, M., and Salvati, N. (2008). Small area estimation: The EBLUP estimator based on spatially correlated random area effects. Statistical Methods and Applications, 17(1), 113-141.
R Core Team (2015). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
Rahman, M.S., Mohiuddin, H., Kafy, A.-A., Sheel, P.K. and Di, L. (2019). Classification of cities in Bangladesh based on remote sensing derived spatial characteristics. Journal of Urban Management, 8(2), 206-224.
Rao, J.N.K., and Molina, I. (2015). Small Area Estimation. Wiley-Interscience.
Rao, J.N.K, and Yu, M. (1994). Small area estimation by combining time series and cross-sectional data. The Canadian Journal of Statistics, 22, 511-528.
Rue, H., and Held, L. (2005). Gaussian Markov Random Fields: Theory and Applications. Chapman and Hall/CRC, 56.
Sakia, R.M. (1992). The Box-Cox transformation technique: A review. The Statistician, 41, 169-178.
Spiegelhalter, D.J., Best, N.G. Carlin, B.P. and van der Linde, A. (2002). Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society B, 64(4), 583-639.
Tibshirani, R. (1996). Regression shrinkage and selection via the Lasso. Journal of the Royal Statistical Society. Series B (Methodological), 267-288.
Watanabe, S. (2010). Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning Research, 11, 3571-3594.
Watanabe, S. (2013). A widely applicable Bayesian information criterion. Journal of Machine Learning Research, 14, 867-897.
West, M. (1984). Outlier models and prior distributions in bayesian linear regression. Journal of the Royal Statistical Society, Series B (Methodological), 431-439.
WHO, UNICEF and Others (2014). Trends in Maternal Mortality: 1990 to 2013: Estimates by WHO, UNICEF, UNFPA, the World Bank and the United Nations Population Division. Technical report, World Health Organization.
Wolter, K. (2007). Introduction to Variance Estimation. Springer.
You, Y. (2008). An integrated modeling approach to unemployment rate estimation for sub-provincial areas of Canada. Survey Methodology, 34, 1, 19-27. Paper available at https://www150.statcan.gc.ca/n1/en/pub/12-001-x/2008001/article/10614-eng.pdf.
You, Y., and Rao, J.N.K. (2000). Hierarchical Bayes estimation of small area means using multi-level models. Survey Methodology, 26, 2, 173-181. Paper available at https://www150.statcan.gc.ca/n1/en/pub/12-001-x/2000002/article/5537-eng.pdf.
You, Y., Rao, J.N.K. and Gambino, J. (2003). Model-based unemployment rate estimation for the Canadian Labour Force Survey: A hierarchical Bayes approach. Survey Methodology, 29, 1, 25-32. Paper available at https://www150.statcan.gc.ca/n1/en/pub/12-001-x/2003001/article/6602-eng.pdf.
- Date modified: