Application of single-level and multi-level modeling approach to examine geographic and socioeconomic variation in underweight, overweight and obesity in Nepal: findings from NDHS 2016

Nepal’s dual burden of undernutrition and over nutrition warrants further exploration of the population level differences in nutritional status. The study aimed to explore, for the first time in Nepal, potential geographic and socioeconomic variation in underweight and overweight and/or obesity prevalence in the country, adjusted for cluster and sample weight. Data came from 14,937 participants, including 6,172 men and 8,765 women, 15 years or older who participated in the 2016 Nepal Demography and Health Survey (NDHS). Single-level and multilevel multi-nominal logistic regression models and Lorenz curves were used to explore the inequalities in weight status. Urban residents had higher odds of being overweight and/or obese (OR: 1.89, 95% CI: 1.62–2.20) and lower odds of being underweight (OR: 0.81, 95% CI: 0.70–0.93) than rural residents. Participants from Provinces 2, and 7 were less likely to be overweight/obese and more likely to be underweight (referent: province-1). Participants from higher wealth quintile households were associated with higher odds of being overweight and/or obese (P-trend < 0.001) and lower odds of being underweight (P-trend < 0.001). Urban females at the highest wealth quintile were more vulnerable to overweight and/or obesity as 49% of them were overweight and/or obese and nearly 39% at the lowest wealth quintile were underweight.


