Health burden due to depressive disorders. The raw data for disability-adjusted life years (DALY) was used to represent the country-level disease burden due to depressive disorders. The data was provided by the World Health Organization (WHO) from the Global Burden of Disease (GBD) study database. Annual country-level estimations were available for the four follow-up years (i.e., 2000, 2010, 2015, and 2016) (http://www.who.int/healthinfo/global_burden_disease/estimates/en/index1.html). DALY, a summary metric of population health, included years lived with disability (YLD) and years of life lost due to premature mortality (YLL). The DALY data is available for the entire population and is grouped by various factors, such as sex and age (e.g., 0–4, 5–14, 15–29, 30–49, 50–59, 60–69, and ≥ 70 years). This study calculated the number of DALY per 100,000 persons for each country. Regarding depressive disorders, the International Classification of Diseases 10th revision codes (F32-F33, F34.1) for non-communicable diseases was applied to identified cases with depressive disorders. In total, 183 WHO member countries across five continents with depressive disorders data were analyzed. Figure S1 shows the spatial distribution of DALY lost due to depressive disorders after adjusting for population size.
Assessment of greenness. The Normalized Difference Vegetation Index (NDVI) measured by a Terra Moderate Resolution Imaging Spectroradiometer (Terra-MODIS) sensor with 1x1 km spatial resolution was used to estimate the presence of greenness in each country9. The data provided by the National Aeronautics and Space Administration (NASA) includes the monitoring and measuring of vegetation, plants, and biomass production, as well as components of greenness including chlorophyll, canopy structure, and leaf10. The monthly NDVI used in this study was MOD13A3 version 6 and for a given pixel has a range of values from − 1.0 to + 1.0. Positive values represent greener vegetation and negative values indicate limit vegetation11. Since recent studies have indicated an association between negative NDVI values and the proximity to water12, pixels with negative values were excluded to avoid misclassification bias due to the effects of water. In this study, satellite images with an acquisition date closer to mid-season were collected for January, April, July, and October. In total, 292 MODIS NDVI images were used to assess levels of green coverage of 183 selected countries. For image integration, a monthly global green coverage map was generated by combining 292 images. Next, we established similar procedures to assess greenness for the four selected months. Finally, monthly greenness concentrations were calculated to estimate the annual average values of greenness for each country. In all, this process integrated a total of 4672 images across the four follow-up years (i.e., 2000, 2010, 2015, and 2016). The spatial distribution of greenness in each country based on the NDVI estimations is shown in Fig. 2.
Covariates. Based on the findings of previous studies, several potential variables were considered for model adjustments including demographics, education, economic status, urbanization level, religion, divorce rate, alcohol consumption, smoking, chronic diseases/comorbidity, PM2.5, and continent (Table S1). Demographic covariates included population size, age, and sex because these factors were found to be associate with depressive disorders13. Moreover, the percentage of the educated population at a country-level was controlled because of its protective effect against depressive disorders14. Given economic status and urbanization level influence accessibility and affordability of medical care, these two factors were also included in the model15,16. Religion, divorce rate, blood pressure, and lifestyle behaviors including smoking and alcohol consumption were treated as covariates because they have been demonstrated to link with psychological health17–21. Given chronic diseases were found positively related to depressive disorders22, this study included malignant neoplasms, diabetes mellitus, cardiovascular diseases, and respiratory diseases as comorbid conditions. Exposure to air pollution such as fine particulate matter (PM2.5) was also assessed due to its negative effect on mental health23. Positive relationships between greenness exposure and mental health have been revealed for multiple continents, including both Americas24, Africa25, Asia26,27, Europe6,28 and Oceania29. In addition, better mental health may lower the risk of experiencing a depressive disorder30; thus, five continents were included so as to corroborate whether different regions have different impacts on depressive disorders resulting from degree of greenness exposure.
To assess the relationship between comorbidities and depressive disorders, Spearman correlation analysis was performed (Table S2). The reason to examine the correlation between comorbidities and depressive disorders was to ensure whether or not the selected comorbidities were eligible for covariates.
Statistical analysis. A series of statistical analyses were conducted to examine the effects of greenness on depressive disorders and the robustness of results, including descriptive analysis, GLMM, sensitivity test, stratified analysis, positive-negative exposure, and outcome controls. Descriptive analysis was performed to present country-level characteristics of all covariates examined in this study, including DALY of depressive disorders, demographic factors (population size, sex, and age), the prevalence rate of education level, economic status, urbanization level, the prevalence rates of populations without religion, divorce rate, lifestyle behaviors (alcohol consumption and smoking), blood pressure, and health burden due to chronic diseases, environmental exposures (greenness and PM2.5), as well as by continent.
The main model with adjustment for the above-listed covariates was developed by using the Generalized Linear Mixed Model (GLMM) with a penalized quasi-likelihood (PQL) algorithm to determine the relationship between exposure to green coverage and depressive disorders. The GLMM accounts for both fixed and random effects and provides a flexible approach for analyzing health outcomes. Given that the high number of variables adjusted for might affect the clustered spatial patterns31, an additional random effect of the country was considered in the GLMMPQL calculation to minimize biases due to the spatial autocorrelation effect. A residual test was employed by using Spatial Autocorrelation Global Moran's I to confirm whether or not statistical significance was affected by spatial autocorrelation. Furthermore, generalized variance-inflation factors (GVIFs) were applied to examine multicollinearity across covariates. The GVIFs value was less than 4 for all covariates. Hence, all variables were included in the model adjustment (Table S3).
To evaluate the robustness of findings, a sensitivity test was applied. Different covariates in ten separate models were adjusted to discern the change in coefficient estimation and significance. Model 1 adjusted covariates, including population size, sex, age, and year. Covariates were added gradually from model 2 to model 6. From model 7 to model 10, different comorbidity was added. Moreover, a stratified analysis was conducted based on the level of greenness in quartiles to assess differences in the impact of greenness on depressive disorder. Subgroup analysis was then applied to assess the association between greenness and a reduction in the burden of depressive disorders, stratified by sex, age groups, economic status, urbanization levels, and continents. Lastly, positive-negative exposure and outcome controls were performed. In the analysis of positive exposure control, the same outcome variable was used but exposure to greenness was replaced. Given a prior study found CO2 to be highly correlated with an increased risk of depression32. The relationship between CO2 exposure (i.e., positive exposure control) and risks of having depressive disorders were checked. For the negative exposure control, the assumption that no relationship exists between wind speed and depressive disorder was examined. Previous studies have found that greenness exposure is negatively associated with prevalence of cardiovascular diseases (CVD)33. Because of this, the associations between exposure to greenness and the burden of disease due to CVD, the positive outcome control, and Human Immunodeficiency Virus (HIV), the negative outcome control, were included in the analyses of positive-negative outcome controls.
All of the spatial and statistical analyses were performed using ArcGIS 10.7.1 (Esri Inc., Redlands, California, United States) and R version 3.6.3 (The R packages Foundation for Statistical Computing, Vienna, Austria). Coefficient estimates with 95% confidence intervals (CI) were reported and p-values < 0.05 were considered statistically significant.