Estimation of child undernutrition for spatiodemographic small domains in Bangladesh An application of multilevel Bayesian model


 Micro-level statistics on child undernutrition are highly prioritized by stakeholders for measuring and monitoring progress on the sustainable development goals. In this regard district-representative data were collected in the Bangladesh Multiple Indicator Cluster Survey 2019 for identifying localised disparities. However, district-level estimates of undernutrition indicators - stunting, wasting and underweight - remain largely unexplored. This study aims to estimate district-level prevalence of these indicators as well as to explore their disparities at sub-national (division) and district level spatio-demographic domains cross-classified by children sex, age-groups, and place of residence. Bayesian multilevel models are developed at the sex-age-residence-district level, accounting for cross-sectional, spatial and spatio-demographic variations. The detailed domain-level predictions are aggregated to higher aggregation levels, which results in numerically consistent and reasonable estimates when compared to the design-based direct estimates. Spatio-demographic distributions of undernutrition indicators indicate south-western districts have lower vulnerability to undernutrition than north-eastern districts, and indicate significant inequalities within and between administrative hierarchies, attributable to child age and place of residence. These disparities in undernutrition at both aggregated and disaggregated spatio-demographic domains can aid policymakers in the social inclusion of the most vulnerable to meet the sustainable development goals by 2030.


Introduction
Malnutrition remains a critical public health problem among children under the age of five years worldwide. Malnutrition is caused by multifarious interlinked factors and has both short and long term detrimental effects (Pelletier & Frongillo, 2003;Pelletier, Frongillo Jr, Schroeder, & Habicht, 1995). It affects the cognitive and physical development of children, increases the risk of infections and significantly contributes adversely to mortality and morbidity (Rice, Sacco, Hyder, & Black, 2000;WHO, 2020). There are three internationally recommended indicators to measure malnutrition and monitor child growth: stunting (low height-for-age), underweight (low weight-for-age), and wasting (low weight-for-height) (De Onis & Blössner, 2003;De Onis et al., 2019;WHO, 2013). While stunting reflects a failure to reach linear growth trajectory, and is attributable to suboptimal health and/or nutritional conditions, being underweight is caused by a child's low body mass relative to their chronological age and as such is influenced by both height and weight. By definition, underweight cannot distinguish between a child who has low weight relative to his/her height or a child who is short for in stature relative his/her age, but may be normal in weight-for-height. Finally, wasting in children is an indicator of acute undernutrition and is a result of more food deprivation or illness. Therefore, on the one hand, stunting and wasting indicate chronic and acute malnutrition respectively. On the other hand, underweight is a composite indicator and includes both acute (wasting) and chronic (stunting) malnutrition.
Malnutrition significantly contributes to the global burden of several diseases (Black et al., 2013;Victora et al., 2008). From the most recent global figures, undernutrition accounts for at least half of all the deaths annually in children under five (Lancet, 2018). The vast majority of these deaths occur in low-and middle-income countries, particularly in Africa and Asia. In the Asian region, Bangladesh has one of highest prevalence of malnutrition over the last two-decades. The most recent Demographic and Health Survey (BDHS) 2017-18 and Multiple Indicator Cluster Survey (MICS) 2019 findings show that near about one-third of preschool-age children are stunted, more than one-fifth are underweight and around one-tenth are wasted (BBS & UNICEF, 2019;NIPORT and ICF, 2020). Due to the current COVID-19 pandemic, the prevalence of wasting is expected to increase to 1.9 million children in 2020 from 1.7 million in 2019 (UNICEF, 2020). Consequently, the impact of this epidemic might obstruct the progress to the Sustainable Development Goal (SDG) target 2.2 of reducing the wasting level below 5% by 2025 (BBS & UNICEF, 2019) as well as the stunting level below 25% by 2022 (NIPORT and ICF, 2020). Though the stunting level has declined considerably over the last two-decades (from 51% in 2004 to 31% in 2017) in Bangladesh (NIPORT and ICF, 2020), the average annual rate of reduction is still below the global recommended rate of 3.9% for stunting (De Onis et al., 2013).
The achievement in the reduction of child undernutrition significantly vary at the sub-national level (Mohsena, Goto, & Mascie-Taylor, 2015) as well as at various socio-economic and demographic levels (Das, Hossain, & Nesa, 2009;Das & Rahman, 2011;Islam et al., 2020). Over the period of 1997-2017, the stunting level remained consistently higher in Sylhet division (above 40%) and lower in the Khulna division (around 25% in 2017) (NIPORT and ICF, 2020). The vulnerability of child stunting and underweight still remain higher among the rural children (NIPORT and ICF, 2020). There is still a huge gap in the undernutrition vulnerability between infants (about 20% in 2017) and other older age-groups (for example, more than 30% for 12-23 age-group). Though undernutrition levels do not significantly vary by sex at the national level, female children are more prone to undernutrition (Haq et al., 2021;Henry et al., 1993). A considerable amount of spatial variation in the prevalence of undernutrition indicators at the lower administrative units such as districts and sub-districts has also been observed in several studies (Haslett & Jones, 2004;Nesa, 2019). However, a proper understanding of the spatio-demographic variation (say, the interaction of administrative units, place of residence, children sex and age) in the prevalence of undernutrition has remained unknown. An expansive range of literature has been devoted to understanding malnutrition and the socio-demographic and spatial determinants through multilevel modelling (for example, Alom, Quddus, and Islam (2012); Amegbor, Zhang, Dalgaard, and Sabel (2020); M.R.K. Chowdhury et al. (2016); Sultana, Rahman, and Akter (2019); Yang et al. (2018)).
To eradicate all forms of malnutrition in line with the SDG targets, the inequalities in the undernutrition vulnerability at the spatial, demographic and spatio-demographic level should be reduced as per the SDG 10 goal of reducing inequalities within country. Proper allocation of resources and policy-making by focusing on the disaggregated administrative units can only help to reduce this inequality. In this regards, reliable statistical information at detailed administrative units as well as their cross-classified demographic domains are required. Since nationally representative surveys fail to provide such precise disaggregated level statistics, small area estimation (Rao & Molina, 2015) has been widely used in poverty estimation (Elbers, Lanjouw, & Lanjouw, 2003;Molina & Rao, 2010;Tzavidis, Salvati, Pratesi, & Chambers, 2008), labour force (Chambers, Salvati, & Tzavidis, 2016), food insecurity (Haslett, Jones, Isidro, & Sefton, 2014;Hossain, Das, Chandra, & Islam, 2020), and child undernutrition (Haslett & Jones, 2004;Haslett, Jones, & Sefton, 2013). In this work, we use a multilevel framework to incorporate contextual factors and increase the precision of estimates under small area estimation for enhanced inference of malnutrition indicators.
The recent MICS conducted in 2019 covered all 64 districts in the survey design for identifying localised disparities as well as the most vulnerable for social inclusion (BBS & UNICEF, 2019). However, district level estimates of three undernutrition indicators have remained largely unexplored and widely unpublished. Availability of the large dataset has given scope to investigate the disparities in the prevalence of children undernutrition at district level as well as at the corresponding spatio-demographic detailed level domains. So this study aims to estimate the prevalence of stunting, wasting and underweight for 1280 small domains which are cross-classified by children's sex (female and male), five age-groups (0-11, 12-23, 24-35, 36-47, and 48-59 months), place of residence (rural and urban) and 64 districts by applying an extension of small area estimation method using multilevel modelling. This approach has explored spatio-demographic variation at both division and district levels (first and second administrative units respectively). These disaggregated level prevalence of children undernutrition indicators as well as the observed spatio-demographic inequalities might help policy makers to understand the spatial disparity in the undernutrition vulnerability and how the vulnerability is distributed unevenly throughout Bangladesh. To the best of our knowledge, such investigation on the dynamics of spatio-demographic disparities in child undernutriton has not been yet done in Bangladesh. So we expect this understanding will provide immense benefit for policy-making in the race of meeting SDG goals by 2030.
In this study, we use a multilevel modelling extension of small area estimation method to estimate stunting, wasting and underweight for spatio-demographic domains of division and districts using MICS 2019 and Bangladesh Population and Housing census 2011 data. The considered multilevel models expressed within a hierarchical Bayesian framework are developed at the detailed level cross-classified domains of sixty-four districts. These multilevel models are known to reduce uncertainty by including variability in the outcome variables at various aggregation levels and hence provide consistent estimates for the domains with very small samples in the survey data (Sugasawa & Kubokawa, 2020). This is accomplished through borrowing strength by accounting for cross-sectional correlation from similar domains, and spatial correlation from the neighbouring domains (Boonstra & van den Brakel, 2018;Boonstra, van den Brakel, & Das, 2021). Spatio-demographic district level estimates are obtained directly from the developed multilevel models and then these disaggregated estimates are aggregated to obtain division and national level estimates, which provides numerically consistent estimates for the disaggregated level domains. Finally, the spatial distribution of children undernutrition indicators as well as their consistency are plotted in bivariate interactive maps.
The remainder of the paper is structured as follows: Section 2 describes preparation of data inputs for multilevel model from the MICS 2019 microdata and the Bayesian multilevel modelling techniques; Section 3 explains development of multilevel models and their assessments; Section 4 illustrates prevalence of three undernutrition indicators at various aggregation levels and district level spatial distribution through interactive bivariate maps; and finally the paper ends with discussion of the results including some policy implications in Section 5.

