6 Analyzing your data

By: Jarvis Chen, ScD

6.1 Overview of Methods

In the years since we initially presented the Public Health Disparities Geocoding Training in 2004, there has been a huge increase in conceptual and methodological work regarding use of ABSMs to document and analyze health inequities, and the computing tools and statistical methods to conduct these analyses. In particular, the availability of software that facilitates data access, mapping, visualization, and fitting of multilevel and spatial models has greatly enhanced the accessibility of these analytic methods for public health scientists and advocates interested in advancing health equity. Our goal in presenting these updated methods is to illustrate their applicability to monitoring, reporting, and investigating health inequities. In doing so, we highlight many of the assumptions underlying the use of these methods and show how they can be applied to real-life data that arise in public health surveillance and health equity monitoring. We comment on some of the pitfalls that may arise when working with imperfect data in real time and discuss analytic and interpretive strategies keeping in mind the overarching goal of accurately documenting health inequities for accountability. Where appropriate, we also comment on ongoing methodologic work to improve analytic approaches.

6.2 Motivating Questions

As the choice of methods for any data analysis depends fundamentally on the goals of the analysis, we will focus our presentation on the following motivating questions:

1. “What is the geographic distribution of area-based social metrics across the study area?” This is of particular interest for exploratory analyses in order to understand who in what areas are most affected by social advantage and disadvantage as captured by area-based social metrics. Mapping and visualization tools can be used to communicate geographic patterns and to elicit local knowledge in order to spur hypothesis generation and prioritization of research questions.

2. “What are the social inequities in health associated with area-based social metrics across the study area?” We highlight the importance of descriptive epidemiology as a first step in understanding inequities in the health of populations living in areas characterized by area-based social metrics. That is, accurate description of where the burden of social inequities falls is a necessary pre-requisite to more causally-oriented investigations or evaluations of interventions. We discuss aggregated analyses, which have the advantages of being straightforward to implement and report using tabulation methods and avoiding the problems of small area estimation, as well as non-spatial regression models. We focus on health equity settings in which understanding the joint patterning of inequities by racialized group and area-based social metric are of interest, and comment on how analyses may need to consider heterogeneity across age strata.

3. “How can geographic variation in health outcomes be modeled in order to facilitate mapping of small-area disease estimates and estimation of social inequities?” We review concepts from the extensive statistical literature on small-area estimation and disease mapping with the goal of visualizing geographic patterns in health outcomes for studies of health equity. While these methods entail more complex modeling frameworks and computational details, we focus on their application to health equity and the interpretation and use of model outputs in the context of monitoring and reporting of health disparities. We consider settings in which multilevel modeling approaches and spatial modeling approaches may be preferred and discuss considerations contributing to the choice of modeling framework.

6.3 Choice of geographic level

• A pragmatic question and approach: at what level(s) are data available?

• For small area level analyses, do the number of events support small area-level analyses?

• census tract (CT) vs. city/town, e.g. breast cancer mortality analysis

• ZIP Codes and US Census defined ZIP Code Tabulation Areas (ZCTAs)

• Modifiable Areal Unit Problem¹

• note that the smallest level of Census geography with intercensal population estimates is the county. The US Census Population Estimates Program (PEP) estimates rely on demographic modeling.

• CT level population estimates are in the US Census American Community Survey (ACS), but ACS specifically cautions against using these as population denominator estimates in small area estimation. However, we generally do not have any other sources of population denominators on which to rely. Five year average estimates may help, but there is still year to year variability. Margin of error estimates are available, but how to incorporate denominator uncertainty into analyses is still an area of active research.

• For a review of issues connected to numerator/ denominator mismatch, please visit Section 4.7 in the “Getting your Data” chapter.

6.3.1 Example: Age-specific all-cause mortality by racialized group in Massachusetts, 2013-2017

We compared the use of the aggregation method to a non-spatial Poisson model for analyzing deaths from all causes by age and racialized group. Using the aggregation method, deaths and population person-time at risk are aggregated into strata by age (0-4, 5-9, 10-14, 15-19, 20-24, 25-29, 30-34, 35-44, 45-54, 55-64, 65-74, 75-84, 85+) and racialized group (White non-Hispanic and Black), and mortality rates per 100,000 person-years are computed using the formulas listed in Section 6.5.1. For the non-spatial Poisson model, the inputs are age- and racialized-group specific deaths for each census tract in the study area, the log(population person-time at risk) as an offset, and a set of indicators for each stratum of age category by racialized group.

When the log(population person-time at risk) variable is zero, however, this presents a problem for the Poisson model-fitting algorithm as the rate is undefined. Census tracts where there are zero deaths and zero person-time at risk in a particular strata can be deleted from the dataset as these observations contribute no information. But situations where the count of deaths in a stratum is greater than zero but the population person-time at risk is zero continue to be an issue (numerator/denominator mismatch). It may be tempting to delete these observations from the dataset as well, reasoning that as these usually involve just one death occurring in a stratum where the population estimate is zero, deleting these observations will not unduly affect the analysis.

Unfortunately, it turns out that numerator/denominator mismatch occurs more frequently for Black observations in the dataset, due to patterns of racialized residential segregation and that fact that the Black population accounts for a much smaller proportion of the population compared with the White non-Hispanic population (10.1% vs. 72.9% in Massachusetts according to the ACS population estimates we used for our population denominators). That is, at the census tract level, there are more likely to be age strata where there are deaths reported by zero population person-time at risk for the Black population compared with the White non-Hispanic population. As a result, when we delete these observations, we end up deleting a larger proportion of the Black deaths across the whole state, which has an impact on our estimates of rates and rate ratios.

To see this, consider the following table where we show the number of deaths and population by age and racialized group aggregated across census tracts. In columns 2-8, we show the deaths, population, and estimated rates per 100,000 person-years for White Non-Hispanics and Blacks in the full dataset, while in columns 9-15 we show the deaths and population after deleting observations where the numerator is greater than zero but the denominator is zero. In columns 16-22 we show the percent bias by strata comparing the dataset with deleted observations to the full dataset. As is evident in column 19, the effect of deleting these observations is to reduce the Black death count in age strata (between 4-30% across strata) to a much greater degree than among the White non-Hispanic age strata (between 0-2%). As a result, the age-specific incidence rate ratios can be depressed by as much as 30% relative to the IRRs calculated using the full aggregated data.

	Full dataset							Adjusted dataset							% Bias
Age	NHW deaths	NHW pop	NHW Rate	Black deaths	Black pop	Black Rate	IRR	NHW deaths	NHW pop	NHW Rate	Black deaths	Black pop	Black Rate	IRR	% bias NHW deaths	% bias NHW pop	% bias NHW Rate	% bias Black deaths	% bias Black pop	% bias Black Rate	% bias IRR
0-4	769	1068730	5.0	316	177968	12.3	2.47	766	1068730	5.0	277	177968	10.8	2.17	0.996	1	0.996	0.877	1	0.877	0.880
5-9	90	1147677	0.6	14	172883	0.6	1.03	88	1147677	0.6	12	172883	0.5	0.91	0.978	1	0.978	0.857	1	0.857	0.877
10-14	114	1301178	0.6	29	176597	1.2	1.87	114	1301178	0.6	25	176597	1.0	1.62	1.000	1	1.000	0.862	1	0.862	0.862
15-19	371	1533896	1.7	100	198678	3.6	2.08	371	1533896	1.7	88	198678	3.2	1.83	1.000	1	1.000	0.880	1	0.880	0.880
20-24	1103	1611714	4.5	163	224987	4.8	1.06	1102	1611714	4.5	155	224987	4.6	1.01	0.999	1	0.999	0.951	1	0.951	0.952
25-29	1788	1615532	7.1	244	212885	7.4	1.04	1787	1615532	7.1	218	212885	6.6	0.93	0.999	1	0.999	0.893	1	0.893	0.894
30-34	2079	1518953	9.7	213	183138	8.3	0.85	2078	1518953	9.7	194	183138	7.5	0.77	1.000	1	1.000	0.911	1	0.911	0.911
35-44	4846	2873290	27.4	517	334861	25.1	0.92	4846	2873290	27.4	496	334861	24.1	0.88	1.000	1	1.000	0.959	1	0.959	0.959
45-54	12386	3756899	44.5	1165	325584	48.2	1.09	12386	3756899	44.5	1113	325584	46.1	1.04	1.000	1	1.000	0.955	1	0.955	0.955
55-64	25334	3735786	59.2	1997	262045	66.5	1.12	25326	3735786	59.1	1899	262045	63.2	1.07	1.000	1	1.000	0.951	1	0.951	0.951
65-74	38765	2527556	101.3	2282	139490	108.0	1.07	38752	2527556	101.2	2081	139490	98.5	0.97	1.000	1	1.000	0.912	1	0.912	0.912
75-84	58096	1335819	195.0	2374	67943	156.7	0.80	58057	1335819	194.9	1974	67943	130.3	0.67	0.999	1	0.999	0.832	1	0.832	0.832
85+	101687	711624	221.6	2571	23810	167.5	0.76	101585	711624	221.4	1787	23810	116.4	0.53	0.999	1	0.999	0.695	1	0.695	0.696