Prevalence of underweight and overweight and/or obesity by wealth quintiles and geography.
The heat map (Fig. 1, Supplementary Fig. 1) for the prevalence of overweight and/or obesity and underweight by age categories and wealth quintiles suggests that middle-aged adults (35-54 years) in the highest wealth quintile (quintile 5) are more vulnerable to overweight and/or obesity; the prevalence for the age category in the given wealth quintile was greater than 50% based on both global classification of BMI (BMI >24.9 kg/m2) and the classification for Asians (BMI >22.9 kg/m2). The heat map for undernutrition shows that older adults, aged 65 and above, in the lower-and mid-wealth quintiles (quintiles 1, 2, 3) were vulnerable to underweight, where more than one-third of the older adults in the given quintiles were underweight (Fig. 1, Supplementary Fig. 1).
The geospatial analyses in Fig. 2 show the prevalence of overweight and/or obesity and underweight by wealth quintiles. The prevalence of underweight was higher among the residents of Province-2 who belonged to lower wealth quintiles (quintile 1, 2, and 3). Overweight and/or obesity was more prevalent (>40% by global cutoff and >50% by Asian cutoff) among the residents of Province-3, 4, and 6 who belonged to higher wealth quintiles (quintile 4 and 5) (Fig. 2).
The geospatial analyses presented in Fig. 3 shows the prevalence of overweight and/or obesity and underweight by gender and residence. Undernutrition, particularly among females and rural residents, is comparatively higher in Province-2; prevalence ranged between 30-40%. Overweight and/or obesity was higher (≥30%) in Provinces-3 and 4, specifically among females and urban residents (Fig. 3).
Single-level multinomial regression. In the single-level multinomial analyses (Tables 2 and 3), to assess the association between participants' weight and their socio-economic characteristics, participant's socio-economic and locational characteristics were associated with their both overweight and/or obesity and underweight status (Tables 2 and 3). After adjusting for age, sex, education, marital status, wealth quintiles, residency, ecological region, and provinces, compared to participants with no or only preschool education, those with any level of formal education had higher odds of being overweight and/or obese, and lower odds of being underweight ( Table 2). Compared to rural residents, urban residents had higher odds of being overweight and/or obese (OR:1.89, 95% CI: 1.62-2.20) or lower odds of being underweight (OR:0.81, 95% CI: 0.70-0.93). Provinces wise, compared to the Province 1, participants from the Province 2, and 7 were less likely to be overweight and/ or obese and more likely to be underweight. Ecologically, residents of Terai region had lower odds of being overweight and/or obese (OR: 0.81, 95% CI: 0.67-0.98) or higher odds of being underweight (OR: 1.74, 95% CI: 1.33-2.27) ( Table 2).
Participants with higher wealth quintiles were associated with higher odds of being overweight and/or obese and lower odds of being underweight (Table 3).
Multinomial regression. The null model, containing no explanatory variables, i.e., only the outcome and a constant, was used to illustrate the total variance in overweight and/or obesity associated with individual and locational (and district) characteristics. The intercept and the Intra-class Correlation Coefficient (ICC) for provinces and districts were significantly different than zero in all the models, suggesting that the log-odds varied significantly between the provinces and the districts. For the null model, the ICC for level 2 and level 3 was 6.8 and 7.3, respectively. This means that about 15% of the variance in overweight and/or obesity is explained by population differences (6.8% due to the systematic differences between provinces and 7.3% due to districts) and more than 85% of the variance is attributable to individual differences. In the fully adjusted model 4, the ICC reduced to 4.1% intra-provincial correlation and 2.9% intra-district correlation suggesting that only a small variation in overweight and/or obesity was explained by random differences between the provinces and the districts. The value of Akaike's information criteria (AIC) and Schwarz Bayesian information criteria (SBIC) were reduced with the addition of covariates across the four models suggesting that the subsequent model improved over the previous model. The model 4 was the best fitting model as suggested by the smallest value of AIC and SBIC compared to the preceding model (Table 4). Compared to the null model, the proportional change in the variance of overweight and/or obesity, due to addition of age, sex, wealth quintiles, education, marital status, and residency in the model 4, was 39.7% across provinces and 60.3% across districts which indicated that 39.7% of the provincial variance and 60.3% of the variance by districts in the empty model was attributable to differences in age, sex, wealth quintiles, education, marital status, and residency.
In the best fitting multilevel model, older ages compared to the youngest (15-25 years), female sex, any level of school education compared to illiterate or with only preschool education, residing in urban areas, and higher wealth quintiles compared to the first quintile was associated with higher odds of being overweight and/or obese. Never married participants were less likely to be overweight and/or obese than those married (OR: 0.34, 95% CI: 0.27-0.42) ( Table 4).
In the multilevel analyses for underweight, the intra-provincial correlation and intra-district correlation for underweight was 3.5% and 2.1%, respectively which slightly increased to 4.6% and 3.8%, respectively in the final model (model 4). Thus, only a small variation in underweight was explained by random differences between the provinces and the districts, and the individual differences explained most of the variance. As suggested by the smaller value of AIC and SBIC, model 4 was the best fitting model. Compared to the null model, there was a decrease in the proportion of provincial (31.4%) and district level (81.0%) variance with the addition of individual-level predictors in the final model.
The best fitting multilevel model for underweight suggests that compared to the youngest age group (15-25 years), adults in the middle ages (25-55 years) had lower odds of being underweight but the older adults (>65 years) had higher odds of being underweight (OR:1.85, 95% CI: 1.50-2.29). Compared to illiterate or with only preschool education, participants with any level of school education were less likely to be underweight. Being married appeared to be protective against underweight compared to being unmarried (OR:2.30, 95% CI:1.99-2.67), and previously married (OR:1.20, 95% CI:1.02-1.42) participants were more likely to be underweight than their married counterparts. Urban residents were less likely to be underweight compared to rural (OR:0.87, 95% CI: 0.76-0.99) but the statistical significance was borderline. The financial privilege was protective against underweight, as participants in the higher wealth quintiles had lower odds of being underweight compared to the first quintile (Table 5).
Socio-economic inequalities in weight status: Lorenz curves analyses. Analyses using the Lorenz curves ( Fig. 4, Supplementary Fig. 2), used for understanding inequality between richer and poorer and disaggregated by sex and residency, further confirmed the earlier findings in multivariable regression that showed a positive association between weight status and wealth quintiles. The Lorenz curves for the prevalence of underweight ( Fig. 4A) shows that underweight status was most heavily concentrated among poor males and females in the urban area and among rich males and females in the rural area. The concentration curve graph shows that, www.nature.com/scientificreports www.nature.com/scientificreports/ urban female at the bottom 20% were more vulnerable as nearly 39% were underweight and at the top 20% nearly 19% were underweight (Fig. 4A).
Overweight and/or obesity prevalence was concentrated among the rich as shown by all the curves below the line of equality. Particularly, urban females at the highest wealth quintile were more vulnerable as 49% were overweight and/or obese (Fig. 4B). This is not very different from the overweight and/or obesity gradient based on Asian BMI cut off ( Supplementary Fig. 2).

