Green space structures and schizophrenia incidence in Taiwan: is there an association?

A growing body of research indicates overall greenness offers potential psychological benefits. However, few studies have explored green space structures and their potential association with mental disorders. The aim of this study was to determine the existence of such an association. From two million randomly sampled people in Taiwan’s National Health Insurance Research Database, we selected 3823 patients that received a first-time diagnosis of schizophrenia from 2005–2016. Moreover, we used a geographic information system and a landscape index to quantify three characteristics of green space structures including area and edge, shape, and proximity. Additionally, we collected the normalized difference vegetation index and enhanced vegetation index data to reconfirm the association between overall greenness and schizophrenia incidence. We used the indices to determine individuals’ exposure according to their residential township. Spearman’s correlation analysis was conducted to select variables by considering their collinearity. Cox proportional-hazards models were applied to assess the relationship between green space exposure and schizophrenia incidence following adjustment for potential confounders, such as air pollution (NO2, SO2, ozone, and PM10), temperature, precipitation, and socioeconomic status, which are risk factors. We found a negative association between most green space structures indices and schizophrenia incidence. Our findings suggest that for green spaces, a larger mean patch area and edge density, higher complex (higher perimeter–area ratio), and greater proximity (higher contiguity index, aggregation index, and contagion index), may reduce the risk of schizophrenia. A sensitivity test and subgroup analysis revealed similar results.