If instead of deleting those strata where deaths>0 and denominator=0, we replace the denominator with the number of deaths, we end up with very slightly larger population counts overall, but the effect on the rates is to bring them more in line with the aggregated analysis. Here is the comparable table comparing the full aggregated dataset (columns 2-8) with a dataset in which the denominators are adjusted to increase the person time at risk to equal the number of deaths the affected strata.

	Full dataset							Adjusted dataset							% Bias
Age	NHW deaths	NHW pop	NHW Rate	Black deaths	Black pop	Black Rate	IRR	NHW deaths	NHW pop	NHW Rate	Black deaths	Black pop	Black Rate	IRR	% bias NHW deaths	% bias NHW pop	% bias NHW Rate	% bias Black deaths	% bias Black pop	% bias Black Rate	% bias IRR
0-4	769	1068730	5.0	316	177968	12.3	2.47	769	1068733	5.0	316	178007	12.3	2.47	1	1	1	1	1.000	1.000	1.000
5-9	90	1147677	0.6	14	172883	0.6	1.03	90	1147679	0.6	14	172885	0.6	1.03	1	1	1	1	1.000	1.000	1.000
10-14	114	1301178	0.6	29	176597	1.2	1.87	114	1301178	0.6	29	176601	1.2	1.87	1	1	1	1	1.000	1.000	1.000
15-19	371	1533896	1.7	100	198678	3.6	2.08	371	1533896	1.7	100	198690	3.6	2.08	1	1	1	1	1.000	1.000	1.000
20-24	1103	1611714	4.5	163	224987	4.8	1.06	1103	1611715	4.5	163	224995	4.8	1.06	1	1	1	1	1.000	1.000	1.000
25-29	1788	1615532	7.1	244	212885	7.4	1.04	1788	1615533	7.1	244	212911	7.4	1.04	1	1	1	1	1.000	1.000	1.000
30-34	2079	1518953	9.7	213	183138	8.3	0.85	2079	1518954	9.7	213	183157	8.3	0.85	1	1	1	1	1.000	1.000	1.000
35-44	4846	2873290	27.4	517	334861	25.1	0.92	4846	2873290	27.4	517	334882	25.1	0.92	1	1	1	1	1.000	1.000	1.000
45-54	12386	3756899	44.5	1165	325584	48.2	1.09	12386	3756899	44.5	1165	325636	48.2	1.09	1	1	1	1	1.000	1.000	1.000
55-64	25334	3735786	59.2	1997	262045	66.5	1.12	25334	3735794	59.2	1997	262143	66.5	1.12	1	1	1	1	1.000	1.000	1.000
65-74	38765	2527556	101.3	2282	139490	108.0	1.07	38765	2527569	101.3	2282	139691	107.9	1.07	1	1	1	1	1.001	0.999	0.999
75-84	58096	1335819	195.0	2374	67943	156.7	0.80	58096	1335858	195.0	2374	68343	155.8	0.80	1	1	1	1	1.006	0.994	0.994
85+	101687	711624	221.6	2571	23810	167.5	0.76	101687	711726	221.6	2571	24594	162.1	0.73	1	1	1	1	1.033	0.968	0.968

An even better solution, and the one we ultimately recommend, is to add a small number (e.g. 0.001) to all of the population person-time denominators to allow the Poisson model fitting algorithm to incorporate these observations into the analysis. The resulting effect on the overall population denominators is negligible, but this allows for estimation of rates and rate ratios of interest. Note that this is only needed when fitting non-spatial regression models to the data. The aggregation method does not incur this problem as small discrepancies are averaged out over areas, and the spatial models we discuss in Section 6.6 address the problem of zero or infinite rates by smoothing.

To visualize the effect of these adjustments on estimates of the age-specific disparity by racialized group, we plot estimates from aggregated analyses and Poisson models with deletion and denominator adjustment in Figure 1.1. The plot confirms that either (a) adjusting denominators to match the number of observed deaths in problematic strata in census tracts or (b) adding a small number to all population person time estimates (preferred) yields comparable estimates to the aggregation method, whereas removing problematic observations results in bias.

Figure 6.1: Comparison of estimated age-specific Black/White mortality rate ratios using different methods to address numerator/denominator mismatch

6.4 Aggregation Method

As we presented in the original Public Health Disparities Geocoding Project Monograph (Krieger et al, 2004), one of the most straightforward ways to incorporate area-based social metrics into analyses of health inequities is to use what we term the Aggregation Method. This method entails using geocodes to append area-based social metrics to health surveillance records, to stratify these records into discrete categories based on ABSM, and to aggregate numerators and denominators over areas, within levels defined by ABSM. Rates, including age-standardized rates, and measures of association (e.g. rate differences or rate ratios) can be easily computed using tabulation methods and formulae that are taught in all basic epidemiologic textbooks. The analyses are straightforward and do not require specialized software: they can even be done in Excel or other spreadsheet programs. One of the key advantages is that aggregation avoids problems with numerator/denominator mismatch since small discrepancies in case counts and population at risk tend to get averaged out in aggregating numerators and denominators over small areas within ABSM strata.

We can describe inequities by categories of ABSMs by aggregating deaths and population at risk from census tracts that share the same values of those ABSMs. This is analogous to what health departments typically do when reporting, for example, statewide cancer mortality rates by age or gender.

This method avoids the problem of unstable rates arising from small areas by assuming that cases and population denominators from areas with similar socioeconomic characteristics can be validly combined into the same strata.

The following steps are used to generate age-standardized disease rates stratified by area-based social metrics once the case data have been geocoded and appropriate ABSMs have been generated from census data.

1. Aggregate the case data into numerators (age cells within areas/geocodes).
2. Aggregate population denominator data into age cells within areas/geocodes.
3. Merge the numerators and denominators with ABSMs, by area/geocode.
4. Aggregate over areas into strata defined by categorical ABSM and age category.
5. Generate age-standardized rates and other summary measures.

Advantages:

• easy to implement using tabulation methods; can even be done in MS Excel

• formulae are taught in basic epidemiologic textbooks

• presentation and interpretation of social inequities is straightforward

• aggregation avoids problems with numerator-denominator mismatch

• age-standardization using the direct method has coherent interpretation even in the presence of effect heterogeneity by age

• Confidence intervals may be wider particularly when age strata with relatively less information get upweighted

Disadvantages:

• precludes identification of small areas with unusually high or low disease rates

• may preclude further adjustment for compositional differences or area-level covariates that vary within ABSM strata

Throughout this section, we will be using data from the Greater Boston Area for the years 2013-2017, drawing on Massachusetts Mortality data and the American Community Survey (ACS) population data.

6.4.1 Direct Age Standardization

• rate estimation
• confidence intervals on the rates
• measures of disparity: incidence rate difference and incidence rate ratio
• population attributable fraction
• a comment on life expectancy

CT % below poverty	Age category	Deaths	Person-years	Weight	Rate per 100,000 p-y
0-4.9%	0-4	27	33162	0.069	5.63
0-4.9%	5-9	3	31593	0.073	0.69
0-4.9%	10-14	1	33800	0.073	0.22
0-4.9%	15-19	7	38060	0.072	1.33
0-4.9%	20-24	12	30408	0.066	2.62
0-4.9%	25-29	20	45937	0.065	2.81
0-4.9%	30-34	31	43138	0.071	5.11
0-4.9%	35-44	51	77292	0.163	10.73
0-4.9%	45-54	145	85527	0.135	22.86
0-4.9%	55-64	330	81124	0.087	35.49
20-100%	0-4	156	107867	0.069	10.00
20-100%	5-9	15	94231	0.073	1.15
20-100%	10-14	11	94853	0.073	0.85
20-100%	15-19	64	209056	0.072	2.21
20-100%	20-24	115	308641	0.066	2.48
20-100%	25-29	168	254985	0.065	4.25
20-100%	30-34	176	182451	0.071	6.85
20-100%	35-44	422	232327	0.163	29.54
20-100%	45-54	867	208522	0.135	56.06
20-100%	55-64	1580	177590	0.087	77.62

CT % below poverty	Deaths	Person-years	Standardized Rate per 100,000 p-y	Var(Rate)	Sum of weights	Sum of weights^2	std_rate_lo95	std_rate_up95
0-4.9%	627	500041	100	0	0.873614	0.086	100.1346	100.1
20-100%	3574	1870523	219	0	0.873614	0.086	218.6464	218.6

6.5 Non-Spatial Regression Methods

6.5.1 Poisson regression

