Risk factors of concurrent malnutrition among children in Ethiopia: a bivariate spatial modeling approach

Abstract

Malnutrition among children remains a challenge to public health in Africa, and it is responsible for the increased infant mortality rate in Ethiopia. This study aims to identify the joint determinant factors of malnutrition among under-age five children in Ethiopia and quantify the regional prevalence across the country. A Bayesian hierarchical linear mixed model, with a bivariate conditional autoregressive model, was adopted to account for the joint spatial prevalence pattern of malnutrition among under-age five children in Ethiopia. The data were provided by the Demography and Health Surveys program. The study revealed that gender, wealth index, mothers' education, toilet system, fever, birth order, birth interval, contraceptive use, and diarrhea, are significant risk factors of child malnutrition. The three malnutrition indicators are most prevalent in the Afar region, while stunting and underweight are most prevalent in Amhara, Beneshangul Gumuz, and Tigray, and the wasting and underweight are most prevalent in the Gambela and Somali. Findings revealed that a stunted child is more likely to be underweight than being wasted, and a wasted child is more likely to be underweight than being stunted. The findings also revealed significant changes in the prevalence of malnutrition for different mothers and children age groups. For cost-effective malnutrition intervention, programs for parents, such as awareness creation about the modern types of contraceptives, appropriate birth spacing, the benefits of antenatal and postnatal care, for example, are a propitious method to mitigate malnutrition in high prevalent regions than direct intervention to the children.

Abbreviations: BICAR, Bivariate Conditional Autoregressive; DHS, Demographic and Health Survey; EDHS, Ethiopia Demographic and Health Survey; EA, Enumeration Area; GMRF, Gaussian Markov Random Field; HAZ, Height for Age Z-sscore; INLA, Integrated Nexted Laplace Approximation; MCAR, Multivariate Conditional Autoregressive; MCMC, Markov Chain Monte Carlo; SNNPR, Southern Nations, Nationalities, and People's Region; WAZ, Weight for Age Z-score; WHZ, Weight for Height Z-score; WHO, World Health Organization

Keywords:

• Malnutrition prevalence
• bivariate conditional autoregressive
• Bayesian spatial model
• hierarchical linear mixed model
• INLAMSM

Introduction

Malnutrition among children remains a challenge to public health in sub-Saharan Africa countries. The rationale is that childhood malnutrition makes children more susceptible to illness due to the weakened immune system caused by malnutrition. Such illnesses include malaria, measles, anemia, and acute respiratory infections (França et al. 2009; Walson and Berkley 2018). For this reason, children suffer negative consequences throughout their lifespan, resulting in shorter height, school dropout, reduced adult income, and decreased offspring birth weight (Victora et al. 2008; Magalhaes and Clements 2011).

Malnutrition among children is classified into stunting, wasting, and underweight, which are calculated by the weight, height, and age indexes of children. Specifically, as stated by WHO et al. (2013), children are stunted if their height-for-age z-score (HAZ) is lesser than negative-two standard deviations ($-2SD$) from the World Health Organization (WHO) Child Growth Standards (WCGSM) median. Underweight occur if weight-for-age z-score (WAZ) is lesser than $-2SD$ from the reference median, while wasting occur if the weight-for-height z-score (WHZ) is lesser than $-2SD$ from the reference median.

In the 2018 UNICEF et al. (2018) report, 149 million children below five years were classified as stunted, which represented an immense number that did not reach expected genetic heights. Statistics revealed that there were 50 million wasted and 40 million underweight children worldwide. Among Eastern and Southern African children under age five, 29 million were stunted, 5.40 million were wasted, and 3.60 million were underweight. Between the years 2000 and 2018, the prevalence in East Africa dropped from 45.5% to 35.2% for stunting and from 7.1% to 2.2 % for wasting. Despite there being evidence of these prevalence declines in East Africa, Ethiopia is still classified among high burden countries (WHO et al. 2020). According to UNICEF et al. (2018), Ethiopia had a 14% decline in stunting, 5% decline in wasting, and 12% decline in underweight between the years 2005 and 2019. Despite these improvements, the prevalence level is high in comparison with WHO thresholds (WHO et al. 2020).

Several kinds of research have been undertaken for decades to understand better the spatial distribution of undernutrition among under-age five children in highly prevalent sub-Saharan African regions. Kandala et al. (2011) used a geo-additive semi-parametric mixed model to measure the spatial distributions of childhood malnutrition in the Democratic Republic of Congo (DRC). Their findings emphasized the importance of accounting for spatial heterogeneity in malnutrition studies. Khatab (2010) adopted Geoadditive models to determine the risk factors of childhood malnutrition in Egypt. Unlike most studies in the developing countries, the authors found that maternal education had little impact on the prevalence of stunting and wasting. Adeyemi et al. (2019) adopted a generalized model to assess the combined spatial distributions of anemia and malnutrition among children in Burkina Faso and Mozambique. The authors found that there is evidence of the comorbidity of anemia and malnutrition. Gayawan et al. (2019) used Bayesian quantile regression to determine the spatial patterns of undernutrition quantiles among under-five children in Nigeria. Recently, Egbon et al. (2021a) used a generalized linear mixed model to quantify the effect of Carbon (IV) Oxide concentration on undernutrition among under-five children in Nigeria. Their finding showed a high correlation of ${\mathrm{C}\mathrm{O}}_2$ concentration to a higher risk of undernutrition in Nigeria. The aforementioned studies identified similar risk factors that amplify the prevalence of malnutrition among children. These include maternal education, age, and body mass index, household wealth quantile, child gender and age, area of residence, among others.

In Ethiopia, studies to better understand the risk factors and spatial distributions of undernutrition among children have attracted many researchers recently. For example, Zewdie and Abebaw (2013) applied logistic regression to examine the key factors of childhood undernutrition in Eastern Hararghe zone, Oromia, and found that farm size and morbidity incidence are among the key risk factors strengthening the prevalence of malnutrition in the study region. Takele and Taye (2014) utilized a semi-parametric model within a Bayesian framework to determine the risk factors of undernutrition in the country, and Gelu et al. (2018) conducted a cross-sectional study and adopted bivariate and multivariate logistic models to determine the risk factors of undernutrition among children in slum areas of Gondar city, northwest Ethiopia. The mentioned studies identified similar demographic and socioeconomic risk factors of malnutrition, such as maternal education and age, household wealth quantile, child gender and age, father's education status, diarrhea, fever, paternal household resource control. Further, Hagos et al. (2017) identified spatial clustering in Meskane Mareko District, Ethiopia, and explore the impact of socioeconomic and demographic factors on the prevalence of stunting using multilevel logistic regression. In a different study, Hailu et al. (2020) investigated the presence of spatial clustering and adopted a multilevel logistic regression to identify risk factors and spatial areas with elevated prevalence of stunting in Ethiopia. In both studies, the authors found evidence of spatial clustering and significant risk factors, which include maternal education, child age, anemia status, multiple births, maternal undernutrition, and household food security. Shebeshi and Asfaw (2020) studied the importance of location to under-five malnutrition in Ethiopia. The authors utilized a quantile interval model to determine the factors causing undernutrition. Takele et al. (2020) adopted a multivariate logistic model to explain the combined spatial effect of malnutrition indicators and found similar risk factors. Egbon et al. (2021b) investigated the risk factors of malnutrition and anemia concurrence using a Bayesian hierarchical linear mixed model, with stochastic partial differential equations. The authors found that under-five children below 20 months are more vulnerable to concurrence.

Findings have shown that there is a concurrence of stunting, wasting, and underweight among under age five children (Khara et al. 2018). However, previous studies failed to account for the spatial concurrence of the classes of malnutrition in Ethiopia. Beside Takele et al. (2020), the studies reviewed for Ethiopia focused on the modeling of a single malnutrition indicator, and contributes little or no information for the understanding of concurrence of malnutrition indicators. Although there is evidence in the literature that revealed an excess risk of underage five mortality due to concurrence of stunting, wasting, and underweight (McDonald et al. 2013), the factors responsible for this state of ailment among under age five children are poorly understood. It is not clear whether there is a direct causal relationship among malnutrition indicators or that the concurrence prevalence varies across regions and clusters. Thus, an adequate statistical methodology must be adopted, especially in Ethiopia with a limited study on the concurrence of malnutrition indicators. Moreover, it is known that observed information is underutilized when subjected to discretization, as the one conducted by Takele et al. (2020) who discretized malnutrition outcome variables, to suit a predefined modeling framework. Such inference could be improved by allowing the outcome variable to decide the most appropriate model to harness its content to determine the unknown process generating the data. In this regard, this study incorporated joint modeling of malnutrition by accounting for all uncertainty in the observed data.