Discussion
Summary of evidence. This study reports the prevalence of underweight and overweight and/or obesity, estimated by measures of height and weight, in a nationally representative sample of Nepalese adults. The findings of this study provide evidence of the dual burden of underweight and overweight and/or obesity, demonstrating that this is an urgent priority for the health of Nepal. Our study demonstrates that the prevalence of overweight and/or obesity was higher in younger adults and undernutrition was higher among the older adults. Based on wealth indicators, individuals from richer households were least likely to be underweight and more likely to overweight and/or obese. The multilevel analyses suggest that only a small variation in overweight and/or obesity and underweight can be explained by random differences between the provinces and the districts, and the individual differences explained most of the variance.
Comparison to existing studies. Though the prevalence and factors affecting underweight and overweight and/or obesity from NDHS 2016 were reported in earlier studies 13,14 , none of these studies examined the geographical and socio-economic variation of weight status with a detailed spatial mapping of BMI status. And, both studies have some methodological limitations. For example, Al Kibria 13 used the same weight for the underweight vs normal weight comparison, and overweight and/or obesity vs normal weight comparisons, though they were analysed separately using multivariable logistic regression analysis. Though, a detailed account of modeling process was lacking in these studies 13,14 . These shortcomings have been addressed in our study. We used a hierarchical multi-nominal regression analysis design where normal weight used as referent category and compared against underweight and overweight and/or obesity under the same model. Further, we used the multilevel analysis to examine the population level differences in BMI status, and Lorenz curves to study socio-economic inequalities in BMI status. Both multi-nominal and multilevel analysis explored the same set of predictor variables and the findings have high congruence between two modeling approaches. Multilevel analysis is accounted for the presence of data hierarchies and the random residual components at each level in the study 17 . It is superior to the traditional logistic regression modeling 17 , however, we did not aim to compare the results between single level Dual burden of underweight and overweight and/or obesity. A notable proportion of participants were found to be underweight. The finding was not surprising given that historically undernutrition has been a significant public health problem in Nepal 18 . Underweight status appeared to be a major concern among the older adults in our study which is in line with previous evidence that the nutritional well-being of older adults is particularly neglected in Nepal 19,20 .
On the other hand, recently overweight and/or obesity has emerged as a notable public health concern. It is particularly of concern among young adults, given that being obese increases the vulnerability to sequelae of cardio-metabolic diseases. The national STEPS survey on non-communicable diseases risk factors, conducted in 2013, also provides evidence that Nepalese adults are becoming overweight and/or obese (obesity: 21%) 9 . Two aspects of Nepalese society, the traditional norms of being obese and rapid urbanization, may be attributed to the rise in obesity. Being overweight is culturally desirable in Nepalese society and is even used as a proxy of family's ability to afford food as well as a measure of their relative wealth status 21 . This traditional phenomenon has been further complicated by the rapid urbanization which has increased the accessibility and affordability of energy-rich, sugary commercialized foods, and simultaneously, a more sedentary lifestyle 22,23,24 . Socio economic status (SES) and nutritional status. The relationship between weight status and SES was expected, yet interesting. We found a high prevalence of overweight and/or obesity among wealthier individuals, measured in terms of wealth index, and a high prevalence of underweight among poorer individuals. This is not unexpected, and consistent with other data from Nepal and globally 15,25 . The poorest quintile were two times more likely to be underweight than the richest quintile, whereas the richest quintile were 7.2 times more likely to be overweight and/or obese than the poorest quintile. However, the disparity in underweight across socio-economic groups is much more evident within urban areas compared to rural areas although there was no statistically significant association of residence with underweight. This large disparity between rich and poor subgroups in Nepal might be associated with ongoing rapid urbanization that will further widen income and social inequities. These findings corroborate to those reported by previous studies conducted in other LMICs 5,15,26 .
A higher burden of underweight status among the lower socio economic group in urban regions could be attributable to rapid growth of cities along with growth of urban poor 27 . Globally, South Asia has the largest proportion of the population living in slums and is urbanizing at the fastest rate 28 . The people living in these www.nature.com/scientificreports www.nature.com/scientificreports/ slums are often beyond the remit of public services and are forced to endure disparities in essential health services and health and nutrition indicators 29 . Policymakers cannot ignore the sheer scale and uncontrollable urban growth and as more and more of the populations in LMICs become urbanized, so the need for special attention for this population subgroup is growing 30 . In addition, the urbanization also incites inflation which makes consumption of fruits and vegetables expensive and increasingly difficult for urban poor. Thus the consumption of www.nature.com/scientificreports www.nature.com/scientificreports/ unhealthy foods along with the transition to sedentary occupations from subsistence agriculture might be the reasons behind the shift in the burden of obesity to lower-SES groups in LMICs 31 . With rapid urbanization, the urban poor subgroup in Nepal is also at risk for double burden of undernutrition and overweight and/or obesity.

