Can Changes in the Distributions of and Associations Between Education and Income Bias Temporal Comparisons of Health Disparities? An Exploration With Causal Graphs and Simulations

Biased estimate of association at a single time point

We begin by developing simple causal directed acyclic graphs (DAGs) to formalize the relationships between observed and unobserved variables of interest at a single time point. Causal DAGs are a simple, flexible device for demonstrating the associations implied by a given set of assumptions about the causal structure relating variables (31–36), and they can help us answer whether a particular statistical analysis will yield an unbiased estimate of a parameter of interest. The DAGs we present are not intended to encode all possible causal structures relating multiple dimensions of socioeconomic position to health, but to focus attention on the structural relationships between observed and unobserved socioeconomic indicators and the outcome.

Consider the situation where we have observed data on a socioeconomic indicator X (e.g., education) that may have a causal effect on an outcome Y (e.g., body mass index). We are interested in the potential effect on estimates of the X − Y association if there is an unobserved socioeconomic indicator Z (e.g., income) that is correlated with X and also has a direct causal effect on Y. Figure 1 shows 4 causal structures that might account for the correlation between X and Z, including the possible effect of an unmeasured common cause of X and Z, W (e.g., parental income when the individual was a child), or conditioning on a common effect, S (e.g., mortality). Under any of these scenarios, absence of observed data on Z (and W and S) precludes statistical adjustment for backdoor paths between X and Y, making estimates of the X − Y association only interpretable as descriptive (rather than causal) (31–36).

Comparing descriptive associations across time

Now consider that we have observed data on X and Y from 2 cross-sectional surveys at different times, and the objective is to determine whether there has been a change in the socioeconomic gradient that is not attributable to shifting associations between the observed and unobserved socioeconomic indicators. Because a change in the strength of a relationship between variables constitutes effect modification by time, we represent this by drawing separate DAGs for each time and indicate qualitative change in effect magnitude by drawing a thicker arrow in the second DAG. Given the potential causal structures leading to correlation between X and Z in Figure 1, Figure 2 shows possible changes in the strength of relationships between variables that can lead to a change in the descriptive association of X and Y over time. We consider scenarios A (X → Y changes) and B (Z → Y changes) to reflect a change in socioeconomic gradient that meets our criterion, albeit noting that scenario C arguably still reflects a change in the socioeconomic gradient captured by the total effect of X on Y (direct effect X → Y and mediated effect X → Z → Y). Even with no other unmeasured confounders, it remains impossible to distinguish changes in the direct effects of X on Y or Z on Y in the absence of observed data on Z. Our simulations explore how the structures depicted in Figure 2 and modeling techniques frequently used by epidemiologists can influence temporal comparisons of the association of a single observed indicator X with Y.

Simulation framework specified as an errors-in-variables model

