Characterizing the Effects of Environmental Exposures on Social Mobility: Bayesian Semi-parametrics for Principal Stratification

Principal stratification provides a robust causal inference framework for the adjustment of post-treatment variables when comparing the effects of a treatment (Frangakis and Rubin, 2002). The post-treatment variables act as intermediate variables between treatment and outcome, providing valuable insights into causal mechanisms.

This framework focuses on classifying subjects into principal strata based on the joint values of the potential outcome for the post-treatment variable at different treatment levels (Antonelli et al., 2023). These strata represent inherent latent characteristics that remain unaffected by treatment, enabling causal analysis through the definition of principal causal effects (Rubin, 1974a). The literature distinguishes between the dissociative effect—representing the treatment effect for the stratum composed of subjects with the post-treatment variable that is not affected by the treatment—and the associative effects, which captures the treatment effect for strata where the post-treatment variable exhibits significant changes between the different treatment levels (Frangakis and Rubin, 2002; Mealli and Mattei, 2012).

Most studies focus on binary treatments and post-treatment variables, simplifying analyses to four combinations (Angrist, Imbens and Rubin, 1996). The extension to continuous post-treatment variables poses interpretational and inferential challenges due to infinite possible strata. Solutions include dichotomizing post-treatment variables (Sjölander et al., 2009; Jiang, Yang and Ding, 2022), parametric or semiparametric modeling (Jin and Rubin, 2008; Conlon, Taylor and Elliott, 2014; Lu, Jiang and Ding, 2023), and nonparametric approaches (Schwartz, Li and Mealli, 2011; Zorzetto et al., 2024a; Kim and Zigler, 2024). Recently, Antonelli et al. (2023) introduced a Bayesian methodology to tackle the continuous treatment issue, advancing nonparametric modeling in this domain.

In this paper, we define a novel Bayesian approach for principal stratification that leverages the flexibility of the Dirichlet process (Ferguson, 1973) as prior. Specifically, we exploit the dependent Dirichlet process (DDP) (Quintana et al., 2022) to define a highly adaptable distribution for the potential outcome for the final outcome, that depends on the confounders, the potential outcome for the post-treatment variable, and the treatment variable.

The current state-of-the-art for Bayesian nonparametric (BNP) in the context of principal stratification with a binary treatment can be found in the works by Schwartz, Li and Mealli (2011) and Zorzetto et al. (2024a). These models, as well as the one we propose in this paper, are all designed for a principal stratification framework where the post-treatment variable and the outcome are continuous, and the treatment is a binary random variable, however they differ in the definition of the nonparametric priors.

First, our proposed model introduces the DDP as prior for the parameters in the potential outcome for final outcome model distribution, while the Schwartz, Li and Mealli (2011) and Zorzetto et al. (2024a)’s models have a DDP prior for the parameters in the post-treatment variable model distribution. Second, the dependency on the covariates in the DDP differs in the three models: in the Schwartz, Li and Mealli (2011)’s model the dependence involves only the atoms, on the model proposed by Zorzetto et al. (2024a) the confounders are included in the weights of the mixture to improve the imputation of missing data. In our proposed model the dependence on the confounders is both on the atoms and on the weights.

The double dependence on confounders in our proposed model is motivated by several key reasons: (i) the dependence in the weights enables fair and informed imputation of missing data, as demonstrated by Zorzetto et al. (2024b) and Zorzetto et al. (2024a) in various causal inference settings; (ii) incorporating confounders and the potential outcomes of the post-treatment variable within the mixture regression enhances the model’s adaptability; and (iii) the flexible structure of the mixture model for the potential outcomes for the final outcome relaxes the rigid assumptions previously imposed on the distribution of the post-treatment variable’s potential outcomes—e.g., the linear regression model—, that induces a mixture distribution for the predictive distribution used in imputing missing data for post-treatment variable.

Following this introduction, Section 2 introduces the motivating application and the causal real-data questions that we aim to address. In Section 3, we introduce the notation, assumptions that we employ throughout the paper, and the principal causal effects. We conduct a simulation study to assess the performance of our proposed model in different scenarios in Section 5. The description of the dataset and the results of the application of are reported in Section 6. Section 7 concludes the paper with a discussion of the proposed model and further research.