Geographic drivers.
In the geospatial analyses, underweight status was more prevalent among the residents of Province 2, Province 6 and Province 7 whereas overweight and/or obesity was more prevalent among the residents of Province 3 and 4. Overweight and/or obesity is more prevalent in provinces and districts with higher affluence and underweight in less affluent provinces. For example, Province 1 and 3, with human development index (HDI) of 0.507 and 0.506, are the most developed provinces in terms of HDI as well as other indicators of economic growth, followed by Province 4 (HDI: 0.493). At the lower end, Province 6 (HDI: 0.39), Province 7 (HDI: 0.416), and Province 2 (HDI: 0.422) have the lowest HDI, with the majority of the population residing in the village and semi-urban area and engaged in agrarian activities as a primary source of income. When comparing the districts, the top three districts in terms of HDI are from Province-3: Kathmandu (HDI: 0.632), Lalitpur (0.601), Kaski (0.576), and the bottom three are from Province 7: Bajura (0.364), Bajhang (0.365), and Kalikot (0.374) 32 . At the individual level-underweight showed a negative trend in dose-response manner with household SES and overweight showed a positive trend. Low SES groups from low HDI districts and high SES groups from high HDI districts are particularly vulnerable to poor nutritional status and needs vigilant attention. In line with our findings, geographic disparitities in overweight and/or obeseity has been reported in previous studies from Nigeria and neighbouring India where overweight and/or obesity was more prevalent in Southern India 33 and south-eastern Nigeria 34 compared to other regions of these countries.
Nepal is in a transitional phase and has recently transformed into a federal setup with seven administrative provinces. The several decades of political instability and ten years of armed conflict that ended in 2006, plunged Nepal into poverty. The situation is particularly grave for the Province 2, Province 6 and Province 7 with approximately half of their population being poor 35 . Though Province 2, the plains of Nepal has been the major region for crop production, food insecurity in this province is on the rise due to annual floods, crop failures and feudalistic distribution of land in agriculture. The difficult terrain in Province 6 and Province 7, along with poor crop harvest and poor income from agriculture and livestock has been the major reasons behind highest poverty incidence in these provinces. Other reasons are out-migration of youth into Persian Gulf countries and in-migration of richer households from rural to urban areas leaving the poorer behind as they are less likely to migrate because of lack of affordability 36 . Along with that, these provinces were the worst affected region during armed conflict. This is also reflected in higher underweight prevalence in these provinces and remains a challenge for Nepal.

Strengths and limitations. Our study adds to the growing literature surrounding nutritional status in
Nepal. To date, nationally representative data are limited in Nepal, and the prevalence, estimated in the earlier studies, are based on small sample sizes and among a homogenized population in a limited geography. This study is based on a large nationally representative sample consisting of both urban and rural populations in Nepal. This is the first study exploring the socio-economic and geographical disparities in the prevalence of underweight and overweight and/or obesity in the Nepalese context. Further, we defined overweight and/or obesity based on two different criteria, the global criteria and the criteria recommended for Asians. Our geospatial analyses may be helpful to the new provincial governments, especially in the context of recent federalization of the country, to design the targeted interventions based on local needs.
Our study also has several limitations. This study is cross-sectional in nature, and thus limits any causal inferences. Likewise, BMI and SES were measured only at one point, but such measurements are likely to change over time. Given the methods of our data sample, we were unable to assess the impact of such changes at different stages of the life course. Additionally, the lack of information regarding dietary habits, alcohol intake, or physical activity in the NDHS database hindered our ability to scrutinize some important determinants of nutritional status and the possibility of residual confounding cannot be ruled out. The geospatial mapping was entirely limited to data-visualization, and further studies needs to be carried out to identify the area-level hotspots of endemic underweight and overweight and/or obesity in Nepal. Despite the limitations, our study provides new insights into socio-economic and geographic determinants of weight status in Nepal.

conclusion
This study reaffirmed the earlier evidence that the dual burden of underweight and overweight and/or obesity is a significant public health problem in Nepal. The BMI status varied by indicators of socio-economic well-being which further confirms the role of the household's economic well-being in the current epidemiological transition. Overweight and/or obesity among the younger adults and undernutrition among the elderly should receive special attention.