Introduction
Schizophrenia is a chronic, severe, and disabling disorder (National Institute of Mental Health 2018) and is one of the 15 leading causes of disability worldwide (Vos et al 2017). Schizophrenia symptoms can be divided into three broad categories: positive (hallucinations, delusions, thought disorders, and movement disorders), negative (flat affect, reduced feelings of pleasure in everyday life, difficulty beginning and sustaining activities, and reduced speaking), and cognitive (poor executive functioning, trouble focusing or paying attention, and problems with working memory) symptoms (National Institute of Mental In particular, recent studies have observed the beneficial associations of green spaces on schizophrenia (Engemann et al 2018, Chang et al 2019. Such studies have used vegetation indices and green space density to assess overall greenness. Further, the pattern and allocation for green space, the indices of green space structures, are available for predicting the associations on human health (Shen and Lung 2016, Tsai et al 2018. The indices of green space structure are commonly used in landscape ecology, a term coined by the German geographer Carl Troll in 1939(Troll 1939) to refer to the study of spatially heterogeneous green spaces characterized by diverse interacting patches or ecosystems, including relatively natural terrestrial systems and human-dominated environments (Rocchini et al 2015, d'Acampora et al 2018. However, the following characteristics of green space structures have been used to determine the association with human health in few studies: the area and edge, shape, and proximity of green spaces. These characteristics are based on theories invoking green space exposure-related psychological mechanisms, such as stress recovery theory (SRT), emphasizing the role of nature in relieving physiological stress (Ulrich 1983, Ulrich et al 1991; attention restoration theory (ART), highlighting the role of nature in relieving mental fatigue (Kaplan and Talbot 1983, Kaplan and Kaplan 1989, Kaplan 1995; and prospect and refuge theory, emphasizing environments that provide security often provide people with the capacity to observe (prospect) without feeling threatened (refuge) (Appleton 1996). However, only one study has examined the association between indices for green space structures and mental health (Tsai et al 2018). Longitudinal studies with appropriate designs should be performed to examine the causal associations between different green spaces and mental health and should consider additional risk factors. Accordingly, the present study analyzed the association between indices for green space structures and the incidence of schizophrenia in Taiwan using data from the Longitudinal Health Insurance Database (LHID) for the period from 2005 to 2016. Based on previous theories and studies, the hypothesis is for green spaces, a larger area and edge density, higher complex, and greater proximity may reduce the risk of schizophrenia. This is the first study to examine the association between indices of green space structures and mental health by retrospectively analyzing a cohort with a large number of individuals (2 million). The findings from this study could provide information to policymakers for future planning allocation of green spaces.

Study area
Taiwan is an island located on the Tropic of Cancer; climatically, it is located in the northern and central subtropical and southern tropical regions. Due to its various climatic regions, Taiwan is a habitat for numerous flora and fauna. Taiwan also lies in the circum-Pacific seismic zone at the intersection of the Yangtze Plate, Okinawa Plate, and Philippine Mobile Belt. This location is characterized by a continual increase in the height above sea level as well as earthquakes. The total area of Taiwan is 36 197 km 2 , and its topography is characterized by rugged mountains running from north to south in the eastern two-thirds of the island and flat plains in the western one-third of the island (Ministry of the Interior 2018). According to the administration of resources, economies, and population, Taiwan can be divided into four regions. The descriptive summary of the township area and population density is depicted in table S1 (available online at (stacks.iop.org/ERL/15/094058/mmedia)). The average of areas is 103.57 km 2 (Q1, Q3: 32.42, 94.58) and population density is 1 611.70 persons/km 2 (Q1, Q3: 293.37, 1 806.57) among 349 townships. The areas of townships in the east are much larger than those in the west, including in northern, central, and southern Taiwan. The country experienced rapid economic growth and industrialization in the second half of the 20th century, resulting in this period being denoted as the period of the 'Taiwan Economic Miracle' (Clough 1991). Between 1952 and 1982, the average rate of economic growth was 8.7% per year; between 1983 and 1986, it was 6.9% per year. The gross national product increased by 360% between 1965 and 1986. However, along with economic development, the prevalence of common mental disorders doubled from 11.5% in 1990 to 23.8% in 2010 (Fu et al 2013).

LHID2000
We used the LHID2000, a subset of the National Health Insurance Research Database (NHIRD), from Health and Welfare Data Science Center, Taiwan, as our main data source to assess the association between indices of green space structures and schizophrenia incidence. Taiwan's National Health Insurance program was established in 1995 and covers more than 99.6% of residents of Taiwan. The NHI is a single-payer and compulsory social insurance system regulating the disbursement of health care funds and guaranteeing equal access to health care for all citizens. By the end of 2014, 99.9% of the citizens in Taiwan (approximately 23 340 000 citizens) had been enrolled in this program. The LHID2000 contains the medical information of two million individuals who had been insured in 2000 randomly sampled from the NHIRD by age, sex, and the registry of regions from the full database population. This database contains original claims data and a registry of the beneficiaries, including demographic characteristics, inpatient and ambulatory care, diagnoses, catastrophic illness certificates, medical expenditures, operations, prescriptions, investigational items, examinations, procedures, treatments, hospital levels, and medical divisions (Lee 2017, Lin et al 2018. Individuals included in the LHID2000 and those in the original Taiwan NHI were demonstrated to not significantly differ, particularly in terms of age, sex, mean monthly health insurance premium, or distribution in townships (National Health Research Institutes 2003). The identification of schizophrenia in the LHID2000 was based on International Classification of Diseases, Ninth Revision (ICD-9) code 295 and ICD, Tenth Revision (ICD-10) code F20. Accordingly, individuals who applied for at least one severe claim in the registry for catastrophic illness during 2005-2016 were selected in this study, and the first application date was treated as the index date. The period from 2000 to 2004 was selected as a wash-out period; patients who had schizophrenia during this period were masked from the analysis. Individuals residing on off-shore islands such as Kinmen and Matsu were also excluded. Finally, for our analysis, we recruited a total of 1 918 501 individuals without a history of schizophrenia. Instead of providing the residential addresses of the patients, the LHID provides their hospital details for confidentiality reasons. Therefore, we assumed all patients lived near the hospital they most often visited for the common cold (ICD-9 code 460 or 465 and ICD-10 code J00-06 or J30-39) during the study period. Figure 1(a) displays the incidence of schizophrenia from 2005 to 2016 at the township level. We also assessed demographic and socioeconomic data, namely age, sex, and monthly health insurance premium (a proxy for personal financial status). Our study protocol was reviewed and approved by the Governance Framework for Human Research Ethics (A-EX-108-009), and we conducted the study in accordance with the Declaration of Helsinki.

Greenness databases
This study used two databases to calculate the indices of green space exposure: the national land survey database and a remotely sensed greenness database. The national land survey database was established by the National Land Surveying and Mapping Center (NLSC). The first national land survey was conducted from 1993 to 1995. To follow the change of land use, the NLSC continually conducts surveys and updates the database. One survey is conducted every 5-10 years and can represent the situation of land use within a survey period. Because the first national land survey was earlier than the study period, we applied the second national land survey database established in 2006 for the primary analysis in the present study.
In addition, the database of the third national land survey starting in 2014 was used to verify the consistency of the result. The data are vector data and the digital data layer of land use is based on the 0.25 m resolution digital orthoimagery in Taiwan. The land use data contains nine major categories, including agriculture, forests, transportation, water, building, public facilities, recreation, mining/salt production, and others (National Land Surveying and Mapping Center 2020). We selected forests as forest spaces. We selected the subcategories from recreation, including a university, park, and grassland as recreational green spaces. Moreover, we combined forest spaces with recreational green spaces to be all green spaces. All of these green spaces were used to quantify indices of green space structures and all green space densities at the township level. We used the normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI) as indices of surrounding greenness throughout our analysis. NDVI and EVI data were extracted using the global moderate-resolution imaging spectroradiometer maintained by the National Aeronautics and Space Administration and United States Geological Survey. The NDVI and EVI are spectrum-based greenness indices for measuring and monitoring plant growth (vigor), vegetation cover, and biomass production using multispectral satellite data. In such data, healthy vegetation exhibits low red-light reflectance and high near-infrared reflectance; thus, healthy vegetation shows high index values. Unhealthy or sparse vegetation reflects more visible light and less near-infrared light. NDVI and EVI values range from −1.0 to 1.0. The difference between the two indices is the EVI is corrected for some distortions in the reflected light caused by the particles in the air as well as the ground cover below the vegetation. This spectroradiometer provides NDVI and EVI data every 16 d; the data are gridded at a 250 m spatial resolution in the sinusoidal projection. Figure  S1(a) illustrates the spatial distribution of the cumulative mean NDVI and EVI around Taiwan for the period from 2000 to 2004. To represent exposure during the study period, the cumulative mean concentration of mean NDVI and EVI was calculated for the period from 2000 to 2004. The indicators of greenspace exposure for individuals at the township level were determined according to the hospital they most often visited for the common cold. The overall greenness assessment is displayed in table S2.

Model adjustment databases
According to previous studies, a total of three databases should be employed for model adjustment, namely meteorological, air pollution, and township urbanization level databases. We obtained monthly weather records from 333 monitoring stations of Taiwan's Central Weather Bureau. We estimated the spatial average of monthly mean temperature and total precipitation at the township level using the ordinary kriging interpolation method through the spherical model. Records of air pollutants, such as records of nitrogen dioxide (NO 2 ), particulate matter with a diameter of ≤10 µm (PM 10 ), ozone (O 3 ), and sulfur dioxide (SO 2 ), were acquired from 76 air quality monitoring stations of the Environmental Protection Administration, Executive Yuan, Republic of China (Taiwan). The spatial average of air pollutants was also estimated using the ordinary kriging interpolation method through the spherical model.
Township urbanization level was determined based on 2000 census and 2004 population statistics. On the basis of the Taiwan Social Change Survey database, the 349 townships and boroughs in mainland Taiwan were reorganized into six strata: metropolis, industrial or commercial district, emerging town or city, traditional industrial town or city, less developed town or city, and aging or remote town or city (Hou et al 2008). To represent exposure during the study period, the cumulative mean concentration of mean temperature, total precipitation, and air pollutants was calculated for the period from 2000 to 2004. The indicators for meteorological variables, air pollution, and township urbanization level were again chosen according to the hospital individuals most frequently visited for the common cold.

Quantification of green space structures
Landscape indices are typically grouped according to eight characteristics of landscape pattern measured. It contained area and edge metrics, shape, core area, contrast, proximity, subdivision, isolation, and diversity (Mcgarigal et al 2002, Mcgarigal 2017). The candidate characteristics and related indices were based on previously conducted studies (Shen and Lung 2016, Tsai et al 2018. Accordingly, we consider the following characteristics: area and edge, shape, and proximity and selected one index for each characteristic, including mean patch area (area and edge), contiguity index (shape), and aggregation index (proximity) as the indices in this study by using FRAGSTATS 4 (Mcgarigal et al 2002, Mcgarigal 2017. We also calculated an edge contrast index as a negative control; this is because it belonged to the contrast characteristic, not to the three characteristics supported by theories, and was revealed to have a random and no significant association with frequent mental distress in a previous study (Tsai et al 2018). The mean patch area shows the mean of the total area of vegetative patches; thus, high index values mean more green spaces in this area. Meanwhile, the contiguity index was used to assess the spatial contiguity or connectedness of cells within a grid-cell patch. The value increases as the vegetative patches become more contiguous. Thus, the aggregation index was used to measure the ratio of the number of like adjacencies to the maximum possible of like adjacencies involving the corresponding class and multiplied by 100 to convert to a percentage. The value increases as the vegetative patches become more physically connected. Finally, edge contrast index represents the percentage of maximum contrast for the average unit of the edge between two patches; thus, high index values mean more difference between adjacent vegetative patches (Mcgarigal et al 2002, Mcgarigal 2017).

Statistical analysis
We hypothesized exposure to the indices of green space structures positively affects schizophrenia incidence. Considering our retrospective cohort design, we used a Cox proportional hazards model to estimate the association of indices of green space structures with the risk of schizophrenia (Cox 2018). We first used the selected indices of green space structures to calculate the hazard ratio (HR) for schizophrenia incidence. To evaluate the nonlinearity between overall greenness or indices of green space structures and schizophrenia incidence, we classified the green space indices into deciles according to procedures utilized in recent studies (Engemann et al 2018(Engemann et al , 2019. Additionally, we used general greenness measures such as the cumulative mean NDVI and EVI from 2000 to 2004 and all green spaces density to reconfirm the association of exposure to overall greenness on schizophrenia incidence. Table S3 presents the cutoff value for each index. Further, the HR for schizophrenia incidence was estimated within the various deciles compared with the HR in the lowest one. We avoided collinearity by employing Spearman's correlation analysis. The inclusion criterion for independent variables was a correlation coefficient of <0.7. To avoid the collinearity issue, Spearman's correlation coefficient <0.7 was used as a criterion for the inclusion of independent variables. Because of the high correlation between the mean patch area and aggregation index (R = 0.9), models were built for each green space structure separately instead of considering all three green space structures in a single model. Moreover, correlations among air pollutants are less than 0.6 (table S4); therefore, NO 2 , PM 10 , O 3 , and SO 2 were adjusted in developed models. Demographic and socioeconomic adjustments were made using representations of a range of well-established factors obtained from the LHID2000, including age (0-30 years and >30 years); sex; and monthly health insurance premium, which was stratified into two groups according to a previous study (Chien et al 2004): less than 20 000 New Taiwan Dollars (NTD) and more than 20 000 NTD. The monthly health insurance premium was adjusted to the nominal exchange rate of 1 USD = 32.214 NTD on December 31, 2016. The groups were transformed into <621 USD and >621 USD. The meteorological and air pollution-related variables and township urbanization level were considered for model adjustment. A P value of 0.05 was considered statistically significant. ArcGIS (version 10.4.1; Environmental Systems Research Institute Inc. Redlands, CA, USA) and SAS 9.4 (Cary, NC, USA) were used for all analyses.

Sensitivity testing and subgroup analysis
To examine the robustness of the indices in determining the association between the indices of green spaces and schizophrenia incidence, we conducted a sensitivity test using six approaches. First, we verified the association between overall greenness and schizophrenia incidence by using 1 and 3 year cumulative mean NDVI and EVI. Second, we used continuous values for the indices of green spaces for robustness verification. Third, we restricted our analysis to individuals in six special municipalities: Kaohsiung, New Taipei City, Taichung, Tainan, Taipei, and Taoyuan. Fourth, we selected other indices based on the three characteristics of green space structures, namely edge density, perimeter-area ratio, and contagion index, for the assessment. Edge density equals the sum of the lengths of all green space edge segments per hectare; thus, more green spaces show high index values. Perimeter-area ratio equals the ratio of the perimeter of green space per square meter; thus, complex green spaces show high index values. Contagion index is minus the sum of the proportional abundance of each type of green spaces multiplied by the proportion of adjacencies between cells of each green space and others, multiplied by the logarithm of the same quantity, summed over each unique adjacency of green spaces; this is divided by 2 times the logarithm of the number of type of green spaces, which is further multiplied by 100 (to convert to a percentage); thus, aggregative green spaces show high index values (Mcgarigal et al 2002, Mcgarigal 2017. Fifth, we used green space structures quantified from the national land survey 2014 for robustness verification. Sixth, we followed the process reported in a recent study (Banay et al 2019) to categorize indices into quintiles for verifying the consistency in the association between greenness exposures and schizophrenia incidence. For our subgroup analysis, three approaches were used. Individuals were divided into the following groups based on sex, age, and monthly health insurance premium: male and female; 0-30 years old and >30 years old; poor (<621 USD) and rich (>621 USD). A continuous exposure was used for subgroup analyses. Western Taiwan consists of major cities, whereas eastern Taiwan consists of natural landscapes. Therefore, we stratified townships into eastern and western Taiwan to assess the effects of green space structures in the two regions. Because a large area of eastern Taiwan is covered with forest, resulting in a narrow range of green index variation, a continuous format of green indices was used for the analysis instead of categorical classifications. We calculated the P value for the HR to determine the statistical significance of the relationship between the indices of green spaces and stratified variables of schizophrenia risk. Table S2 presents descriptive statistics for the indices of green space structures. Since more than 90% of green spaces were forest spaces, we observed a similar trend between forest spaces and all green spaces. The mean patch area of all green spaces, forest spaces, and recreational green spaces were 42.22, 48.95, and 0.85 ha, respectively, and the corresponding contiguity indices were 0.48, 0.48, and 0.47; therefore, all green spaces were the most contiguous. The aggregation indices for all green spaces, forest spaces, and recreational green spaces were 90.40%, 89.40%, and 84.35%, respectively, indicating all green spaces or forest spaces exhibited greater proximity than did recreational green spaces. Figure 1(b) and figures S1(b)-(c) provide a graphical representation of the statistics in table S2.

Descriptive statistics
We included a total of 1 918 501 individuals for this study. During 2005-2016, 3 823 new cases of schizophrenia were identified. We stratified the NDVI values into the following percentiles and evaluated the distribution of the individuals' demographic characteristics in these percentiles: <10th, 10th-20th, 20th-30th, 30th-40th, 40th-50th, 50th-60th, 60th-70th, 70th-80th, 80th-90th, and >90th in table S5. The mean NDVI values in these percentiles were 0.18, 0.23, 0.29, 0.34, 0.38, 0.41, 0.45, 0.50, 0.55, and 0.64. The highest mean monthly health insurance premium was 988 USD and was observed in the 10th-20th percentile (NDVI from 0.20 to 0.26); the lowest was 840 USD and was observed in the >90th percentile (NDVI >0.58). A higher proportion of men lived in greener places. The highest proportion of individuals were aged 0-30 years, the peak age of onset of schizophrenia, was 22.15% and was observed in the 60th-70th percentile (NDVI from 0.43 to 0.48); the lowest was 18.15% and was observed in the <10th percentile (NDVI <0.20). The highest proportion of individuals with monthly health insurance premiums of <621 USD was 22.99% and was observed in the <10th percentile; the lowest was 16.74% and was observed in the 50th-60th percentile (NDVI from 0.39 to 0.43). Although the <10th percentile had the lowest values regarding sociodemographic characteristics, it did not include the highest incidence of schizophrenia (0.21%), implying environmental factors may play a role in the incidence of schizophrenia. Figure 2 illustrates the association between overall exposure to greenness or indices of green space structures and schizophrenia incidence in Taiwan. An HR of <1 was considered to indicate the beneficial associations of green space exposure on a population. An HR of >1 was considered to signify green space exposure was associated with a higher risk of schizophrenia. After adjusting for demographic, socioeconomic, meteorological, and air pollution factors, the coefficients of overall greenness were associated with lower HRs for schizophrenia incidence of almost percentiles compared with the reference level (the lowest decile). In addition, the association between the indices of green space structures and schizophrenia incidence was supported by our hypothesis. A larger mean patch area, higher contiguity index, and higher aggregation index were associated with lower HRs for schizophrenia incidence. It implied that green spaces, a larger area, and greater proximity, may reduce the risk of schizophrenia. Although no significant doseresponse relationship was observed for all green space indices, almost all HRs were <1 compared with the reference level. This demonstrates the beneficial associations of green spaces among all green spaces, forest spaces, and recreational green spaces. In all green spaces, the most favorable HR [95% confidence interval (CI)] for mean patch area was 0.77 (0.66, 0.89) and was observed in the sixth percentile; that for the contiguity index was 0.80 (0.69, 0.94) and was observed in the seventh percentile, and that for the aggregation index was 0.85 (0.73, 0.99) and was observed in the eighth percentile relative to the reference level.

Construction of the main model
In addition, the edge contrast index exhibited a random association, no consistently positive and negative associations, and the results seemed corroborated.

Sensitivity testing and subgroup analysis
The six sensitivity analysis approaches yielded similar results. Figure S2 depicts the HRs derived for 1 and 3 years cumulative mean NDVI and EVI. Table 1 presents the continuous values of green space exposure. The assessments were restricted to six special municipalities, and the corresponding results are illustrated in figures S3(a) and (b). Figure S3(c) illustrates the HRs derived for other green space indices such as edge density, perimeter-area ratio, and contagion index. Figure S4 depicts the associations between green space structures quantified from the national land survey 2014 and schizophrenia incidence. Figure S5 presents the assessment results for green space indices categorized into quintiles. As presented in figure S2, HRs for all categorization of overall greenness were <1, indicating protective associations between overall greenness and schizophrenia incidence regardless of cumulative years. In addition, the HRs for overall greenness is also <1 in table 1. In particular, the contiguity index and aggregation index were significant in all green spaces (HR = 0.57, 95% CI = 0.33, 0.96; HR = 0.99, 95% CI = 0.98, 0.99) and aggregation index was significant in recreational green spaces (HR = 0.99, 95% CI = 0.98, 0.99). In general, all classes remained associated with low HRs of schizophrenia incidence in the six municipalities ( figure S3(a)). In particular, the mean patch area and aggregation index results were significant ( figure  S3(b)). Furthermore, consistent results were observed for other green space indices ( figure S3(c)). The most favorable HR (95% CI) for edge density was 0.82 (0.69, 0.97) and was observed in the fourth percentile, that for perimeter-area ratio was 0.79 (0.67, 0.93) and was observed in the fourth percentile, and that for the contagion index was 0.75 (0.65, 0.8) and was observed in the fourth percentile relative to the reference level. Negative associations were found between Association between HRs (estimates with 95% CIs) for schizophrenia and indices of (a) overall greenness and (b) green space structures. All models were adjusted for age, sex, health insurance premium, township urbanization level, temperature, precipitation, NO2, PM10, O3, and SO2. green space structures and schizophrenia incidence based on data from national land survey 2014, and results are illustrated in figure S4. Additionally, consistent results were found for green space indices categorized into quintiles and are depicted in figure S5. Moreover, the variables of sex, age, and monthly health insurance premium were used for stratification. The results are shown in figure 3 and figures S6 and S7. Overall greenness indices were negatively associated with schizophrenia incidence. Similar associations tended to be observed for male individuals, individuals with higher monthly health insurance premiums, and individuals aged <30 years. In addition, the mean patch area index was associated with lower HRs for schizophrenia incidence on individuals aged <30 years, and the contiguity index was associated with lower HRs for schizophrenia incidence on female individuals or those with a health insurance premium of <621 USD month −1 . Similar findings were observed from subgroup analyses regardless of age, sex, and monthly health insurance premium (tables S6-S8). The results of a subgroup analysis for the eastern and western regions of Taiwan are presented in table S9. In general, consistent findings were observed in the two regions.

Discussion
This study used NDVI, EVI, and all green spaces density to assess overall greenness and determine the beneficial associations of overall greenness on schizophrenia incidence. The beneficial associations of overall greenness from this study are in line with those of previous studies (Engemann et al 2018, Chang et al 2019. In addition, compared with our previous study (Chang et al 2019), the present study improved the time scale of exposure variables by using time-dependent average variables. Newly diagnosed schizophrenia was identified using a catastrophic illness database rather than outpatient records and was adjusted for air pollution (NO 2 , SO 2 , O 3 , and PM 10 ). The possible pathophysiology of schizophrenia includes an imbalance in the complex and interrelated chemical reactions of the brain involving neurotransmitters, dopamine, glutamate, energy metabolism of mitochondrial abnormalities, and abnormal inflammatory processes (National Institute of Mental Health 2016, Radhakrishnan et al 2017, Guest 2019, Marder and Cannon 2019. Previous research found green space exposure could reduce stress (Frumkin et al 2017); reduce prefrontal cortex activity during negative-emotion processing (Tost et al 2019); inhibit inflammatory cytokine production, such as interleukin (IL)-1β, IL-6, and nuclear factor-κB (Cho et al 2017); and improve the human innate immune system (Li 2010). In addition, green spaces could mitigate the negative effects of stress and environmental pollution in modern society to improve mental health indirectly (Van den Berg et al 2015, Fong et al 2018. The present study observed green spaces tended to have distinct protective associations on male individuals, individuals with a relatively high monthly health insurance premium, and individuals aged <30 years. The possible reason for this is men are more likely than women to use green spaces or recreational areas for physical activity (Cohen et al 2007) and women may not access green spaces if they perceive them as unsafe (Virden and Walker 1999). Recent studies indicated socioeconomic status could influence the structure of the developing brain (Brito andNoble 2014, Farah 2017). Lower family income and lower parental education have been reported to be associated with reduced volumes of cortical grey matter, the hippocampus, and the amygdala (Luby et al 2013) as well as reduced cortical surface area and cortical thickness in children and adolescents (Mackey et al 2015, Noble et al 2015. In addition, the brain develops steadily during prenatal and early postnatal periods, which are considered the most vulnerable windows for the effects of environmental exposures (Grandjean and Landrigan 2014); therefore, exposure to greenness early in life may beneficially influence brain development and cognitive function (Dadvand et al 2018). This study found beneficial associations of different indices of green space structures. We noted that green spaces, a larger mean patch area and edge density, higher contiguity index, higher perimeterarea ratio, and higher aggregation index and contagion index were associated with a lower schizophrenia incidence. This suggests a larger area and edge density, higher complex, and greater proximity for green space, which may reduce the risk of schizophrenia. We observed similar results in three types of green spaces and throughout Taiwan or metropolitan areas. In addition to SRT (Ulrich 1983, Ulrich et al 1991, ART (Kaplan and Talbot 1983, Kaplan and Kaplan 1989, Kaplan 1995, and prospect and refuge theory (Appleton 1996), our findings are supported by previous studies. For example, greater edge density and larger green spaces have been reported to be more beneficial for recreational and physical activities (Tsai et al 2016). Similarly, contiguous or complex green spaces were reported to be positively correlated with neighborhood satisfaction (Lee et al 2008). More access and opportunities are created for people to enjoy green spaces where may improve mental health. The mean patch area index and contiguity index had protective associations against schizophrenia incidence in individuals aged <30 years old   Figure 3. Subgroup analysis of HRs (estimates with 95% CIs) for the association between schizophrenia incidence and indices of (a) overall greenness and (b) green space structures according to sex. All models were adjusted for age, health insurance premium, township urbanization level, temperature, precipitation, NO2, PM10, O3, and SO2. and in female individuals or individuals with a health insurance premium of <621 USD month −1 , respectively. The possible reason is younger individuals were more likely to access green spaces or recreational areas for physical activity (Cohen et al 2007). Although the present study examined associations between green space structures and schizophrenia incidence in different stratifications, evidence explaining linkages among sex, age, and health insurance premiums and green space structures remains limited. The results should be confirmed in future studies. Although we found an association between indices of green space structures and schizophrenia incidence, our study has some limitations. First, given the privacy policy, the residential address of each individual was masked by the database providers; instead, the location of the hospital visited for the common cold (because people usually visit the nearest hospital for minor illnesses) was used to determine the individuals' exposure to surrounding green spaces. This approach could have induced some bias. For example, in rural townships with limited medical services, the most frequently visited hospital or clinic might not be close to the individuals' residence. Any relocation events or trips abroad were not considered; therefore, the actual exposure during the period may not be very precise and would possibly lead to random misclassification and under-estimation. Furthermore, the residential location information was not provided by database providers; therefore, the individuals' exposure data can only be estimated based on the location of hospital the individuals most frequently visited for common cold. It may cause misclassification of exposure estimation and may not be fully representative of the individuals' exposure. Table S1 presents a descriptive summary of the township area and population density. The average township area and population density are 83.38 km 2 (Q1, Q3: 31.22, 78.81) and 1 891.43 persons/km 2 (Q1, Q3: 328.82, 2 380.92) for the western region and 264.01 km 2 (Q1, Q3: 75.13, 320.55) and 119.82 persons/km 2 (Q1, Q3: 68.38, 215.66) for the eastern region, respectively. In addition, populations in Taiwan are unequally distributed. The western region has 95.33% of the Taiwan population, and the eastern region has only 4.67%. Although dominance in population distribution was found in the western region, the results of our findings might be predominantly influenced by people living in the western region. Moreover, smaller township areas coupled with a high vehicle density of approximately 93.8 vehicles per 100 people (Ministry of Transportation and Communications 2020) could imply better accessibility within the township in the western region, reducing the uncertainty in the exposure assessment. Additionally, similar associations between the indices of green spaces and schizophrenia were observed in eastern and western Taiwan  (table S9), which may imply the limited uncertainty in the present study. Second, the choice of residential location is potentially affected by genetic confounding, which cannot be completely ignored. According to a previous study, people with higher genetic loading for schizophrenia lived in denser urban areas, as determined on the basis of polygenic risk scores (Colodro-Conde et al 2018); however, another study found the correlation between the risk of schizophrenia and the urban city was not explained by genetic liability (Paksarian et al 2018). Therefore, the effect of genetically determined choices of residential location was unlikely to drive the association between exposure to the indices of green spaces and schizophrenia incidence entirely. Further, the individuals' gene information and family history could not be determined in this study. However, a previous study observed genetic variation accounted for a much smaller proportion of variance than did twin-based heritability estimates; the study suggested the rate of genetic involvement was 20%-25% in schizophrenia (Uher 2014). Another study suggested the effects of single-nucleotide polymorphisms were considerably less relevant than environmental factors (Lencz et al 2013). According to the above findings, missing genetic information and family history would not bias or limit our conclusions regarding the relationship between exposure to the indices of green spaces and schizophrenia incidence. This study used data from the national land survey database acquired for only one period to calculate the indices of green spaces. Additionally, we assumed there was no considerable change in green spaces in our study period. Therefore, we recommend future studies use satellite images to assess the distribution of land use in different periods and substantiate these analyses after considering the aforementioned limitations.

Conclusions
This is the first study to use a national retrospective cohort database to examine the association between indices of green space structures and schizophrenia incidence. The results reveal that green spaces, a larger mean patch area and edge density, higher contiguity index, higher perimeter-area ratio, and higher aggregation index and contagion index were associated with a lower schizophrenia incidence. Considering the social and medical burden of schizophrenia and the quality of long-term care, indices of green space structures are crucial elements that should be considered by policymakers when planning allocation measures for green spaces in urban areas. Further research is required in additional locations to verify the associations of indices of green space structures on different ethnic groups. We suggest future studies duplicate these analyses in different countries to assess the associations of geographical differences on the benefits of the indices of green space structures.