Asthma is a chronic respiratory disease with variable and recurring symptoms, airflow obstruction, bronchial hyperresponsiveness, and inflammation [3]. Its prevalence rose by 15% in the last 10 years [4]. Based on a Wisconsin Department of Health Services asthma surveillance report, approximately 14% of adults and 10% of children have been diagnosed with asthma in Wisconsin [5]. In 2009, 5,300 people were hospitalized and 21,000 went to an emergency department with a principal diagnosis of asthma. Eleven percent of adults with asthma had an emergency department visit and 20% had urgent care visits for symptoms [5].

Asthma onset is associated with multiple, complex factors. While some are non-modifiable such as sex and age [6], many others are associated with the environment and residential location. These include educational attainment, household income, health insurance, smoking, physical activity, and obesity [6]. Medical conditions influenced by the environment and associated with asthma include atopy [7], allergic reactions [8], airway hyperreactivity [9], and airway responsiveness [10]. Over 370 outdoor and indoor environmental factors have been associated with asthma including substances from building materials, cleaning products, personal care products, central heating systems, maintenance, and humidification devices [11].

Geographic information system (GIS) analyses have been used to study geographic environmental variables associated with asthma. The most studied variable was air pollution [12], which has been measured via passive measurement, direct measurement, proximity to roadways, and traffic carbon emissions. Besides air pollution, asthma was associated with climate differences [13], latitude [14], and socioeconomic status [15]. Socioeconomic status, specifically male employment, was positively associated with asthma in a Southern California study, where access to care and the hygiene hypothesis—the idea that limited exposure to bacterial and viral pathogens during childhood result in a predisposition to allergy [16,17]—were proposed as explanations.

Fewer asthma studies have incorporated local environmental variables aggregated at the level of census tracts or block groups. Census tracts and block groups are geographic areas developed by the United State Census Bureau and contain 1,500–8,000 and 600–3,000 people, respectively. Using census tract data, asthma diagnosis was correlated with houses facing highway intersection [18] and sociodemographic characteristics of race, sex, and education [19]. Fewer studies have used block group level variables. Socioeconomic status was associated with asthma diagnosis using block group level data [15]. Many of these analyses used questionnaire data to determine asthma diagnosis, which may be limited by self-report bias [20]. These analyses involved less than 5700 participants, 10 environmental variables, and census geographic regions from only a portion of a state.

Environmental variables and built environments have been studied using EHR data. For example, nitrogen oxides were tested for association with diseases including asthma diagnoses obtained from EHR datasets in primary care [21]. Body mass index (BMI) calculated from EHR data was positively associated with the number of fast food restaurants near a person’s home [22].

Schwartz et al. [23] used an EHR dataset, environmental community-level variables, and multilevel statistical analysis to demonstrate that lower BMI was associated with higher socioeconomic status and areas with more venues for physical activity.

Spatial statistics offer methods to incorporate geographic location to identify risk factors associated with disease [24]. The spatial statistics utilized in this study included a generalized additive model. Generalized additive models [25] are generalized linear models with predictors that involve a linear sum of smooth functions.

Previous health studies that utilized spatial generalized additive models investigated the association of air pollution and mortality, tuberculosis drug resistance patterns in Peru [26], and geographic distribution of heart disease [27].

Spatial statistics, specifically additive models, have been combined with sparsity. COSSO [28] and SpAM [29] extended the lasso estimator [30] while another approach created a new sparsity-smoothness penalty [31].

Our goal was to identify an interpretable set of environmental risk factors of asthma distributed geographically. Other studies have combined environmental variables and EHR data, spatial statistics and disease, and spatial statistics and sparsity. Our main contribution is the addition of sparsity to spatial statistics. As applied to geographically distributed EHR and environmental datasets, we describe this methodology as Sparse Spatial Environmental Analysis (SASEA).

