The patterns of Non-communicable disease Multimorbidity in Iran: A Multilevel Analysis

The prevalence of non-communicable diseases is increasing worldwide. Multimorbidity and long-term medical conditions is common among these patients. This study aimed to investigate the patterns of non-communicable disease multimorbidity and their risk factors at the individual and aggregated level. Data was inquired from the nationwide survey performed in 2011, according to the WHO stepwise approach on NCD risk factors. A latent class analysis on multimorbidity components (11 chronic diseases) was performed and the association of some individual and aggregated risk factors (urbanization) with the latent subclasses was accessed using multilevel multinomial logistic regression. Latent class analysis revealed four distinct subclasses of multimorbidity among the Iranian population (10069 participants). Musculoskeletal diseases and asthma classes were seen in both genders. In males, the odds of membership in the diabetes class was 41% less by increasing physical activity; but with increased BMI, the odds of membership in the diabetes class was 1.90 times higher. Tobacco smoking increased the odds of membership in the musculoskeletal diseases class, 1.37 and 2.30 times for males and females, respectively. Increased BMI and low education increased the chances of females’ membership in all subclasses of multimorbidity. At the province level, with increase in urbanization, the odds of membership in the diabetes class was 1.28 times higher among males (P = 0.027). Increased age, higher BMI, tobacco smoking and low education are the most important risk factors associated with NCD multimorbidity among Iranians. Interventions and policies should be implemented to control these risk factors.

Statistical analysis. The analytic strategy in the current study was in two steps. First, we used a latent class analysis to identify homogeneous subgroups based on multimorbidity with NCDs (mentioned above), which have a common pattern of co-morbidity within a heterogeneous population. Second, we examined the relations between some risk factors and class memberships using multilevel multinomial logistic regression. The goal of this second step of the analysis was to determine which individual-or province level characteristic was associated with multimorbidity latent class membership among the nationwide population, while controlling for the nested nature of the data (population within different provinces).
LCA was fit as an exploratory and iterative process with an increasing number of classes. The final number of classes was determined by considering model fitness indices, empirical evidence and interpretability.
Model fitness was assessed using the Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), Sample Size Adjusted Bayesian Information Criterion (SSABIC), the Lo-Mendell-Rubin likelihood ratio (LMR LR) test, Vuong-Lo-Mendell-Rubin likelihood ratio test (VLMR), and bootstrap likelihood ratio test (BLRT). Smaller values of AIC, BIC, and SSA BIC indicated a better model fit. The LMR is a likelihood ratio test and compares the estimated model to a model with one less class (k-1). If the LMR is not significant, model fitness will not improve by including an additional class into the model 30 . The VLMR and BLRT P-values are interpreted in the same way as the LMR LR test. Entropy was used to assess the quality of member classification; a value closer to 1.0 indicates better classification 31 .
Two-level multinomial logistic regression analysis with random intercept was conducted to explore the association of both individual and province-level predictors with non-communicable disease multimorbidity, in males and females, separately. The "low-risk of disease" category was used as the reference category. Each class was compared with the low-risk class in the same gender.
Initially an empty model with no predictors was built, and then individual-level and province-level variables were added. The variance components (random effects) of outcomes in both models for males and females were estimated as the province-level variance.
Data preparation was done in Stata14; latent class analysis and regression models were fit in Mplus 7.4. In this study the full-information maximum likelihood (FIML) estimation was used, which is recommended for handling missing data under missing at random (MAR) or missing completely at random (MCAR) situations. This method uses the robust maximum likelihood estimator (MLR) which uses model-based methods 32 .
Initially, an empty model with no predictors was fitted to examine the variation of the NCD multimorbidity classes across provinces, for males and females; separately. The variance component and the standard error of provinces in the empty model for males was 0.083 and 0.035 respectively; and the ICC was 2.4%. These values were 0.021, and 0.035 and the ICC was 0.6% in females. These low ICCs in both male and females show that the variation of the NCD multimorbidity classes at province-level is low.
Although, authors suggest that ICCs smaller than 0.05 may benefit less from multi-level analysis, Muthén states that regardless of the ICC, design effects greater than 2 need to take into account multilevel-analysis 33 . In the current study, the design effect was 4.36 for females and 14.44 for males according to this equation: 1 + (average cluster size -1)*intraclass correlation Thus, a two-level multinomial logistic regression analysis (random intercept) was conducted to explore the concurrent association of individual and province-level predictors with NCDs, for males and females, separately.
In the first part of the study which aimed at extracting multimorbidity patterns in the Iranian population, there was no missing data. In the second part, which examined the association between multi-morbidity patterns; 20% of the physical activity variable and 28% of the fruit and vegetable consumption variable were missing.

