Methods
Overall approach
The Diabetes LEAD Network is a Centers for Disease Control and Prevention (CDC)-funded research collaboration of the following academic centers: Drexel University, Geisinger and Johns Hopkins University, New York University Grossman School of Medicine, and the University of Alabama at Birmingham.28 The goal of the Network is to identify modifiable community-level determinants of T2D and cardiometabolic conditions using electronic health records (EHRs) and survey data from across the USA. Led by the Drexel University Data Coordinating Center, the Diabetes LEAD Network partners conducted analyses to reduce, where possible, between-study heterogeneity and aid interpretation of results. First, we created harmonized measurements of key neighborhood-level variables, including SES, PA, land use, racial/ethnic composition, and community type. Second, we developed complementary analysis plans to estimate mediation effects within three independent study populations, described in detail further, with different populations, study designs, and geographic contexts. Baseline address data were geocoded and used to link individuals to census tract-level community features.
Veterans Administration Diabetes Risk (VADR) cohort study
The VADR cohort study enrolled veterans in the national Veterans Health Administration EHR for primary care with sufficient residential address information within 2 years after cohort entry (1 January 2008).29 Participants were either (1) T2D-free at cohort entry with at least two T2D-free primary care visits at least 30 days apart within any 5 year period since 1 January 2003 or (2) subsequently enrolled and T2D free at cohort entry through 31 December 2016. Prevalent and incident T2D diagnoses were defined using encounter diagnoses, medication orders, and laboratory test results prior to and during the study period (online supplemental table S1). Subjects were followed from the date of cohort entry until censoring, defined as date of incident T2D, death, no encounters with the EHR for 2 years, or the end of the study (31 December 2018).
REasons for Geographic and Racial Differences in Stroke (REGARDS) cohort study
The REGARDS study enrolled adults age 45 years and older at baseline (2003–2007) from the contiguous USA with oversampling in the Southeast.30 This analysis included 11 208 study participants without prevalent T2D at baseline and who completed a follow-up exam in 2013–2016. Incident T2D was defined as a fasting glucose ≥126 mg/dL or random glucose ≥200 mg/dL or use of T2D medication at the follow-up exam (online supplemental table S1). Follow-up time was measured as the time between baseline and follow-up home visits.
Geisinger EHR nested case–control study
Geisinger is an integrated health system that provides primary, specialty, urgent, and emergency healthcare services at community practice clinics and hospitals in central and northeastern Pennsylvania. The design of the nested case–control study of new onset T2D in the Geisinger EHR has been described previously.31 Between 2008 and 2016, individuals with T2D (n=15 888) were identified using T2D encounter diagnoses, medication orders, and laboratory test results (online supplemental table S1). We required at least two encounters on different days with a primary care provider prior to ensure that we could detect T2D if present. To exclude prevalent T2D, we required individuals to have at least one encounter with the health system without evidence of T2D at least 2 years prior to the T2D onset date. We randomly selected with replacement five control encounters for each case (n=79 435, with 65 084 unique persons) from individuals who never met any of the criteria used to define T2D cases and frequency matched to cases on age, sex, and year of encounter.
Harmonized community type
The distributions of the primary exposure of interest, area-level neighborhood socioeconomic disadvantage, and the mediator variables of interest, density of PA facilities and distance to parks, vary considerably across urban, suburban, and rural community contexts. We assigned census tracts to one of four community type categories (higher density urban, lower density urban, suburban/small town, and rural) using a previously described modification of the US Department of Agriculture Rural-Urban Commuting Area methodology to better differentiate urban cores from surrounding non-rural areas.32 We stratified all analyses by community type due to concerns about census tract-level residual confounding and non-positivity, which can lead to bias when exposure variables do not overlap within strata of confounding variables.
Harmonized neighborhood socioeconomic environment (NSEE)
As part of the Diabetes LEAD Network, we developed a harmonized definition of area-level socioeconomic disadvantage relevant across the US, hereafter referred to as neighborhood socioeconomic environment (NSEE). We defined NSEE as a z-score sum of six census tract variables (percentage of persons with less than a high-school education, persons unemployed, households earning less than $30 000/year, population with income below poverty level, households on public assistance, and occupied housing units with no cars), based on the work of Xiao and colleagues,33 from the 2000 decennial Census and 2010 5-year American Community Survey. Census data from the year 2000 were converted to 2010 tract boundaries using the interpolation tool from Brown Universities Longitudinal Tract Data Base.34 In VADR, participants who entered the cohort before 2010 were assigned NSEE based on the 2000 Census data; otherwise, they were assigned NSEE based on 2006–2010 ACS data. In REGARDS, participants were assigned NSEE based on the 2000 Census data. In Geisinger, cases and matched controls in years 2008–2012 were assigned NSEE based on 2000 Census, while cases and matched controls in years 2013–2016 were assigned NSEE based on 2006–2010 ACS data. Higher NSEE z-scores indicated greater disadvantage. To improve interpretability of comparisons within a community type and across studies, we rescaled NSEE to range from 0 to 100 and categorized NSEE separately within each community type.
Harmonized density of PA facilities within street network buffers
We defined PA facility density as the 5-year average (ending the year prior to study entry or selection) of the number of commercial facilities for fitness or recreation per square kilometer within a street network buffer centered on the population-weighted centroid of a census tract. We obtained annual data (1998–2014) on PA facilities from the Retail Environment and Cardiovascular Disease (RECVD) study,35 which classified health-related neighborhood amenities using the National Establishment Time Series (NETS) Database. Details on classification methods and the NETS data, which was based on information collected by Dun and Bradstreet (D&B, Short Hills, New Jersey, USA), has been described elsewhere.35 The RECVD team re-geocoded the data to improve locational accuracy and recategorized facilities using Standard Industrial Classification codes and word searches or name searches to enhance classification. Using ESRI ArcGIS9.3, different street network buffers were calculated depending on the LEAD community type: 1 mile (1.6 km) walking buffer in high density urban, 2 mile (3.2 km) driving buffer in lower density urban, and 6 mile (9.7 km) driving buffers in suburban/small town and rural. Buffer sizes were selected based on prior literature measuring the commercial PA facility environment, although few past studies could inform appropriate buffer sizes in rural areas.36–40
Harmonized population-weighted distance to seven closest parks
We defined neighborhood spatial access to parks for each census tract using a population-weighted measure of the closest seven parks, as described previously.25 41 Briefly, this measure accounts for distance (miles) and size of the seven closest parks while adjusting for the heterogeneous population distribution within a census tract. Data on national, state, and local parks in 2010, provided by the CDC Division of Population Health, were derived from two sources: The Homeland Security Infrastructure Program Gold database and ESRI ArcGIS 10.1 Data DVD.
Harmonized census tract-level covariates
We measured percent Hispanic and percent non-Hispanic Black and land use environment within each census tract as potential confounders. We obtained race/ethnicity data from the 2000 and 2010 decennial census. To measure the land use environment, we used factor scores from a multiple group confirmatory factor analysis, stratified by community type, which identified a single latent variable from seven components of the built environment: average block length, average block size, intersection density, street connectivity, establishment density, percent developed land, and household density. Data were derived from ESRI 2009 Vintage Street and computed via ArcGIS Pro 2.3. Details on the land use measurement model have been described previously.42 The z-score sums within each community type were scaled to range from 0 to 100, with a higher value indicating a more walkable neighborhood.
Individual-level covariates
Individual-level covariates, including age, sex, race and ethnicity, smoking status, and available proxies for individual or family SES, were defined within each study. Age was defined at baseline in VADR and REGARDS or at the time of case or control selection in Geisinger. We considered race/ethnicity as a proxy for race-based discrimination, a social construct that is correlated with social determinants of health. Categories of race/ethnicity were defined within each study (VADR: non-Hispanic White, non-Hispanic Black, Hispanic, Asian, and other/unknown; REGARDS: non-Hispanic White vs non-Hispanic Black; Geisinger: non-Hispanic Black, non-Hispanic White, or other race/ethnicity). SES was defined by baseline data on individual disability and income in VADR (disabled, low-income/non-disabled, or neither), as baseline household income in REGARDS (<$20 000, $20 000–$34 999, $35 000–$74 999, ≥$75 000, or refused to answer) and by receipt of medical assistance prior to case or control selection (0% vs >0% of time prior to case or control selection) in Geisinger.
Statistical analysis
Using a counterfactual framework, we estimated the average natural direct effect of NSEE on risk of T2D, not operating through the two hypothesized mediators, and the average natural indirect effect of NSEE on risk of T2D operating through the LTPA environment, adjusted a priori for confounding variables at the individual-level and census tract-level (online supplemental figure S1). When interpretable (ie, when direct and indirect effects were in the same direction), we calculated the proportion mediated (ie, the proportion of the total effect operating through the indirect pathway of interest).43 No sample size calculations were conducted.
In the VADR cohort study, we first examined exposure–outcome and mediator–outcome associations using piecewise exponential survival models44 with 2-year intervals and county random effects fit in R (V.4.04) using the ‘lme4’ R package.45 County random effects were used because of convergence issues when clustering by census tract. We used generalized linear mixed effects Poisson regression models with a log link function and an offset of logarithm of time-at-risk during each interval to estimate HRs assuming a constant hazard function within intervals over time. We estimated total, direct, and indirect effects46 as differences in 2-year T2D incidence rates, with quasi-Bayesian approximation 95% CIs, at each NSEE quartile compared with the first quartile, using the ‘mediation’ R package.47 Differences in incidence rates were multiplied by 100 for interpretability.
In the REGARDS cohort study, we examined exposure–outcome associations using Poisson mixed models with robust variance estimation to account for correlation of participants within census tracts, and exposure–mediator associations using linear mixed models with robust variance estimation in R (V.4.04) using the ‘lme4’ package.45 We estimated total, direct, and indirect effects46 as differences in T2D incidence rates, with bootstrapped 95% CIs, at each NSEE quartile compared with the first quartile, using the ‘mediation’ R package.47 Differences in incidence rates were multiplied by 100 for interpretability.
In the Geisinger EHR case–control study, we estimated total, direct, and indirect effects using logistic regression models as ORs, with bootstrapped 95% CIs,48 at each NSEE quartile compared with the first quartile, in R (VV.4.04) using the ‘medflex’” R package49 in order to fit models appropriate for our study design and to use tract-level bootstrap resampling to account for clustering within census tracts.
Final models were adjusted for Network-harmonized census tract-level race/ethnicity and land use environment (factor score) and study-specific age, sex, race/ethnicity, smoking status (REGARDS and Geisinger only), and SES. We treated age as a continuous variable with linear and quadratic terms. All VADR models included linear and quadratic age terms, while REGARDS and Geisinger models included quadratic age if the quadratic term was statistically significant (p<0.05). In VADR, we did not adjust for smoking status because of the large proportion of missingness (66%) at cohort entry.
Sensitivity analyses
We conducted sensitivity analyses to evaluate the robustness of the main effects. First, we refitted mediation models for PA facility density with larger street network buffers (2 mile driving buffers in higher density urban, 6 mile driving buffers in lower density urban, and 10 mile driving buffers in suburban/small town and rural communities). Second, we removed land use environment from models because this variable could be on the causal pathway. Third, we took advantage of the national VADR cohort to restrict VADR models to census tracts shared with either REGARDS or Geisinger. Although the effects across studies cannot be quantitatively compared due to differences in study design and analytic approaches, we qualitatively evaluated the consistency (eg, presence vs absence) of mediation results. We hypothesized that comparisons between primary mediation results from REGARDS or Geisinger and the geographically restricted VADR models, which compared different populations within the same census tracts, could suggest whether the observed discrepancies in primary findings across sites could partly be attributed to differences in geographic context or in population composition. Fourth, we conducted a post hoc analysis in the REGARDS rural community model in which we treated the PA facility density mediator as a binary variable, dichotomized at the median, to investigate whether the observed associations in REGARDS rural communities of mediation via PA facility density could be due to a misspecification of the PA facility density mediator as a continuous variable.