While the aggregation method is attractive for its conceptual simplicity and ease of implementation, some may prefer to take a regression approach to analysis of health inequities by ABSM. As we will describe below, this makes it possible to relax the strong assumption of a homogenous Poisson process within category of ABSM. We might therefore choose to model the data (deaths and population at risk by age stratum within census tracts) using a Poisson loglinear model – that is, a generalized linear model with a Poisson error distribution and a log link.

Let \(Y_{ij}\) be the count of deaths in age-stratum \(j\) in census tract \(i\) and \(n_{ij}\) be the corresponding person-time at risk. We fit a Poisson log linear model to the data and include dummy variables for ABSM categories, e.g. with \(x_1\), \(x_2\), \(x_3\), and \(x_4\) coded 1 if the census tract is in the corresponding category of CT % below poverty and 0 otherwise. We model age categories as a set of dummy variables. We also include \(\log(n_{ij})\) as an offset in the model for \(\log(\mu_{ij})\). \[ \begin{aligned} Y_{ij} & \sim \mbox{Poisson}(\mu_i) \\ \log(\mu_{ij}) & = \beta_0 + \beta_1 I(pov=5-9.9\%) + \beta_2 I(pov=10-19.9\%) + \beta_3 I(pov=20-100\%) \\&+\sum_{j=2}^{J} \alpha_j I(age=age_j) + \log(n_{ij}) \end{aligned} \]

6.5.2 Quasi-poisson regression

A key assumption of the Poisson distribution is that the variance of the expected count equals the mean, \[ \mbox{Var}(Y_{ij}) = \mathbb{E}(Y_{ij}) = \mu_{ij} \] In real life, however, we often find data where the empirical variance is larger than the mean. This is known as overdispersion.

To diagnose potential overdispersion, we can look at the standardized residuals, defined as \[ \begin{align*} z_{i} & = \frac{y_{i} - \hat{y}_{i}}{sd(\hat{y}_{i})} \end{align*} \] If the Poisson model is true, then the \(\{z_{i}\}\)s should be approximately independent, each with mean 0 and standard deviation 1. If there is overdispersion, we would expect the \(\{z_{i}\}\)s to be larger in absolute value, reflecting the extra variation beyond what is predicted under the Poisson model.

We can test for overdispersion by computing the sum of squares of the standardized residuals and comparing this to the \(\chi^{2}_{n-k}\) distribution where \(n\) is the sample size and \(k\) is the number of parameters in the Poisson log linear model.

> df_overdispersion <- df_forModel 
> df_overdispersion$yhat <-  predict(poisson_apINDPOV,
+                                          newdata=df_forModel %>%
+                                          dplyr::select(denominator, 
+                                                        apINDPOV_2, 
+                                                        apINDPOV_3, 
+                                                        apINDPOV_4,
+                                                        apINDPOV_NA,
+                                                        agecat), type="response")
> 
> df_overdispersion$z <- (df_overdispersion$numerator - df_overdispersion$yhat)/sqrt(df_overdispersion$yhat)
> n.obs <- length(df_overdispersion$numerator)
> k <- 4
> cat("overdispersion ratio is ",sum(df_overdispersion$z^2, na.rm=T)/(n.obs-k),"\n")
overdispersion ratio is  2.144269 
> cat("p-value of overdispersion test is ", pchisq(sum(df_overdispersion$z^2, na.rm=T), n.obs-k, lower.tail=FALSE), "\n")
p-value of overdispersion test is  0 

We can also visualize the predicted counts vs. the standardized residuals.

Figure 6.2: Predicted counts vs. standardized residuals. The standardized residuals should have mean 0 and standard deviation 1 (hence the dashed lines at \(\pm2\) indicating approximate 95% error bounds). The variance of the standardized residuals is larger than 1, indicating in this case a moderate amount of overdispersion.

6.5.3 Negative binomial regression

• aggregation method and Poisson regression coincide in the analysis of crude rates
• Poisson regression assumes that the variance is equal to the mean, but real-life data often exhibit overdispersion (actually more complicated with age-specific data, where some age strata may be underdispersed and others may be overdispersed)
• Quasi-poisson regression estimates an extra scale parameter: results in wider confidence limits when there is overdispersion. ABSM estimates are identical to the Poisson model fit.
• Negative binomial model assumes a different distribution to allow for overdispersion. Can be conceived of as a mixture of Poisson distributions where the latent variable is gamma distributed
• Negative binomial model yields ABSM estimates that can vary from the Poisson and Quasi-Poisson fits

How do these models handle numerator/denominator mismatch? Zero numerator/zero denominator observations need to be deleted. Non-zero numerator/zero denominator areas present a problem. If we delete these areas, the undercount of cases can be differential by racialized group; a better approach is to add a small number (e.g., 0.1) (the effect on the total denominators is negligible, but this yields rates and ABSM estimates more similar to the aggregation method).

6.5.4 Comparison of non-spatial regression estimates

Figure 6.3: Comparison of estimates of social inequities by ABSM from the aggregated method and from non-spatial regression models

6.6 Small Area Estimation

6.6.1 Indirect Age Standardization

We observe \(y_{ij}\), the count of deaths in area \(i\) in age stratum \(j\). We define \(O_{i}\) as the observed number of deaths in area \(i\), summed up over the age strata: \(O_{i} = \sum_{j} y_{ij}\). The expected number of deaths, \(E_{i}\), is computed by applying age-specific mortality rates from a reference population to the age-specific population counts in each area: \[ E_{i} = \sum_{j} N_{ij} \times R_{j} \] where \(N_{ij}\) is the person-time at risk in age group \(j\) in area \(i\) (computed by taking the population in age-group \(j\) in area \(i \times\) number of years) and \(R_{j}\) is the mortality rate age group \(j\) of a suitable reference population. The ratio \(O_{i}/E_{i}\) is known as the Standardized Incidence Ratio (SIR), or, in the case of mortality outcomes, the Standardized Mortality Ratio (SMR). This ratio is an estimate of the within each area, \(\theta_{i}\). For rare, non-communicable diseases, the standard statistical model for \(O_{i}\) is that of the Poisson distribution: \[ O_{i} \sim \mathrm{Poisson} (\theta_{i}E_{i}). \] The maximum likelihood estimator of \(\theta_{i}\) is \[\begin{align*} \hat{\theta}_{i} = \mathrm{SMR}_{i} & = \frac{O_{i}}{E_{i}} \end{align*}\] with \(\mathrm{Var}(\mathrm{SMR}_{i}) = \theta_{i}/E_{i}\), estimated by \(O_{i}/E_{i}^2\). Note that \(\mathrm{Var}(\mathrm{SMR}_{i})\) is inversely proportional to \(E_{i}\). To translate the SMR into an age-standardized rate, we can multiply the SMR by the overall mortality rate from the reference population. In these examples, however, we will present the small area estimates on the SMR scale, since we are interested in the pattern of increased or decreased relative risk across small areas.

Indirect age standardization and SMRs are particularly useful when the age-specific counts of deaths are not available for each area, since all that is required are age-specific population estimates for each area, a set of age-specific reference rates, and the total cases in each area. Interestingly, this method also produces rate estimates with smaller asymptotic variance than the corresponding direct standardization method (Pickle and White, 1995).

As summary measures, SMRs also have certain drawbacks. They are based on ratio estimators, and are thus sensitive to small changes in \(E_{i}\). In particular, when \(E_{i}\) is close to zero, the SMR will be very large for any positive count. As the estimate of \(\mathrm{Var}(\mathrm{SMR}_{i})\) is proportional to \(1/E_{i}\), SMRs of zero do not distinguish variation in expected counts. Most importantly, interpretability and comparability of SMRs based on indirectly age standardized data across areas depends on the assumption of independent area and age effects with respect to the standard population. This is known as the proportionality assumption, and assumes that \(r_{ij} = \theta_{i} \times \alpha_{j}\) where \(\alpha_{j}\) is an age effect that does not vary over areas.

To determine whether an \(SMR>1\) is inconsistent with a chance occurrence, we can carry out a statistical hypothesis test: \[ \begin{align*} \mbox{Null hypothesis} & : \theta_{i}=1 \\ \mbox{Alternative hypothesis} & : \theta_{i}>1 \end{align*} \] We test this hypothesis as follows: Under the null hypothesis, the mean of the Poisson distribution is \(E_{i}\) \[ O_{i} \sim \mathrm{Poisson}(E_{i}) \]

Calculate the probability of obtaining at least \(O_{i}\) cases by chance from a Poisson distribution with mean \(E_{i}\): \[ \begin{align*} \Pr(X \geq O_{i} | \mathrm{Exp}=E_{i}) & = 1 - \Pr(X < O_{i} | \mathrm{Exp}=E_{i}) \\ & = 1 - \sum^{O_{i}-1}_{X=0} \frac{E_{i}^{X} e^{-E_{i}}}{X!} \\ & = p \end{align*} \] This probability \(p\) is the (1-sided) p-value. If \(p \leq 0.05\), we usually reject the null hypothesis and accept the alternative. We say that the excess risk of disease in area \(i\) is statistically significant at the 5% level. Such a finding might warrant further epidemiological investigation.