Our research group developed the University of Wisconsin Electronic Health Record - Public Health Information Exchange (UW eHealth-PHINEX), an EHR data exchange between University of Wisconsin (UW) Departments of Family Medicine, Internal Medicine, and Pediatrics and the Wisconsin Division of Public Health. Further details have been described previously [32]. Briefly, the database contains clinical care variables such as disease diagnoses, medications, and laboratory test results. Patient home addresses from year 2012 were geocoded to year 2000 block groups, the smallest geographic area the US Census Bureau publishes. Block groups were linked to detailed demographic and environmental data from the ESRI Business Analyst database [33]. The data exchange is a HIPAA Privacy Rule compliant-limited dataset, and the Wisconsin Division of Public Health is blinded to patient/provider specific information. All patient identifiers were removed from the data except birth month and year, ZIP code, and census block group of the patient’s address. Random accession numbers were used for patients, primary care providers, and clinics. This study was approved by the UW Institutional Review Board protocol M2009-1273 and UW Health with data use agreements.

UW Departments of Family Medicine, Internal Medicine and Pediatrics provide care in 42 clinics throughout Wisconsin, but most are located in southcentral Wisconsin. Patients represent various environmental and socioeconomic strata in rural and urban regions.

The dataset study period was from 2007–2009. Patients were identified as asthma cases when an asthma ICD-9 code of 493.xx was associated with a Current Procedural Terminology (CPT) codes for hospital discharges (CPT codes 99238 and 99239) or office visits (CPT codes 99201–99205 and 99211–99215). Patients were identified as controls if they did not have a hospital discharge or office visit associated with an asthma ICD-9 code over the study period, but were seen at least once in the UW Departments of Family Medicine, Internal Medicine, or Pediatrics. Participants in the study were restricted to be 5 to 50 years of age. There were no additional exclusion criteria.

This study included 199,220 participants [32]. There were 103,690 patients living in 2,186 block groups with sufficient data also linked with ESRI data to perform the analysis described in section 2.3.

The ESRI Business Analyst environmental database [33] consisted of 1,117 variables, which included demographics (age, income, education), living conditions (household members, rental property, pets, rural living), behaviors (food consumption, transportation, smoking, television), health (drug prescriptions), and businesses (types of employees and employers). Most variables (992 of the 1,117) represented data from year 2010 while the remaining variables represented data from the year 2000 (please see Appendix Table 1). Variables were normalized to the number of participants or number of households when appropriate and standardized to (0,1).

The large-scale spatial variation of asthma was estimated using a Logistic generalized additive model with a thin plate regression spline smoothing term [34], which we refer to as a Logistic thin plate regression spline model. As described in the Introduction, generalized additive models [25] are generalized linear models with predictors that involve a linear sum of smooth functions. Smooth functions allow a more flexible model specification that can account for the spatial location of variables. A thin plate regression spline is considered an optimal smooth function as it was developed for optimal smoothness and data fitting using a more computationally feasible low rank approximation [34]. Thin plate regression splines do not require user-specified locations of knots and are multivariate, penalized low rank approximations of a smooth function with optimal data fitting and smoothness [34]. Tensor product smooths were not used as both longitude and latitude were scaled similarly. The geographic area of Wisconsin was small and did not necessitate pseudosplines on a sphere [35]. The thin plate regression spline was represented by a bivariate smooth term with the longitude and latitude of the block group centroid.

ArcGIS software [36] was used to map the total number of patients, prevalence, and Logistic thin plate regression spline modeling predicted prevalence per block group. Block groups with ≤ 20 total participants (asthmatic and non-asthmatic) were mapped with a different coloring scheme than block groups with > 20 total participants.

The Logistic thin plate regression spline model with covariates was: (1)log(asthmai,j)=f(xi,yi)+αjblockj+β1agei+β2sexi+β3racei+β4ethnicityi+β5(BMI)i+β6(encounterdays)i+β7(distance)i where i is a participant and j is the block group participant i’s home address is geocoded to. The thin plate regression spline is f(xi,yi)=∑k=1qck(xi,yi)ζk where c_k(x_i, y_i) is the k^th basis function, ζ_k is an unknown parameter, and x_i and y_i are the latitude and longitude for the centroid of the block group participant’s geocoded home address. α_jblock_j is the block group random effect allowing for hierarchical0 structuring of the model. The basis dimension, q, was chosen to be 80, which was twice the estimated degrees of freedom to allow for appropriate smoothness. BMI was the body mass index at first encounter. The encounter days covariate was defined as the number of days between a patient’s first and last encounter in the EHR dataset. Encounter days controlled for the differences between patients who utilized the University of Wisconsin’s hospitals and clinics over a short amount of time (e.g., those who had one visit to the emergency department) versus patients who utilized the hospitals and clinics over a longer amount of time (e.g., those who received the majority of their medical care at the University of Wisconsin). The distance covariate was defined as the Euclidean distance between a patient’s home address and the address of the primary care office with the most frequent visits.