This study aims to determine risk factors of malnutrition prevalence among under-five children in Ethiopia and quantify the joint regional prevalence and its concurrence across the country using a Bayesian hierarchical linear mixed model, jointly accounting for the spatial heterogeneity. It is essential to understand the spatial distribution of the concurrence of malnutrition to identify the most susceptible regions and their consequential health impact on public health for sustainable development. Spatial contiguity is an essential ingredient to incorporate regional heterogeneity into the model, such that, closer regions share similar malnutrition risk compared to farther ones. The spatial contiguity is incorporated into a hierarchical model using a random field model. This study adopted a Bivariate Gaussian Markov random field in the form of a Bivariate Intrinsic Conditional Autoregressive (BICAR) model to enable the pairwise estimation of the conditional spatial correlation between malnutrition indicators, such as stunting and underweight (HAZ and WAZ), stunting and wasting (HAZ and WHZ), and underweight and wasting (WAZ and WHZ). It provides a flexible framework to account for the joint spatial distribution between two latent Gaussian fields. Moreover, the probability of high malnutrition prevalence was estimated.

The remaining sections are organized as follows. Section Materials and method discusses the sampling scheme of the data set, study area, the statistical model adopted, and the estimation procedure. The results of analysis is presented in Section Results, while Section Discussion discuss the findings. The conclusion is in Section Conclusion.

Materials and method

Study area and data collection

The study scope is the administrative division of Ethiopia, with longitude ${38}^{\circ }$00 E and latitude $8^{\circ }00$ N, in the eastern part of Africa. Figure 1 represents the 11 administrative division of Ethiopia.The data used in this study were provided and granted access by Demography and Health Surveys (DHS) program upon request. Specifically, Ethiopia Demography and Health Service Program EDHS (2016) data provide valuable information about the structure of populations and health indicators over time. The sampling scope was the enumeration areas (EAs) prepared for the country's Census. The survey was designed to provide valuable information about key health and demographic indicators at the country and regional levels. For administrative purposes, Ethiopia has nine regions and two administrative cities, adding up to 11 regions (Figure 1).

Figure 1. Map of Ethiopia.

A two-stage cluster design sampling scheme was utilized. At first stage, 645 (202-urban, 443-rural) EA were selected with sampling probability proportionate to EA size. In the second stage, a systematic sampling with equal probability was adopted to select 28 households in each EA.

In all, 18,008 (5659 from urban and 12,349 from rural) households were selected, of which, 16,650 households were interviewed. That is 5232 from urban and 11,418 from rural settlements with 98% of response rate. All women aged 15–45 present in the household at the enumeration time were eligible for interview. In total, data for 11,144 (1219 from urban and 9925 from rural) women were observed, which amount to 95% response rate EDHS (2016). After data cleaning, the total number of underage five children included in the study is 8940.

The covariates included in this study were chosen based on the global burden of disease and risk factors for children (Lopez et al. 2006). The choices were guided by the previous studies (Hailu et al. 2020; Shebeshi and Asfaw 2020; Takele et al. 2020) and the outcome variables of interest considered were the malnutrition indicators: HAZ, WAZ, and WHZ, which, respectively, measures stunting, underweight, and wasting.

Table 1 presents all the dependent and independent variables considered in this study.

Table 1. Study variables.

Variable	Description	Summary
Dependent/Outcome variable
Stunting (HAZ)	Standardized Height-for-Age.	Continuous: $\mathbf{M}\mathbf{i}\mathbf{n}\mathbf{i}\mathbf{m}\mathbf{u}\mathbf{m}=-5.98$, $\mathbf{M}\mathbf{e}\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{n}=-1.45$, $\mathbf{M}\mathbf{a}\mathbf{x}\mathbf{i}\mathbf{m}\mathbf{u}\mathbf{m}=5.98$
Underweight (WAZ)	Standardized Weight-for-Age.	Continuous: $\mathbf{M}\mathbf{i}\mathbf{n}\mathbf{i}\mathbf{m}\mathbf{u}\mathbf{m}=-5.81$, $\mathbf{M}\mathbf{e}\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{n}=-1.19$, $\mathbf{M}\mathbf{a}\mathbf{x}\mathbf{i}\mathbf{m}\mathbf{u}\mathbf{m}=4.14$
Wasting (WHZ)	Standardized Weight-for-Height.	Continuous: $\mathbf{M}\mathbf{i}\mathbf{n}\mathbf{i}\mathbf{m}\mathbf{u}\mathbf{m}=-4.97,\mathbf{M}\mathbf{e}\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{n}=-0.57,\mathbf{M}\mathbf{a}\mathbf{x}\mathbf{i}\mathbf{m}\mathbf{u}\mathbf{m}=4.91$
Fixed Effect Covariate
Education	Mothers' highest education level.	Categorical: 0- No Education, 1- Primary, 2- Secondary, 3 - Higher
Wealth	Family wealth index.	Categorical: 0- Poorest, 1- Poorer, 2- Middle, 3- Rich, 4- Richer, 5- Richest
Settlement	Place of settlement.	Binary: 0- Rural, 1- Urban
Gender	Child's gender.	Binary: 0- Female, 1- Male
Water source	House holed water source.	Binary: 0- Unimproved System, 1- Improved System
Toilet type	House hold type of toilet system.	Binary: 0- Unimproved, 1- Improved
Presence of cough	Child have cough in the last two weeks.	Binary: 0- No, 1- Yes
Presence of diarrhea	Child have diarrhea in the last two weeks.	Binary: 0- No, 1- Yes
Presence of fever	Child have fever in the last two weeks.	Binary: 0- No, 1- Yes
Breastfed	Child is being breast fed.	Binary: 0- No, 1- Yes
Mass media	Parents or caregiver listens to radio, TV, or reads newspapers.	Binary: 0- No, 1- Yes
Random Effect Covariate
Child's age	The age of the child (in months) on the day of survey.	$\mathbf{M}\mathbf{i}\mathbf{n}=0,\mathbf{M}\mathbf{e}\mathbf{d}=28,\mathbf{M}\mathbf{a}\mathbf{x}=59$
Mother's age	The age of the mother or caregiver (in years) on the day of survey.	$\mathbf{M}\mathbf{i}\mathbf{n}=15,\mathbf{M}\mathbf{e}\mathbf{d}=28,\mathbf{M}\mathbf{a}\mathbf{x}=49$
Spatial Effect Covariate
Region	The region identity among the eleven regions of the country.	Categorical $\in \{1,2,3,\dots ,11\}$
Antenatal and Postnatal care visit	Frequency of child visits to antenatal/postnatal care	Categorical $\in \{0,1,2,\dots 12\}$
Contraceptive	Type of contraceptive used	Categorical: 0- No use, 1- Traditional, 2- Modern
Birth order	A child birth position in the household	Categorical: 1- First-born, 2- Second-born, 3- Third born,…, 14- Fourteen-born
Birth interval	Birth spacing: the interval between successive births	Categorical: 1- less than a year, 2- less than two years,…, 20-less than twenty years
Cluster	Cluster identification number	Numeric: $1,2,3,\dots ,635$

Statistical method

The aim is to quantify the regional prevalence of malnutrition and identify the factors causing slow recovery while accounting for the joint spatial correlation between two anthropocentric indicators. The generalized hierarchical joint spatial model is a powerful statistical tool that accounts for the joint spatial heterogeneity of malnutrition indicators through a Gaussian Random Field (GRF).