Consider, however, that we have 306 census tracts in the Greater Boston mortality example. The statistical test described above was just for the comparison of one census tract’s SMR relative to the null hypothesis. If we want to identify all of the census tracts in our study area with significantly elevated rates, we would have to repeat this tests many times, which would create a multiple hypothesis testing situation. Moreover, we could imagine that there would be dependence of tests for nearby areas if neighboring areas shared some risk factor which increased the SMR for these clusters of census tracts.

To overcome this variability, hierarchical models can be used to “smooth” the raw rates. When faced with the problem of making inferences on many parameters \(\{\theta_{i}\}=\theta_{1}, \ldots, \theta_{n}\), measured on \(n\) areas, one can imagine two possible extreme assumptions:

• We could assume that all of the \(\{\theta_{i}\}\) are identical, in which case all the data can be pooled, and the individual units ignored. This is what we typically do when presenting summary rates and rate ratios over the whole study area.

• At the other extreme, we could assume that all the \(\{\theta_{i}\}\) are independent and entirely unrelated. In this case, the SMR from each area would have to be estimated independent of the data for other areas. As we have saw in the SMR example above, this leads to statistical instability if the numbers are small.

A third possible assumption lies somewhere between these two extremes. One could assume that the \(\{\theta_{i}\}\) are “similar” in the sense that the area labels convey no additional information. This is known as exchangeability, and is equivalent to assuming that \(\{\theta_{i}\}\) are drawn from a common prior distribution with unknown parameters.

6.6.2 Poisson gamma model

A classic example of this approach is presented by Clayton and Kaldor (1987), who developed a Bayesian analysis of a Poisson likelihood model. Their model is a useful introduction to the idea of hierarchical modelling of disease rates, as the second stage distribution for the area variability is analytically tractable and helps to build intuition about how smoothing works to stabilize the SMR estimates.

In the first stage of the hierarchy, we assume that the observed death counts for each area are Poisson distributed: \[ O_{i} \sim \mathrm{Poisson}(\theta_{i} E_{i}). \] In the second stage, a hierarchical prior is placed on \(\theta_{i}\):
\[ \theta_{i} \sim \mathrm{Gamma}(\nu, \alpha) . \] Recalling that the gamma distribution with parameters \(\nu\) and \(\alpha\) has mean \(\nu/\alpha\), this simply states that we expect the distribution of \(\{ \theta_{i}\}\) to follow a gamma distribution with mean \(\nu/\alpha\) and variance \(\nu/\alpha^{2}\). Since the gamma distribution is the conjugate prior of the Poisson, the posterior distribution of \(p(\theta_{i} | O_{i},E_{i})\) also follows a gamma distribution: \[ \mathrm{Gamma}(\nu + O_{i}, \alpha + E_{i}) \] with mean given by \[ \mathbb{E}(\theta_{i}|O_{i},\nu, \alpha) = \frac{\nu+O_{i}}{\alpha + E_{i}} = w_{i}\mbox{SMR}_{i} + (1-w_{i}) \frac{\nu}{\alpha} \] where \[w_{i} = \frac{E_{i}}{\alpha+E_{i}}. \] The expression for \(\mathbb{E}(\theta_{i}|O_{i},\nu,\alpha)\) shows that the posterior mean of the relative risk for the \(i\)th area is a weighted average of the observed SMR for the \(i\)th area and the average relative risk (\(\nu/\alpha\)) over all areas. The weight is inversely proportional to the variance of the SMR. Accordingly, when \(E_{i}\) is small (for rare diseases or small population counts), the variance is large, so the weight \(w_{i}\) is small and the posterior mean is dominated by the prior mean, \(\nu/\alpha\). In areas with abundant data, the posterior mean is close to the observed \(SMR_i = O_{i}/E_{i}\). This feature, whereby the amount of smoothing is proportional to the amount of information available for a particular area, is known as precision weighting. It has an intuitive appeal in that, when one does not observe a lot of information about an area (because the sample size is small and the risk estimate is unstable), one’s “best guess” concerning that area’s mortality risk should be weighted towards what little is known from prior knowledge, i.e. that the risk is, on average, \(\nu/\alpha\). In contrast, if one observes a lot of information for an area (e.g. because the sample size is large), one is more likely to believe what the data say about mortality risk in that particular area, and thus the “best guess” would reasonably be weighted towards the observed SMR for that specific area.

In the Empirical Bayes approach developed by Clayton and Kaldor (1987), \(\nu\) and \(\alpha\) are replaced by their estimates, \(\hat{\nu}\) and \(\hat{\alpha}\), which can be calculated by means of an iterative procedure using the following two equations:

\[ \begin{align*} \frac{\hat{\nu}}{\hat{\alpha}} & = \frac{1}{n} \sum_{i}^{n} \frac{O_i + \hat{\nu}}{E_i + \hat{\alpha}} = \frac{1}{n} \sum_i \hat{\theta}_i \\ \frac{\hat{\nu}}{\hat{\alpha}^2} & = \frac{1}{n-1} \sum_{i}^{n} \left(1 + \frac{\hat{\alpha}}{E_i} \right) \left(\hat{\theta}_i - \frac{\hat{\nu}}{\hat{\alpha}}\right)^2 \end{align*} \] where \(\{\hat{\theta}_i\}\) are the empirical Bayes estimates. Together, these two equations can be used recursively to compute \(\hat{\nu}\) and \(\hat{\alpha}\). At each stage of the iteration, the \(\{\hat{\theta}_i\}\) are calculated from the current estimates of \(\nu\) and \(\alpha\), and then the right hand sides of the two equations are used to provide new estimates of \(\nu\) and \(\alpha\) (Clayton and Kaldor, 1987).

There is an interesting connection here to negative binomial regression in that the marginal posterior distribution of \(O_i\) (unconditional on \(\theta_i\)) is negative binomial with size \(\nu\) and probability \(\alpha/(E_i + \alpha)\).

We can apply this algorithm to the observed and expected death counts in the census tracts of our study area to obtain empirical Bayes estimates of the census tract level SMRs, which we compare to the raw \(SMR_i = O_i/E_i\):

# Use DCluster::empbaysmooth to fit the Poisson Gamma model as proposed by Clayton and Kaldor (1987)
poisson_gamma <- DCluster::empbaysmooth(df_indirect_ordered$O, df_indirect_ordered$E)

# Append these estimates to the dataset and also calculate naive empirical Bayes credible intervals
# using the gamma distribution.
df_eb <- df_indirect_ordered %>%
         mutate(nu1 = O + poisson_gamma$nu,
                alpha1 = E + poisson_gamma$alpha,
                empbayes_SMR = poisson_gamma$smthrr,
                empbayes_CI95low = qgamma(0.025, nu1, alpha1),
                empbayes_CI95up = qgamma(0.975, nu1, alpha1),
                eb_sig = factor(case_when(
                  empbayes_CI95low>1 & empbayes_CI95up>1 ~ 1,
                  empbayes_CI95low<1 & empbayes_CI95up<1 ~ -1,
                  TRUE ~ 0)),
                raw_SMR = O/E, # recenter rawSMRs by adding intercept from intercept only model
                raw_SMR_CI95low = pois.exact(x=O, pt=E, conf.level=0.95)[,4],
                raw_SMR_CI95up = pois.exact(x=O, pt=E, conf.level=0.94)[,5],
                raw_sig = factor(case_when(
                  raw_SMR_CI95low>1 & raw_SMR_CI95up>1 ~ 1,
                  raw_SMR_CI95low<1 & raw_SMR_CI95up<1 ~ -1,
                  TRUE ~ 0))) %>%
          mutate(empbayes_SMR = ifelse(empbayes_SMR==0, nu1/alpha1, empbayes_SMR))

# caterpillar plot of the raw SMRs
rawSMR_caterpillar_eb <- ggplot(df_eb %>% arrange(raw_SMR) %>% mutate(id_order = row_number()),
       aes(x=id_order, y=raw_SMR, ymin=raw_SMR_CI95low, ymax=raw_SMR_CI95up, col=raw_sig)) +
  geom_point() +
  geom_errorbar() +
  geom_hline(yintercept = 1, col="red", linetype="dotted") +
  ylim(c(0,6)) +
  labs(x="", y="raw SMR") +
  scale_color_manual(name="",
                       values=c("#018571", "grey50", "#A6611A"),
                       labels=c("CI below 1",
                                "CI crosses 1",
                                "CI above 1")) +
  theme(legend.position="bottom")