Recent empirical evidence suggests a concerning trend in US social mobility, with studies documenting stagnation or decline in intergenerational economic mobility (Piketty, Saez and Zucman, 2018; Song et al., 2019). In a very influential work, Chetty et al. (2017) found that the fraction of children who earn more than their parents has fallen from approximately 90% for children born in 1940 to 50% for children born in the 1980s.

The decline in intergenerational mobility raises significant societal concerns, as reduced economic mobility can weaken social cohesion. In fact, as argued by Stiglitz (2012), the growing inequality often equates to diminishing equality of opportunity and exacerbating economic and social inequalities.

Recent literature has focused on searching for the root causes of the decline in economic mobility. While multiple factors—such as cultural (Platt, 2019), demographic (Van Bavel et al., 2011; Salvanes, 2023), geographical (Connor and Storper, 2020; Salvanes, 2023), labor market (Choi, Kim and Kim, 2023), welfare state (Heckman and Landersø, 2022) as well as macro-economic secular trends (Piketty, Saez and Zucman, 2018)—play an important role, recent works have highlighted that environmental and educational factors are key determinants of the decline in social mobility.

On the one hand, in a recent contribution Lee, Merlo and Dominici (2024) identified a direct causal link between pm_2.5 exposure and social mobility, finding that a $1\mu g/m\textthreesuperior$ increase in childhood exposure to fine particulate matter (pm_2.5 ) leads to a $1.146\%$ decrease in absolute upward mobility. This seminal study highlights the potential for pm_2.5 to directly hinder socioeconomic advancement. However, this study falls short in investigating the potential causal pathways through which exposure to higher levels of air pollution hinders social mobility.

On the other hand, educational attainment is a factor influencing both economic and social mobility. Empirical evidence consistently demonstrates its significant returns in terms of earnings and intergenerational mobility. For instance, studies have shown that higher levels of education correlate with increased earnings potential, as highlighted by Card (2001), who discusses the effects of education on earnings. Blanden et al. (2014), Pfeffer and Hertel (2015) and Brown and James (2020) provide evidence on how educational attainment can enhance intergenerational mobility, allowing individuals from lower socio-economic backgrounds to improve their status across generations.

As highlighted in the previous section, the relationship between social mobility, education, and air pollution has been the focus of a growing body of research. However, the literature has so far neglected to study the interconnectedness of environmental and educational factors and social mobility.

Air pollution significantly impacts educational outcomes, both in the short and long term. Recent studies have revealed that long-term exposure to pollutants earlier in life can lead to decreased academic performance and higher rates of cognitive disabilities and mental health problems (Zhang, Chen and Zhang, 2018; Braithwaite et al., 2019). These findings are further supported by studies indicating that even short-term air pollution exposure can impair cognitive function, reduce concentration, and increase absenteeism among students (Sunyer et al., 2015; Shehab and Pope, 2019).

The effects of pollution on education are not limited to individual performance. Air pollution may also contribute to the exacerbation of existing inequalities, since schools in more polluted areas often serve disadvantaged populations (Mohai et al., 2011). Furthermore, there is critically important evidence that suggests that the economic impact of pollution-related educational deficits can be substantial, with significant losses in lifetime earnings for affected students (Currie, Neidell and Schmieder, 2009), which could, in turn, substantially hinder social mobility.

To disentangle the complex causal relationships linking air pollution (as an exposure), education (as a post-treatment), and social mobility (as a final outcome), we harmonize multiple publicly available datasets and develop a methodology to analyze such data within a coherent principal stratification framework.

We assume to observe $n$ independent and identically distributed (iid) units. For each unit $i\in\{1,\dots,n\}$, we let ${\bf X}_{i}\in\mathcal{X}\subseteq{\mathds{R}}^{P}$ be the vector of observed covariates, $T_{i}\in\{0,1\}$ be the observed binary treatment variable, let $P_{i}\in{\mathds{R}}$ be a post-treatment (or intermediate) variable, and $Y_{i}\in{\mathds{R}}$ be a final outcome of interest. Capital Roman letters denote random variables, realized values are denoted in lower case—e.g., $t$ are the observed values for the random treatment $T$, while the bold letters denote vectors—such as the $p$ dimensional vector for the confounder random variables ${\bf X}$.