Results
participant characteristics. The descriptive statistics of all variables are shown in Table 1.

Multimorbidity patterns.
A series of latent class models were ran in which the number of classes ranged from 2 to 6 for males and 2 to 7 for females. Table 2 shows the goodness-of-fit indices for each model. Although the VLMR and LMR suggested 5-class models for males and females, we selected the 4-class model for both males and females, because it had a smaller BIC, SSABIC and a high entropy value.
The response probabilities for each of the classes in males and females are displayed in Table 3.
In males the first class (labeled "Diabetes"), constituted of 6.5% of the male sample, and was characterized by a high probability of diabetes (0.54) and a low probability of other diseases. This group also had the highest probability of Coronary Artery Disease (0.20). Class two (2.2% of the male sample, labeled "Asthma and Wheezing") was characterized by a high probability of Asthma (0.55) and Wheezing (0.73). The third class (labeled "Musculoskeletal diseases"), comprised 9.2% of the male sample and was characterized by a high probability of Back Pain (0.88), Neck Pain (0.51), and Knee Pain (0.81) and a low probability of other disease. This group also had the highest probability of Osteoporosis (0.16). The fourth class (labeled "Low risk of diseases"), comprised 82.1% of the male sample and was characterized by a low probability of all diseases.
For females the first class (labeled "Asthma and Wheezing"), constituted 3.8% of the female sample, and was characterized by a high probability of Asthma (0.58) and Wheezing (0.56). Class two (25.0% of the female sample, labeled "pre-skeletal diseases") was characterized by a moderate probability of Osteoporosis (0.25), Back Pain (0.38), Neck Pain (0.17), and Knee Pain (0.59). The third class (labeled "Musculoskeletal diseases"), which comprised 11.7% of the female sample was characterized by a high probability of Back Pain (0.99), Neck Pain (0.83), and Knee Pain (0.96) and a low probability of other disease except Osteoporosis. This group also had the highest probability of Osteoporosis (0.49). The fourth class (labeled "Low risk of diseases"), comprised 59.6% of the female sample and was characterized by a low probability of all diseases. The probability of Knee Pain was high in class 1 to 3. Colorectal Cancer and Cerebrovascular disease had a low probability in all male and female classes. (2020) 10:3034 | https://doi.org/10.1038/s41598-020-59668-y www.nature.com/scientificreports www.nature.com/scientificreports/ Multimorbidity and multilevel analysis. The variance component and the standard error of provinces in the empty model for males was 0.083 and 0.035 respectively; and the ICC was 2.4%. These values were 0.021, and 0.035 and the ICC was 0.6% in females. These low ICCs in both male and females show that the variation of the NCD multimorbidity classes at province-level is low. In the current study, the design effect was 4.36 for females and 14.44 for males. Table 4 shows the association between the predictors and the NCD multimorbidity classes in males. At individual-level, the odds of membership in all classes (class 1 to class 3) compared to the low-risk of disease class were higher as age increased. The hypertensive individuals had a higher odds for membership in class 1 (Diabetes) compared to the low-risk class (OR = 1.59; 95% CI: 1.15-2.21). The odds of membership in the diabetes class compared to the low-risk class was 41% less by increasing physical activity from one level to the other (levels: low, moderate, and high). With one level increase in BMI (levels: normal, over-weight, and obese), the odds  www.nature.com/scientificreports www.nature.com/scientificreports/ of membership in the diabetes class was 1.90 (95% CI: 1.45-2.49) times higher. Comparing the tobacco-using individuals to the non-using individuals, the odds of membership in the musculoskeletal diseases class versus the low-risk class was 1.37 times higher.
Illiterate individuals compared with those with academic education had a higher odd for membership in the musculoskeletal diseases class (OR = 1.96; 95% CI: 1. 16-3.31). Also, individuals with elementary-level education had higher odds for membership in the diabetes class compared with individuals who had academic education (OR = 2.00; 95% CI: 1.31-3.05). Students and soldiers compared with governmental officers had a 58% less odds  www.nature.com/scientificreports www.nature.com/scientificreports/ for membership in the skeletal-diseases class versus the low-risk class. At province-level, with one unit increase in the urbanization score, the odds of membership in the diabetes class were 1.28 times higher. Table 5 shows the associations between the predictors and the NCD multimorbidity classes in females. At individual-level, the odds of membership in all classes (1 to 3) were higher as age increased. Hypertensive individuals had a lower odds for membership in class 2 (OR = 0.77; 95% CI: 0.61-0.98). With one level increase in BMI (levels: normal, over-weight, and obese), the odds of membership in all classes versus the low-risk class increased. Among tobacco-using individuals, the odds of membership in the Musculoskeletal diseases class was 2.30 (95% CI; 1.35-3.93) times higher. Illiterate individuals had higher odds for membership in all classes compared with those who had academic education. Also, the individuals with elementary education had higher odds for membership in class 1 and 2 compared with those who had academic education. But the individuals with a high school diploma had higher odds for membership only in class 2 compared with those who had academic education. Retirees or un-employed people had a 54% less odds for membership in the musculoskeletal diseases class compared to governmental officers. Private employees had very low odds for membership in the asthma and wheezing class compared to governmental employees. At province-level, urbanization score was not associated with NCD multimorbidity in females.
The variance in the non-communicable multimorbidity classes in both males and females at the province-level did not increase significantly after including individual-level and province-level in the model.
Although the Charlson Comorbidity Index has also been used in similar research, but it was not appropriate for this study, because it includes 17 morbidities and this study did not have data for all of them. On the other hand, the LCA approach, takes into account response measurement errors and provides several goodness of fit indexes, and therefore provides more realistic co-morbidity patterns for the same population, at the same time, which is an advantage. Results are more reliable in LCA 22 .
The world population is aging and policymakers will have to face multimorbidity as a global challenge in the near future, especially in developing countries. This indicates that health care systems need to deliver early preventive and treatment programs 7 . NCDs are one of the biggest challenges for health care providers and policymakers in Eastern Mediterranean countries with a low socio-economic status, or involved in economic, demographic, geographic, or epidemiological transition. Longer-lived individuals are more likely to acquire multimorbidity, because of the increasing trend of chronic diseases in this area 4 . Therefore, studies should be conducted for the promotion of healthy lifestyles to prevent or manage multimorbidity 34
This study provides the first nationally representative estimate of NCDs multimorbidity patterns among Iranian adults. Other studies on multimorbidity conducted in other countries have included only young or middle-aged adults 35,36 , but our study included people aged from 20 to 70. In this study, the prevalence of most risk factors, such as hypertension, low physical activity, obesity, and low fruit and vegetable intake were higher in females than males.
In Ghana, India, Mexico, South Africa and the Russian Federation, the prevalence of hypertension was higher in females as well 37 . Whereas, a systematic review about blood pressure in several world countries showed that males had a higher average blood pressure than females in all world regions, expect West Africa 38 .
Higher life expectancy among females and more exposure to NCD risk factors can be one of the reasons for the higher prevalence of NCDs among females 13,39 .
In this study, the prevalence of most NCDs such as Diabetes, Coronary Artery Disease, Osteoporosis, Back Pain, Neck Pain and Knee Pain was significantly higher in females. Multimorbidity was more common among females in some previous studies 40,41 , but not is studies from Finland, Poland, and Spain 42 . A study from Kuwait also showed that the prevalence of NCDs was higher among males 43 . A potential cause for the higher prevalence of NCDs among females could be their relatively higher tendency to mention their condition in self-reports 41 and another reason may be gender inequality in access to healthcare 44 .
In our analysis, LCA revealed four distinct subclasses of multimorbidity among the Iranian population, for males and for females. Overall, it seemed liked heterogeneity and class separation was better performed in males than females, although both genders had the same number of latent classes.
Apparently in Iran, females are more prone to develop skeletal disorders and osteoporosis, and the reason might be vitamin D deficiency and low sun exposure. Also, the reason might be that females are more susceptible, report their pain more often or are more inclined to somaticize their feelings 45,46 . In a study conducted over a 28-year period in Iran and 6 other countries, low back pain was the leading cause of YLDs (Year Lost due to Disease) in both sexes 47 . The 2016 GBD reported that the burden of disease associated with musculoskeletal health increased 19.6%, from 2006 to 2016 48 . During 1990-2013, the total DALYs of musculoskeletal disorders increased by 105.2% in the EMR, although it increased only 58.0% in the rest of the world; and low back pain and neck pain had the highest burden in EMR countries 49 . Our findings also confirm that prevention and control programmers for musculoskeletal disorders should be implemented in different levels of national health care programs.
In this study, in both genders, the chance of membership in subclasses 1 through 3 increased with age. It is evident that chronic conditions and multimorbidity typically increase with age 50 . Authors state that age is related to the duration of exposure to risk factors 37 . However, a study reported that in South Africa, multimorbidity prevalence did not increase with age and a gradual decrease after 60 years was observed 51 . The reason for this pattern might be shorter life time or mortality from communicable diseases in Africa.  Table 5. Odds Ratios and 95% Confidence Intervals of the Association Between Non-communicable Multimorbidity Classes and Predictors in Individual-and provinces-Level Using Two-Level Multinomial Logistic Regression Analysis for females. *P-value < 0.05, **P-value < 0.01, ***P < 0.001.

