SPATIAL CLUSTERING OF STUNTING CASES IN INDONESIA: A BAYESIAN APPROACH

,


INTRODUCTION
Stunting is a problem of chronic malnutrition in children which is influenced by various factors.
The height of a stunted child is shorter than the normal height of children of their age. Likewise, the brain of a stunted child will likely never develop to normal cognitive potential. Globally, approximately 144 million children under five are stunted. It is stated by the World Health Organization (WHO) that the average prevalence of stunting in Indonesia from 2005-2017 is 36.4% which exceeds the WHO threshold (30%) [1]. In 2017, the prevalence of stunted toddlers aged 0-59 months in Indonesia was 29.6%. Indonesian Nutritional Status Study (SSGI) in 2021 stated that the percentage of stunting in Indonesia is 24.4%. One indicator to determine whether a child is stunted or normal is to look at height for age. The number of toddlers 0-59 months whose height was measured in Indonesia in 2021 was 24958590 with the number of stunting cases being Some research on spatial modelling of stunting has been done [2][3][4][5][6][7], but the use of spatial analysis using the Bayesian method is still rare. The risk of stunting in South Sulawesi Province has been investigated by using the Bayesian spatial Conditional Autoregressive (BSCAR) Leroux model and found that the best model for modeling the relative risk (RR) cases of stunting under five in 2019 is the BSCAR Leroux model with hyperprior Inverse Gamma IG (0.5; 0.0005). They found that Toraja district and Pare-Pare have the highest RR (2.02) of stunting followed by Enrekang district (1.94). In contrast, Gowa district has the lowest RR (0.46) of stunting followed by Makassar City (0.52), and Pinrang district (0.67) [8]. Another research on modelling stunting cases used the BSCAR Leroux model in South Sulawesi Province in 2020 involving some covariates in the model [9]. They found that factors that significantly influence the incidence of stunting are the percentage of poverty and undernourished children aged 0-59 months. These both factors have a positive effect on stunting incidence. However, these studies have not implemented the BSCAR localised model in modelling the RR and focused only on one province in Indonesia. This study aims to provide the most suitable BSCAR localised (clustering) model and identify the RR of stunting in each province in Indonesia.

Data
Data on the number of toddlers 0-59 months whose height is measured and the number of stunted toddlers in each of the 34 provinces in Indonesia in 2021 were used. These data were gathered from Badan Pusat Statistik [10].

Models
A BSCAR localised model [11] was used to estimate the risk of stunting and examine the clusters of stunting cases. This model is similar to BSCAR Leroux [12] model but there is a clustering component ( ) in the BSCAR localised model which allow adjacency different over areas.
The number of stunting cases ( ) are modelled by using a Poisson log-linear model, commonly applied for disease mapping. This model is explained as follows: ~Poisson( ) for i = 1, 2, 3, …, 34 provinces and are the number of expected cases and the RR in the th area, respectively. β 0 is the overall level of RR. The spatial random effect ( ) is modelled using an intrinsic CAR prior as follows: is the spatial adjacency matrix based on queen contiguity. is defined by using the simplest form of weights (binary spatial weight) where = 1 if locations i and j are neighbours. Otherwise, where 0 = −∞ and +1 = +∞ A variable assigns the allocation of the th area to a cluster, 4 ANDI ASMAWATI AZIS, ASWI ASWI G is suggested to be a small and odd number [11]. We choose different BSCAR localised models with a different maximum number of clusters: two clusters (G=2), three clusters (G=3), and five clusters (G=5) for each hyperprior.
All models were fit using the CARBayes package version 5.3 [13] in R software version 4.2.0 [14]. To check the convergence of the MCMC sample, trace plots, as well as density plots, were used. MCMC samples were generated based on 106,000 iterations. The burn-in of 6,000 samples was discarded. As a result, 100,000 MCMC samples were collected. The goodness of fit of the model formulation was compared using some criteria: Watanabe Akaike Information Criterion (WAIC) [15], the Deviance Information Criterion (DIC) [16], and Modified Moran's I (MMI) [17,18] for the residuals. A model with a smaller value of DIC and WAIC fits the model well. Furthermore, a closer to zero of MMI residual signifies a better model fit.
Moran's I statistics on raw data were used to check the spatial autocorrelation and it is calculated as follows: is the number of areas (provinces), and are the observed data in the certain province i and province j, ̅ is the average of all the values over the provinces, is the spatial adjacency matrix. Moran's I statistics on residuals can also be used to detect model goodness of fit.
However, it is recommended to use MMI for a small number of areas. An MMI was expanded [17] for detecting spatial dependence which works even for a small number of areas. Research on MMI is available in some pieces of literature [17,18]. MMI statistic is calculated as follows:

Descriptive Analysis
The

Bayesian Spatial CAR localised Model
A BSCAR localised model was used with two (G =2), three (G=3), and five (G =5) clusters with different five hyperpriors. The DIC, WAIC, and MMI values for the residuals and the number of areas (provinces) included in each cluster are given in Table 1. Table 1 indicates that the BSCAR localised model with G=2 in each hyperprior has a lower DIC value than the model with G=3 and G=5 except for hyperprior IG(0.5, 0.0005) (M15).
Localised models with G=2 also have lower WAIC values compared to models with G=3 and G=5 for the five hyperpriors. The models with G=2 (M1=M4=M7=M10=M13) and G=3 (M2=M5=M8=M11=M14) have the same clustering structure for five different hyperpriors. Meanwhile, the model with G=5 has a different grouping structure for each hyperprior (M3=M6≠M9≠M12≠M15). The BSCAR localised model with G=3 with hyperprior IG(1, 0.1) (M5) has an MMI value for the residual that is the closest to zero (-0.32) and it is relatively similar to Model M10 (-0.41). Generally, based on the criteria used in this research, the preferred model for estimating the RR of stunting cases is BSCAR localised with hyperprior IG (0.5, 0.05) (M10) with two clusters (G=2).
The localised structure (LS), as well as RR values for each province for G=2 with hyperprior IG (0.5, 0.05) (M10), were given in   Figure 1 and  The RR maps of stunting cases based on the BSCAR localised model with G=2 (M10) are given in Figure 3. Out of ten provinces on the Sumatera island, only two provinces (Aceh, and Sumatera Barat) have a relative risk above the average (RR>1). The RR values of Aceh province (ID=1) and Sumatera Barat province (ID=31) are 1.42 and 1.28, respectively.

CONCLUSIONS
The results indicated that the BSCAR localised with hyperprior IG (0.5, 0.05) and IG (1, 0.1) are preferred for two and three clusters, respectively. Our results identified the high-risk areas for stunting. Approximately 56% of provinces in Indonesia are at a high risk of stunting. Sulawesi Barat has the highest RR for stunting followed by Nusa Tenggara Timur dan Papua Barat. In contrast, Jakarta has the lowest RR of stunting followed by Sulawesi Utara and Sumatera Selatan.
Government should pay more attention to areas that are most at high risk of stunting. Including some covariates in the model could be possible for future work.