In the case of our application, the units are the counties in the continental United States, the treatment variable is the pm_2.5 exposure dichotomized with respect to the median of the observed values as the threshold, the post-treatment variable is a continuous variable that indicates the education level, the outcome of the social mobility, and the confounders are socio-economic and meteorological information. Further detail on variables definition are provided in Section 6.

Following to Rubin’s Causal Model (Rubin, 1974b; Holland, 1986), we postulate the existence of two potential outcomes for the final outcome, $\{Y_{i}(0),Y_{i}(1)\}\in{\mathds{R}}^{2}$, for each $i=1,\dots,n$, represent the collection of the final outcomes when the unit $i$ is assigned to the control group—i.e., when $t=0$—or the treatment group—i.e., when $t=1$—, respectively. Similarly, following the principal stratification framework (Frangakis and Rubin, 2002), we also assume the existence of two potential outcomes for the post-treatment variable, $\{P_{i}(0),P_{i}(1)\}\in{\mathds{R}}^{2}$, which represent a collection of outcomes for the post-treatment variable when the unit $i$ is assigned to the control or treatment group, respectively.

We need to invoke the stable unit treatment value assumption (SUTVA) to properly define the potential outcomes for the final outcome of interest and the post-treatment variables.

Specifically, SUTVA is a combination of two assumptions: no interference between units—i.e., the potential values of the final outcome and post-treatment variable of the unit $i$ do not depend on the treatment applied to other units—and consistency—i.e., no different versions of the treatment levels assigned to each unit (Rubin, 1986).

In the principal stratification framework (Frangakis and Rubin, 2002; VanderWeele, 2011; Mealli and Mattei, 2012), we assume that the treatment affects not only the final outcome but also the post-treatment variable. Therefore, we want to study the causal effect of the treatment on the final outcome through specific values assumed by the potential outcome of the post-treatment variable. Specifically, we are interested in estimating the causal effects given the principal strata, where the principal strata are defined by groups of units with similar (in some literature also identical) combinations of the potential outcome for post-treatment variables (Frangakis and Rubin, 2002; Mealli and Mattei, 2012; Feller, Mealli and Miratrix, 2017; VanderWeele, 2011). Due to its ability to establish—in an elegant mathematical manner—the causal pathways through which the treatment affects the final outcome via the post-treatment variable, this approach has been extensively adopted in applied statistics (see, e.g., VanderWeele, 2011; Mealli and Mattei, 2012; Feller, Mealli and Miratrix, 2017; Ding et al., 2017; Lu, Jiang and Ding, 2023; Zorzetto et al., 2024a).

However, the continuous nature of the potential outcome for post-treatment variables makes the definition of the strata challenging due to the infinite combination of the values for $(P(0),P(1))$. Thus, following the seminal contribution of Zigler, Choirat and Dominici (2018), which also proposes an application of principal stratification for environmental exposures, we define the units allocation in the principal strata introducing a threshold $\xi$, such that each unit $i$ belongs to he associative positive stratum if $P_{i}(1)-P_{i}(0)\geq\xi$, to the associative negative stratum if $P_{i}(1)-P_{i}(0)\leq-\xi$, or to the dissociative stratum otherwise.

According with that, the causal estimands proposed by Zigler, Choirat and Dominici (2018) are respectively the Expected Associated Effect for the positive stratum (EAE₊) and for the negative stratum (EAE_-) and the Expected Dissociative Effect (EDE):

The underlying motivation of causal estimands definition is their relevance and interpretability in the real-data application. In particular, the EDE defines the causal effect of the pm_2.5 exposure on social mobility given that the counties where the pm_2.5 exposure does not affect the education level. Conversely, the EAE_- (EAE₊) quantifies the causal effect of the pm_2.5 exposure on social mobility given that the counties where the pm_2.5 exposure decreases (improves) the education level. These estimands are particularly important because they allow us to quantify and understand the heterogeneous pathways through which pm_2.5 exposure affects social mobility. For example, we can evaluate if exposure to a higher level of air pollution can degrade, improve, or not change the social mobility, for three important strata (in this peculiar case, corresponding to groups of counties): where the education level is not modified by the air pollution exposure, where the percentage of the people obtaining the diploma/degree is decreasing (or increasing) if the students are exposed to higher level of air pollution.