# caterpillar plot of the empirical Bayes SMRs
ebSMR_caterpillar_eb <- ggplot(df_eb %>% arrange(empbayes_SMR) %>% mutate(id_order = row_number()),
       aes(x=id_order, y=empbayes_SMR, ymin=empbayes_CI95low, ymax=empbayes_CI95up, col=eb_sig)) +
  geom_point() +
  geom_errorbar() +
  geom_hline(yintercept = 1, col="red", linetype="dotted") +
  ylim(c(0,6)) +
  labs(x="", y="empirical Bayes SMR") +
  scale_color_manual(name="",
                       values=c("#018571", "grey50", "#A6611A"),
                       labels=c("CI below 1",
                                "CI crosses 1",
                                "CI above 1")) +
  theme(legend.position="bottom")

rawSMR_caterpillar_eb + ebSMR_caterpillar_eb

From the caterpillar plots, we observe that the empirial Bayes estimates are substantially smoothed relative to the raw SMRs. The observations at the extremes of the distribution with huge confidence intervals are smoothed towards the mean, and the spread of smoothed SMRs is generally narrower than the raw SMRs.

Figure 6.4: Caterpillar plots of raw SMRs with 95% CIs and smoothed SMRs from a Poisson gamma model with 95% credible intervals

This is also evident if we compare maps of the raw SMRs to the empirical Bayes SMRs, where the more extreme low or high SMRs are smoothed towards the middle values.

Figure 6.5: Maps of raw SMRs (left) and smoothed SMRs from a Poisson gamma model

6.6.3 Poisson log normal model

While a gamma prior for \(\{\theta_{i}\}\) is mathematically convenient, it can be restrictive in that it is difficult to incorporate observed covariates into the model. A more flexible model makes use of a normal prior for \(\log(\theta_{i})\) (Wakefield et al., 2000; Lawson et al., 2003): \[ \begin{align*} O_{i} & \sim \mathrm{Poisson}(\theta_{i} E_{i}) \\ \log(\theta_{i}) & = \alpha + v_{i} \\ v_{i} & \sim \mathrm{Normal}(0, \sigma^{2}_{v}) \end{align*} \] where \(\alpha\) is an intercept term representing the overall log relative risk of disease in the whole study region compared to the reference rate and \(v_{i}\) is the residual log relative risk in area \(i\) compared with the average over the study region. Here, the \(\{v_{i}\}\) are assumed to arise from a normal distribution with mean 0 and variance \(\sigma^{2}_{v}\). You may recognize this model as a generalized linear mixed model, where “mixed” refers to the fact that the model contains both “random effects” (the \(\{v_{i}\}\)), as well as accommodating “fixed” covariate effects, e.g. \[ \log(\theta_{i}) = \alpha + \beta_{1} x_{1} + \ldots + \beta_{p} x_{p} + v_{i} . \] In this case the \(\{v_{i}\}\) are interpretable as *residual area-specific effects conditional on the fixed covariates. Similarly, the \(\beta\)s are interpretable as covariate effects conditional on the area random effects.

Fixed vs Random Effects: in standard regression models, fixed and random effects refer to the type of statistical model being used. When using fixed effects analysis of variance (ANOVA), the assumptions being made are about the independent variable and the error distributions for the variable. Fixed effects are estimated with maximum likelihood (the traditional beta estimates we see in regression models). They are most appropriate when trying to generalize results to values for the fixed variables used in the study (fixed variables are “assumed to be measured without error … and assumed that the values of a fixed variable in one study are the same as values of the fixed variable in another study”) (Newsom, 2019). If however, the researcher is seeking to make inferences beyond particular values of the independent variable, a random effects model is used. Random effects are estimated with shrinkage (partial pooling or linear unbiased predictions). Random effects models are accounting for “additional expected random variation on the independent variable” and instead of the value of the variable itself being of interest, the random variables “are assumed to be values that are drawn from a larger population of values and thus will represent them” (Newsom, 2019; Gelman, 2005)

It should be noted here that, as pointed out by Wolpert and Ickstadt (1998), the Poisson log normal model does not aggregate consistently. That is, if one specifies a log normal distribution for each of the relative risks and then combines two areas and specifies a log normal distribution for the relative risk of the combined area, then these distributions are inconsistent (because the sum of log normal distributions is not log normal). This can be understood as a form of aggregation bias whereby risk relationships do not remain constant across levels of aggregation (Wakefield et al., 2000). Nevertheless, a normal second-stage distribution has been observed empirically to provide a good model for log relative risks over a range of aggregations, and does present advantages with respect to model flexibility and ease of computation.

Unlike the Poisson gamma model, there is no analytically tractable closed-form solution for the posterior distributions. Instead, we will use INLA to fit this model.

# Use INLA to fit Poisson gamma model
model_form <- O ~ 1 + f(id_order, model="iid")
model_iid <- inla(model_form, family="poisson", 
                data=df_indirect_ordered, E=E, # E points to the expected count field
                control.predictor=list(compute=TRUE), # computes transformed posterior marginals
                control.compute=list(dic=TRUE, waic=TRUE)) # computes DIC for model fit
summary(model_iid)

We can see from the caterpillar plots of the raw SMRs vs. the smoothed SMRs from the Poisson lognormal model that, similar to the Poisson gamma model, the Poisson lognormal model smooths extreme SMRs towards the overall mean.

Figure 6.6: Caterpillar plots of raw SMRs with 95% CIs and smoothed SMRs from a Poisson lognormal model with 95% credible intervals

Figure 6.7: Maps of raw SMRs (left) and smoothed SMRs from a Poisson lognormal model

The estimate of the SMR for area \(i\) relative to the reference mortality rate (from the indirect age standardization) is \(\exp(\alpha + v_{i})\), while \(\exp(v_{i})\) is the residual relative risk for area \(i\) relative to the mean over the study area. Note that if an internal set of reference rates for indirect age standardization has been used, i.e. based on the observed rates over the whole study area itself, then \(\alpha\) will generally be close to zero.

6.6.4 Poisson multilevel model

In both the Poisson gamma and Poisson log normal models, the smoothing is global: the random effects are treated as exchangeable and come from the same common distribution. As a result, all of the area specific relative risks are shrunk to the same overall mean. This kind of global smoothing does not allow for spatial correlation between risks in nearby areas, as might be expected if there is local clustering in the spatial pattern of risks. Local clustering of risks may be due to risk factors that are shared within and across census tracts: individuals who share these spatially varying risk factors would be expected to have similar outcomes. If such local clustering exists, this is a violation of the assumption of exchangeability. If these risk factors are known and observed, the most straightforward solution would be to include them as covariates in the model. However, in most observational settings, one can rarely measure, or even know about, all the relevant risk factors. If one suspects that these risk factors vary by area, it will be necessary to consider ways to allow for local clustering in models for \(\{\theta_{i}\}\).

In order to accommodate local smoothing, we can specify a Poisson multilevel model with random effects at the census tract and neighborhood levels. In the Greater Boston example, we will treat the neighborhoods of Boston as roughly equivalent as a level to the city/town designation in the communities outside of Boston that make up the greater Boston study area.

\[ \begin{align*} O_{ij} & \sim \mathrm{Poisson}(\theta_{ij} E_{ij}) \\ \log(\theta_{ij}) & = \alpha + v_{ij} + \eta_j \\ v_{ij} & \sim \mathrm{Normal}(0, \sigma^{2}_{v})\\ \eta_{j} & \sim \mathrm{Normal}(0, \sigma^{2}_{\eta}) \end{align*} \]

This can be in INLA using the following code:

model_form <- O ~ 1 + f(id_order, model="iid") + f(cityid, model="iid")
model_multi <- inla(model_form, family="poisson", 
                data=df_indirect_ordered, E=E, # E points to the expected count field
                control.predictor=list(compute=TRUE), # computes transformed posterior marginals
                control.compute=list(dic=TRUE, waic=TRUE)) # computes DIC for model fit
summary(model_multi)

Figure 6.8: Maps of CT raw SMRs (left) and smoothed SMRs from a Poisson multilevel model

Figure 6.9: Map of the neighborhood level variation in relative risk (Level 2 of the Poisson multilevel model)

6.6.5 Poisson BYM model

In the multilevel Poisson log normal model with independent census tracts and neighborhoods, the neighborhood units were specified a priori. In a sense, the neighborhood boundaries are fixed, and crossing a boundary means that one’s relative risk may jump by quite a lot. What if one wants to account for spatial correlation in a smoother manner? One way to do this is to specify a multivariate normal prior distribution for all the area parameters with a spatially-structured covariance matrix. Many different ways have been proposed to specify spatially structured multivariate normal distributions for \(\log (\theta_{i})\) (Wakefield et al., 2000; Banerjee et al., 2004).