Individual-Level
www.nature.com/scientificreports www.nature.com/scientificreports/ In this study, in males, latent class 1 (Diabetes) was strongly associated with hypertension, physical inactivity and BMI. Hypertension increased the chance of male membership in the diabetic subclass, but in women reduced the chance of membership in the Pre-skeletal subclass. A study conducted in the population over the age of 20 years in Korea also showed an inverse relation between high blood pressure and low back pain 52 . Also a long-term cohort reported that high systolic or diastolic blood pressure was directly associated with low back pain occurrence. This association may be due to using anti-hypertensive medication and diminished pain sensitivity among patients with hypertension 53 .
In the present study, with increase in physical activity, the chance of male's membership in the diabetic subclass decreased, but no significant association was observed in females. This may be because in Iran, males have more opportunity for doing physical activity outdoors or in their workplaces. Low physical activity is estimated to be the cause of roughly 27% of diabetes in Iran 27 . The KORA-Age Augsburg study, performed on people 65 to 94 years, also showed an inverse association between physical activity and multimorbidity among males. Multi-morbidity referred to having at least two non-communicable diseases that were self-reported. In this study, 13 non-communicable diseases were defined according to the Charlson comorbidity index 54 . Similar to these findings Alimohammadian et al. showed an association between decreased physical activity and multimorbidity in adults, which was more evident in males 41 . Another study from Quebec, reported no association between multimorbidity and physical activity in both genders 55 .
In the present study, increased BMI significantly increased the chance of female membership in all disease subclasses, while in males; this association was seen only in the diabetic group. Similar findings have been previously reported in several studies in which obesity was associated with multimorbidity, and obesity increased the risk diabetes and other NCDs 13,56 .
In our analysis, tobacco consumption increased the chance of both genders membership in the musculoskeletal disorders subclass, which indicates that smoking in populations might be associated with musculoskeletal pain. Studies have shown that smoking increases the risk of multimorbidity 9,13 . Studies have suggested that both smoking and female gender were independently associated with a higher odds of reporting moderate to severe pain. Also, smoking was more strongly related to increased pain among women 46 . Also a review showed that tobacco smoking has negative effects on the musculoskeletal system 57 , and cigarette smoking may be associated with decreased muscle strength. However, more research about the mechanism that smoking affects the musculoskeletal system is needed. Compared with other non-communicable diseases, impaired musculoskeletal health is responsible for the greatest loss of productive years in the workforce 48 , and requires special attention.
In this study, males with lower education were more likely to be in the diabetic and musculoskeletal disorders subclass; while in females, they were more likely to be in the asthma and peri-skeletal diseases groups. Similar findings have been previously reported from several nations, and have shown that the prevalence of multimorbidity of chronic illnesses is significantly higher among people who were older, retired and less educated 43,58 . Also in West Asia, with increased levels of education, a decrease was seen in the prevalence of multimorbidity 41 . However, in Qatar, with increasing levels of education, multimorbidity increased 59 . In this present study, at province-level, the urbanization score was not associated with non-communicable multimorbidity subclasses in females. But, urbanization increased the odds of membership in the diabetes class among males. According to GBD2017, diabetes ranked fourth as the leading cause of age-standardized YLDs in the world 47 . Adekanmbi et al. conducted a multilevel analysis study in Namibia and showed increased odds of diabetes among individuals with the highest urbanization and development level 60 .
In this study, similar to previous studies 41, 61 , lower socioeconomic groups were associated with multimorbidity. This might partially be explained by the risk factors that the lower socioeconomic class faces such as unplanned urbanization, living in slums, physical inactivity, and consuming unhealthy food 41 . However, studies from China have shown that high income was associated with increased multimorbidity, which may be due to the high cost of health care in China 34 which make poor families less inclined to check their health.
WHO data from low-and middle-income countries revealed a reciprocal relation between the prevalence of multimorbidity and development 62 . In high-income countries, multimorbidity has been directly associated with development 42 . Differences in the prevalence of non-communicable diseases in urban areas may be due to the impact of health services and health literacy in these areas 63 .