Data source and input estimates
The microdata of children nutrition status have been extracted from the Bangladesh MICS 2019. The sampling plan of MICS 2019 was designed to produce consistent estimates of population, health and nutrition indicators at division (first administrative geography) and district (second administrative geography) levels. In this regard, rural and urban parts of each district are considered as sampling strata. Then a two-stage stratified sampling procedure is followed: clusters, which are census enumeration areas in the 2011 Census of Bangladesh, are selected from each stratum systematically with probability proportional to size. Subsequently, 20 households were selected from each of the selected cluster. Finally, a total of 64,400 households (in 3,220 clusters) were enumerated in the MICS 2019 sample.
Weights and heights for a total of 24,686 children under 5 years of age were measured in MICS 2019 to calculate their anthropometric indices for heightfor-age, weight-for-age and weight-for-height expressed in standard deviation units (z-scores) from the median of the reference population. According to MICS 2019 report, about 2.8%, 4.5%, and 4.7% children were excluded from calculations of the weight-for-age, height-for-age, and weight-for-height indicators, respectively due to either missing weight and height measurements or measurements were out of normal range (BBS & UNICEF, 2019). Finally, anthropometric information of 22,106, 22,484, and 22,063 children were utilized to calculate height-for-age z-scores, weight-for-age z-scores and weight-forheight z-score respectively. A child is considered as stunted, underweight and wasted when his/her z-scores are below -2.0 standard deviation units respectively (De Onis et al., 2019).
Though the MICS 2019 provides a large data set representative at the district level, the sample size of the targeted cross-classified domains are comparatively very small to estimate undernutrition indicators with enough precision. The sample size for the rural sampled domains varies over about 11-60 with a mean of about 28, while for the urban domains it ranges between 1-66 with a mean of 6.5. Notably, there are about 14 domains with no observation for all the three indicators. The comparison of domain-specific sample size with their respective population size obtained from Census 2011 (shown in Figure S.1 in Supplementary file) also indicates sample size for most of the urban domains are comparatively very small than those of rural domains having under five children below 10,000. Consequently, it is expected the estimates of child undernutrition will be comparatively highly inconsistent for the urban domains.
Standard errors of the estimated prevalence of stunting, wasting and underweight for each of the target domains are calculated by accounting for the sampling design and sampling weights through the Taylor series linearisation method (BBS & UNICEF, 2019). The standard errors (SEs) for the urban domains are comparatively very higher and noisy than those of rural domains due to very small sample as expected. Since these variance estimates are treated as fixed and known in the hierarchical Bayesian model as the Fay-Herriot area-level model (Fay & Herriot, 1979), they are smoothed using Generalized Variance Functions (GVF) method (Wolter (2007), Chapter 7) following the approach described in Boonstra et al. (2021). Figure 1 shows that the direct SEs are positively correlated with the estimates for all indicators and also the direct input estimates are very noisy for urban domains. As such, some transformations were examined to reduce the high dependencies among where X is a design matrix of D × p, β is a p-vector of fixed effects, Z (i) are D × q (i) design matrices for q (i) -dimensional random effect vectors v (i) , and e = (e 1 , . . . , e D ) ′ is the vector of sampling errors. The first part of the right hand side refers to fixed effects and the second part refers to inclusion of several random effects terms defined at different levels (e.g., district and spatio-demographic domains) simultaneously. The sampling errors are assumed normally distributed as e ∼ N (0, Σ), where the covariance matrix Σ can be a full or diagonal matrix of dimension D × D. In this study, Σ is considered as a diagonal matrix with the above mentioned smoothed error variances. Assuming normality of the sampling errors e, the likelihood function conditional to the fixed and random effects parameters can be defined as where η = Xβ + α Z (i) v (i) is the linear predictor. Each of the random effect vectors are formed based on two factor variables. For example, five levels of the Age variable are assumed to vary over 64 levels of the District variable. As such, the former factor variable (Age) with d = 5 levels and the later variable (District) with l = 64 levels constitute a random effect vector of length q = d × l = 5 × 64 with variance covariance matrix A ⊗ V where V and A are respectively d × d and l × l covariance matrices. The covariance matrix A is assumed known and its precision matrix Q A = A −1 (instead of A) is used because of computational efficiency (Rue & Held, 2005). The covariance matrix V for the d varying effects can be parameterized as a full covariance matrix (or diagonal matrix) with unequal (or equal) diagonal elements. The fixed effects parameters are assumed to follow a weakly informative prior distributions as p(β) = N (0, 100I). The scaled-inverse Wishart and half-Cauchy priors are used for the standard deviation parameters when 7 unstructured (Gelman & Hill, 2007) and diagonal (Gelman, 2006) covariance structure are assumed respectively.
A number of models structured as (1) were fitted using Markov Chain Monte Carlo simulation method for each of the undernutrition indicators in R (R Core Team, 2015) using package mcmcsae (Boonstra, 2018). A detail of our model formulation can also found in Boonstra and van den Brakel (2018) and Boonstra et al. (2021). Two model performance measures, namely the Widely Applicable Information Criterion or Watanabe-Akaike Information Criterion (WAIC) (Watanabe, 2010(Watanabe, , 2013 and the Deviance Information Criterion (DIC) (Spiegelhalter, Best, Carlin, & van der Linde, 2002) are used to compare models developed with same input estimates. The selected best multilevel model for each of the indicators was run finally with 2000 iterations so that longer simulations provide the Gelman-Rubin potential scale reduction factor below the recommended 1.10 value for all model parameters and model predictions at convergence (Gelman & Rubin, 1992).
3 Model development and model assessment

Model development
For each undernutrition indicator, initially the main effects of children demographic variables (Age, Sex, Residence) along with the regional variable Division are used as fixed effect components. As a random effect (RE) component (denoted by W N DIST ), domain specific random effects are assumed to follow the normal distribution with a common variance. This initial model with W N DIST component is a perfect example of Fay-Herriot model (Fay & Herriot, 1979). Since this model failed to account for variation by place of residence, interaction effects of Residence with Age, Sex and Division are examined as fixed effects, initially. However, it was found that inclusion of these interaction effects as random effects improved the model further in terms of the model selection criteria, DIC and WAIC. In this regard, a number of random effect components defined at various aggregation levels shown in Table 1 were included by specifying separate variance parameters for rural and urban domains. Three residence specific random effects components are defined at District (denoted by RES DIST ), District × Age (denoted by RES DIST AGE) and District × Sex (denoted by RES DIST SEX) level. The random effects under RES DIST component can be interpreted as random intercepts for the domains cross-classified by District and Residence following normal distributions with different variance components for rural and urban domains.
In a similar way, age-specific differentials at district level are accounted by an age-specific random effect component (denoted by AGE DIST ) with diagonal variance structure (See Table 1). The full covariance structure of V for the AGE DIST component has an additional 10 correlation parameters. However, inclusion of a spatial component (denoted by SP DIST ) in the model 8 Undernutrition prevalence in spatio-demographic domains to include spatial correlations among district through a conditional autoregressive (CAR) model (Besag & Kooperberg, 1995) reduced the importance of these additional correlation parameters of AGE DIST component. So, the diagonal variance structure is used for AGE DIST to borrow strength over space by children age-groups. The component SP DIST indicates spatial random effects for each of 64 districts conditional on the others, and are assumed to vary with equal variance. This spatial component mainly accounts for the spatial variation of undernutrition vulnerability over sub-regions. A similar approach was taken in the specification of the spatial random effect terms in Besag and Kooperberg (1995) and Rue and Held (2005). Finally, like the W N DIST component, a division level random effects component (denoted by W N DIV ) is added to specify random intercept for division level crossclassified domains of Age, Sex and Residence. The notation of the variance parameters related to the random effect components are shown in Table 2. After extensive examination, the main effects of the cross-classification and division variables are considered as fixed effects in the models of all three indicators. The posterior mean and the corresponding t-value for the fixed effects parameters are shown in Table 3. Table 3 shows that the cross-classified domains of children aged 12-23 months have higher stunting and underweight levels but lower wasting level compared to the cross-classified domains of infant children (0-11 months). While stunting and underweight levels are found highest among 24-35 months old children. Domains belong to the Sylhet division also experience higher levels of stunting and underweight. Though at the national level there was no significant difference by residence and sex, urban domains seem to experience lower stunting and underweight levels, while male children were observed to have higher wasting levels compared to females.
All the random effects specified in Table 1 had a significant contribution only in the selected underweight model (see Table 2). However, the variance parameters of RES DIST component in stunting model and the variance parameter of SP DIST in wasting model were not significant and so they are removed from the final models. According to the t-statistics, the SP DIST component seems to have higher contribution in the stunting model than the underweight model. The variance parameters of RES DIST AGE and RES DIST SEX components indicate that random effects of urban domains have higher variance than those of rural domains for all the three indicators. The AGE DIST component indicates random effects for the domains relating to infant (0-11) and toddler (12-23) children have comparatively higher variances. Variance estimates of W N DIV and W N DIST components reveal that variations at the disaggregated level domains are still significant for all the three indicators after accounting for the variations at the district level crossclassified domains, which are two-and three-order interactions of districts with children age, sex and residence.  Table 1 Summary of the random effect components for the multilevel models for stunting, wasting and underweight. The second and third columns refer to the varying effects with covariance matrix V , whereas the fourth column refers to the factor variable associated with A. The last two columns contain the total number of parameters and random effects for each random effects component. 0.026 2.327 ** 0.044 4.235 ** 0.036 3.208 * * p-value < 0.01, * * p-value < 0.05 Table 2 Posterior mean of random effect parameters and their t-values in the multilevel models of stunting, wasting and underweight.

Model assessment
Information criteria WAIC and DIC were utilized to find the best model among the models developed with same input estimates for each of the three undernutrition indicators. As model diagnostics, however, four discrepancy measures (i) relative bias (RB), (ii) absolute relative bias (ARB), (iii) relative reduction of the standard errors (RRSE), and the ratio of coefficients of variation (CV) denoted by CVR are calculated by comparing the model-based estimates denoted byθ d with the corresponding direct estimatesŶ d for domain d to evaluate how the multilevel models perform. Hereafter model-based and design-based estimates are denoted by SAE and DIR respectively. The RB and ARB defined by  Table 3 Posterior mean of fixed effect parameters and their t-values in the multilevel models of stunting, wasting and underweight.
at the higher aggregation level due to a greater chance of small (and zero) samples. The fourth measure CVR helps to examine the relative SEs of SAE estimates, their CVs are compared to those of DIR estimates. The values of CVR less than 100% indicates that the SAE estimates provide better relative accuracy than the DIR estimates. All these measures are expressed in percentage in Table 4. The measures are calculated at different (i.e., higher and lower) aggregation levels such as division, district, division×age, district×age, district×age×residence, and district×age×residence×sex.
For all the indicators, it is observed that, on one hand, the values of RB, ARB and RRSE increase with disaggregation levels, while the CV ratios decrease. These are expected since the direct estimates are more reliable and consistent at the higher aggregation level than the disaggregation level. On the other hand, the SAE estimates are as consistent as the DIR estimates at the higher aggregation level (there is little difference between the measures when compared to the direct estimates). More importantly, however, at the disaggregated level SAE estimates are more accurate than the DIR estimates. The SAE estimates appear to provide more gains for the cross-classified domains of district. The highest gains, in terms of RRSE and CVR, are observed at the most detailed level at which the multilevel models are developed. Figure S.2 in Supplementary file also supports that the developed multilevel models provide approximately unbiased and consistent estimates with reasonably better CVs at the most detailed level, particularly for the urban domains. This does make intuitive sense since it highlights the importance of the model estimation having more gains at the detailed disaggregated level.
At higher aggregation levels, the values of RB and ARB are found slightly higher for stunting and underweight compared to wasting, which indicates that the SAE estimates provide slightly higher estimates than the DIR estimates (only about 3% of the DIR estimates). The opposite pattern of ARB, at the detailed level, indicates that the SAE estimates are slightly overestimated for wasting than stunting and underweight. The purported reason for this is that the DIR estimates of wasting for many detailed level urban domains are too volatile shown in Figure S.2. Figure S.2 also shows how the developed multilevel models provide reasonable estimates for the domains with zero and one DIR estimates as well as zero sample size. For these domains, the huge uncertainty mean that the SAE estimates greatly improve upon the DIR estimates through using cross-sectional and spatial information from other domains.
At the detailed level, the multilevel model-based SAE estimator is expected to outperform the DIR estimator in terms of RRSE and CVR shown in Table 4. To show how the models provide reliable estimates by bringing strength from similar domains and spatial location, the SAE estimates of stunting (with 95% CI) for the detailed level domains under eight districts (Barguna, Bandarban, Dhaka, Sherpur, Khulna, Sirajganj, Dinajpur, and Sunamganj) as the representatives of eight divisions as well as some data characteristics such as domains with zero (and one) estimates, and the degree of undernutrition vulnerability. Model-based estimates with their 95% CI (coloured dots and error-bar lines) are plotted along with the DIR estimates (black circular dots) in Figure 2. All the considered districts except Dhaka have domains with zero estimates, where SAE estimators provide reasonable estimates. Some districts (Sherpur, Satkhira, Sirajganj, and Sunamganj) have one or more urban domains with the DIR estimates of one (sample size ranges 1-2 for most cases). The developed multilevel model improves on these inconsistent estimates and provided reasonable estimates by bringing strength from the relevant domains, particularly in the highly vulnerable age-groups 12-23 and 24-35 months. The SAE estimates follow the trend of comparatively more consistent DIR estimates, as for example most of the urban domains of Dhaka district. The performance of the developed multilevel models for wasting and underweight are very similar and so the similar comparisons for all the indicators are shown by division in Supplementary Figures from S.26 to S.51.

Results
Multilevel models are developed on square-root scale transformation, and so model predictions are back-transformed to original scale with a small bias correction following Boonstra et al. (2021) to obtain the prevalence of stunting, wasting and underweight along with their SEs. The prevalence of stunting, wasting and underweight at the national level are estimated as about 29%, 10% and 23% based on the developed multilevel models (Tables 5). These estimates are very similar to the direct estimates (28.0%, 9.8% and 22.6% respectively) documented in the Bangladesh MICS 2019 report (BBS & UNICEF, 2019).

12
Undernutrition prevalence in spatio-demographic domains  Table 4 Mean of relative bias (MRB, in %), mean of relative absolute bias (MARB, in %), mean of relative reduction of the standard errors (MRRSE, in %), and mean of CV ratios (CVR, in %) at different aggregation levels for the multilevel models of stunting, wasting and underweight.
The prevalence of these undernutrition indicators does not vary much by sex, however it varies considerably by children age-groups. The levels of wasting (10.6% in rural and 8.5% in urban) and underweight (24.7% in rural and 19.7% in urban) vary by place of residence compared to their magnitude at the national level. The levels of both stunting and underweight follow an increasing trend with the children age-groups up to 24-35 months (35.9% for stunting and 26.4% for underweight) and then reduce, but the reduction rate is higher for stunting than that of underweight. Though the children belonging to 24-35 months age-group are highly vulnerable to both stunting and underweight, they seem less vulnerable to wasting. Children belonging to 0-11, 12-23 and 48-59 months age-groups are more vulnerable to wasting. At division level, the prevalence of stunting and underweight are higher in Sylhet (38.5% and 33.1% respectively) and Mymensingh (34.4% and 26.3% respectively) divisions and lower in Khulna division (22.3% and 20.2% respectively). Though the prevalence of underweight (20.2%) and wasting (9.0%) are lower in Dhaka division, stunting is considerably higher (29.1%) than in Khulna division. The variations in the prevalence of stunting, wasting and underweight by children age, sex, and place of residence are expected to increase for their cross-classified domains with division and district. These spatio-demographic variations are explored in the following two sub-sections, along with the spatial distribution of child undernurition vulnerability at the district level through mapping the prevalence levels.

Spatio-demographic variation at division level
The prevalence levels of stunting, wasting and underweight for the crossclassified domains of division with children sex, age and place of residence are shown in Table 5. Like the national level, there is no sex difference in the division level prevalence of child undernutrition. However, rural-urban difference increases at the division level for all the indicators. The prevalence of stunting is found higher in rural parts of all divisions except that it is reversed in Dhaka division, with 31.1% of urban children are stunted compared to 28.2% of rural children. The rural-urban difference has been doubled even in Khulna division (23.0% vs 18.6%) compared to the difference at national level (29.7% vs 27.4% for stunting). The highest rural-urban difference is observed in Mymensingh division for both stunting (35.3% vs 28.8%) and underweight (27.2% vs 20.3%).
The increasing trends of stunting and underweight with children age-groups are also visible for all the divisions. The stunting and underweight levels have picked up to 45.2% and 38.2% respectively in Sylhet division followed by 42.4% and 30.4% in Barishal division among the 24-35 months old children ( Table  5). The concerning issue is that the highest level of stunting and underweight for all age-groups are both found in Sylhet division. However, infants (0 -11 months) of Khulna division have experienced lowest level of stunting (14.6%) and underweight (14.2%). For the most vulnerable 24-35 months children, the level of stunting and underweight increased to only 27.9% and 23.0% respectively in Khulna division. In the case of wasting, no age-specific pattern is obvious by division (Table 5). However, the wasting level in Sylhet and Barishal divisions are more than 12% for the younger two age-groups 0-11 and 12-23 months.
Undernutrition vulnerability does not vary considerably by sex and residence, however it varies considerably for the cross-classified domains of children age-group, sex and place of residence shown in Figure 3. Undernutrition levels are expected to be lower in urban domains, however for the infants, stunting level is found to be slightly higher for the urban domains than the rural domains. Though rural-urban difference is not obvious for stunting, the difference is significant for both wasting and underweight. For underweight, rural-urban differences have increased with the child age-groups and this difference is higher for the female children. The prevalence of stunting and underweight has increased until 24-35 months and then they decline for stunting but remain relatively stable for underweight with different slope by sex. As for example, the underweight levels sharply increased with children age-groups for female children compared to the male children. For wasting, there are no specific trends but the prevalence remains high for the 0-11, 12-23 and 48-59 months age-groups for both male and female children. Figure 3 also indicates that the confidence bands for urban domains are wider than the rural domains even at this higher aggregation levels. Figure 4 shows age-sex-residence level differences by division for all the undernutrition indicators. At a glance, Figure 4 shows considerable rural-urban differences at different degrees by children age-groups and division. For example, the difference is considerably higher for the female children belonging to 12-23 and 24-35 months age-groups of Rangpur division for all indicators. For underweight, the rural-urban differences seem higher among the female older age-groups (more specifically 24+ months age-groups) in Dhaka and Rajshahi divisions. More interestingly, this rural-urban difference in the prevalence of underweight increases with the age-groups (which pattern is also observed in Figure 3). This difference occurs due to stable estimates for rural domains compared to the declining trends in urban domains. In Sylhet and Barishal divisions, the prevalence of stunting increases steadily up to age-group 36-47 months in the urban domains, while the trends of rural domains follow the overall trend of increasing up to 24-35 months and then declining.

Spatio-demographic variation at district level
Since the MICS 2019 considered district in the sampling design, it is expected that the direct estimates of district level stunting, wasting and underweight are relatively consistent. However, the model-based estimates along with their accuracy for the cross-classified domains constructed by interaction of district, residence, children sex, and children age are expected to vary with the order of interaction. The distribution of estimated stunting, wasting and underweight based on the developed multilevel models are shown in Figure 5 (box-plots of top row) along their CVs (box-plots of bottom row) for districts (denoted by D) and their cross-classified domains with place of residence (R), children sex (S), and age-group (A). The spread of the estimates (which also represents inequality in the domains) increase with the order of interactions as expected, while spread of CVs (which represents relative accuracy of the estimates) increase exponentially. The distribution of CVs indicates that the estimates are more accurate for stunting and underweight than those for wasting. For stunting and underweight, the values of CV are below 20% for most of the domains constructed from second-order interactions of districts, while the range increase to 30% for third-order interactions and to 40% for the detailed level domains. Figure 5 shows that the values of CV are comparatively higher for wasting than those for stunting and underweight even at the district level. At the detailed level, the values of CV considerably vary over the ranges of 8-53%, 20-61% and 11-52% for stunting, wasting and underweight respectively (please see supplementary Table S.4). These differences among CVs are expected according to the mean level of stunting (29%), wasting (10%) and underweight (23%) as the denominator of CV.
District level prevalence of stunting, wasting and underweight are mapped in Figure 6, 7 and 8 as univariate and bivariate choropleth maps. The bivariate choropleth maps are created through the "Vizumap" R package (Lucchesi, Kuhnert, & Wikle, 2021), where the estimates and their accuracy in terms of SEs are shown simultaneously with single colour by encoding them in a bivariate grid of terciles. For example, the grey color of the first diagonal grid of Figure 6 indicates lower estimates of stunting with higher accuracy (grid of 1st tercile). The constructed bivariate choropleth maps help local decision makers to ensure that their decisions are made based on the statistical estimates with their uncertainty.
The maps in Figures 6 and 8    and underweight than the other regions of the country. For wasting, the spatial distribution shows that districts surrounding the capital city Dhaka have lower prevalence of wasting (Figure 7). A higher wasting level represented by yellowish colour is scattered around the country. For example, Maulvibazar in Sylhet division, Sherpur in Mymensingh division, Lakxmipur in Chittagong division and Pirojpur in Barishal division. Spatial maps of stunting and underweight indicates high positive correlation among the north-eastern (higher prevalence) and the south-western (lower prevalence) districts, however the correlation seems moderate in the districts situated in central and south-eastern regions where vulnerability due to underweight is less compared to that for stunting. These estimates are also more precise, under model-based small area estimation. The bivariate map in Figure 6 shows a higher prevalence of stunting (≥ 30%) with higher SEs (≥ 2%) in north-eastern region, hilly districts of southeastern region and in some southern districts. For underweight, there is higher prevalence of underweight (≥ 24.9%) with corresponding higher SEs (≥ 2%) in a number of districts belonging to the north, north-eastern and southern regions (Figure 8). For both stunting and underweight, lower prevalence (≤ 25.9% for stunting and ≤ 21.2% for underweight) with higher accuracy is observed in most of the districts of south-western districts (Khulna division). The capital city Dhaka has higher stunting but lower underweight levels with higher accuracy and the port city of Chittagong (now known as Chattogram) has lower stunting and underweight level with higher accuracy, since both the districts have considerably larger sample sizes in the survey data (being large populous regions in Bangladesh). In case of wasting, Figure 7 shows that many districts in the central region have lower prevalence (≤ 7%) with higher accuracy.
District level prevalence of children undernutrition indicators by place of residence are mapped in Figure 9. The maps explicitly reveal higher ruralurban difference in wasting and underweight compared to stunting and higher heterogeneity in the spatial distribution among the urban domains than the rural domains as expected according to the estimated random effects variance parameters shown in 2. For all the indicators, urban domains of south-western districts seem the least vulnerable groups. To investigate in detail, the estimated prevalence by place of residence are plotted in Figure 10 for Khulna, Dhaka, and Sunamganj districts, which belong to south-western, central and north-eastern regions and have different degrees of child undernutrition vulnerability. Figure 10 shows that prevalence of stunting is higher in the urban parts of Dhaka district (highly dense urban areas) though the prevalence is higher in most of the rural parts of Bangladesh. Interestingly, the precision for rural parts of Dhaka district is comparatively lower than that of urban parts only in Dhaka district. These patterns of estimate and precision are reversed in case of Khulna and Sunamganj districts, who have the lowest and highest prevalence of stunting respectively. In highly populated districts, the rural-urban difference in either stunting or underweight is found minimal. Some examples are the second highest populated district (Chittagong) and two other densely populated districts (Gazipur and Narayanganj) adjacent to Dhaka district (please see Figures S.4 Figure  10 also shows that stunting and underweight levels are highly underestimated by the direct estimator for Khulna district. Model-based estimators overcome this underestimation which, in turn, reduces the rural-urban differences. Ruralurban difference in the prevalence of child undernutrition for other districts are shown in supplementary figures S.3-S.9. Figure 11 shows that the trend of prevalence with the increase of children age remain same in the district level with significant differences by districts. Stunting level is found lowest in the oldest age-group instead of the infant agegroup in Dhaka district, While stunting level for the older age-groups (more than 40%) is significantly far way from the level of infants (about 25%) in Sunamganj district. The underweight level decreases with the age-groups in Dhaka while the levels in Khulna and Sunamganj district follow the usual trend of increasing up to 24-35 months age group. In case of wasting, no significant difference is observed among the districts. Age-group specific prevalence level of stunting, wasting and underweight for other districts are shown in supplementary Figures S.11-S.17.

and S.5). The comparison of DIR and SAE estimates in
District level prevalence of undernutrition indicators by children age-groups are also mapped in Figure S.10. The maps show that the spatial pattern varies with the increase of children's age for both stunting and underweight. On the one hand, for stunting, the map's colour tends to move from deep purple for younger children to orange and light yellow colours for older children (specifically those aged 24-35 months). For underweight, this pattern change from purple to orange-yellow colours remain the same with increasing children age. On the other hand, for wasting, no spatial pattern is visible except some differences for the 12-24 age-group. Instead of mapping the district level estimates of stunting, wasting and underweight for the corresponding cross-classified domains of children agegroups and place of residence, prevalence levels for three districts Dhaka, Khulna and Sunamganj are plotted in Figure 12, to demonstrate the benefits of multilevel modelling. This figure reveals higher stunting levels for the urban children aged 12-23 and 24-35 months and negligible rural-urban difference for other age-groups lead to higher prevalence of stunting in urban parts of Dhaka district. However, considerable rural-urban differences are observed in case of underweight and wasting for all age-groups in Dhaka district. In Khulna district, considerable rural-urban differences exist for most of the age-groups (except 48-59 group for wasting), however these difference seem negligible in the case of Sunamganj district. Confidence bands of the estimates indicate that estimates for rural domains of Dhaka district and urban domains of Khulna and Sunamganj districts are less precise. These findings are also supported by the fact that the direct estimates (black circle), particularly for many rural domains of Dhaka district and urban domains of Khulna and Sunamganj districts, are unreasonably very high or low compared to the expected level with respect to the corresponding age-group and residence (lying outside the SAE confidence bands). The differences in the prevalence level of stunting, wasting and underweight by children age-group and place of residence for other districts are shown in supplementary Figures S.18-S.25.
The estimated stunting level for the detailed level domains of districtresidence-age-sex shown in Figure 2 shows some sex differences particularly for the 12-23 and 24-35 months age-groups. These differences mainly arises in the DIR estimates for sample size differences particularly for the urban domains, and SAE estimators overcome these differences. Figure 2 shows considerable gender-inequality in the prevalence of stunting for the domains of Sunamganj district. These rural-urban differences also exhibit by gender in some districts of Rajshahi (see Joypurhat and Naogaon in Figure

Discussion
In this study, multilevel models expressed in a Bayesian framework have been developed for child malnutrition indicators of stunting, wasting and underweight in Bangladesh for spatio-demographic cross-classified small domains of Undernutrition prevalence in spatio-demographic domains sixty-four districts, two places of residence (urban and rural), sex and their five age-groups. Data from the recent Bangladesh Multiple Indicator Cluster Survey 2019, in which all the 64 districts are considered in the sampling design is used to demonstrate this. Survey-weighted direct estimates of each undernutrition indicator and their standard errors for the cross-classified domains calculated from the MICS 2019 microdata have been used as inputs in developing the multilevel models. The estimated standard errors for the domains with very small observation are very noisy (with a significant number of zero estimates). Hence, the direct standard error estimates are smoothed employing a generalized variance function approach before utilizing them in the models as inputs. The prevalence of stunting, wasting and underweight for the target spatio-demographic small domains are predicted from the respective multilevel models, and then the predictions are aggregated to estimate the prevalence at higher aggregation levels. Doing this aggregation procedure from bottom to top ensures that numerically consistent estimates are produced at various aggregated administrative units and their cross-classified demographic domains, for instance, we can provide unbiased estimates of rural female infants (aged 0-11 months) in the Dhaka division. District level prevalence and their standard errors are plotted simultaneously as bivariate choropleth maps, which help endusers and policy makers to understand the accuracy of vulnerability through visual presentation of the derived model-based estimates.
The developed multilevel models have allowed for cross-sectional and spatial correlations that increase the effective sample size of the detailed level domains through specification of various fixed effects and random effects of the cross-classified variables of district (district), place of residence (residence), sex of child (sex) and age-group (age) along with the sub-national variable (division). All the three finally selected multilevel models have common fixed effects components, which consist of only the main effects of sex, age, residence and division. The interaction effects of the demographic variables with administrative units are allocated as random effects in the model (to account for the unobserved differences that exist at the district and division levels). After accounting for the district and division level random effects, the fixed effects components indicate that children aged 24-35 months, living in Sylhet division and rural domains have higher possibility of being stunted and underweight. These patterns are expected, based on existing studies which have explored risk factors of children undernutrition in Bangladesh (Bhowmik & Das, 2019;M.R.K. Chowdhury et al., 2016;Das & Rahman, 2011). For wasting, most of the fixed effects are found insignificant except for sex (male domains have slightly higher prevalence) and children age-group (12-23 months age-groups has lower prevalence than infants). In some studies (T.R. Chowdhury et al., 2020;S.J. Rahman et al., 2021), male children were found to have higher risk of being wasted. However, infants were found to have lower risk of wasting than other children in older age-groups. This is reflected in findings in most studies. Due to having inconsistent direct estimates for the urban domains compared to rural domains, the rural-urban variations are found significant at various disaggregation levels, such as district-residence, district-residence-sex, and district-residence-age. These variations are accounted for by specifying random effects for the respective detailed level domains with unequal variance components for rural and urban domains. In a similar way, to account for age-specific variation at the district level, random regression coefficients for children age-groups are assumed to vary over districts with diagonal variance structure (i.e. with zero correlation between age-groups). The correlation between the districts are mainly captured by the spatial random effects component, which borrow spatial strength in the model from the neighbouring districts.
District level variation for 0-11 and 12-23 months age-groups are found significant for all models, while the variation for 48-59 months age-group is found significant for stunting and underweight. The spatial component had less contribution for wasting and so it was removed from the finally selected model of wasting. The spatial distribution of stunting, wasting and underweight at the district (Figures 6,7,8) and district-residence ( Figure 9) levels also suggest that spatial correlation among the districts are stronger for stunting and underweight compared to wasting. In addition to this random effects components, two white noise components are added finally to account for unexplained (residual) variation at the detailed level domains of district (district-residence-age-sex) and division (division-residence-age-sex).
The direct estimates at the detailed level were too much volatile and inconsistent, particularly for wasting there were many domains with zero (and also one) estimates as well as zero standard errors. The developed multilevel models provide reasonable estimates with better accuracy for these domains with zero and one estimates. Model assessment through four model discrepancy measures (relative bias, absolute relative bias, relative reduction of standard errors and ratio of coefficient of variations) confirms that multilevel model-based estimator provides approximately unbiased and more consistent estimate for most of the cross-classified domains of division and district. The aggregation from the very detailed level to the national level also confirms that the model-based estimates are consistent to the direct estimates at the higher aggregation levels such as national, division, residence, children sex and age-class.
The national-level prevalence of stunting and wasting have been declined in 2019, however these levels are still within the label of "high" prevalence (20-<30% for stunting and 10-<15% for wasting) according to the recent WHO threshold of stunting and wasting (De Onis et al., 2019). The level of underweight is also high (20-<30%) according to old referenced WHO threshold of underweight (WHO, 2010). Division level prevalence indicate that only the Khulna division reached to the SDG target of around 20% stunting level, while wasting level is almost double of 5% for all the divisions. Children of Sylhet division are highly vulnerable to all the undernutrition indicators, followed by the children of Mymensingh and Barishal divisions. Vulnerability scenarios differ only in Dhaka division, where children are more stunted compared to the underweight and wasting levels. This different scenario is attributed to the densely populated capital district Dhaka where stunting level is higher among the urban children than the rural children. Though the rural-urban difference at the national level is lower for all the indicators, the difference increases with the disaggregation for many domains.
Children aged more than one year particularly 12-23 and 24-35 months old are highly vulnerable to stunting and underweight compared to infants. However, more than one-quarter infants of Sylhet and Mymensingh divisions are stunted. On the other hand, infants are prone to loss their weight very quickly, however the wasting level does not vary too much by children agegroups, but their cross-classification with sex and place of residence indicate consistent rural-urban difference over all age-group and male children have higher wasting level. For underweight, the rural-urban difference increase with children age-groups and rural females are seemed to be more vulnerable than rural males. These age-sex-residence specific differences vary considerably with different degrees at the division level. As for example, urban children of Dhaka division have slightly higher stunting level for most of the age-groups, which is opposite to other divisions. On the other hand, rural-urban differences in underweight and wasting are considerably higher in less populated divisions, particularly for the female children.
Spatial distributions of district level stunting and underweight indicate that children of north and north-eastern districts are highly vulnerable to stunting and underweight compared to those districts in south-western districts, while children in the central region are more vulnerable to stunting than underweight and wasting. The spatial patterns of stunting also reveal that the level started to increase with the move from the south-western region to the north-eastern region through the central and north regions of the country. According to the bivariate choropleth maps, the estimates of stunting and underweight are more accurate in the south-western and central regions. It may happen due to two reasons either the sample size or the magnitude of the estimates. The estimates for the districts under Sylhet division (smaller sample size with higher estimates) are found with higher standard errors compared to those of in Khulna (smaller sample size with lower estimate) and Dhaka (bigger sample size with moderate estimate) divisions. The choropleth map of wasting suggests districts with higher wasting level with higher standard errors are sparsely distributed across the country.
District level spatial distributions by place of residence are expected to have higher heterogeneity among the urban parts of districts compared to their rural counterparts, since urban-specific variance components were considerably higher than those for rural parts in the model development. All the undernutrition indicators are higher in most of rural parts of districts with higher precision except in the Dhaka district. Only in Dhaka district, stunting level for urban part is slightly higher with higher precision than rural part. The ruralurban difference become wider in the district level (Kushtia, Narail, Khulna, Netrokona, Nawabganj, and Panchagarh are some example) compared to the division level as expected, however when urban domains have higher prevalence, then the rural-urban difference is minimal particularly in case of stunting and underweight (Comilla, Dhaka, Patuakhali, Kishoreganj, Manikganj, Sherpur, and Joypurhat are some example). In the highly urbanized districts (like Dhaka, Chittagong, Gazipur, and Narayanganj), the rural-urban difference are lower, which may be due to having many slums in the densely populated city areas (Angeles et al., 2009;M.S. Rahman, Mohiuddin, Kafy, Sheel, & Di, 2019).
District level spatial distributions by children age indicate stunting and underweight levels are about 10-15% higher among the older age-groups (particularly 12-23 and 24-35 months) than the infants, and these scenario is more prevalent in most of the districts belonging to vulnerable Sylhet, Mymensingh and Barishal divisions. In the Khulna division, the differences in the prevalence of stunting and underweight between infants and younger age-groups are also persistent but the levels for the infants are around 15% for most of the districts. So significant inequalities in the prevalence of stunting and underweight are mainly due to the infants and the other older age-groups. Since age-class has lower influence in the prevalence of wasting, modeling wasting is highly complex than the stunting and underweight.
In case of cross-classified domains of districts with children age-group and residence, the rural-urban difference has increased as expected but the patterns of rural-urban differences are very similar to those at the district-by-residence level for each of the age-groups. These patterns remain similar when the disaggregation down to district-residence-age-sex. The direct estimates at these detailed cross-classified domains are too volatile and model-based estimators have reduced this volatility in most of the cases. At the detailed level, the urban domains with unreliable estimates of one are adjusted by the modelbased estimates but the rural-urban differences remain comparatively larger for some domains. Some models are developed by ignoring these domains to examine whether these point estimates and their smoothed standard errors have any influence in such higher level predictions. The underlying differences were negligible and so all the domains were considered in the final multilevel models.
The main contribution of this study is the prediction of child undernutriton indicators at various disaggregation levels of division and district by developing models at the detailed cross-classified domains of children sex, age, place of residence and district by accounting variation in the target outcomes at the detailed as well as higher aggregation levels. The considered multilevel modelling approach helped to provide reliable and consistent official statistics as well as to explore the administrative hierarchies and demographic factors for which inequalities in the undernutrition vulnerability are still persistent in Bangladesh. Additionally, the bivariate choropleth map might help the policymakers to take their decisions based on the estimates as well as their quality. Undernutrition prevalence in spatio-demographic domains Concurrence of stunting, wasting and underweight may lead to correlation in between their prevalence level at various aggregation levels (Khara, Mwangome, Ngari, & Dolan, 2018). Developing univariate models for each of the undernutrition indicators by ignoring these correlations is one of the limitations of this study. A multivariate small area model would bring additional strength from the cross-correlation among the three indicators (Benavent & Morales, 2016;Nesa, 2019). However, the multivariate multilevel model at the considered detailed level would require survey-weighted co-variance matrix of dimension 3 × 3 for each of the detailed level 1280 domains. Due to having small sample size, estimation of consistent covariances would be more complex for many domains. In addition, the unavailability of Z-scores for the three indicators for all children would lead to loss of information for a significant proportion of domains. Most importantly, the estimation of a multivariate multilevel model would be more complex, time and computationally intensive at the considered detailed level. Furthermore, smoothing a variance-covariance matrix (compared to the vector of sampling error variance) might be more complex through generalized variance function approach. Nonetheless, these issues related to multivariate multilevel model are current areas of future research. UNICEF (2020). An additional 3.9 million children under 5 could suffer from wasting in South Asia this year due to COVID-19. https://www.unicef.org/bangladesh/en/press-releases/additional-39 -million-children-under-5-could-suffer-wasting-south-asia-year-due. (Accessed: 2021-08-25) Victora, C.G., Adair, L., Fall, C., Hallal, P.C., Martorell, R., Richter, L., . . . others (2008