In order to identify these causal estimands, a second assumption needs to be introduced.

The strongly ignorable treatment assignment declares that the potential outcome for the final outcome and the post-treatment variable are independent of the treatment conditional on the set of covariates and all units have a positive chance of receiving the treatment. In our application, this means that the potential values of the social mobility and the education level under the two levels of exposure are independent of the (underlying) mechanism that controls the effective level of pm_2.5 in each county, conditional to the socio-economic and weather characteristics. Moreover, the values of socioeconomic and weather characteristics for any county cannot preclude the possibility of observing one of the two levels of pm_2.5 exposure.

Leveraging the Assumptions 1 and 2, we can rewrite the principal causal estimand defined in (1) as the following statistical estimand:

Indicating with $f(P_{i}(1),P_{i}(0))$ the general definition of the three strata, then the inner expectations in each statistical estimand $\mathds{E}\big{[}Y_{i}\mid T_{i}=t,{\bf X}_{i}={\bf x},f(P_{i}(1),P_{i}(0))\big{]}$ is estimated with the outcome model $\{Y_{i}\mid T_{i},{\bf X}_{i},P_{i}(1),P_{i}(0)\}$ defined by (4) in Section 4. Similarly, the probability can be rewritten as follows:

where $\Pr(P_{i}(1)-P_{i}(0)\geq\xi\mid{\bf X}_{i}={\bf x})$ is estimated with the potential post-treatment variables model defined in (8), while the $\Pr({\bf X}_{i})$ is observed.

Following the Bayesian paradigm, the joint distribution of confounders ${\bf X}$, treatment $T$, potential outcome for post-treatment variable $(P(0),P(1))$, and potential outcome for the final outcome $(Y(0),Y(1))$ can be rewritten as follows:

where $p(\theta)$ is the prior distribution for all the involved parameters $\theta$ that take values in the parametric space $\Theta$. The inner probability distribution can be factorized into:

The SUTVA and strong ignorability assumptions, defined in Section 3, allow us to simplify the conditional probability for the treatment variable as dependent only to the confounders, such as $\Pr(T_{i}\mid Y_{i}(0),Y_{i}(1),P_{i}(0),P_{i}(1),{\bf X}_{i},\theta)=\Pr(T_{i}\mid{\bf X}_{i},\theta)$. Following Schwartz, Li and Mealli (2011), the treatment and covariate distributions—accounting for the empirical distribution then $\Pr({\bf X}_{i}\mid\theta)=\Pr({\bf X}_{i})$—are directly observed, so that they do not require to be modeled. While the distribution of the potential outcome of the post-treatment variable, conditioned on confounders and treatment, and the distribution of potential outcome for final outcome, conditioned on potential outcome for post-treatment variable, confounders and treatment, require to be modeled.

Here, motivated by our data application, we want to define a flexible model that allow us to capture the heterogeneity in the data and impute correctly the missing post-treatment variables and the potential outcome for final outcomes, such as it can allow us to estimate the principal causal effects defined in the Section 3.2. To achieve this goal we take advantage of the Bayesian nonparametric (BNP) priors and in particular the dependent Dirichlet process (DDP) (MacEachern, 2000; Barrientos, Jara and Quintana, 2012; Quintana et al., 2022). These prior are well-known for they flexibility to capture complex heterogeneity in the data and adaptability in different real-data context. Moreover, in the Bayesian framework the imputation of missing data is straightforward.

Specifically for principal stratification, few works leverage Dirichlet process as prior. The seminal contribution of Schwartz, Li and Mealli (2011) defines a Dirichlet process mixture model for the potential outcome for post-treatment variable. More recently, Zorzetto et al. (2024a) define a DDP prior that share information between the two potential outcome for post-treatment variable distributions. However, both works define the nonparametric priors for the potential outcome for post-treatment variable distribution, while assuming a linear regression model for the potential outcome for final outcome distribution.

Reminding that the main purpose on principal stratification framework is to estimate the principal causal effects—defined in (2), as difference of potential outcome for final outcomes given the strata—, our focus is into potential outcome for final outcome distribution. Therefore, we model the potential outcome for final outcome probability distribution for the final outcome, by exploiting the DDP as prior, such that for each unit $i=1,\dots n$:

where ${\mathcal{Q}}(\cdot;{\bf x},p(0),p(1),\theta)$ represents a continuous density function for every $\theta\in\Theta$ which depends on both the confounders ${\bf x}$ and on the potential outcome for post-treatment variables $\big{(}p(0),p(1)\big{)}$. For each treatment level $t\in\{0,1\}$, also the the discrete random probability measure $H_{{\bf x},p(0),p(1)}^{(t)}$ depends on both the confounders ${\bf x}$ and the outcomes for post-treatment variables $(p(0),p(1))$. The characterization of the random probability measure allows us to rewrite it as an infinite mixture as follows:

where the infinite sequences of the weights $\{\lambda^{(t,m)}(\cdot)\}_{m\geq 1}$ and of the atoms $\{\mathbb{\theta}^{(t,m)}(\cdot)\}_{m\geq 1}$ are stochastic processes. Both stochastic processes depend on the confounders ${\bf x}$, and the atoms depend also on the post-treatment variables. This definition allow us to make the potential outcome for final outcome distribution dependent to both confounders and post-treatment variable, essential to distinguish the principal strata, while the confounders in the infinite sequences of the weights allow to characterize the heterogeneity in the potential outcome for final outcomes and improve the imputation of the missing data.

Following the stick-breaking representation (Sethuraman, 1994) and, in particular, exploiting the Probit stick breaking process (PSBP) introduced by (Rodriguez and Dunson, 2011), the random weights $\{\lambda^{(t,m)}({\bf x}_{i})\}_{m\geq 1}$, for each unit $i$, can be defined as following:

where $\Phi(\cdot)$ describes the cumulative distribution function of a standard Normal distribution and $\{\gamma^{(t,m)}({\bf x})\}_{m\geq 1}$ has Gaussian distributions with mean a linear combination of the confounders ${\bf x}$, such that $\epsilon^{(t,m)}$ are the regression parameters, including the intercept, for each cluster $m$ and treatment level $t$.

The weights $\{\lambda^{(t,m)}({\bf x})\}_{m\geq 1}$ implicitly describe the probability to belong to a specific cluster $m$, under treatment $t$ and given confounders ${\bf x}$. Therefore we can introduce a latent categorical variable $V_{i}^{(t)}$, for each unit $i\in\{1,\dots,n\}$ and each treatment level $t\in\{0,1\}$, such that

We can rewrite (4), for the final outcome $Y_{i}$ given the cluster allocation $V_{i}^{(t)}=m$ as:

Furthermore, assuming the kernel $\mathcal{Q}(\cdot;x_{i},p_{i}(0),p_{i}(1),\theta)$ be normally distributed, we define the distribution in (6) for the two treatment levels as follows:

where $\{\eta^{(t,m)}\}_{m\geq 1}$ are the regression parameters for the mean of the Gaussian distribution and $\{{\sigma_{y}^{(t,m)}}\}_{m\geq 1}$ are the scale parameters, such that $\theta_{m}=(\eta^{(t,m)},{\sigma_{y}^{(t,m)}}^{2})$ in (6).

Moreover, we assume for the potential outcome for post-treatment variables distribution a linear regression with normal errors, such that the post-treatment variables are independent given the treatment and dependent to the confounders ${\bf X}$, such that:

The simulation study is designed to allow us to test our proposed model’s performance to (i) impute the missing post-treatment variables, (ii) impute the potential outcomes for the final outcome, and (iii) correctly estimate the causal effects. To evaluate that we estimate the bias for the Sample Average Treatment Effect (SATE) for the potential outcome for post-treatment variable and for the potential outcome for the final outcome, respectively:

where $n_{s}$ is the sample size of each generated sample. Please note that, since we are in a simulated setting, SATE is observed as we have knowledge of both potential outcome for final outcome under treatment and control.

We evaluate our proposed model under three data generating processes, where we generate different levels of heterogeneity in post-treatment variables and final outcomes. Specifically, the third Scenario mimics the characteristics of the dataset analyzed in the data application in Section 6. The results are compared with the semiparametric model for principal stratification introduced by Schwartz, Li and Mealli (2011).

