SAE methods can be broadly divided into design-based methods (estimates constructed from the sampling design) and model-based methods (estimates relying on a specified model) [12]. However, design-based methods fail for undercovered areas, and for these areas, only a model-based method can be utilized. Multilevel regression and poststratification (MRP) is a model-based SAE approach that was first used to estimate state-level indicators from national polls and has since been extended to track an array of health outcomes from traditional surveys, notably by the CDC’s PLACES project [11]. Additionally, the model-based nature of MRP can adjust the estimation from a nonrepresentative survey as a result of recruiting difficulties [13]. MRPs do not need to be restricted to just survey data — respresentative or not — but can also be applied to EHR data, which can be viewed as a highly nonrepresentative survey. We apply MRP techniques to five health conditions — asthma, diabetes, hypertension, obesity, and smoking — and to evaluate their performance against BRFSS’ 500 Cities estimates for Massachusetts cities and towns. We will first introduce the data sources employed in this study, and then describe the details of our MRP.
Data sources
MDPHnet: MDPHnet [5] is a distributed health data network in Massachusetts that permits authorized public health officials to submit detailed queries against the EHR data of three large multispecialty group practices that serve a combined patient population of approximately 1.1 million people (about 19% of the state population as of 2020). Each practice group uses open-source software (Electronic Medical Record Support for Public Health) [14] to create, host, and maintain a data repository behind its firewall in accordance with a common data specification. MDPHnet’s current practice partners include Atrius Health, a large multisite, multispecialty ambulatory group serving primarily a well-insured population (approximately 800,000 patients, 29 clinical locations), the Massachusetts League of Community Health Centers (MLCHC), a network of community health centers targeting underserved populations (approximately 500,000 patients, 15 clinical locations), and Cambridge Health Alliance (CHA), a combination safety net–general practice hospital and ambulatory group (approximately 200,000 patients, 15 clinical locations). The overall patient return rate across all 3 practice groups is 82% [15]. Analysis for this study was reviewed and approved by the institutional review board of Harvard Pilgrim Health Care Institute.
We analyzed 5 test conditions — diabetes, asthma, smoking, hypertension, and obesity — using EHR-data derived from MDPHnet (disease definitions are listed in Additional file 1). To facilitate comparisons against BRFSS 500 Cities estimates conducted in 2016, we restricted our analysis to anonymized, individual-level monthly data from MDPHnet from 2016 on individuals aged ≥20 with at least one outpatient visit of any kind in the health care system in the 2 years prior to each index month in 2016. The individual-level characteristics available within MDPHnet are sex, race, and age group. Because MDPHnet has measurements taken each month, we aggregated all 12 months of data for 2016 and included a time trend variable to account for temporal changes. Note that each provider will carry over mostly the same set of patients each month, although some patients may change providers and patient IDs were not provided for this analysis, making it impossible to account for within-person correlations over the months. Therefore, while we utilize all 12 months of data for estimation purposes, we inflated standard errors as if we had 1 month of data to provide conservative confidence intervals.
Community-Level Data: In addition to sex, race, and age compositions, the following seven sociodemographic data were derived from the ACS 2012–2016 5-year estimates for each of the 537 Zip Code Tabulation Areas (ZCTA) of Massachusetts: % never married, % single householder, % with bachelor’s degree, % English spoken at home, % unemployed, % receiving food stamps, and per capita income.
Hospital Discharge Data (“Case Mix”): MA hospital discharge data due to asthma, diabetes (with or without complications), and hypertension (including essential and secondary) for those aged 18 and older were obtained from Massachusetts’ case mix database. Because hospitalizations due to a particular chronic illness are far rarer than prevalent within the population, we developed a calibration model for the discharge data and blended the resulting estimates with our SAE estimators.
Behavioral Risk Factor Surveillance System Data: Our gold standard for validating our model estimates was the CDC’s Behavioral Risk Factor Surveillance System (BRFSS) MA statewide and 500 Cities estimates. For many years, the BRFSS has been a mainstay of public health surveillance for chronic diseases [16]. This telephone survey generates national and state-specific estimates of the prevalence of the major chronic diseases and risk factors that state and local public health departments rely upon to monitor health, plan interventions, and monitor their impact. Recent innovations in BRFSS methods in response to evolving surveillance needs include adding mobile telephone numbers to the sample, imputing measures for the 500 largest cities, incorporating area-level poverty as a predictor, and fielding follow-up surveys to gather clinical care information for conditions such as asthma, diabetes, hypertension, obesity, and smoking. Despite this, BRFSS is disadvantaged by the limited number of questions asked, the self-reported nature of the data, the lack of clinical data, declining response rates [17], and differential nonresponse over sex, age, race/ethnicity, income, and rurality [18].
MRP estimation procedure
We present comparisons of 2016 BRFSS statewide estimates from five candidate models:
-
M0: Crude estimates with no adjustments
-
M1: Adjustment for sex, race, age
-
M2: M1 + American Community Survey
-
\({\mathrm{M}}_1^{\mathrm{C}}\): M1 + calibrated Case Mix data
-
\({\mathrm{M}}_2^{\mathrm{C}}\): M2 + calibrated Case Mix data
To test the influence of separate practice groups’ influence on the estimators, we fit M0, M1, M2 individually for each of the three practice groups and in aggregate; the aggregate model pooled the three separate fits by weighting the estimates by relative coverage by provider per ZCTA. This relative coverage was computed as the coverage of each provider divided by the coverage of all three providers and therefore the three relative coverages sum up to 100% (see Step 3 below for more details). We included this relative coverage to account for the large differences in disease prevalence rates amongst different providers; therefore, provider affiliation was used as a proxy for additional behavioral and sociodemographic characteristics associated with disease risk. We note that for models M1 and M2, we also weighted each race, age, and sex strata by the corresponding census weight in a procedure known as poststratification [9]. We did not fit \({\mathrm{M}}_1^{\mathrm{C}}\) nor \({\mathrm{M}}_2^{\mathrm{C}}\) for each provider because the calibration procedure, as described shortly, needs sufficient coverage per provider over the available ZCTAs which was not the situation. The estimation steps are described briefly below, and more mathematical details provided in Additional File 2. All analyses were performed using R version 3.6.1.
-
Step 1 (Crude estimation of ZCTA-level disease prevalence per practice group): Direct prevalence estimates of the disease outcomes were computed based only off the available patients within each practice group in each ZCTA; many ZCTAs were inestimable due to lack of coverage from any of the three providers. This step forms M0 for Atrius, CHA, and MLCHC separately.
-
Step 2 (Adjusted estimation of ZCTA-level disease prevalence per practice group): A generalized linear mixed model (GLMM) was used to estimate ZCTA-level prevalence for each outcome within each practice group as a function of predictors. Model M1 included sex, race, and age while M2 included the full set of predictors listed under Community Level Data.
-
Step 3 (Pooling of crude and adjusted estimates from practice groups): Models M0, M1, M2 were then averaged over providers according to a relative coverage weighting, the weights were computed in each ZCTA as the number of patients subscribed to each provider divided by the total number subscribed to all three providers in that particular ZCTA. Stratum-specific weights were computed in a similar fashion. If data were not available for a stratum within a ZCTA, we rolled-up and replaced with the relative coverage for the overall ZCTA. If no data existed for the ZCTA, we rolled-up once again and used the overall relative coverage of the three data sources. This step forms M0, M1, M2 over all three providers.
-
Step 4 (Incorporation of case mix data with pooled): Hospitalization data was available for three of the five disease outcomes (asthma, diabetes, hypertension). These data contain rates of hospitalization due to a specified complication as measured by number of hospitalizations with that specified complication divided by the population of each ZCTA over the years 2011–2015. Hospitalizations do not equate to general prevalences; we therefore calibrated these crude hospitalization proportions such that the overall proportions for MA matched the overall predictions from M1 and M2 and the variation of proportions over ZCTAs matched the variation from M1 and M2 predictions. This formed \({\mathrm{M}}_1^{\mathrm{C}}\) and \({\mathrm{M}}_2^{\mathrm{C}}\), our calibrated case mix estimates. For the other two outcomes without available case mix (obesity and smoking), their model values are deferred to M1 and M2.
-
Step 5 (Roll-up from ZCTAs to municipalities and statewide estimates): The 537 ZCTAs of MA were aggregated to the 351 municipalities of MA and to the entirety of MA.
Ranking municipalities
Ranking was determined by the following two-step procedure. First, the interquartile range (IQR) outlier detection test retains municipalities if their prevalence is greater than Q3 + 1.5(Q3 − Q1), where Q1 is the first quartile (25 percentile) and Q3 is the third quartile (75 percentile). Next, this shortlist is ordered based off the lower bound in a 95% confidence interval for the true proportion in that municipality; this system is known as Wilson score ranking and accounts for high prevalences which may be due to small sample sizes.