To develop our simulations, we put a parametric form on the qualitative relationships depicted in Figures 1 and 2 and associate values with the parameters governing these relationships (refer to the simplified data-generating graph in Figure 3). For continuous X, Y, and Z, we specify a series of linear regression models. Here, it is convenient to induce correlation between X and Z by specifying a data-generating model for X given Z: This has a connection to the classical errors-in-variables model (29) when there is no direct effect of X on Y, noting that equivalent data structures can be generated by assuming a model for Z given X, or generating X and Z as bivariate normal given a vector of means and a variance-covariance matrix (Web Appendix 1 available at http://aje.oxfordjournals.org/).

We generate Z at time t as a normally distributed continuous variable with mean μ_z and variance  (noting that our results do not require that Z be normally distributed). At a given time, the ith subject's outcome y_i depends on the variables z_i via the linear regression

(1)

where ɛ_i is a normally distributed error term with mean zero and variance . The observed x_i at a given time depends on z_i via

(2)

where η_i is a normally distributed error term with mean zero and variance . We assume that z_i, ɛ_i, and η_i are independent.

At a given time, if the y_i's are regressed on the observed x_i's, then the estimator for the slope coefficient from this regression as n → ∞ is the following:

(3)

where β is the true effect of Z on Y, and β* is the observed association between the observed X and Y (29, 30). Details of this derivation are presented in Web Appendix 1. As this expression clarifies, changes in any 1 of the 4 parameters will influence the association between X and Y over time.

As shown in Web Appendix 2, we considered 6 different scenarios in our simulations, variously involving changes in the distribution and relationships of the variables Y, Z, and X. Values for the various parameters are summarized in Table 1. We calibrated the parameters to yield distributions of Y, Z, and X roughly similar to those observed in NHANES I (1971–1975) and NHANES 2005–2008 (19) for body mass index, log(adjusted family income), and education (modeled as years of education, to keep all values continuous). For each of scenarios 1–6, we simulated 1,000 data sets of 1,000 subjects at each of the 2 time points.

Table 1.

Simulation Parameters for the 6 Different Simulation Scenarios, as Specified in Web Appendix 2

Scenario 1			Scenario 2a			Scenario 3			Scenario 4			Scenario 5			Scenario 6
Parameters	Time Point		Parameters	Time Point		Parameters	Time Point		Parameters	Time Point		Parameters	Time Point		Parameters	Time Point
	t = 1	t = 2		t = 1	t = 2		t = 1	t = 2		t = 1	t = 2		t = 1	t = 2		t = 1	t = 2
μz	10.0	10.0	μz	10.0	10.0	μz	10.0	10.0	μz	10.0	10.0	μz	10.0	10.0	μz	10.0	10.0
	1.0	1.0		1.0	1.0		1.0	1.0	b	1.0	4.0		1.0	1.0	b	1.0	4.0
α0b	−2.0	0.0	α0b	−18.0	0	α0	−2.0	−2.0	α0	−2.0	−2.0	α0	−2.0	−2.0	α0b	−2.0	−12.0
α1	1.5	1.5	α1b	3.0	1.5	α1	1.5	1.5	α1	1.5	1.5	α1	1.5	1.5	α1b	1.5	3.0
	0.3	0.3		0.3	0.3	b	0.3	1.0		0.3	0.3		0.3	0.3	b	0.3	0.6
β0	29.0	29.0	β0	29.0	29.0	β0	29.0	29.0	β0	29.0	29.0	β0	29.0	29.0	β0	29.0	29.0
β1	−0.5	−0.5	β1	−0.5	−0.5	β1	−0.5	−0.5	β1	−0.5	−0.5	β1b	−0.5	−0.2	β1b	−0.5	−0.3
	4.0	4.0		4.0	4.0		4.0	4.0		4.0	4.0		4.0	4.0		4.0	4.0

Open in a new tab

^a In order to keep X in an appropriate range for years of education, in our simulation model we varied both α₀ and α₁.

^b Indicates that parameter changes between time 1 and 2.

Modeling methods

To analyze the simulated data, we employed 6 different modeling approaches (designated I–VI), also explicated in Web Appendix 2: These approaches include different regression models using (unobserved) Z and (observed) X, an imputation model, and the use of deciles, SII, and the relative distribution method, as we explain below. For each scenario and modeling approach, we estimated socioeconomic gradients in the outcome, calculated the difference in the observed estimates between time periods, , and tested this against the null hypothesis of .

Relative index of inequality, slope index of inequality, and quantile rescaling

The relative index of inequality (RII) is a widely used measure for assessing trends in health disparities that is intended to address the changing composition of socioeconomic categories over time (8, 17, 20, 37, 38). First employed in health research by Pamuk in 1985 to analyze trends in occupational class gradients in mortality in the United Kingdom (spurred by the shrinking size of class V (unskilled manual labor)), the RII was devised as an index “based on the mortality levels of all five social classes, and their proportionate share of the population” (20, p. 18). It provides an estimate of the relative disparity in fitted health effects comparing the best with worst off in the socioeconomic distribution. Recent work also provides estimators for the RII that relax the assumption of linearity (39).

Figure 4 depicts how the RII (and the related slope index of inequality (SII)) are methods of rescaling X as measured to the population distribution of X (39), analogous to taking quantiles (e.g., deciles) of a continuous variable. In the context of this simulation and the additive model for Y given Z, it is the SII that we expect to be robust to changes in the scale of X, not the RII (a relative measure which additionally depends on the observed intercept in the regression of Y on X). We accordingly focus on the SII as method VI.

Although the RII has been discussed primarily as a measure of “total population impact” (37, 38), we note that it involves a deliberate data transformation that rescales X at each time onto the population distribution scale. The RII accordingly implies that the health effects at levels of X (as measured on its usual scale) at the 2 different times are not comparable, and that rescaling onto the population distribution scale is required to allow comparisons that are unaffected by temporal differences in X unrelated to the true socioeconomic gradient in Y. This is equivalent to positing a latent version of X* on the population distribution scale that gives rise to observed X and assuming that only the measurement relationship between X and X* changes over time, while the “true” relationship between X* and Y is unchanged.

In spite of the widespread use of the RII/SII for comparing socioeconomic gradients over time, we are aware of no systematic investigation of how effectively the RII/SII protects against the multiple sources of bias identified in equation 3.

Relative distribution method

New work in economics and sociology on trends in income inequality has also proposed a less extreme method of rescaling, termed the relative distribution approach (40–42). Here, the distribution at one time is characterized in relation to quantile cutpoints defined by a reference distribution (e.g., the distribution at the first time), thereby preserving marked differences in the shape of the second distribution relative to the first (Figure 5). To our knowledge, the relative distribution approach has not yet been used to study trends in socioeconomic inequalities in health and may be worth considering.

Imputation model

If we knew in advance that Z was an important socioeconomic predictor of Y, we would ideally measure it and include it in the analysis at each time point. In many settings, however, Z is not available, or if it is available, it may be subject to missingness. If data are available from an external validation source or from a subsample of the survey, an imputation model can be specified for Z given X, thereby permitting comparison of the association between Z and Y across time. We implemented a single imputation version of this as method III in our simulation by assuming that data on Z and X were available for a random sample of 30% of the population, and applied the imputation model to the rest of the population.

Simulation results

We summarize the contrast in effect estimates over time under the different scenarios and models in Figures 6 and 7. The boxplots show the distribution of the point estimates for the temporal contrast over 1,000 simulations, with method I showing the “true” temporal contrast, that is, the contrast that would be obtained if we observed and analyzed Z at both time points. In Table 2, we also report the mean temporal contrast over 1,000 simulations, the 95% simulation intervals, and the proportion of simulations in which the various methods detected a statistically significant temporal contrast assuming a type I error rate of 5%. In scenarios 1, 2, 3, and 4, this reflects the probability of detecting a significant temporal effect given no true change in the socioeconomic gradient, while for scenarios 5 and 6, this is the power to detect a significant change in the socioeconomic gradient given that there was a true change. In scenarios 1–4, the least biased methods are those whose median value is at or close to zero.

In scenario 1 (change in the location of the X distribution—i.e., α₀—with no other changes in the “meaning” of education and no change in the underlying exposure-disease relationship), all methods yielded, on average, the correct characterization of no change in the socioeconomic relationship, with the exception of the relative distribution method (method IV), which detected on average a decrease in the socioeconomic gradient of (d = −0.45) and would have erroneously concluded this was significant 22% of the time.

In scenario 2 (change in α₁, i.e., change in the strength of the relationship between X and Z), method II (regression of Y on Z) yielded, as expected, a biased estimate of the temporal comparison. In this case, as the strength of the relationship between Z and X decreased (from 3.0 to 1.5), this resulted in a decrease in the relationship between X and Y of 0.14, when there was no real change to the underlying relationship between Z and Y. The type I error rate was 87%. Both the quantile (method V) and SII (method VI) approaches were robust to this source of bias over 1,000 simulations: Methods V and VI had type I error rates of only 5% and 4%, respectively. Method IV performed quite poorly, with a mean temporal difference of −1.83 and a type I error rate of 92%.

In scenario 3 (change in , i.e., change in error variance of X), method II also led to biased estimates of the temporal effect, with a type I error rate of 29%. Methods IV, V, and VI also produced simulation means substantially greater than zero, with type I error rates of 19%, 12%, and 10%, respectively.

In scenario 4 (change in ), the bias using method II was relatively small: The simulation mean was −0.03 and the type I error rate was 9%. However, methods IV, V, and VI all perform poorly: Simulation means were substantially less than zero, and type I error rates approached 100%.

In scenario 5 (true change in β₁), all methods detected the temporal change in the socioeconomic relationship with the outcome, suggesting that all methods are sensitive, if not specific. Power to detect a significant temporal change (given that there was a true change of 0.3) was somewhat less with methods IV, V, and VI.

In scenario 6 (change in α₁, , , and β₁), methods I, II, and III did detect a significant change in the socioeconomic gradient over time. The combined effect of the multiple varying parameters on the estimated temporal change in the association of Y and X was to attenuate it relative to what would be estimated if we could regress Y on Z. However, methods IV, V, and VI all performed poorly, with simulation means going in the opposite direction of the true change in socioeconomic gradient and virtually nonexistent power for the decile and SII approaches.

Over all the simulations, method III (validation data) consistently yielded estimates of the temporal effect that were virtually identical to those obtained by the “true” regression of Y on Z. The slightly increased variability over 1,000 simulations shows the additional uncertainty induced by using a 30% subsample as validation data and applying the imputation model to the rest of the data set.