For each of the following three Scenarios, we consider 200 replicates. Each replicate has $500$ units in the first and second Scenario, and $300$ units in the third Scenario—similarly to the analyzed dataset in the Section 6.

Scenario 1: We define five confounders $X_{1:5}$—two Bernoulli random variables and three standard Gaussian random variables—and a Bernoulli treatment variable such us the as $T\sim\mbox{Be}(f_{1}(X_{1:5}))$. The potential outcome for post-treatment variables are sampled from a linear regression normally distributed that differs for the two treatment levels, while the potential outcome for final outcomes are sampled from a mixture of Gaussian distributions. The units are divided into three clusters $m\in\{1,2,3\}$ according to the values of the confounders $X_{1:2}$–see further details in the Supplementary Material—, that determine the allocation of the units to the different components in the mixture distribution of the potential outcome for final outcomes. The variable $M_{i}$ indicates the cluster allocation for each unit $i\in\{1,\dots,n\}$. Specifically, we use the following distributions to sample potential outcome for post-treatment variables and potential outcome for final outcomes for each unit $i\in\{1,\dots,n\}$:

where the function $g_{0,1}(\cdot)$ and $g_{1,1}(\cdot)$ are nonlinear functions, and $\mathbb{I}_{\{\cdot\}}$ is an indicator variable.

Scenario 2: We increase the number of confounders to ten, and the treatment level is defined as $T\sim\mbox{Be}(f_{2}(X_{1:10}))$. The cluster allocation depends on the binary confounders $X_{1:2}$ as in Scenario 1, while the potential outcome for post-treatment variable distribution depends on the remaining confounders $X_{3:10}$. Choosing different confounders for the cluster allocation and to include in the potential outcome for post-treatment variable distribution induces a different heterogeneity scheme compared with Scenario 1. The potential outcome for final outcome distribution is defined similarly to the Scenario 1, but including all the ten confounders in the regression of the means in the Gaussian distributions in the mixture. See further details in the Supplementary Material.

Scenario 3: This Scenario closely mimics the characteristics of the dataset used in the application in the Section 6. We reduced the number of units and increased once again the number of covariates to 14, some of them Bernoulli distributed and some others normally distributed with different variances. The treatment variable, the cluster allocation, and the potential outcome for post-treatment variable distribution are defined similarly to Scenario 2 but include a bigger number of covariates. See the Supplementary Material for details.

The Gibbs sampler—implemented in R and available in GitHub—has the option to define different covariates matrices for the post-treatment variable regression, final outcome regression, and regression on the weights—i.e., the covariates that describe the cluster allocation.

The Table 1 reports the comparison of the bias for $SATE_{P}$ and $SATE_{Y}$ between our proposed model—indicated with Y_BNP—and the Schwartz, Li and Mealli (2011)’s model—indicated with SLM. The median and the interquartile range (IQR) are estimated over the $200$ replicates for each Scenario. See the supplementary materials for further visualization of the results.

The Y_BNP model exhibits superior performance in terms of bias and variability across all Scenarios when compared to the SLM model. Specifically, the Y_BNP model estimates values closer to zero for the median and a smaller IQR compared with the SLM model in each of the three Scenarios and for both the causal quantities $SATE_{P}$ and $SATE_{Y}$.

The results for $SATE_{P}$ show that the Y_BNP model can capture well the heterogeneity of the post-treatment variable distribution and imputed correctly the missing data. Indeed, the median bias has values across all Scenarios, ranging from -0.0129 to -0.0182, and narrow IQRs that include the zeros values for all three Scenarios. This indicates that the Y_BNP model performs well on the post-treatment variable, avoiding substantial over- or underestimation, even though we assume a linear regression model. As explained previously, however, the predictive distribution of the potential outcome for post-treatment variable is not still a Gaussian distribution with linear regression but a mixture of them. In contrast, the SLM model shows larger median bias values, particularly notable in Scenarios 2 (0.4037) and 3 (0.4340), and the IQRs are almost three times those for the Y_BNP model.