An adapted Logistic generalized additive model fitting with subsampling for smoothing spline fitting was used to accommodate the large dataset [37,38]. Subsampling was a technique used for faster computation and did not cause parameter estimate variability. The smoothing splines were first set using a subsample of the data. In each subsequent step of the penalized iteratively re-weighted least squares (PIRLS) algorithm, the weighted model matrix was constructed in blocks with the corresponding QR decomposition so as not to form the entire model matrix. This method is justified for restricted maximum likelihood estimation because of asymptotic multivariate normality of Q’z, where z is the pseudodata. This adapted method was previously implemented in the R package mgcv using the bam function with tp parameter [34].

The 1,117 environmental variables from ESRI were dimensionally reduced using sparse principal component analysis (SPCA) [39] before testing for association with asthma. SPCA is in contrast to principal component analysis (PCA). In PCA, the principal components are a linear combination of the original variables. SPCA uses only a small number of non-zero weighted original variables to create each principal component. By having a small number of the original variables constitute each principal component, we can more easily discuss groupings of variables. The simplest SPCA implementation first identifies principal components with traditional PCA. Each principal component can then be regressed using the original variables with a lasso penalty. We chose twenty as the number of non-zero variables to be included for each sparse principal component for ease of interpretability. The SPCA algorithm determined which environmental variables were chosen. We utilized the spca function in the elasticnet package from R [39].

The sparse principal components were used to determine how environmental variables were associated with asthma. Starting with the first sparse principal component, which represented the greatest variance of the ESRI dataset, sparse principal components were added sequentially to the Logistic thin plate regression spline model with covariates as shown below.

(2)log(asthmai,j)=f(xi,yi)+αjblockj+β1agei+β2sexi+β3racei+β4ethnicityi+β5(BMI)i+β6(encounterdays)i+β7(distance)i+∑m=1rδm(SPC)j,m where r = {1, …, 18} and (SPC)_j,m is the value of sparse principal component m at block group j. The largest model tested included the thin plate regression spline, covariates, and sparse principal components one through eighteen. Bayesian Information Criterion (BIC) was used to compare models without sparse principal components and with r = {1, …, 18}. Eighteen was chosen as the maximum number of sparse principal components we would be willing to investigate, as interpretability of environmental variables was a major goal. As models with increasing parameters can have a greater likelihood, BIC is a score used in model selection that penalizes the likelihood by the number of parameters. BIC = −2 * ln(L) + k * ln(N), where L is the likelihood, k is the number of parameters estimated and N is the number of observations [40]. The model with the lowest BIC is optimal.

We summarize the number of variables used in modeling. There are 1,117 environmental variables. Using sparse principal components analysis, 20 environmental variables were selected to represent each sparse principal component (SPC). By using SPCA, SPCs were ranked by importance based on the variance each SPCs represented from the original environmental variable dataset. To determine which SPCs to add to the model, we added the SPCs in order from rank #1 to rank #18. For example, we tested if the model was best fit if SPC 1 was added; if SPC 1 and 2 were added; if SPC 1, 2, and 3 were added; etc…; and if SPC 1 through 18 were added. The model also included 6 non-environmental covariates to control for variables that likely affect asthma diagnosis.

The change in log odds of asthma diagnosis per unit measure of sparse principal component m, δ_m(SPC)_j,m, was examined for each Wisconsin block group j. As (SPC)j=∑n=120ηn(EV)n where EV is an environmental variable from the ESRI database, the associated effect on the change in log odds of asthma diagnosis for an individual environmental variable could be assessed via the sign of η_n and η_m. All statistical analyses were performed in R.[34,41]

The graphical abstract summarizes the SASEA methods integrated in this study. We began with electronic health record data (covariates, asthma diagnosis as defined above, and the block group participants resided in) and environmental variables from Esri (values represent measurements from a block group). We applied sparse principal component analysis to the environmental variables. We combined the EHR dataset with the sparse principal components from the environmental variable dataset. We ran a logistic thin plate regression spline model on this combined dataset. Bayesian information criterion was used to select the number of sparse principal components added to the model. The odds ratios for variables in the logistic regression model were reported. The change in log odds value was color coded and mapped to block groups.