Suppose that ${\mathbf{Y}}_i=(y_{1i},y_{2i}{)}^T$ and ${\mathbf{Y}}_i\mid {\mathit{\theta }}_i\sim N({\mathit{\theta }}_i,{\tau }^{-1}\mathbf{I})$ is a bivariate standardised malnutrition indicators $y_{1i}$, and $y_{2i}$ of a child in region i, assuming a Gaussian bivariate distribution with mean ${\mathit{\theta }}_i$ and a constant covariance matrix ${\tau }^{-1}\mathbf{I}$, where ${\mathit{\theta }}_i=({\theta }_{1i},{\theta }_{2i}{)}^T$ is the expected value of a child's malnutrition indicator in region i, $\mathbf{I}$ is an identity matrix, $y_{1i},y_{2i}\in \left\{HA{Z}_i,WA{Z}_i,WH{Z}_i\right\}$, $y_{1i}\ne y_{2i}$, and $i=1,2,\dots ,11$. The interest is to account for the factors determining ${\mathit{\theta }}_i$ through a bi-generalized additive linear mixed model given as (1) $\begin{array}{rl}{\theta }_{1i}& =x_i^T\text{β}+{\mathbf{c}}_i^T\mathit{\varphi }+\sum _{k=1}^L{f}_k\left(z_i\right)+u_{1i},\\ {\theta }_{2i}& =x_i^T\text{β}+{\mathbf{c}}_{\mathbf{i}}^{\mathbf{T}}\mathit{\varphi }+\sum _{k=1}^L{f}_k\left(z_i\right)+u_{2i},\end{array}$(1) where ${\mathbf{x}}_i=(x_{i1},x_{i2},\dots ,x_{ip}{)}^T$ is the column vector of a $p\times n$ design matrix $\mathbf{X}$, $\text{β}$ is a $p\times 1$ fixed effect vector, ${\mathbf{c}}_i=(0,0,\dots ,0,1,0,\dots ,0{)}^T$ is a $m\times 1$ sparse column vector whose only qth entry is 1 if the current child in question is from cluster q in region i, $\mathit{\varphi }$ is an $m\times 1$ vector of cluster-specific random intercept, $f_k$ is a nonlinear function of kth random effect modeled in a nonlinear form, and $u_{si},s\in \{1,2\}$ is the spatial random effect of region i for indicator s, to account for spatial heterogeneity. Equation (1(1) $\begin{array}{rl}{\theta }_{1i}& =x_i^T\text{β}+{\mathbf{c}}_i^T\mathit{\varphi }+\sum _{k=1}^L{f}_k\left(z_i\right)+u_{1i},\\ {\theta }_{2i}& =x_i^T\text{β}+{\mathbf{c}}_{\mathbf{i}}^{\mathbf{T}}\mathit{\varphi }+\sum _{k=1}^L{f}_k\left(z_i\right)+u_{2i},\end{array}$(1) ) can be re-written in vector form of the model effects, given as (2) ${\mathit{\theta }}_i={\left[\begin{array}{c}{x}_i\\ x_i\end{array}\right]}^T{\text{β}}^{\ast }+{\left[\begin{array}{c}{c}_i\\ c_i\end{array}\right]}^T{\mathit{\varphi }}^{\ast }+\sum _{k=1}^L{f}_k\left(\begin{array}{c}{z}_i\\ z_i\end{array}\right)+{\mathbf{u}}_{\mathbf{i}},$(2) where ${\text{β}}^{\ast }=(\text{β},\text{β}{)}^T$, ${\mathit{\varphi }}^{\ast }=(\mathit{\varphi },\mathit{\varphi }{)}^T$, and ${\mathbf{u}}_i=(u_{1i},u_{2i}{)}^T$. Specification (2(2) ${\mathit{\theta }}_i={\left[\begin{array}{c}{x}_i\\ x_i\end{array}\right]}^T{\text{β}}^{\ast }+{\left[\begin{array}{c}{c}_i\\ c_i\end{array}\right]}^T{\mathit{\varphi }}^{\ast }+\sum _{k=1}^L{f}_k\left(\begin{array}{c}{z}_i\\ z_i\end{array}\right)+{\mathbf{u}}_{\mathbf{i}},$(2) ) will be independent, as in the case of modeling single malnutrition indicators if we assume that all the effects are independent. In this analysis, we relaxed the assumption that the spatial heterogeneity of malnutrition indicators is independent and considered cases where they are spatially correlated. In the Bayesian perspective, every parameter is assumed to be random, and thus, is assigned a prior distribution. The priors could be non-informative, objective, or it can be elicited from expert opinion, or data in the form of power prior (Ibrahim et al. 2014; Louzada et al. 2021). The bayesian spatial framework provides the possibility to incorporate regional heterogeneity into the model, such that close regions share similar properties compared to farther regions. The regional dependency structure is incorporated into the model through the prior distribution assigned on the spatial effect ${\mathbf{u}}_i$.