Clear differences between the two models are underlined by the results for $SATE_{Y}$. Again, the results for the Y_BNP model are more accurate and show smaller median bias values—with a range between -0.0333 and -0.0896—and a controlled variability—the IQRs take values between 0.0818 and 0.1790. While the results for SLM show that the bias and the variability increase significantly. Scenarios 2 and 3 have a median of the bias more than ten times the median estimated by the Y_BNP model and IQR five times the one estimated with the Y_BNP model. The high values for Scenario 1 show the difficulty of the SLM model to capture a more complex heterogeneity in the data and controlling the variability propagation between post-treatment variable and final outcome. Indeed, the SLM model defines a flexible model for the post-treatment variable, but it is not enough due to the linear model used for the final outcome distribution.

The accuracy and the stable estimation of treatment effects are particularly important in the context of principal stratification analysis, where both the missing post-treatment variable and the potential outcome for final outcomes have to be imputed. Leveraging Bayesian nonparametric prior to the final outcome offers a more robust and reliable framework to capture the heterogeneity in the data and control the variability propagation that can be present in this type of Scenario.

The dataset used in this study is constructed from multiple publicly available sources, mimicking and building of the dataset used by Lee, Merlo and Dominici (2024). The socio-economic and demographic data were obtained from the U.S. Census (1990–2000), while meteorological variables were sourced from the Daymet dataset (1982) (Castro, Yazdi and Schwartz, 2024). Air pollution data, specifically pm_2.5 exposure levels, were derived from high-resolution satellite estimates at a $0.9km\times 1.1km$ grid scale (Colmer, Voorheis and Williams, 2023), which was then aggregated to the census tract level. Social mobility data, used as the final outcome, and different education levels measured as a rate over the considered population, used as post-treatment variables, were taken from the Opportunity Atlas (Chetty et al., 2018). These combined sources provided the comprehensive dataset used for our analysis.

The study was conducted at the county level in the continental U.S., encompassing 3,009 counties. Although most variables were initially available at the census tract level, all data was aggregated to the county level to align with education levels, which are only available at this scale.

The final dataset for analysis retains several key variables listed in 2, where their descriptive statistics and data sources are reported. Census data includes the percentage of college-educated in the year 2000, the median household income in 1990, the population density in the year 2000 per square mile, the share of people who live below poverty levels, and the share of black, white Hispanic, and Asian. Meteorological data included average daily minimum and maximum temperatures for winter (December to February) and summer (June to September), along with seasonal precipitation for both seasons. Spatial controls included a four-level categorical variable based on US. Census regions (Northeast, South, Midwest, and West). The treatment variable was pm_2.5 exposure levels for 1982, measured in $\mu g/m^{3}$ at the census tract level and aggregated to the county level. Specifically, we observe the mean of $17.03\mu g/m^{3}$ and the median of $17.47\mu g/m^{3}$. Binarization of pm_2.5 —which helps reducing the complexity of the exposure which is continuous in nature—has been previously explored and justified in the literature (Lee, Small and Dominici, 2021; Bargagli-Stoffi et al., 2020; Zorzetto et al., 2024b). The post-treatment variables include the community college graduation rate, the high school graduation rate, and the college graduation rate, chosen for their relevance to the literature on social mobility and data availability. The final outcome is the Absolute Upward Mobility (AUM), which is defined as the mean income percentile at adulthood of individuals born between 1978 and 1983 into families at the $25^{th}$ percentile of the national parent income distribution. Income rank is measured in 2015 (ages 31–37).

A preliminary analysis shows that the community college graduation rate has the strongest correlation with AUM ($66.45\%$) and a significant negative correlation with pm_2.5 exposure ($-38.19\%$). High school and college graduation rates also exhibited significant correlations with AUM ($55.55\%$ and $53.47\%$, respectively) and negative correlations with pm_2.5 exposure ($-36.32\%$ and $-28.86\%$, respectively).

We applied our proposed model in Section 4.1 to the constructed dataset. Specifically, we replicate the analysis for the three education levels: high school, community college, and college. The three analyses share the confounders—i.e., socio-economic characteristics, meteorological information as well as the spatial confounder—; the treatment, that is defined as exposure to high or lower level of pm_2.5—the lower (high) level when pm_2.5 is below (above) the threshold of $17.47\mu g/m^{3}$—; and the final outcome represented by social mobility which is measured by AUM.