Methods
Data sources. This study is based on a secondary analysis of the 2016 NDHS, a cross-sectional and nationally representative data. The 2016 NDHS uses multi-stage (two stages in rural areas and three stages in urban areas) stratified cluster sampling technique. A total of 383 primary sampling units (PSU) or clusters were selected nationally, and 30 households were selected from each of these PSU with an equal probability proportion to size criterion yielding a final sample size of 11,040 households. From these households, 14,937 individuals, including 6,172 males and 8,765 females, provided height and weight measurements. Data includes individuals from all the 77 districts in the country, across all seven provinces. The full report on the 2016 NDHS study design and sampling technique is available elsewhere 37 www.nature.com/scientificreports www.nature.com/scientificreports/ Anthropometry measurements. The trained field research staffs measured the height and weight of the participants. BMI was calculated as weight in kilograms divided by the square of the height in meters (kg/m 2 ). We followed two standard classifications of BMI. First, we used the global standard to categories BMI as underweight (<18.5 kg/m 2 ), normal (18.5-24.9 kg/m 2 ), overweight (25-29.9 kg/m 2 ), and obese (>30 kg/m 2 ); the latter two categories were combined to form a single category overweight and/or obese (≥25 kg/m 2 ) 39 . Further, to provide a comparative context in the figures (heat map, geospatial maps, and Lorenz curves), we also calculated and presented the prevalence of overweight and/or obese (≥23 kg/m 2 ) following the recommendation by WHO's expert consultation on BMI classifications for Asians 40 .
Explanatory variables. The explanatory variables that were selected for analysis can be grouped into socio-demographic, and geographic variables. We considered age of the participants, sex, educational status and wealth index under socio-demographic variables. Educational status was grouped into four categories: no formal schooling, primary schooling, secondary schooling and those with higher education (includes high school, college/university, postgraduate and above). The principal component analysis was used to determine wealth index which included information on number and kinds of consumer goods such as bicycle or car owned by household and housing characteristics (such as the source of drinking water, availability of toilet facilities 41 . Place of residence (rural/urban), Province and Ecological zone (Mountain/Hill/Terai) were included under geographic variables. The province was categorized into seven administrative divisions, according to the current administrative structure of Nepal, and compared against Province-1. Province 1 was chosen a reference category as it is most advantaged in political and economic means Nepal 42 . Wealth quintiles were computed by dividing the distribution into five equal categories, each comprising 20% of the population. Data analysis. Multivariable analysis. Analyses were adjusted for cluster and sample weight to ensure representativeness of provinces and to the urban and rural areas. Rao-Scott chi-square tests were used to assess the association between weight status and the explanatory variables in the bivariate analyses. Univariate multinomial logistic regression was performed to assess the association of weight status with various socio-demographic factors, as listed in Table 1, with normal weight as the reference category. We built three hierarchical models under multinomial analyses. Model 1 was adjusted for age and sex; model 2 was adjusted for variables from the model 1 plus education (E), marital status (M), wealth quintiles (SES), and residency (R). Model 3 was adjusted by adding ecological region (Ec) and provinces (P) to the model 2. Model 4 contained all the variables from the model 3 excluding the wealth quintiles. The fully adjusted model is summarized as:  www.nature.com/scientificreports www.nature.com/scientificreports/ Multilevel analysis. The hierarchical nature of NDHS data allowed us to perform a hierarchical generalized linear modeling using PROC GLIMMIX procedure in SAS 9.4 45 . Such dataset is ideal for exploring socio-economic and geographic differences in BMI status at individual and population level. We performed three-level analyses with individuals at level 1 (a), districts at level 2 (b), and provinces at level 3 (c), using the same set of predictor variable used in multivariable analysis. In the multi-level analyses, the model 1 is an empty model. The detailed model building process for three-level multilevel analysis is described by Kim et al. 46 using BMI as a continuous outcome variable, and is adapted for categorical BMI as a categorical variable. The study found that the residuals are independently and identically distributed at the individual, and population level in NDHS 2016 46    www.nature.com/scientificreports www.nature.com/scientificreports/ The intraclass correlation coefficient (ICC) (a measure of the amount of variation in outcome due to a given level) is calculated using the formula below where τ 00 is the covariance parameter estimate. The proportional change in variance (a measure of change in the area level variance between the empty model and a given successive model) was calculated based on the earlier established procedure in SAS v9.4 45 . The value of AIC and SBIC were used as model fit statistics.
Geographic analysis and Lorenz curves. The geospatial analyses of provinces of Nepal were performed in SAS v9.4 using GMAP procedure 47 . The Lorenz curve assessed the socio-economic inequality in underweight and overweight and/or obesity among rural and urban residents. This curve has been used to assess the inequality gradient related to socio-economic status (SES) in various health indicators [48][49][50] . We plotted the cumulative percentage of the sample, ranked by SES on X-axis and the corresponding cumulative percentage of weight status, by sex and urban-rural residence, was plotted on the Y-axis. A curve above the line of equality indicates a greater concentration of the outcome among the poor and a curve below the line indicates a greater concentration of outcome among the rich.
Ethical approval and consent to participate. The 2016 Nepal Demographic and Health Survey protocol received ethical approval by the Nepal Health Research (NHRC) -Ethical Review Board and ICF Institutional Review Board. The survey was conducted in accordance with relevant guidelines/regulations. An informed written consent was obtained from all the participants.