Let ${\mathbf{u}}^s=(u_{s1},u_{s2},\dots ,u_{sI}{)}^T$, for $s\in \{1,2\}$ be a GRF of malnutrition indicator s with positive definite precision matrix ${\mathbf{Q}}_s$ corresponding to a graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where $\mathcal{E}$ is a set of the edges connecting the regions, and $\mathcal{V}$ is a set of the nodes, each denotes the eleven administrative regions in the country. For computational speed, ${\mathbf{Q}}_s$ is considered to be sparse. Let $\mathbf{W}=\left(w_{ij}\right)$ be an $11\times 11$ matrix such that $w_{ij}=1$ if node i and j are neighbors and $w_{ij}=0$ otherwise, and diagonal matrix $\mathbf{D}$ with entries equals to the numbers of neighbors for each nodes. Therefore, (3) $Q_s={\tau }_{{\mathbf{u}}^{\mathbf{s}}}\mathbf{D}(\mathbf{I}-\mathbf{W}),$(3) where ${\tau }_{{\mathbf{u}}^{\mathbf{s}}}$ is the precision associated to indicator s. Note that by this specification, ${\mathbf{Q}}_s$ is sparse, consequently ${\mathbf{u}}^s$ is said to be a Gaussian Markov Random Field (GMRF) according to $\mathcal{G}$, such that $u_{si}$ is independent of $u_{sj}$ given ${\mathbf{u}}_{-sij}$ ($u_{si}\perp \perp u_{sj}|{\mathbf{u}}_{-sij}$) if and only if the edge $\{i,j\}\notin \mathcal{E}$, where ${\mathbf{u}}_{-sij}$ indicates ${\mathbf{u}}^s$ excluding nodes i and j. Further, define $\mathbf{u}=({\mathbf{u}}^1,{\mathbf{u}}^2{)}^T$, where ${\mathbf{u}}^1$ and ${\mathbf{u}}^2$ are considered correlated GMRF with respect to $\mathbf{G}$, and define ${\mathrm{\Lambda }}^{-1}$ as (4) ${\mathbf{\Lambda }}^{-1}=\left[\begin{array}{cc}{\displaystyle \frac{1}{{\tau }_{{\mathbf{u}}^{\mathbf{1}}}}}& {\displaystyle \frac{\rho }{\sqrt{{\tau }_{{\mathbf{u}}^{\mathbf{1}}}{\tau }_{{\mathbf{u}}^{\mathbf{2}}}}}}\\ {\displaystyle \frac{\rho }{\sqrt{{\tau }_{{\mathbf{u}}^{\mathbf{1}}}{\tau }_{{\mathbf{u}}^{\mathbf{2}}}}}}& {\displaystyle \frac{1}{{\tau }_{{\mathbf{u}}^{\mathbf{2}}}}}\end{array}\right].$(4) Here, ρ is the correlation coefficient that measures the strength and direction of the relationship between the GMRFs. The precision and <roman > log < /roman > prior distribution on $\mathbf{u}$ is then given respectively as (5) $\begin{array}{rl}\mathbf{Q}& =({\mathrm{\Lambda }}^{-1}\otimes D(\mathbf{I}-\mathbf{W})\\ \pi \left(\mathbf{u}\right)& \propto \frac{2n}{2}\log \left(|Q|\right)(-\frac{1}{2}{\mathbf{u}}^{T}Q\mathbf{u}),\end{array}$(5) The above specification is a Bivariate Conditional Autoregressive (BCAR) model, which is a special case of the Multivariate Conditional Autoregressive MCAR (Ippoliti et al. 2015; Palmi-Perales et al. 2019). The precision matrix ${\mathrm{\Lambda }}^{-1}$ is assigned a Wishart prior distribution with r degrees of freedom and a fixed $2\times 2$ positive definite matrix ${\mathbf{R}}^{-1}$.

A second-order random walk (RW2) prior was assigned for the nonlinear effects on child's and mother's age. Let $f_k=(f_{k1},f_{k2},\dots ,f_{kT}{)}^T$, where $f_{kt}$ is the model effect on the knot t of the nonlinear function $f_k$. Assuming independent second-order increment and Gaussian error, then ${\delta }_{f_{kt}}^2=f_{kt}-2f_{k(t+1)}+f_{k(t+2)}\sim N(\mathbf{0},{\sigma }_{f_k}^2),$and thus, the distribution of $f_k$ follows as (6) $\pi \left(f_k\right)\propto {\left(\frac{1}{{\sigma }_{f_k}^2}\right)}^{(n-2)/2}\exp \{-\frac{1}{2{\sigma }_{f_k}^2}\sum _{t=1}^{n-2}({\delta }_{f_{kt}}^2{)}^2\}.$(6) Here $f_1$ accounts for the effect of child's age, and $f_2$ accounts for the effects of mother's age among under age five malnutrition.

Equation (6(6) $\pi \left(f_k\right)\propto {\left(\frac{1}{{\sigma }_{f_k}^2}\right)}^{(n-2)/2}\exp \{-\frac{1}{2{\sigma }_{f_k}^2}\sum _{t=1}^{n-2}({\delta }_{f_{kt}}^2{)}^2\}.$(6) ) is equivalent to the one given by Lindgren and Rue (2008) for a regularly ordered label t. The child's age was ordered at regular intervals from month 0 to 59, likewise the mother's age between 15 and 49 years. In each case, $t\in \{0,1,\dots ,59\}$ for child's age, and $t\in \{15,19,\dots ,49\}$ for mother's age. Thus, $f_{kt},k=1,2$ is the effect at knot t. The RW2 model is flexible due to its invariant property under addition, and also computationally convenient since it exhibits the Markov property (Lindgren and Rue 2008). The prior (6(6) $\pi \left(f_k\right)\propto {\left(\frac{1}{{\sigma }_{f_k}^2}\right)}^{(n-2)/2}\exp \{-\frac{1}{2{\sigma }_{f_k}^2}\sum _{t=1}^{n-2}({\delta }_{f_{kt}}^2{)}^2\}.$(6) ) is proportional to a normal distribution, which facilitates sampling from the posterior distribution during computation. An independent and identically distributed Gaussian prior, centered at zero with a variance ${\sigma }_{\varphi }^2$ was assigned to each component of $\mathit{\varphi }$. A multivariate Gaussian prior was assumed for $\text{β}\sim N(\mathbf{0},{\sigma }_{\text{β}}^2\mathbf{1})$, and an Gamma hyper-prior was assigned for the precision parameter, $\tau \sim Gamma(a,b)$. The Penalizing Complexity (PC) prior was assigned for $({\sigma }_{f_k}^2{)}^{1/2}$, with parameters $({\sigma }_0$, $\alpha )$, such that $p\left(\right({\sigma }_{f_k}^2{)}^{1/2}>{\sigma }_0)=\alpha$ and the distribution follows as (7) $\begin{array}{rl}\pi \left({\sigma }_{f_k}^2\right)& \propto \frac{\lambda }{2}({\sigma }_{f_k}^2{)}^{3/2}\exp (-\lambda \left({\sigma }_{f_k}^2{)}^{1/2}\right),\\ {\sigma }_{f_k}^2& >0,k=1,2,\end{array}$(7) where $\lambda =-\log \left(\alpha \right)/{\sigma }_0$. The same prior is also assigned for $({\sigma }_{\varphi }^2{)}^{1/2}$ with similar parameterization.

The joint posterior distribution of the parameters of interest is expressed as (8) $\begin{array}{rl}& \pi (\text{β},\mathit{u},f_1,f_2,\mathit{\varphi }|\mathbf{Y},\mathbf{X})\propto L(\text{β},\mathit{u},f_1,f_2,\mathit{\varphi }\mid \mathbf{Y},\mathbf{X})\\ & \pi \left(\text{β}\right)\pi \left(\mathit{\varphi }\right)\pi \left(\mathbf{u}|\mathrm{\Lambda }\right)\pi (f_1,f_2|{\sigma }_{f_1}^2,{\sigma }_{f_2}^2)\pi \left(\mathbf{M}\right),\end{array}$(8) where $\mathbf{M}=\{\mathrm{\Lambda },\tau ,{\sigma }_{f_1}^2,{\sigma }_{f_2}^2,{\sigma }_{\text{β}}^2,{\sigma }_{\varphi }^2\}$.

Estimation

Different estimation procedures exist in the literature, and the most commonly used method is the Markov Chain Monte Carlo (MCMC) algorithm. However, this study utilized the Integrated Nested Laplace Approximation (INLA) package in the R (R Core Team 2019) through the INLAMSM (Palmi-Perales et al. 2019) package to estimate the parameters of interest. Unlike the MCMC approach, INLA estimates the posterior marginal distributions rather than the joint posterior distribution. It has been shown in simulation studies that the INLA performs considerably faster than the MCMC estimation procedure, and efficiently manages computation of GMRF precision matrix. The hyperparameters used the for the Pc-priors are ${\sigma }_0=1$, and $\alpha =0.001$, and INLA default priors were assigned for the hyperparameter of Λ, and a default prior for the remaining hyperparameters, $a,b,r,{\sigma }_{\text{β}}^2,{\sigma }_{\varphi }^2$ and $R=\mathbf{I}$. The INLA estimation procedure can be found in Ruiz-Cárdenas et al. (2012), Martins et al. (2013), Blangiardo et al. (2013), Palmi-Perales et al. (2019), Egbon et al. (2021a).

The spatial-specific posterior probability, in which low values of the estimate is associated with high prevalence of malnutrition for each region, compared with the whole country was obtained from the spatial effect posterior marginal density, which follows as, (9) $p(u_{si}<0|D)={\int }_{-\mathrm{\infty }}^0\pi \left(u_{si}|D\right)\mathrm{d}{u}_{si},\mathrm{\forall }i,s,$(9) where $\pi \left(u_{si}|D\right)$ is the posterior marginal distribution of region i for indicator s. The threshold 0 was used since it is considered as the country global average of malnutrition. Regions with high probability indicates high prevalence, and consequently high risk of malnutrition. The approximation of (9(9) $p(u_{si}<0|D)={\int }_{-\mathrm{\infty }}^0\pi \left(u_{si}|D\right)\mathrm{d}{u}_{si},\mathrm{\forall }i,s,$(9) ) can be obtained using a Monte Carlo integration by computing the expected value $\mathbf{E}\left(I\right(u_{si}<0\left)|D\right)$ that is (10) $p(u_{si}<0|D)=\frac{\sum _h^{n}I(u_{s{i}_h}<0)}{n},$(10) where $u_{s{i}_h}$ are random draws from the posterior marginal distribution of spatial effect $u_{si}$, and $I(\cdot )$ is an indicator function.

Results

This section presents the results of the bivariate model described in Section Materials and method. That is, the bivariate model of stunting and underweight, stunting and wasting, and underweight and wasting. An independent single analysis of malnutrition indicators was also conducted using a univariate CAR model, and the results are presented in the Appendix. Note that low estimated coefficients of the risk factors indicate that such risk factors exacerbate the prevalence of malnutrition. That is, the lower the expected value of the coefficients, the lower the predicted values of HAZ, WAZ, and WHZ. We used the word ‘prevalence’ here to denote the risk factor's impact on the fraction of children whose predicted values of HAZ, WAZ, and WHZ scores are less than $-2.0$. The inclusion and non-inclusion of 0 in the confidence intervals were used to determine the significance of the estimated effects. The estimated spatial conditional correlation of stunting and underweight is $\rho =0.22$, stunting and wasting is $\rho =-0.31$, and wasting and underweight is $\rho =0.26$. Findings indicate that stunting and wasting are correlated positively with underweight, and whereas, they are negatively correlated together.

Fixed effect

Table 2 presents the estimates of the fixed effects for each bivariate model. The first panel contains the results of the bivariate estimates of stunting and underweight, the second panel contains the bivariate estimates of stunting and wasting, and the third panel contains the bivariate estimates of underweight and wasting.

Table 2. The fixed effect estimates of the bivariate model.

	Stunting and Underweight				Stunting and Wasting				Underweight and Wasting
	mean	sd	2.50%	97.50%	mean	sd	2.50%	97.50%	mean	sd	2.50%	97.50%
Gender (Ref=Female)
Male*	−0.057	0.023	−0.103	−0.011	−0.066	0.025	−0.115	−0.016	−0.051	0.021	−0.092	−0.010
Highest Education Level (Ref=No Edu)
Primary*	0.113	0.031	0.051	0.174	0.082	0.034	0.016	0.148	0.077	0.028	0.022	0.132
Secondary*	0.325	0.064	0.201	0.450	0.204	0.068	0.070	0.338	0.117	0.056	0.006	0.227
Higher*	0.322	0.088	0.149	0.494	0.248	0.094	0.062	0.433	0.224	0.078	0.071	0.377
Household Wealth Index (Ref=Poorest)
Poorer	−0.021	0.036	−0.093	0.050	−0.007	0.039	−0.084	0.069	0.007	0.032	−0.056	0.071
Middle*	0.207	0.040	0.129	0.285	0.136	0.043	0.052	0.220	0.094	0.035	0.025	0.164
Rich*	0.276	0.044	0.190	0.361	0.201	0.047	0.109	0.293	0.169	0.039	0.093	0.245
Richest*	0.366	0.059	0.250	0.482	0.299	0.064	0.174	0.424	0.303	0.053	0.200	0.406
Improved Toilet System (Ref=No)
Yes*	0.299	0.072	0.158	0.439	0.213	0.077	0.062	0.364	0.170	0.064	0.046	0.295
No. of children under 5 (Ref=0 & 1)
2–3	−0.043	0.038	−0.118	0.032	−0.032	0.041	−0.113	0.048	−0.024	0.034	−0.091	0.042
>4	0.108	0.083	−0.054	0.270	0.039	0.089	−0.135	0.213	−0.025	0.073	−0.169	0.119
Birth order (Ref=1)
2–4*	−1.461	0.055	−1.569	−1.353	−1.091	0.057	−1.203	−0.978	−0.990	0.047	−1.083	−0.898
>5*	−1.568	0.054	−1.674	−1.463	−1.168	0.056	−1.278	−1.059	−1.056	0.046	−1.146	−0.966
Birth Interval (Ref < 2 years)
2-5 years*	0.146	0.030	0.088	0.204	0.107	0.032	0.045	0.169	0.082	0.026	0.030	0.134
>5years*	0.346	0.049	0.249	0.442	0.263	0.053	0.158	0.367	0.217	0.044	0.131	0.303
Antenatal and Postnatal care visit (Ref=0&1)
2–3	0.159*	0.037	0.086	0.233	0.089*	0.040	0.010	0.167	0.021	0.033	−0.044	0.086
>3	0.551*	0.115	0.325	0.776	0.310*	0.123	0.068	0.552	0.079	0.102	−0.122	0.279
Type of contraceptive used (Ref=N0)
Traditional	−0.043	0.181	−0.398	0.311	0.076	0.194	−0.305	0.457	0.194	0.161	−0.122	0.508
Modern*	0.115	0.031	0.054	0.175	0.089	0.033	0.024	0.153	0.085	0.027	0.032	0.139
Settlement (Ref=Rural)
Urban	−0.070	0.053	−0.173	0.034	−0.061	0.057	−0.172	0.050	−0.055	0.047	−0.147	0.037
Mass media(Ref=No)
Yes*	0.111	0.031	0.051	0.172	0.096	0.033	0.031	0.161	0.103	0.027	0.050	0.157
Presence of Diarrhea (Ref=No)
Yes*	−0.100	0.039	−0.177	−0.022	−0.096	0.042	−0.179	−0.013	−0.121	0.035	−0.190	−0.052
Presence of Fever (Ref=No)
Yes*	−0.105	0.042	−0.188	−0.023	−0.104	0.045	−0.193	−0.016	−0.130	0.037	−0.204	−0.057
Presence of Cough(Ref=No)
Yes	0.000	0.039	−0.077	0.077	0.001	0.042	−0.082	0.084	0.004	0.035	−0.064	0.072
Child Being Breatfed (Ref=No)
Yes*	−0.244	0.029	−0.301	−0.186	−0.173	0.031	−0.234	−0.111	−0.127	0.026	−0.178	−0.076

The reference category is denoted as Ref. The significant categories are indicated by *, identified by the noninclusion of zero in the credible interval. A significant lower posterior mean indicates the risk factor contributing relatively larger to the prevalence of malnutrition.

In the first panel, stunting and underweight are significantly more prevalent among male children compared to their female counterparts. Likewise for stunting and wasting, and underweight and wasting in the second and third panels. Concerning the children whose mothers have no education, the prevalence of both stunting and underweight decreases for children whose mothers' attained primary, secondary and higher education level. Similar results were obtained for stunting and wasting and underweight and wasting. It means that children whose mothers' education level is at least primary are less likely to be malnourished compared with children whose mothers have no formal education. These findings indicate that educated mothers are more likely to grow healthier children than uneducated ones since they are more likely to give children fortified food and seek expert opinions. Turning attention to the household wealth index, children whose family has at least a middle wealth index are significantly less prevalent in stunting and underweight compared with the poorest families. Likewise for stunting and wasting and underweight and wasting. Similar to the univariate case (Appendix 2, Table A1), the malnutrition prevalence significantly declines for middle, richer, and richest. Nevertheless, family wealth is inversely related to the prevalence of malnutrition. That is, it is more likely that a child is malnourished in lower-income families than in higher-income families. The effect of the toilet system reveals that the prevalence of stunting and underweight decreases in the presence of an improved toilet system. The same occurs in stunting and wasting, and underweight and wasting. The number of children under five years in the household does not significantly contribute to malnutrition. However, as the birth order increases, the joint prevalence of stunting and underweight significantly increases, and likewise for stunting and wasting, and underweight and wasting. This implies that mothers or caregivers diminish the attention given to their children as they give births to more children. Prevalence of stunting and underweight among children with two years and above birth interval is significantly less likely. Similarly for stunting and wasting, and underweight and wasting. The frequency of visits to antenatal and postnatal care services significantly decreases the joint prevalence of stunting and underweight, stunting and wasting, but not underweight and wasting. A modern type of contraceptive used significantly lowers the joint prevalence of stunting and underweight, stunting and wasting, and underweight and wasting compared to no use of contraceptive, whereas, the traditional type is not significantly different from no use of contraceptive. The prevalence of malnutrition across settlement types is not significant. Prevalence of stunting and underweight among children whose mothers listen to mass media are significantly less likely compared to those mothers who do not listen to mass media. A child having diarrhea, fever is significantly more likely to be stunting and underweight, stunting and wasting and underweight and wasting. That is, children having any of these ailments increases the risk of malnutrition. These results show that episodes of diarrhea and fever are detrimental to child growth.

Spatial effect

The spatial effect estimated from the BCAR model is shown in Table 3, which contains three panels, where each panel is further divided into top and bottom, showing marginal estimates of one of the two anthropometric variables jointly modeled. Moreover, the left panel contains the marginal results of stunting from the joint modeling of stunting and underweight (top) and stunting and wasting (bottom). The middle panel contains the marginal results of wasting from the joint modeling of stunting and wasting (top) and underweight and wasting (bottom). Further, the right panel contains the marginal results of underweight from the joint modeling of stunting and underweight (top) and of underweight and wasting (bottom). Figure 2 presents the probability of excess prevalence across the regions. This was calculated using samples drawn from the corresponding posterior distribution to evaluate Equation (10(10) $p(u_{si}<0|D)=\frac{\sum _h^{n}I(u_{s{i}_h}<0)}{n},$(10) ) for each region.

Figure 2. The probability of excess prevalence of (a) stunting in the presence of underweight, (b) wasting in the presence of stunting, (c) underweight in the presence of stunting, (d) stunting in the presence of wasting, (e) wasting in the presence of underweight, and (f) underweight in the presence of wasting. High probability value signifies strong evidence of malnutrition prevalence that exceeds the country's average prevalence, and consequently, increase in susceptibility of children in the region. (a) Stunting${}_{-\left(\mathrm{s}\mathrm{t}\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}\right)}$, (b) Wasting${}_{-\left(\mathrm{s}\mathrm{t}\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{w}\mathrm{a}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{g}\right)}$, (c) Underweight${}_{-\left(\mathrm{s}\mathrm{t}\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}\right)}$, (d) Stunting${}_{-\left(\mathrm{s}\mathrm{t}\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{w}\mathrm{a}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{g}\right)}$, (e) Wasting${}_{-\left(\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{w}\mathrm{a}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{g}\right)}$ and (f) Underweight${}_{-\left(\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{w}\mathrm{a}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{g}\right)}$.

Table 3. Spatial effect estimated from the BCAR model.

	Joint model of Stunting and Underweight				Joint model of Stunting and Wasting				Joint model of Stunting and Underweight
	Stunting				Wasting				Underweight
Region	mean	sd	2.5%	97.5%	mean	sd	2.5%	97.5%	mean	sd	2.5%	97.5%
Addis Ababa	0.258*	0.092	0.078	0.438	−0.206*	0.098	−0.398	−0.014	0.253*	0.090	0.077	0.430
Afar	−0.165*	0.054	−0.271	−0.058	−0.196*	0.058	−0.310	−0.082	−0.284*	0.053	−0.388	−0.180
Amhara	−0.457*	0.051	−0.557	−0.357	0.162*	0.055	0.055	0.269	−0.198*	0.050	−0.296	−0.100
Beneshangul Gumuz	−0.356*	0.055	−0.464	−0.249	0.104	0.059	−0.011	0.219	−0.160*	0.054	−0.266	−0.054
Dire Dawa	0.021	0.073	−0.123	0.164	−0.373*	0.078	−0.527	−0.220	−0.169*	0.072	−0.310	−0.029
Gambela	0.391*	0.064	0.266	0.516	−0.174*	0.068	−0.308	−0.041	0.128*	0.063	0.005	0.251
Harari	0.114	0.070	−0.025	0.251	−0.136	0.075	−0.284	0.012	0.065	0.070	−0.072	0.201
Oromia	0.013	0.042	−0.069	0.095	0.357*	0.045	0.269	0.445	0.162*	0.040	0.083	0.241
SNNPR	−0.085	0.047	−0.176	0.007	0.518*	0.050	0.420	0.616	0.230*	0.046	0.139	0.320
Somali	0.573*	0.046	0.483	0.663	−0.149*	0.049	−0.246	−0.053	0.127*	0.045	0.038	0.216
Tigray	−0.307*	0.050	−0.406	−0.208	0.093	0.054	−0.013	0.199	−0.154*	0.050	−0.252	−0.057
	Joint model of Stunting and Wasting				Joint model of Underweight and Wasting				Joint model of Underweight and Wasting
	Stunting				Wasting				Underweight
Region	mean	sd	2.5%	97.5%	mean	sd	2.5%	97.5%	mean	sd	2.5%	97.5%
Addis Ababa	0.729*	0.097	0.538	0.921	−0.080	0.081	−0.239	0.080	0.848*	0.081	0.689	1.007
Afar	−0.251*	0.057	−0.364	−0.139	−0.214*	0.048	−0.309	−0.120	−0.399*	0.048	−0.493	−0.305
Amhara	−0.554*	0.054	−0.660	−0.448	0.131*	0.045	0.042	0.219	−0.325*	0.045	−0.413	−0.237
Beneshangul Gumuz	−0.428*	0.058	−0.543	−0.313	0.075	0.049	−0.020	0.170	−0.260*	0.048	−0.355	−0.166
Dire Dawa	0.172*	0.078	0.020	0.324	−0.357*	0.065	−0.484	−0.230	−0.019	0.064	−0.144	0.107
Gambela	0.429*	0.068	0.296	0.562	−0.149*	0.056	−0.260	−0.038	0.185*	0.056	0.076	0.295
Harari	0.264*	0.075	0.117	0.412	−0.096	0.062	−0.219	0.026	0.253*	0.062	0.131	0.375
Oromia	−0.157*	0.044	−0.243	−0.071	0.308*	0.037	0.236	0.381	−0.037	0.036	−0.108	0.035
SNNPR	−0.221*	0.050	−0.319	−0.123	0.479*	0.041	0.398	0.561	0.061	0.041	−0.020	0.141
Somali	0.409*	0.049	0.313	0.505	−0.170*	0.041	−0.250	−0.090	−0.051	0.040	−0.131	0.028
Tigray	−0.392*	0.054	−0.498	−0.287	0.072	0.045	−0.015	0.160	−0.257*	0.044	−0.344	−0.170

The spatial effects significantly different from zero are indicated by *, and identified by the noninclusion of zero in the credible interval. Relatively to the country's average prevalence, a lower posterior mean signifies the region with more aggravated prevalence of malnutrition.

From the top first panel, region Amhara and Beneshangul-Gumuz are the highest prevalent regions of stunting, in the presence of underweight, compared to the whole country. These are followed orderly by the regions Tigray, and Afar. This conclusion is evident in the probability plot in Figure 2(a). The black-colored regions are those regions with the highest probability of stunting in the presence of underweight. The results indicate that children who reside in these regions are more likely to be stunted in the presence of underweight compared to the whole country. However, the prevalence at Dire Dawa, Harari, Oromia and SNNPR is not significant. Addis Ababa, Gambela, and Somali are the regions with a lesser prevalence of stunting in the presence of underweight. Turning attention to the bottom first panel, the prevalence of stunting in the existence of wasting increased for Afar, Amhara, Beneshangul-Gumuz, Oromia, Tigray, and SNNPR compared to stunting in the presence of underweight (top panel). This result is evident in the probability plot in Figure 2(d), where these listed regions represent the black-colored regions, which are the regions with the highest probabilities of stunting in the presence of wasting. Observed that Oromia and SNNPR regions attained significance in contrast to the results of stunting in the existence of underweight. It implies that under-five children residing in these regions are more likely to be stunted in the presence of wasting compared to the presence of underweight.

From the top middle panel, Dire Dawa has the highest prevalence of wasting in the presence of stunting. Dire Dawa is followed orderly by Addis Ababa, Afar, Gambela, and Somali. These regions are represented in black color in the probability plot, indicating the regions with higher probabilities of wasting in the presence of stunting (Figure 2(b)). The prevalence of wasting in the presence of stunting in Harari, Beneshangul-Gumuz, and Tigray regions is not significant. They represent the grey portion in Figure 2(b), which indicates moderate prevalent probability. However, Addis Ababa, Amhara, Oromia, and SNNPR are regions with a lesser prevalence of wasting in the presence of stunting. Turning attention to the bottom of the middle panel, a similar prevalence pattern was obtained for wasting in the presence of underweight and in the presence of stunting. However, the prevalence of wasting in the presence of underweight increases for Somali, Afar, Amhara, Oromia, and SNNPR and decreases for Dire Dawa and Gambela compared to the prevalence of wasting in the presence of stunting. Moreover, Dire Dawa has the highest prevalence of wasting in the presence of underweight. It is followed orderly by Afar, Somali, and Gambela, representing the black-colored regions in the probability plot (Figure 2(e)), which indicates high probabilities of wasting in the presence of underweight. Similarly, the prevalence of wasting in the presence of underweight is not significant in Harari, Beneshangul-Gumuz, and Tigray regions.

From the top right panel, Afar has the highest prevalence of a child being underweight in the presence of stunting. Afar is followed orderly by Amhara, Dire Dawa, Beneshamgul-Gumuz, and Tigray. These regions are represented in black color (Figure 2(c)), which indicates the high prevalent regions of underweight in the presence of stunting. However, Addis Ababa, Gambela, Oromia, Somali, and SNNPR are regions with the lowest prevalence of underweight in the presence of stunting. Moreover, Harari are not significantly different from the country's average prevalence, and these regions are represented in grey in Figure 2(c). Turning attention to the bottom right panel, the prevalence of underweight in the presence of wasting further increased for Afar, Amhara, Beneshamgul-Gumuz, and Tigray, compared to the presence of stunting. These regions are represented bythe black-colored portion of Figure 2(f). This implies that children with wasting are more likely to be underweight compared to children with stunting. Moreover, the prevalence further decreased for Addis Ababa, Harari, Gambela, and Dire Dawa, which represent the lightest regions. This indicates that children residing in these regions are less likely to be underweight in the presence of wasting compared to the presence of stunting.

These results indicate that the prevalence of stunting in the presence of wasting and the prevalence of wasting in the presence of stunting is not equivalent across Ethiopian regions. Likewise for stunting and underweight, underweight and stunting, wasting and underweight and underweight and wasting across the country. Compared to the univariate models (Appendix 2, Figure A5), the bivariate model revealed high-prevalent regions that are ordinarily hidden in the univariate case, and thus, introducing a type two error when investigated in a univariate case. Regions where stunting is prevalent, wasting is not, and verse versa Figure 2(a,b,d,e). That is, the spatial relationship between stunting and wasting, concerning the Ethiopian graph is inversely proportionate. This characteristic was inherited from the univariate case. However, the prevalence of underweight partially shares high-prevalent regions between stunting and wasting.

Nonlinear effect

Figure 3 presents the relative change in the prevalence of malnutrition for different mothers and children age groups. Figure 3(a–c) is the joint prevalence of malnutrition to a child's age group, while Figure 3(d,e) is the joint prevalence of malnutrition indicators to mothers' age group. Figure 3(a) presents the prevalence of stunting and underweight, Figure 3(b) presents the prevalence of stunting and wasting, while Figure 3(c) presents the prevalence of wasting and underweight according to a child's age distribution. Figure 3(d) presents the prevalence of stunting and underweight, Figure 3(e) presents the prevalence of stunting and wasting, while Figure 3(f) presents the prevalence of wasting and underweight according to a mother's age distribution. It is important to note that the lower the posterior estimate, the higher the prevalence, since, lower malnutrition indicator score indicates the severity of malnutrition.

Figure 3. The relative change in the prevalence of (a) stunting and underweight, (b) stunting and wasting, (c) wasting and underweight with respect to a child's age, and (d) stunting and underweight, (e) stunting and wasting, (f) wasting and underweight with respect to a mothers' age. A declining curve signifies increasing prevalence, and a rising curve signifies a lowering prevalence.

The joint prevalence of stunting and underweight increases consistently until a child attains age 20 months, where it decreases shortly for an additional three months. The prevalence further increased until age 33 months, and then decrease shortly between age 33 months and 40 months, and nearly remain constant afterward until age 59 months (Figure 3(a)). Compared to the univariate analysis of stunting and underweight, the joint prevalence of stunting and underweight appears more like the prevalence of stunting (Appendix 2, Figure A4). This indicates that a child is more likely to be stunted and underweight with age. Moreover, in Figure 3(b), the prevalence of stunting and wasting increases in a similar way as that of stunting and underweight. The prevalence increases consistently from age 0 to age 20 months and further decreases and increases slightly until age 30 months and then decreases and levels off until age 59 months. In Figure 3(c), the prevalence of wasting and underweight increases gradually until age 30 months, and then maintains a constant pattern until age 40, which then increases gradually until age 59 months. The prevalence pattern of malnutrition to age is as expected because children are more exposed to poor environmental conditions as they grow older.

Considering mothers age, the prevalence of stunting and underweight decreases as the mother's age increases until age 22 years, then increase slowly until age 25, and gradually decreases until age 49 (Figure 3(d)). A similar shape was obtained for the joint prevalence of stunting and wasting (Figure 3(e)), and wasting and underweight (Figure 3(f)). It implies that children are less likely to be malnourished as mothers get older. The prevalent curve is as expected because, as the mothers grow older, their parental care improves due to gathered experience, maternal or formal education attainment.

The estimates of cluster-specific random effects for 635 clusters across the country, grouped by regions, are presented in Appendix 1. Each Figure contains 11 panels. The total number of clusters adopted for each region are as follows: Addis Ababa (55), Afar (53), Amhara (71), Beneshangul Gumuz (49), Dire Dawa (43), Gambela (49), Harari (41), Oromia (73), SNNPR (71), Somali (67), Tigray (63). The Oromia region had the highest number of clusters, whereas, Harari had the least number of clusters. Figure A1 presents the cluster-specific random effects of stunting and underweight. The result shows that there are cluster variations within the regions. While the prevalence of stunting and underweight is high for some clusters within a specific region, others are moderate and less prevalent. A similar pattern is also shown for stunting and wasting (Figure A2), and wasting and underweight (Figure A3).

Discussion

This study used a Bayesian hierarchical linear mixed model to determine the key factors causing malnutrition among children under age five in Ethiopia. This method was employed to determine the joint regional prevalence and its concurrence across the country. A bivariate intrinsic conditional autoregressive model was adopted to measure the joint spatial prevalence of malnutrition across Ethiopia. Findings reveal that place of residence, mothers' educational status, child gender, family wealth index, mother and child age, presence of diarrhea, fever, coughs, breastfeeding, and toilet type are significant risk factors of malnutrition.

Mothers' education level significantly determines whether a child is malnourished, implying that children whose mothers attain higher education are less likely to suffer from malnutrition. This finding conform with one reported in Ethiopia (Darsene et al. 2017; Shebeshi and Asfaw 2020; Takele et al. 2020) and other developing countries such as Nigeria (Gayawan et al. 2019), Democratic Republic of the Congo Kandala et al. (2011) and Bangladesh (Islam et al. 2013). However, based on UNICEF reports in 2018, only 53% of Ethiopia women complete primary education and only 25% of the aged girls attend up to secondary school, which consequentially accumulate to a high rate of uneducated mothers (UNICEF Ethiopia 2018). Culturally, not until recently, most ethnic groups in Ethiopia believe that women are for household activities, and as such, they engage in farming, household chores, and early marriage. A positive step to achieve sustainable development goal is to promote formal and informal maternal education, especially in the rural region of the country.

Further, children whose family wealth index is low are more likely to suffer concurrence of malnutrition. This may be a result of the poverty level and food insecurity, where children are denied adequate diet due to insufficient family income. On contrary, the wealthy family is more likely to provide fortified nutritious intake for their children, which helps the body and mental development and prevention of malnutrition (Yadav et al. 2015).

Male children are more likely to be malnourished compared with female children. This finding confirms the biological result that male children are more exposed to ailments compared with female children (Costa et al. 2017; Egbon et al. 2021b). This could also be attributed to attitude of parents or caregivers towards male children. For instance, parents wrongly assume that by nature, male children are more immune than their female counterparts. This assumption could consequently lead to exposure of male children to environmental contaminants causing malnutrition. This finding is consistent with study conducted in Ethiopia (Darsene et al. 2017), Nigeria (Gayawan et al. 2019), and Bangladesh (Islam et al. 2013). On contrary, a study conducted in Tanzania (Mbago and Namfua 1992), and Kenya (Gewa and Yandell 2012) reported that female children are more malnourished compared to male children.

Finding also revealed that diarrhea and fever significantly contributes to the concurrence of stunting, wasting and underweight among children, preventing rapid recovery and causing child mortality. Episodes of diarrhea, majorly caused by poor hydration hygiene among children and could weaken their immune system. This fragility creates room for infections and diseases that escalates concurrences of malnutrition. These findings coincide with previous studies in Nigeria (Gayawan et al. 2019) and Kenya (Gewa and Yandell 2012). Findings indicate that children under breastfeeding are more likely to suffer from joint occurrences of malnutrition. This finding is contrary to scientific knowledge that breastfeeds milk supplies adequate nutrients for child development (Motee and Jeewon 2014). This could be the causality effect, where mothers engage in prolong breastfeeding to make better the health of malnourished children (Gewa and Yandell 2012). As expected, unimproved toilet system such as pit toilets, open deification in household backyards, woods, and streams significantly increases the prevalence of the concurrence of malnutrition. These practices expose children to common infections and diseases and preventing recovery. This finding is consistent with the results obtained in Ethiopia (Takele et al. 2020) and Egypt (Khatab 2010) and contrary to the study reported in Nigeria (Egbon et al. 2021a), where the toilet system was not significant. Place of residence and the attitude of mothers and caregivers not listening to any form of mass media are not significant. This implies that children whose family resides in rural and urban are faced with similar risk of malnutrition. Findings revealed that the higher the birth order the higher the prevalence of malnutrition. This could be associated with strain on family resources (Horton 1988). Similarly, birth interval, contraceptive use, antenatal, postnatal care service, and mass media significantly influence child malnutrition. That is, low birth spacing, low frequency of antenatal and postnatal care service, use of no contraceptive or traditional contraceptive and no mass media significantly increases the prevalence of malnutrition.

As expected, an increase in child's age significantly contributes to the concurrence of malnutrition in children. It implies, younger children are less likely to be malnourished compared with older ones. It could be attributed to inappropriate feeding practices and insufficient nutritional content in solid food as they grow older (Hall et al. 2020). It could also be a consequence of children's exposure to environmental risk factors, such as ambient and household pollution, infections, airborne diseases, and more (Sinharoy et al. 2020). The extent of a child's exposure to environmental factors increases with age, and the mother's or caregiver's control over the children decreases, placing them at a higher risk of malnutrition. The concurrence prevalence of child malnutrition decreases with an increase in maternal age. It is expected that older-age mothers are more experienced in child caregiving than younger-aged mothers. That is, mothers at higher ages have gained more experience compared with younger ones. This result coincides with the study conducted in Ghana (Wemakor et al. 2018).

Findings revealed that the concurrences of malnutrition indicators vary across regions in the country. This may be due to variation in cultural practices and ethnic structure, variation in agricultural practices, religion, traditional living habits, education sector, and health facilities across the country. This regional diversity and resources are associated with the differences in the spatial distribution of malnutrition prevalence across the country.

Findings revealed that the three malnutrition indicators are most severe in the Afar region. The stunting and under-weight were consistently severe in Amhara, Beneshangul Gumuz, and Tigray, which occupy the northern part of the country. This finding is consistent with the results obtained in Ethiopia (Hagos et al. 2017; Takele et al. 2020). This could be associated with low socio-economic status in the regions. Findings show that more than 30% of inhabitants living in Tigray and Amhara are below the food poverty line (OCHA 2019). And around 25% of inhabitants living in Afar and Benishangul Gumuz are below the food poverty line (OCHA 2019). This leads children to a high degree of food deprivation. Also, this could be associated with unfavorable weather conditions for farming, inadequate water supplies, poor irrigation system, decades of locust invasion, which might have limited development gain in these regions (USAID 2021). It is common among Ethiopians to believe that the Oromia, SNNPR have good weather conditions and water supplies. However, the concurrence probability of stunting and being underweight is severe. This could be associated with political disruption and social unrest in the regions. For instance, the political movement (Abebe 2020) could affect emotional behavior in children and prevents access to health care services. Finding shows that wasting and underweight are highly prevalent in Somali (Shebeshi and Asfaw 2020; Takele et al. 2020), Dire Dawa and Gambela (Tusa et al. 2020) regions. This could be due to ethnic structure, traditional believes, and the living habit of household members in the regions. It is surprising that the prevalence of wasting for a stunted child in Addis Ababa is higher than the countries overall prevalence and could require a closer investigation. Other findings indicate that a relationship between stunting and underweight, wasting and underweight were spatially positive, while stunting and wasting have a negative spatial correlation. A stunted child is more likely to be underweight than being wasted, and a wasted child is more likely to be underweight than being stunted. Findings revealed cluster heterogeneity among regions in the country. While Oromia had the highest number of clusters, Harari had the least number of clusters. The differences in the cluster-specific estimates could be associated with the extent of ethnic diversity and disparity in cultural believes in these regions, which could have a considerable impact on parents' and community's perception of malnutrition among children.

An important finding revealed that the prevalent factors associated with parents and caregivers are more likely to determine the malnourished status of a child. This is obvious in the coefficients of these factors in the univariate results and similar to the multivariate results. That is, the positive impact of educating a mother from no education to secondary education is equivalent to the impact of making the poorest family wealth index to the richest category. Also, the impact of having a good toilet system is equivalent to making the poorest family wealth index to the richer category. The impact on the provision of mass media for educative purposes would generate a higher impact compared to common diseases treatment supplements. These findings are essential for building a cost-effective intervention for short and long terms programs.

Conclusion

This study identified the risk factors of the concurrences of malnutrition among Ethiopian children under the age of five. These significant risk factors were identified to be the child's gender, mother's age, mothers' educational status, family wealth index, presence of diarrhea, fever, coughs, breastfeeding, and toilet type. Findings revealed that there is a spatial correlation in the concurrence of malnutrition indicators, and identified the Afar region to have a consistently higher prevalence of all malnutrition indicators compared to other regions. Based on a clinical trial (Humphrey et al. 2019), improved water, sanitation, and hygiene are not significant in reducing malnutrition, and therefore, focus should be directed towards parents and caregivers. In general, cost-effective childhood malnutrition intervention programs for parents and caregivers is a propitious method to mitigate malnutrition in high prevalent regions, which should be considered with utmost priority than direct methods on the children. For example, programs that enlighten women on the modern types of contraceptives, appropriate birth spacing, the benefits of antenatal and postnatal care, and educate mothers or caregivers on the negative impacts of their food preferences in child malnutrition. It is cost-effective, in the sense that, consistent enlightenment could be cheaply conducted through telecommunication programs. The generated maps pinpointed the locations with severe malnutrition that could be adopted to design a localized-based malnutrition intervention program, a randomized clinical trial, a localized framework for public health program evaluation, and an intervention causality effect. Also, the map could be overlaid with the spatial variability maps of another indicator, such as anemia or diarrhea prevalence to allow a better understanding of the spatial distribution of these indicators in Ethiopia for effective public health interventions. This study recommends nutrition and maternal education in the interior residence of the hotspot regions, supportive health care services, and poverty alleviation programs among less privileged households in the country.

One of the main limitations of this study is the cross-sectional nature of the data, which do not allow for causal inferences. Another limitation is that the developed model did not account for sampling weight nor adjust for sampling design on variance estimates. Also, the study did not include some regional aggregate factors which are directly relevant for under-five nutritional status, including food security and nutritional program coverage, which were not captured by the survey data. Despite this, the use of a spatial and temporal analytical approach has not provided sufficient contribution to the under-five malnutrition literature.

Future research can investigate changes in undernutrition patterns over time using point estimation at smaller spatial scales or a continuous geostatistical model in Ethiopia. Also, future work could account for the impact of the sampling scheme on the variance estimates, and extend the idea developed in this work to small area estimation techniques.

Funding 

No funding was received.

Acknowledgments

The authors would like to acknowledge the DHS program for granting access to the data used in this study.

Availability of data

The data utilised in this study is available publicly from the Demographic and Health Survey Program data base (EDHS) at https://www.dhsprogram.com/, reference ET_2016_DHS_02242021_1746_159735. The authors were given the permission to use the data for this study. A fraction of the cleaned and re-coded data, which is very particular to the developed model in this study, is available on http://dx.doi.org/10.17632/y33986sbyj.1.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Authors' contributions

OAE conceived the study idea, AMB sourced the data, OAE performed the computation, MAB validated the model, and MAB interpreted the results. All authors participated in the manuscript writing, and AMB revised it critically for intellectual content. All authors read and accepted the final manuscript revision. All authors agree to be accountable for all aspects of the work.

Appendices

Appendix 1.

Cluster-specific random intercept

Figure A1. Cluster-specific random intercept of stunting and underweight in Ethiopia. The vertical bars are the credible intervals.

Figure A2. Cluster-specific random intercept of stunting and wasting in Ethiopia. The vertical bars are the credible intervals.

Figure A3. Cluster-specific random intercept of underweight and wasting in Ethiopia. The vertical bars are the credible intervals.

Figure A4. The relative change in the prevalence of stunting with respect to a child's age (a) and mothers' age (b), underweight with respect to a child's age (c) and mothers' age (d) and wasting with respect to a child's age (e) and mothers' age (f).

Figure A5. Prevalence of (a) stunting, (d) underweight, (g) wasting, the significance of (b) stunting (e) underweight, (h) wasting, and the probability of high prevalence of (c) stunting, (f) underweight, and (i) wasting. The black portions indicates regions with significant high prevalence and the white portion indicates significant lower prevalence region.

Appendix 2.

Univariate analysis

Table C1. Fixed effects estimates of the univariate model.

	Stunting			Underweight			Wasting
		Credible interval			Credible interval			Credible interval
	mean	0.025quant	0.975quant	mean	0.025quant	0.975quant	mean	0.025quant	0.975quant
(Intercept)	−1.617	−45.497	42.227	−1.394	−45.177	42.354	−0.652	−44.456	43.116
Education
No education	0.000
Primary education	0.117*	0.019	0.216	0.109*	0.037	0.182	0.046	−0.031	0.123
Secondary	0.422*	0.223	0.620	0.246*	0.099	0.393	−0.002	−0.157	0.154
Higher	0.348*	0.073	0.623	0.304*	0.101	0.507	0.161	−0.054	0.376
Wealth index
Poorest	0.000
Poorer	−0.039	−0.152	0.075	−0.002	−0.085	0.082	0.027	−0.062	0.116
Milddle	0.243*	0.118	0.368	0.168*	0.076	0.260	0.027	−0.071	0.124
Rich	0.303*	0.167	0.440	0.246*	0.145	0.346	0.098	−0.008	0.205
Richest	0.362*	0.177	0.547	0.370*	0.233	0.507	0.238*	0.093	0.383
Source of water
Unimprove	0.000
Improve	−0.109*	−0.193	−0.025	−0.076*	−0.138	−0.014	−0.009	−0.075	0.057
Toilet type
Unimprove	0.000
Improve	0.342*	0.118	0.566	0.261*	0.096	0.426	0.088	−0.087	0.263
Settlement
Rural	0.000
Urban	−0.067	−0.232	0.098	−0.060	−0.181	0.062	−0.047	−0.176	0.081
Mass media
No	0.000
Yes	0.102*	0.006	0.199	0.120*	0.048	0.191	0.090*	0.015	0.166
Gender
Female	0.000
Male	−0.071	−0.145	0.002	−0.041	−0.096	0.013	−0.058*	−0.116	−0.001
No of children under 5
$<2$	0.000
$2\mathrm{\&}3$	−0.049	−0.169	0.071	−0.034	−0.123	0.054	−0.011	−0.105	0.082
$>3$	0.177	−0.082	0.436	0.043	−0.148	0.235	−0.100	−0.303	0.102
Birth order
$<2$	0.000
2-4	0.051	−43.829	43.894	0.037	−43.747	43.784	0.018	−43.786	43.785
$>5$	−0.054	−43.933	43.789	−0.052	−43.835	43.695	−0.024	−43.827	43.743
No of children
<2years	0.000
2-5 years	0.171*	0.078	0.264	0.127*	0.059	0.196	0.041	−0.032	0.113
$>5$years	0.397*	0.242	0.551	0.310*	0.196	0.424	0.129*	0.008	0.249
Antenatal and Postnatal care visit
$<2$	0.000
2-3	0.227*	0.110	0.344	0.095*	0.009	0.182	−0.056	−0.147	0.035
$>3$	0.788*	0.428	1.148	0.340*	0.074	0.605	−0.173	−0.454	0.107
Contraceptive
No	0.000
Traditional	−0.159	−0.725	0.406	0.069	−0.349	0.487	0.341	−0.101	0.783
Modern	0.113*	0.016	0.209	0.105*	0.034	0.175	0.065	−0.010	0.140
Presence of diarrhea
No	0.000
Yes	−0.071	−0.195	0.052	−0.122*	−0.213	−0.030	−0.125*	−0.222	−0.029
Presence of Fever
No	0.000
Yes	−0.078	−0.209	0.054	−0.134*	−0.232	−0.037	−0.128*	−0.231	−0.025
Presence of cough
No	0.000
Yes	−0.005	−0.128	0.118	0.003	−0.087	0.094	0.003	−0.093	0.099
Breast feeding
No	0.000
Yes	−0.290*	−0.382	−0.199	−0.202*	−0.270	−0.135	−0.057	−0.129	0.014

The significant factors are indicated by *.

Appendix 3.

R Code