We assume the threshold $\epsilon$ equal to $0.01$ for the principal strata, defined in (1), for all three education rates. This choice is justified by a change in education level that is deemed notable considering the range of the estimated difference of the education rates under treatment and under control. Table 3 shows that the three strata are balanced in all three educational levels. More than $50\%$ of the units fall into the associative negative stratum, indicating that for these cases, the pm_2.5 exposure has a negative effect on the education level. In this stratum, a county exposed to a higher level of pm_2.5 exhibits lower educational levels compared to if it had been exposed to a lower pm_2.5 level. This observation aligns with findings previously reported in the literature, supporting its plausibility (see e.g., Brochu et al., 2011; Bell and Ebisu, 2012; Binelli, Loveless and Whitefield, 2015; Hajat, Hsia and O’Neill, 2015; Grunewald et al., 2017; Yang and Liu, 2018). A percentage from $10.4\%$ and $14.8\%$—in the three education levels—of the analyzed counties are included in the dissociative stratum. I.e., the counties where the pm_2.5 exposure does not affect the education level. While, in all three education levels, approximately $30\%$ of units belong to the associative positive stratum. This stratum includes counties that experience increased education rates when exposed to higher levels of pm_2.5.

Figure 1 reports the posterior distribution of the principal causal effects of pm_2.5 exposure on social mobility (AUM) using as post-treatment variable the education level, respectively high school, community college, and college.

The posterior distributions of the $EAE_{-}$ have a mean of -3.9%, -5.0%, and -4.7%, respectively, for the analysis with high school, community college, and college as post-treatment variables. This means that the counties where the exposure to high levels of pm_2.5 reduces significantly the education rate for community college also induce a reduction in the social mobility of 5%. With similar meanings for high school and college. The credible intervals at the 95% level are $(-5.6\%,-2.1\%)$, $(-6.8\%,-3.1\%)$, and $(-6.5\%,-2.8\%)$, for the three education levels.

Also, the posterior distributions of the $EDE$ confirm the negative effect of the pm_2.5 exposure and social mobility. In fact, we observe for each educational level a posterior mean equal to -2.3%, -1.5%, -2.0% and a credible interval at the 95% level equal to $(-4.4\%,-0.3\%)$, $(-3.6\%,0.7\%)$, and $(-4.0\%,0\%)$. The zero is only included in the credible interval for the community college analysis. Therefore, also for the counties where the education level is not affected by the pm_2.5 exposure—defined as the dissociative stratum—the social mobility is reduced when exposed to higher levels of air pollution.

Lastly, we estimated posterior distributions for $EAE_{+}$ that gravitate around 0. Specifically, the posterior mean is equal to -0.9%, 0.3%, and 0%, with credible intervals at the 95% level that include the zero for all three education variables. These results indicate that for the counties where the education level increases with a higher pm_2.5 level, then the social mobility is not significantly affected by air pollution exposure.

Overall, we can conclude that the pm_2.5 exposure does have a negative effect on social mobility, coherently to the conclusion of Lee, Merlo and Dominici (2024). However, the effects are negative even for the dissociative strata, meaning that even when pm_2.5 does not hinder educational attainments, the effects on social mobility are still negative. The principal causal effect has a bigger range when the intermediate variable is the college rate and a smaller range for the high school graduation rate.

Analyzing the distribution of characteristics between strata provides insights for designing targeted social policies to promote equity within the population. Figure 2 highlights that the primary factors that differ between strata include weather conditions, population density, median household income, and ethnic composition (percentage of white and black populations).

Weather conditions exhibit varying patterns based on education levels: counties with higher temperatures and more rainfall are associated with the positive associative stratum when high school is the mediator of the causal effect of exposure to pm_2.5 on social mobility. In contrast, for college and community college as mediators, these weather features define the negative associative stratum.

A consistent finding across the three levels of education is that counties with higher population density and higher median household income—typically urban areas—are more adversely affected by pm_2.5 exposure, leading to significant reductions in social mobility (negative associative stratum). Furthermore, counties with a higher percentage of white residents are predominantly assigned to the positive associative stratum, aligning with literature that indicates that white populations are less affected by air pollution and are more likely to experience upward social mobility. In contrast, counties with a higher proportion of Black residents are concentrated in the negative associative stratum, emphasizing the need for targeted social policies to mitigate the increased risks faced by black populations, particularly in urban areas, due to air pollution.