One particular form of the multivariate normal distribution that is commonly used is the intrinsic Gaussian conditional autoregressive (CAR) prior suggested by Clayton and Kaldor (1987) and developed by Besag et al. (1991). This is one of the most popular ways of dealing with spatial autocorrelation. The spatial structure is formulated through a set of conditional autoregressions, which uses the fact that if a vector of random variables has a multivariate normal distribution, then the distribution of each element of that vector conditional on all the other elements in the vector is also normal, with mean and variance that depend on the original multivariate mean and covariance matrix.

As with the other models, the first stage model assumes that the observed counts are Poisson distributed, and that an additive model for \(\log(\theta_{i})\) can be specified for accommodating covariate effects: \[\begin{align*} O_{i} & \sim \mathrm{Poisson}(\theta_{i} E_{i}) \\ \log (\theta_{i}) & = \alpha + u_{i} \end{align*}\] Instead of the independent normal prior for the distribution of the census tract effects, one models \[ u_{i}|u_{j, j \neq i} \sim \mathrm{Normal}(\mu_{i},\tau^{2}_{u}/m_{i}) \] where \[ \mu_{i} = \frac{\sum_{j} w_{ij} u_{j}}{\sum_{j} w_{ij}}, \qquad \sigma^{2}_{i} = \frac{\tau^{2}_{u}}{\sum_{j} w_{ij}} \] and \[ w_{ij} = \begin{cases} 1& \text{if $i, j$ are adjacent},\\ 0 & \text{if they are not} \end{cases} \]

and \(m_{i}\) is the number of adjacent areas. To understand this, consider the red highlighted census tract \(i\) in Figure 6.10. If \(\partial_{i}\) is the set of areas adjacent to \(i\) (shaded in yellow in the figure), and one sets \(w_{ij}\) to 1 for areas \(j \in \partial_{i}\) and zero otherwise, then the prior distribution for \(u_{i}\) has conditional mean equal to the average of the neighbouring \(u_{j}\)’s and variance inversely proportional to the number of adjacent neighbours. The effect is to smooth \(u_{i}\) toward the mean risk in the set of neighbouring areas. Note that \(\tau^{2}_{u}\) is the variance (scaled by \(\sum_{j \in \partial_{i}} {w_{ij}} = m_{i}\), i.e. the number of neighbours); to emphasize that it is only interpretable conditionally, we have labeled it \(\tau^{2}_{u}\) instead of \(\sigma^{2}_{u}\).

Figure 6.10: Illustration of the CAR idea

Besag, Yorke and Mollié (1991) recommend combining the CAR prior and the standard normal prior to allow for both spatially unstructured latent covariates spatially correlated latent covariates. Their model, called the BYM model or the convolution model in the literature, is thus

\[ \begin{align*} O_{i} & \sim \mbox{Poisson}(\theta_{i} E_{i}) \\ \mbox{log}(\theta_{i}) & = \beta_0 + u_{i} + v_{i} \\ v_{i} & \sim \mbox{Normal}(0,\sigma^{2}_{v}) \\ u_{i}|u_{j, j \neq i} & \sim \mbox{Normal}(\mu_{i},\sigma^{2}_{u}/m_{i}) \end{align*} \]

Here, the unstructured \(v_{i}\) can be thought of as capturing correlation within areas and the spatially structured \(u_{i}\) capture spatial correlation across areas.

Note, however, that \(\sigma^{2}_{v}\) (unstructured heterogeneity variance) and \(\tau^{2}_{u}\) (spatial variance) are not directly comparable: \(\sigma^{2}_{v}\) reflects the variability of the unstructured random effects between areas, while \(\tau^{2}_{u}\) is the variance of the spatial effect in area \(i\), conditional on values of neighboring spatial effects. No closed-form expression is available for the between-area variance of the spatial effects. However, in the Bayesian approach, the marginal spatial variance \(s^{2}_{u}\) can be estimated empirically from the posterior samples of \(\{u_{i}\}\): \[ s^{2}_{u} = \sum_{i} (u_{i}-\bar{u})^2 / (n-1) \] where \(\bar{u}\) is the average of the \(\{u_{i}\}\).

With an estimate of the marginal between-area variance, one can also characterize the relative contribution of spatial vs. unstructured heterogeneity: \[ \mathrm{frac_{spatial}} = s^{2}_{u} / (s^{2}_{u} + \sigma^{2}_{v}) \] When \(\mathrm{frac_{spatial}}\) is close to 1, spatial heterogeneity dominates. When \(\mathrm{frac_{spatial}}\) is close to 0, unstructured heterogeneity dominates (Best et al. 2005).

To aid with interpretation of the variability in the posterior estimates of \(SMR_i = \exp(u_i + v_i)\), we define a quantity which we will call QR90 (Best et al. 2005) \[ QR_{90} = \mbox{exp}(q_{95\%} - q_{5\%}) = \mbox{ relative risk between top and bottom 5\% of areas} \] where \[ \begin{align*} q_{5\%} & = \mbox{ log relative risk for the area ranked at the 5th percentile} \\ q_{95\%} & = \mbox{ log relative risk for the area ranked at the 95th percentile} \end{align*} \] The advantage of QR90 is that it expresses the variability in the SMRs on the relative risk scale, which facilitates comparison of the disparities between areas to the estimated inequities (also on the relative risk scale) attributable to categories of ABSMs.

Figure 6.11: Depiction of QR90: the relative risk comparing the 95th and 5th quantiles of the random effects distribution

Figure 6.12: Caterpillar plots of raw SMRs with 95% CIs and smoothed SMRs from a Poisson BYM model with 95% credible intervals

Figure 6.13: Comparison of premature mortality raw SMRs and smoothed (BYM) SMRs

6.7 Estimating ABSM effects

6.7.1 Premature mortality

Figure 6.14: ABSMs vs. smoothed SMRs (BYM)

Figure 6.15: Comparison of premature mortality inequities by ABSM categories across all analytic methods and models

Table 6.1: Comparison of estimated inequities in premature mortality by category of CT ABSMs from aggregated analysis and Quasipoisson, Negative Binomial, and BYM models

		Aggregated			Quasipoisson			Negative binomial			BYM
ABSM	Category	IRR	conf.low	conf.high	IRR	conf.low	conf.high	IRR	conf.low	conf.high	IRR	conf.low	conf.high
CT % below poverty
apINDPOV	5-9.9%	1.25	0.94	1.66	1.25	1.09	1.43	1.19	1.04	1.36	1.15	0.96	1.39
apINDPOV	10-19.9%	1.78	1.37	2.31	1.77	1.57	2.01	1.69	1.49	1.92	1.53	1.27	1.83
apINDPOV	20-100%	2.18	1.68	2.84	2.18	1.93	2.48	2.17	1.92	2.47	1.75	1.45	2.11
CT Index of Concentration at the Extremes
qICEwnhinc	[-0.711,0.0557]	2.37	1.85	3.04	2.39	2.13	2.70	2.56	2.29	2.86	2.13	1.84	2.47
qICEwnhinc	(0.0557,0.238]	1.84	1.41	2.40	1.85	1.63	2.10	1.78	1.59	2.01	1.74	1.49	2.03
qICEwnhinc	(0.238,0.358]	1.46	1.08	1.97	1.46	1.27	1.69	1.43	1.25	1.63	1.44	1.22	1.72
qICEwnhinc	(0.358,0.464]	1.33	0.99	1.80	1.35	1.17	1.55	1.34	1.17	1.53	1.35	1.14	1.60
CT Social Vulnerability Index
qsvi	(0.125,0.246]	1.09	0.78	1.51	1.08	0.92	1.27	1.05	0.91	1.21	1.10	0.92	1.30
qsvi	(0.246,0.407]	1.17	0.87	1.58	1.18	1.03	1.37	1.15	1.01	1.32	1.30	1.10	1.52
qsvi	(0.407,0.683]	1.55	1.17	2.05	1.56	1.36	1.79	1.48	1.31	1.68	1.53	1.31	1.78
qsvi	(0.683,1]	2.11	1.61	2.76	2.15	1.89	2.45	2.28	2.03	2.57	2.13	1.84	2.48

Figure 6.16: Premature mortality: comparison of smoothed SMRs (BYM) before and after adjustment for CT % below poverty

Figure 6.17: Premature mortality: comparison of smoothed SMRs (BYM) before and after adjustment for CT Index of Concentration at the Extremes (racialized economic segregation)

Figure 6.18: Premature mortality: comparison of smoothed SMRs (BYM) before and after adjustment for CT Social Vulnerability Index

6.7.2 Lung cancer mortality

Figure 6.19: Lung cancer mortality: CT ABSMs vs. smoothed SMRs (BYM)

Figure 6.20: Comparison of lung cancer mortality inequities by CT ABSM categories across all analytic methods and models