Limitations
The current study had several limitations. First, this study was based on cross-sectional data which cannot show a causal relation between different risk factors and multimorbidity. Second, multimorbidity was assessed according to self-report and there is a possibility of recall bias. However, this approach is used in most national surveys, because of limited time and resources. Using self-reports also avoids potential recording errors that occur in administrative data, such as incomplete documentation, and miscoding 5 . Third, this study in comparison to the Charlson Comorbidity Index included a smaller number of diseases; and this is one reason for the few extracted classes. We suggest that future investigations include other chronic diseases such as mental disorders, kidney or liver diseases as well.
However, our study has several strengths as well. We used a large community-based sample which included adults from all provinces of Iran (3). In this study by using LCA both individual and community-level analyses was performed to describe the economic and social context in which individuals live and experience multimorbidity.
In order to control multimorbidity, health promotion and educational methods should be used to enhance public awareness about modifiable risk factors such as physical activity and smoking. These results highlight that policies, strategies and health programs for non-communicable diseases, must include musculoskeletal health as an integral component, particularly in low socioeconomic settings.
The identification of distinct subtypes of public population is important for describing the etiological pathways of morbidities which often co-occur with different NCDs. Other implications of the results of this study are helping to identify groups at risk; and can be used in developing countries with limited resources to tailor health education and services, to special groups. Also by recognizing these groups, the health system can allocate active services to them; and researchers can conduct etiologic studies by comparing these groups with other populations.

conclusion
The findings of this study showed that behavioral factors such as physical inactivity, tobacco use and urbanization are associated with multimorbidity in Iran; and there is a need to allocate resources for controlling, preventing and managing these risk factors. Prevention and control of multimorbidity requires health promotion programs that increase public awareness about modifiable risk factors, particularly among the at risk populations. A deeper understanding of these patterns may lead to the development of preventive measures to reduce the burden of these diseases, and offer new and comprehensive ways to manage these common conditions.