We attempted to account for multiple comparisons of the many variables and identify a smaller set of interpretable risk factors. The SASEA method performs sparse principal component analysis outside of the regression model as a means to prevent overfitting. Twenty non-zero loading variables for each sparse principal component were chosen to consider small groups of variables. Sparse principal components were sequentially added using BIC for model selection given the greater variance represented by higher ranked components. The sequential addition allowed for further structured and sparse variable evaluation. Although a set of sparse principal components were selected by BIC (four in this study), only some may be significant based on the odds ratio (two in this study). This feature of SASEA enhances sparsity as well.

The integrations of various scalable methods accommodated analysis of the EHR, environmental, and geographical datasets. Use of adaptable statistical model fitting based on well-studied algorithms was an asset that allowed for simple extension to the large number of patients and variables.

Similarly to our study, two additional studies [14,15] investigated community environmental variables associated with asthma. In our study, asthma was defined based on EHRs compared to survey data in the other two studies [14,15]. Many variables overlapped among these three studies. Our study and Krstić’s study [14] used latitude and longitude. We did not use insolation, air temperature or air pollution. Shankardass et al. [15] and our study had the individual variables of age, race, gender, and BMI. Shankardass et al. [15] included more individual variables including freeway distance while our study included more community environmental variables. We did not have male unemployment, which Shankardass et al. [15] found significantly associated with asthma. However we had other variables similar to socioeconomic status such as disposable income and employed civilian population in sparse principal components.

For analyses, Krstić [14] used linear regression, Shankardass et al. [15] used multilevel logistic random effect modeling, and our study used logistic thin plate regression spline modeling. The random effect modeling likely was more applicable to Shankardass et al. [15] as communities were concentrated. In our study the random effect in addition to thin plate regression spline based on latitude and longitude was chosen because of the distribution of patients throughout the state of Wisconsin.

As seen in other studies [6,40], higher asthma prevalence was associated with increased BMI, female sex, and black race, while lower asthma prevalence was associated with Hispanic ethnicity. Age, encounter days in the EHR dataset, and distance to most frequented clinic had little association with asthma diagnosis. Sparse principal component 2 represented by dog ownership and renter occupied housing units in addition to sparse principal component 4 represented by disposable income less than $15,000 were significantly associated with asthma. The individual variables representing sparse principal components likely contributed a small effect size.

Previous studies support the association of asthma and the environmental variables representing the sparse principal components in this study. In this study, dog ownership had a negative association with asthma. Other studies have shown perinatal and early life exposure to dog allergen was associated with reduced allergy and asthma risk later in life [43,44]. Renter occupied housing units were positively associated with asthma in a Brazil study [45]. Rental housing was associated with cold and damp housing, which in turn were associated with increased asthma [46]. Lastly, lower socioeconomic status as reflected by disposable income less than $15,000 was associated with greater asthma. Previously mentioned studies came to similar conclusion [5,6]. However, these results contradicted the positive association of socioeconomic status with asthma found in the Shankardass et al. [15]. Thus, the SASEA method used in this study identified variables that were previously associated with asthma risk, suggesting that these methods may have a role to studying chronic disease.

Mapping the associated change in log odds of asthma for a sparse principal component highlighted the geographic distribution of these sparse principal components and high loading environmental variables. The urban and rural discrepancy seen in differences in renter occupied housing units may be driven by the built environment, the human-made space where people live and work [47].

The use of EHR and block group characteristics merits comparison with traditional forms of health surveys including self-report and public health measured data. Canadian studies suggested census aggregate-level measures of income and education did not approximate individual level measures well [48–51]. There was similarity between self-reported variables and clinically measured variables. Self-reported colon cancer screening was similar to EHR imputed data [52]. Public health measured data were similar to EHR measured data. For example, BMI-based childhood obesity was 18% in both an EHR dataset and the National Health and Nutrition Evaluation Survey [53].