Table 6.2: Comparison of estimated inequities in lung cancer mortality by category of CT ABSMs from aggregated analysis and Quasipoisson, Negative Binomial, and BYM models

		Aggregated			Quasipoisson			Negative binomial			BYM
ABSM	Category	IRR	conf.low	conf.high	IRR	conf.low	conf.high	IRR	conf.low	conf.high	IRR	conf.low	conf.high
CT % below poverty
apINDPOV	5-9.9%	1.04	0.62	1.76	1.07	0.86	1.33	1.05	0.87	1.26	1.01	0.83	1.23
apINDPOV	10-19.9%	1.32	0.81	2.17	1.35	1.10	1.66	1.33	1.12	1.58	1.18	0.97	1.43
apINDPOV	20-100%	1.14	0.68	1.91	1.19	0.96	1.47	1.18	0.99	1.41	1.00	0.81	1.23
CT Index of Concentration at the Extremes
qICEwnhinc	[-0.711,0.0557]	1.29	0.79	2.10	1.33	1.08	1.63	1.32	1.12	1.56	1.12	0.93	1.35
qICEwnhinc	(0.0557,0.238]	1.42	0.86	2.36	1.44	1.17	1.79	1.42	1.19	1.69	1.29	1.06	1.55
qICEwnhinc	(0.238,0.358]	1.37	0.79	2.37	1.37	1.09	1.73	1.34	1.10	1.62	1.28	1.04	1.57
qICEwnhinc	(0.358,0.464]	1.15	0.65	2.05	1.16	0.92	1.48	1.14	0.94	1.39	1.10	0.90	1.36
CT Social Vulnerability Index
qsvi	(0.125,0.246]	0.88	0.47	1.64	0.92	0.71	1.19	0.90	0.73	1.11	0.91	0.73	1.14
qsvi	(0.246,0.407]	0.84	0.48	1.48	0.88	0.70	1.11	0.86	0.71	1.05	0.93	0.76	1.15
qsvi	(0.407,0.683]	1.00	0.58	1.72	1.04	0.83	1.31	1.03	0.86	1.25	1.02	0.84	1.24
qsvi	(0.683,1]	0.99	0.58	1.70	1.05	0.85	1.32	1.04	0.87	1.25	0.95	0.78	1.16

Figure 6.21: Lung cancer: comparison of smoothed SMRs (BYM) before and after adjustment for CT % below poverty

Figure 6.22: Lung cancer: comparison of smoothed SMRs (BYM) before and after adjustment for CT Index of Concentration at the Extremes (racialized economic segregation)

Figure 6.23: Lung cancer: comparison of smoothed SMRs (BYM) before and after adjustment for CT Social Vulnerability Index

6.8 Intersectional analysis of inequities by racialized group and CT ABSMs

In this section, we explore intersectional patterns of inequities in premature mortality by racialized group and CT ABSMs in the census tracts comprising the Greater Boston study area. That is, we focus on the joint patterning of inequities by membership in racialized groups (as captured by individual-level membership in racialized groups) and residence in census tracts characterized by area-based social metrics as indicative of racialized societal inequities. To accomplish this, we operationalize intersectional effects using methods for analyzing statistical interactions (Knol and VanderWeele 2012; Jackson et al. 2016; VanderWeele 2015). In contrast to analyses that focus on racial/ethnic disparities “controlling” for area-based social metrics or, conversely, disparities by ABSM “controlling” for racialized group membership, this framework emphasizes the importance of describing the joint patterning of racialized social inequities.

Drawing on recommendations on reporting interactions (Knol and VanderWeele 2012), we propose that health disparities researchers should be interested in reporting these interactions in three different ways:

• highlighting inequities by racialized group within categories of CT ABSMs

• highlighting inequities by CT ABSM within racialized groups

• highlighting the joint inequities by racialized group and CT ABSM relative to a common reference group

Moreover, a thorough reporting of interaction analyses should present measures of effect measure modification on both the additive and multiplicative scales.

6.8.1 Aggregated analysis

We begin by computing age-standardized premature mortality rates by racialized group and CT ABSMs. Similar to how we computed these age-standardized rates by CT ABSMs alone, we aggregate age-specific death counts and population person-time by racialized group across census tracts within categories of CT ABSM and apply direct age-standardization. We present theses estimated age-standardized premature mortality rates and 95% confidence intervals in Figure 6.24.

Several patterns are apparent from the visualization of the age-standardized rates. Firstly, gradients in inequities by CT ABSM are apparent for all three ABSMs, with premature mortality rates higher in more disadvantaged census tracts. This gradient is more pronounced among the non-Hispanic White population, with steeper gradients and tighter confidence limits. Among the Black population, inequities across CT ABSM categories are more shallow, with wider confidence limits in the more advantaged categories reflecting the smaller population sizes of Black individuals living in these census tracts.

Focusing on the Black/White inequities across CT ABSM categories, there is some suggestion of qualitative interaction. That is, Black/White inequities are largest in the most advantaged census tracts, but in more disadvantaged census tracts, non-Hispanic Whites have higher rates of premature mortality than Black populations living in the same areas, particular by quintile of CT ICE and SVI.

Figure 6.24: Age-standardized premature mortality rates by racialized group and CT ABSM, Greater Boston, 2013-2017, as computed by the aggregated method

To evaluate the magnitude of the CT ABSM inequities within racialized group, we present (a) age-standardized rate differences per 100,000 person-years and (b) age-standardized incidence rate ratios in Figure 6.25. On both the additive (incidence rate difference) and multiplicative scales (incidence rate ratio), ABSM gradients are much steeper among the non-Hispanic White population compared with the Black population. ABSM gradients among the Black population are attenuated towards the null and have much wider confidence limits.

Figure 6.25: ABSM gradients in premature mortality by racialized group (Blacks vs. Non-Hispanic Whites): age-standardized rate differences per 100,000 person-years (top row) and rate ratios (bottom row) computed using the aggregation method, Greater Boston, 2013-2017. In each plot, the reference group is the most advantaged racialized group-specific ABSM category.

Meanwhile, Black/White inequities as depicted in Figure 6.26 confirm the pattern of Black/White crossover suggested in Figure 6.24, with non-Hispanic White premature mortality rates significantly elevated in the most disadvantaged census tracts by CT Index of Concentration at the Extremes and CT Social Vulnerability Index.

Figure 6.26: Black/White inequities in premature mortality by category of CT ABSM: age-standardized rate differences per 100,000 person-years (top row) and rate ratios (bottom row) computed using the aggregation method, Greater Boston, 2013-2017.

We visualize the joint pattern of inequities on the additive (incidence rate difference) and multiplicative (incidence rate ratio) scales relative to a common reference group in Figure 6.27. In general, it gives insights similar to visualization of age-standardized incidence rates by racialized group and CT ABSM in Figure 6.24 above.