Agreement between disease prevalence based on health surveys and disease prevalence based on EHR datasets varies depending on disease. EHR datasets had prevalence similar to that from surveys for test-based conditions (e.g. diabetes) and decreased prevalence for minor conditions (e.g. back pain, headache, skin conditions) [54–56]. Specifically, two Spanish studies showed that the asthma prevalence calculated from an EHR dataset was lower compared with asthma prevalence calculated from population surveys [54,55]. However asthma prevalence based on UW eHealth-PHINEX (8.4%) was similar to the Wisconsin health survey, Behavioral Risk Factor Surveillance System (8.0%) [20]. As there is no single lab test for diagnosis of asthma, ICD-9 codes likely under-identify asthma when compared with “gold standard” manual record review [57] but may be more objective compared with population surveys.

Further analysis to determine the individual variables from sparse principal components that are associated with asthma could be performed using traditional methods such as stepwise model selection with BIC. This analysis could be performed with UW eHealth-PHINEX data from other years (e.g. 2009–2012), a UW eHealth-PHINEX hold out dataset, or a non-UW eHealth-PHINEX EHR dataset in another geographic region.

There are many future directions for this work regarding diseases and methods. Our methods could be applied to asthma control, other chronic diseases, and different communities. The census block groups and ESRI environmental data are already available nationwide. It is foreseeable that with the integration of a national EHR dataset, this type of analysis will be utilized to identify spatial risk factors to allow investigation or evaluation of interventions in any geographic region [58].

Alternative methods could have been used in this study. Traditional Logistic regression without the smoothing term does not account for the unknown orientation of spatial correlation among asthma due to geography, nor does it directly address difficulties in high dimensional data by constructing sparse models. Other spatial models included conditional auto-regressive models ^[59]. As the four sparse principal components accounted for a small percentage of variance from the original dataset, other methods such as traditional principal components analysis or clustering could have been utilized. However, traditional principal component analysis maintained all variables in each principal component preventing sparse interpretation, and clustering environmental variables added complexity. Few variables could have been associated more directly with the Logistic thin plate regression spline model using least absolute shrinkage and selector operator [30] such as COSSO [28]. However, a new set of variables would be identified for different diseases and variables could not be grouped. Allowing regression coefficients to vary over space as in geographically weighted regression [60] could be accomplished with spatial smoothing spline interaction terms.

There were limitations to the study. Although measures were taken to prevent overfitting and accommodate high dimensionality, this was an ecological, data-mining study without a priori variable hypotheses. This additive non-linear model likely does not fully capture the complexity of environmental factors influencing asthma.

The associations noted in the study may be due to confounding factors. One must be cognizant of ecological bias, because results about groups of people do not necessarily translate to the same findings about individuals. However, the neighborhood in which an individual lives in has been associated with health outcomes [61].

Multiple studies have shown the importance of EHR disease phenotype definitions, algorithm development, and validation [62–64]. In this study, asthma cases were defined based on ICD-9 codes. Some have argued this may under-estimate asthma prevalence [57]. Aside from this study’s EHR asthma phenotype definition, which is similar to the validated definition of Gershon et al. [65], other definitions such as the Healthcare Effectiveness Data and Information Set [66], have not been validated. We are currently validating alternative EHR phenotype definitions, which will also be used to segment asthma severity.

The results may be biased as UW eHealth-PHINEX data is not a complete representation of all block groups or persons in the state of Wisconsin. UW Family Medicine, Internal Medicine, and Pediatrics departments are an integrated health care system, but patients can receive care in at least two other major systems in the same catchment area. Because many hospitals and clinics have, or will soon have, an EHR system, sharing data through a statewide information exchange could mitigate this issue.

Another potential limitation is that the analysis included data elements from different years. While the EHR dataset represented years 2007–2009, the patient addresses were from the date of EHR data extraction in year 2012. However, compared with other states, Wisconsin residents tend to move less frequently. Wisconsin is the fifth “stickiest” state, with 68.6% of the current residents having been born in Wisconsin, an indicator of decreased residential mobility [67]. Patient addresses were geocoded to year 2000 block groups to match the ESRI database. ESRI database variables were mostly from year 2010, while some census variables were from year 2000. There is minimal change in block group from year to year and the goal was to identify general trends of larger geographic areas. The ESRI year 2010 variables were closest to the EHR database dates and the census year 2010 variables were not yet available.

Our main contribution is the incorporation of sparsity in spatial modeling. The sequential addition of sparse principal components to Logistic thin plate regression allowed interpretable analysis of geographically distributed EHR and environmental datasets. Understanding spatial disease variation and environmental risk factors using methods such as SASEA can allow better explanation of geographical disease disparity.