Figure 6.27: Intersectional age-standardized incidence rate differences per 100,000 person-years (top row) and age-standardized rate ratios (bottom row) by racialized group and category of CT ABSM computed using the aggregation method, Greater Boston, 2013-2017. Presented with a common reference group (non-Hispanic Whites in the most advantaged categories of CT ABSMs

ABSM gradients among the Black population are attenuated towards the null and have much wider confidence limits. Why do you think this might be the case?

6.8.2 Intersectional inequities as estimated by non-spatial, multilevel, and spatial regression models

Intersectional inequities by racialized group and CT ABSM can also be estimated using regression models.

Figure 6.28: Comparison of inequities in premature mortality by racialized group and CT ABSM as estimated by the aggregated method, Poisson regression with robust variance estimator, multilevel Poisson regression (census tracts in city/town/neighborhoods), and spatial Poisson regression with BYM prior, Greater Boston, 2013-2017. Here, we focus on estimated incidence rate ratios with a common reference group (non-Hispanic Whites in the most advantaged census tracts)

Table 6.3: Comparison of intersectional inequities in premature mortality by racialized group and category of CT ABSMs from aggregated analysis, Poisson with robust variance estimator, multilevel, and BYM models. Note the common reference group (non-Hispanic Whites in the most advantaged CTs)

		Aggregated			Poisson robust			Multilevel			BYM
Racialized group	ABSM	IRR	conf.low	conf.high	IRR	conf.low	conf.high	IRR	conf.low	conf.high	IRR	conf.low	conf.high
CT % below poverty
NHW	0-4.9%	1.00			1.00			1.00			1.00
NHW	5-9.9%	1.28	0.93	1.75	1.28	1.10	1.48	1.20	1.07	1.34	1.25	1.02	1.54
NHW	10-19.9%	1.99	1.47	2.69	1.98	1.73	2.27	1.63	1.45	1.84	1.87	1.52	2.30
NHW	20-100%	2.49	1.81	3.43	2.40	2.05	2.82	1.87	1.64	2.12	2.07	1.66	2.58
Black	0-4.9%	1.29	0.46	3.61	1.29	0.89	1.89	1.19	0.84	1.67	1.14	0.79	1.65
Black	5-9.9%	1.69	0.99	2.87	1.68	1.34	2.10	1.42	1.18	1.71	1.42	1.09	1.83
Black	10-19.9%	1.76	1.23	2.52	1.79	1.49	2.15	1.38	1.20	1.58	1.54	1.23	1.92
Black	20-100%	2.30	1.69	3.15	2.34	2.03	2.70	1.53	1.34	1.76	1.90	1.52	2.37
CT Index of Concentration at the Extremes
NHW	(0.464,0.747]	1.00			1.00			1.00			1.00
NHW	(0.358,0.464]	1.37	0.99	1.91	1.39	1.18	1.63	1.35	1.21	1.51	1.32	1.11	1.58
NHW	(0.238,0.358]	1.55	1.11	2.16	1.56	1.34	1.82	1.58	1.41	1.78	1.57	1.31	1.89
NHW	(0.0557,0.238]	2.09	1.55	2.83	2.10	1.83	2.42	2.01	1.79	2.25	1.98	1.67	2.33
NHW	[-0.711,0.0557]	3.27	2.43	4.39	3.21	2.77	3.71	2.74	2.44	3.07	2.98	2.52	3.52
Black	(0.464,0.747]	2.02	0.65	6.31	1.90	1.26	2.87	2.02	1.41	2.88	1.81	1.28	2.55
Black	(0.358,0.464]	2.35	1.08	5.10	2.43	1.75	3.38	2.36	1.82	3.06	2.32	1.72	3.12
Black	(0.238,0.358]	1.90	0.97	3.73	1.89	1.44	2.47	1.91	1.52	2.39	1.98	1.51	2.59
Black	(0.0557,0.238]	1.95	1.29	2.95	1.99	1.62	2.44	1.92	1.65	2.22	1.92	1.58	2.33
Black	[-0.711,0.0557]	2.37	1.77	3.18	2.41	2.09	2.79	1.91	1.69	2.16	2.19	1.85	2.60
CT Social Vulnerability Index
NHW	[0.0003,0.125]	1.00			1.00			1.00			1.00
NHW	(0.125,0.246]	1.11	0.78	1.60	1.11	0.94	1.32	1.18	1.04	1.33	1.10	0.92	1.33
NHW	(0.246,0.407]	1.23	0.89	1.71	1.26	1.07	1.48	1.48	1.32	1.65	1.32	1.11	1.56
NHW	(0.407,0.683]	1.68	1.23	2.31	1.70	1.46	1.98	1.69	1.51	1.89	1.63	1.38	1.92
NHW	(0.683,1]	3.00	2.20	4.11	3.03	2.58	3.57	2.66	2.36	3.00	3.16	2.67	3.74
Black	[0.0003,0.125]	1.52	0.43	5.34	1.54	0.94	2.51	1.61	1.07	2.43	1.49	1.01	2.21
Black	(0.125,0.246]	1.52	0.69	3.35	1.51	1.10	2.06	1.49	1.14	1.93	1.39	1.03	1.88
Black	(0.246,0.407]	1.72	0.91	3.25	1.73	1.31	2.26	1.95	1.58	2.41	1.74	1.36	2.23
Black	(0.407,0.683]	1.77	1.19	2.63	1.82	1.49	2.21	1.82	1.59	2.10	1.82	1.51	2.19
Black	(0.683,1]	2.05	1.51	2.79	2.09	1.79	2.45	1.72	1.52	1.95	2.08	1.76	2.47

Figure 6.29: Comparison of CT ABSM inequities in premature mortality within racialized group (Black and Non-Hispanic White) as estimated by the aggregated method, Poisson regression with robust variance estimator, multilevel Poisson regression (census tracts in city/town/neighborhoods), and spatial Poisson regression with BYM prior, Greater Boston, 2013-2017. In each plot, the reference group is the most advantaged ABSM category within racialized group.

As seen in Figure 6.29, inequities in premature mortality across CT ABSM categories are much more pronounced for the non-Hispanic White population, whereas the gradient in age-adjusted incidence rate ratios for CT ABSM categories within the Black population is much shallower.

Figure 6.30: Comparison of Black vs. non-Hispanic White inequities in premature mortality within category of CT ABSM as estimated by the aggregated method, Poisson regression with robust variance estimator, multilevel Poisson regression (census tracts in city/town/neighborhoods), and spatial Poisson regression with BYM prior, Greater Boston, 2013-2017

As seen in Figure 6.30, Black-White inequities as measured by incidence rate ratios are largest in the most advantaged census tracts, whereas inequities by racialized group are smaller or even reversed in the most disadvantaged census tracts. That is, Non-Hispanic Whites living in the most disadvantaged census tracts have higher observed age-adjusted premature mortality rates than Black populations living in the same areas. Why do you think this might be the case?

Across methods, we note that the confidence limits on age-standardized rates, incidence rate differences, and incidence rate ratios are much wider for the aggregated method compared with model-based estimates. This reflects the impact of age-standardization, whereby age strata with relatively less information may nevertheless contribute strongly to the age-standardized estimates if those age strata are highly weighted in the standard population. The result is that standard errors and confidence limits from the aggregated method will usually be larger than their corresponding model-based estimates. Thus, the aggregated method, the most commonly used, has less power to detect statistically significant inequities.

6.9 REFERENCES

Banerjee A, Dhillon I, Ghosh J, Merugu S. An information theoretic analysis of maximum likelihood mixture estimation for exponential families. In Proceedings of the twenty-first international conference on Machine learning 2004 Jul 4 (p. 8). https://dl.acm.org/doi/10.1145/1015330.1015431 Accessed June 4, 2022

Besag J, York J, Mollié A. Bayesian image restoration, with two applications in spatial statistics. Annals of the institute of statistical mathematics. 1991 Mar;43(1):1-20.

Best N, Richardson S, Thomson A. A comparison of Bayesian spatial models for disease mapping. Statistical methods in medical research. 2005 Feb;14(1):35-59.

Clayton D, Kaldor J. Empirical Bayes estimates of age-standardized relative risks for use in disease mapping. Biometrics. 1987 Sep 1:671-81.

Gelman, A. Why I don’t use the term “fixed and random effects.” Post in Statistical Modeling, Causal Inference, and Social Science. 01,25,2005. Available at: https://statmodeling.stat.columbia.edu/2005/01/25/why_i_dont_use/

Jackson JW, Williams DR, VanderWeele TJ. Disparities at the intersection of marginalized groups. Soc Psychiatry Psychiatr Epidemiol. 2016 Oct;51(10):1349-1359. doi: 10.1007/s00127-016-1276-6. Epub 2016 Aug 16. PMID: 27531592; PMCID: PMC5350011.

Knol MJ, VanderWeele TJ. Recommendations for presenting analyses of effect modification and interaction. Int J Epidemiol. 2012 Apr;41(2):514-20. doi: 10.1093/ije/dyr218. Epub 2012 Jan 9. PMID: 22253321; PMCID: PMC3324457.

Krieger N, Waterman PD, Chen JT, Rehkopf DH, Subramanian SV. The Public Health Disparities Geocoding Project Monograph. Available as of June 30, 2004 at: http://www.hsph.harvard.edu/thegeocodingproject

Lawson AB, Browne WJ, Rodeiro CL. Disease mapping with WinBUGS and MLwiN. John Wiley & Sons; 2003 Sep 12.

Nethery RC, Rushovich T, Peterson E, Chen JT, Waterman PD, Krieger N, Waller L, Coull BA. Comparing denominator sources for real-time disease incidence modeling: American Community Survey and WorldPop. SSM-Population Health. 2021 Jun 1;14:100786.

Newsom, J. Distinguishing Between Random and Fixed: Variables, Effects, and Coefficients. Psy 526 Lecture Material: Multilevel Regression. Spring 2019. Available at: http://www.web.pdx.edu/~newsomj/mlrclass/ho_randfixd.pdf

Pickle LW, White AA. Effects of the choice of age‐adjustment method on maps of death rates. Statistics in medicine. 1995 Mar 15;14(5‐7):615-27.

United States Census Bureau. Population and Housing Estimate Program (PEP). Updated 2022. Available at: https://www.census.gov/programs-surveys/popest.html

VanderWeele TJ. Explanation in Causal Inference: Methods for Mediation and Interaction. New York, NY: Oxford University Press, 2015.

Wakefield JC, Best NG, and Waller L. Bayesian approaches to diseasemapping, In: Elliott, P., Wakefield, J.C., Best, N.G. and Briggs, D.J. SpatialEpidemiology: Methods and Applications, Oxford University Press, London, UK.pp 104–127, 2000.

Wolpert RL, Ickstadt K. Poisson/gamma random field models for spatial statistics. Biometrika. 1998 Jun 1;85(2):251-67.