Bayesian Nonparametric Trees for Principal Causal Effects

The approach is formalized within a potential outcomes framework (Rubin, 1974). Let $(A_{i},M_{i},Y_{i})$ denote the exposure, intermediate, and outcome variables of interest, and $\bm{X}_{i}$ be a set of confounders for observation $i$, all conceived as fixed pre-treatment features. Throughout the paper, we assume the binary exposure. Let $M_{i}(\bm{A})$ denote the potential value of the intermediate variable that would be observed under the vector of exposures $\bm{A}=(A_{1},\ldots,A_{n})$ for subjects $i=1,\ldots,n$. Similarly, $Y_{i}(\bm{A})$ denotes the potential value of the outcome that would be observed under the vector of exposures $\bm{A}=(A_{1},\ldots,A_{n})$. Under the stable unit treatment value assumption (SUTVA, Rubin (1980)), potential intermediate and outcome values for subject $i$ do not depend on exposures assigned to other subjects (i.e., $M_{i}(\bm{A})=M_{i}(A_{i})$ and $Y_{i}(\bm{A})=Y_{i}(A_{i})$) and are essentially equal to the observed intermediate and outcome values (i.e., $M_{i}=M_{i}(A_{i})$ and $Y_{i}=Y_{i}(A_{i})$). The exposure ($A$) in the motivating air quality study, is the presence of a flue-gas desulfurization scrubber (henceforth, “scrubber”), measured at each coal-fired power plant in 2014; the intermediate variable ($M$) is the amount of SO₂  emissions from the power plant during 2014; and the outcome of interest is the annual ambient average PM_2.5  concentration within a 50km radius of the power plant ($Y$). $M_{i}(A_{i})$ and $Y_{i}(A_{i})$ are the potential SO₂  emissions and surrounding PM_2.5  values based on the state of scrubber installation ($A_{i}$) at the $i$-th power plant.

Principal stratification (Frangakis and Rubin, 2002) partitions the population based on how $A_{i}$ affects the intermediate variable, encoded as membership in a principal stratum defined as $S_{i}=(M_{i}(0),M_{i}(1))$. This principal stratum membership $S_{i}$ is a latent variable that is unaffected by exposure and can be considered the baseline characteristic of the $i$-th subject, encoding features that dictate how the primary treatment impacts the intermediate variable. In the power plant case, $S_{i}$ denotes a combination of the effectiveness of scrubbers for reducing SO₂  as well as other plant characteristics dictating variability in this effect that may also relate to the formation of ambient PM_2.5  such as patterns of operation, responses to regional and seasonal regulatory programs, or atmospheric and weather conditions surrounding the plant. Our goal is to examine causal effects of $A_{i}$ on $Y_{i}$ across principal strata defined by different values of $(M_{i}(0),M_{i}(1))$, known as principal causal effects (PCEs, Frangakis and Rubin (2002)).

PCEs can be described as a surface of $\mathbb{E}[Y(1)-Y(0)|M(0),M(1)]$ across values of $(M(0),M(1))$, which has been termed the “causal effect predictiveness” surface (Gilbert and Hudgens, 2008). Summarizing this surface can be challenging with continuously-scaled $M$, and deriving summary estimands will depend on context. Expected dissociative effects have been proposed and used to summarize average causal effects of $A$ on $Y$ among principal strata for which $A$ has little or no impact on $M$, representing something similar to a “direct effect” of $A$ on $Y$ when $M$ remains unaffected. Expected associative effects can summarize average effects of $A$ on $Y$ across strata where $M(0)\neq M(1)$. Here, we extend previous work in Kim et al. (2020) to define multiple average principal effects based on a pre-determined threshold vectors, $\bm{\epsilon}=(\epsilon_{\text{lower}},\epsilon_{\text{upper}})$ and define expected PCEs as:

where various values of $\bm{\epsilon}$ will define PCEs among strata with different impacts of $A$ on $M$. For example, specifying $|\epsilon_{\text{lower}}|=|\epsilon_{\text{upper}}|\approx 0$ would correspond to an expected dissociative effect among strata of power plants for which scrubbers have little to no effect on SO₂  emissions. If we consider the sample average estimand as the target estimand, the expression for $\Delta_{\text{PCE}}(\bm{\epsilon})$ is as follows:

where $n_{\epsilon}=\mathbb{I}\{\sum_{i=1}^{n}\epsilon_{\text{lower}}<M_{i}(1)-M_{i}(0)<\epsilon_{\text{upper}}\}$ and $\mathbb{I}\{\}$ is an indicator function.

Note that the fundamental problem of causal inference applies in this case to both the values of $Y_{i}(A)$ and $M_{i}(A)$. Thus, the principal stratum membership $S_{i}$ is unobserved and quantities such as $\Delta_{\text{PCE}}(\bm{\epsilon})$ are not identified based only on observed data. Unobserved values $M_{i}(1-A_{i})$ will determine the value of $S_{i}$ and must be estimated along with the unobserved values of $Y_{i}(1-A_{i})$.

To estimate the principal causal effects, we first formulate an assignment mechanism, the distribution of $\bm{A}$ conditional on the potential (intermediate) outcomes $(\bm{Y}(0),\bm{Y}(1),\bm{M}(0),\bm{M}(1)$) and observed covariates where $\bm{Y}(a)=(Y_{1}(a),\ldots,Y_{n}(a))^{\prime}$ and $\bm{M}(a)=(M_{1}(a),\ldots,M_{n}(a))^{\prime}$ for $a=0,1$. Under the assumption that the exposure assignment for each unit is independent of information from other units, the exposure assignment we consider for unit $i$ has the following form: $\text{Pr}(A_{i}|Y_{i}(0),Y_{i}(1),M_{i}(0),M_{i}(1),\bm{X}_{i})$. The ignorability assumption (Rosenbaum and Rubin, 1983) states that with a proper set of covariates $\bm{X}_{i}$, the exposure assignment mechanism can be reduced to a simpler form.

Our target estimands are sample average PCEs for the $n$ sample units. Unlike corresponding population average estimands, the parameters governing associations between $Y_{i}(0)$ and $Y_{i}(1)$ and $M_{i}(0)$ and $M_{i}(1)$ cannot be ignored (Schwartz et al., 2011; Ding and Li, 2018; Kim et al., 2019; Li et al., 2023) and we need to apply the notion of inference used in Jin and Rubin (2008). With the Bayesian causal forest approach described in the next section, we adopt the following conditional independence assumptions to specify two joint conditional distributions of $(Y_{i}(1),Y_{i}(0))$ and $(M_{i}(1),M_{i}(0))$: $Y_{i}(1-A_{i})\perp Y_{i}(A_{i})|A_{i},M_{i}(0),M_{i}(1),\bm{X}_{i}$ and $M_{i}(1-A_{i})\perp M_{i}(A_{i})|A_{i},\bm{X}_{i}$. Then, this assumption also further factorizes the joint conditional distributions of potential (intermediate) outcomes as follows: $Pr(Y_{i}(0),Y_{i}(1)|M_{i}(0),M_{i}(1),\bm{X}_{i},\bm{\phi})=Pr(Y_{i}(0)|M_{i}(0),M_{i}(1),\bm{X}_{i},\bm{\phi})Pr(Y_{i}(1)|M_{i}(0),M_{i}(1),\bm{X}_{i},\bm{\phi})$ and $Pr(M_{i}(0),M_{i}(1)|\bm{X}_{i},\bm{\phi})=Pr(M_{i}(0)|\bm{X}_{i},\bm{\phi})Pr(M_{i}(1)|\bm{X}_{i},\bm{\phi}),$ where $A_{i}$ is omitted in both sides of the equations under the unconfoundedness assumption. This assumption is conservative in that it is stronger than assuming that $(M_{i}(0),M_{i}(1))$, and/or $(Y_{i}(0),Y_{i}(1))$ are conditionally associated.

Among many causal inference methods based on BART, Hahn et al. (2020) proposed the Bayesian causal forest model (BCF) to accurately measure heterogeneous effects with the goal of minimizing confounding. Ignoring for the moment the presence of the intermediate variable, $M$, the BCF method specifies the response surface as

where independent BART priors are specified on the functions $\mu$ and $\tau$. Here, $\hat{\pi}(\bm{x}_{i})$ is an estimated propensity score, which corresponds to the exposure assignment mechanism under Assumption 1.

A key point of the BCF specification is the role of the $\mu$ function, also known as the prognostic function, which is equivalent to the expected outcomes of the unexposed (i.e., $\mu(\bm{x})=\mathbb{E}(Y(0)|\bm{x})$). If the predicted value of the outcome in the absence of exposure relates to the assignment of exposure, the prognostic function will exhibit some dependence on the propensity score $\hat{\pi}(\bm{x}_{i})$. This phenomenon is referred to as targeted selection (Hahn et al., 2020). One issue with targeted selection is that it can lead standard regularization methods - including the regularization implied by BART - to present regularization-induced confounding (RIC), particularly when conditional outcome distributions depend more on covariates in $\bm{X}$ than they do on the treatment, $A$. The inclusion of the estiamted propensity score within the function $\mu$ is designed specifically to mitigate the threat of RIC, and is a key feature of the BCF model in (1). In the power plant setting, we expect targeted selection to arise because scrubbers are more likely to be installed in power plants emitting a significant amount of pollutants (SO₂) and in areas in danger of violating ambient air quality standards, that is, in areas with the higher values of PM_2.5  in the absence of emissions controls.

The other main component of BCF is its inclusion of the function $\tau(\bm{x}_{i})$, which is the feature of the model specifically targeted to a characterization of treatment effect heterogeneity. Specifically, a flexibly-modeled $\tau(\bm{x}_{i})$ would represent a complex function of covariates that modify the relationship between $A_{i}$ and the expected outcome. Estimation in the BCF approach is achieved by specifying an independent BART prior for the $\mu(\bm{x}_{i})$ and $\tau(\bm{x}_{i})$ functions. We provide further details within the context of the proposed model, which incorporates intermediate variables, in a later section.

The observed data models of $Y$ and $M$ can be expressed using the following form of BCF:

where $\epsilon_{i}\sim N(0,\sigma_{y}^{2})$ and $\epsilon^{\prime}_{i}\sim N(0,\sigma_{m}^{2})$. This essentially corresponds to specifying two BCF models; one for the intermediate, $M$, and another for the outcome, $Y$, conditional on both potential intermediate values. These two BCF models correspond to the two joint distributions described as consequences of the conditional independence assumptions described in Section 2.2.: $Pr(Y_{i}(0),Y_{i}(1)|M_{i}(0),M_{i}(1),\bm{X}_{i},\bm{\phi})$, and $Pr(M_{i}(0),M_{i}(1)|\bm{X}_{i},\bm{\phi})$. Here, independent BART priors are specified for the functions $\mu_{M}$, $\mu_{Y}$, $\tau_{M}$ and $\tau_{Y}$, and $\hat{\pi}(\bm{x_{i}})$ represents the estimated propensity score $\hat{\pi}(\bm{x}_{i},\bm{\phi}_{A})=Pr(A_{i}=1\,|\,\bm{x}_{i},\bm{\phi}_{A})$ for each $i$-th subject where $\bm{\phi}_{A}$ is a subvector of the global parameter $\bm{\phi}$. Key points to highlight include: (1) the flexibility of the outcome and intermediate models in mitigating the threat of RIC while effectively capturing the causal effects of $A$ on both $Y$ and $M$; (2) the reliance of the outcome model on potential intermediate variables $M_{i}(1)$ and $M_{i}(0)$ rather than an observed intermediate variable $M_{i}^{\text{obs}}$; and (3) the importance of the function $\tau_{Y}$ in the outcome model, which allows the treatment effect of $A$ on $Y$ to be modified by a complex function of principal stratum membership $S_{i}=(M_{i}(0),M_{i}(1))$ and covariates. The simulation study compares against an approach that excludes $\bm{X}_{i}$ from the specification of $\tau_{Y}$. In either case, estimated PCEs will marginalize over values of $\bm{X}_{i}$ within principal strata defined by $S_{i}$.

By incorporating BART priors, the entire set of observed data models can be reformulated:

where, for each $f\in\{\mu_{Y},\tau_{Y},\mu_{M},\tau_{M}\}$, each of $H_{f}$ distinct tree structures is denoted by $\mathcal{T}_{h}^{f}(h=1,\cdots,H_{f})$ and the parameters at the terminal nodes of the $h$-th tree are denoted by $\mathcal{M}_{h}^{f}=\{\mu_{h,1}^{f},\cdots,\mu_{h,n_{h}}^{f}\}$ where $n_{h}$ is the number of terminal nodes of $\mathcal{T}_{h}^{f}$. And $\mathcal{F}_{f}(\bm{x},\mathcal{T}_{h}^{f},\mathcal{M}_{h}^{f})$ represents $\mu_{t,\eta}^{f}\in\mathcal{M}_{h}^{f}$ if $\bm{x}$ is associated wtih the $\eta$-th terminal node in tree $\mathcal{T}_{h}^{f}$.

As suggested in Hahn et al. (2020), we specify different BART priors on $\mu$ and $\tau$ functions. In general, the default parameter setting in Chipman et al. (2010) is used: (1) a node at depth $d$ splits with the probability $\alpha(1+d)^{-\beta}$ ($\alpha=0.95$ and $\beta=2$); (2) for $\tau_{Y}$ in the outcome model, independent priors $\mu_{h,j}^{\tau_{Y}}\sim N(0,\sigma_{0}^{2})$ are assigned to the leaf parameters where $\sigma_{0}=\sigma_{y}/\sqrt{H_{\tau_{Y}}}$; and (3) for $\tau_{M}$ in the intermediate model, independent priors $\mu_{h,j}^{\tau_{M}}\sim N(0,\sigma_{0}^{2})$ are assigned to the leaf parameters where $\sigma_{0}=\sigma_{m}/\sqrt{H_{\tau_{M}}}$. For $\mu_{Y}$ and $\mu_{M}$, however, we specify a weakly informative half-Cauchy prior on the scale parameter of the leaf parameters. We set the prior third quartile to twice the standard deviation of the observed $Y$ (or $M$).

The detailed Markov chain Monte Carlo (MCMC) approach for posterior computation process is illustrated in the pseudo-code (Algorithm 1) in the Appendix. As shown in Algorithm 1 in the Appendix, a total of four BART priors (two for $M$ model and two for $Y$ model) were updated using “Bayesian backfitting” (Hastie and Tibshirani, 2000) as a Metropolis-within-Gibbs sampler in which each tree is fit sequentially through the unexplained responses. The number of trees included in each BART prior is set differently, as suggested by Hahn et al. (2020). The BART priors assigned to the prognostic functions $\mu_{Y}$ and $\mu_{M}$ that adjust for confounders consist of a large number of trees (e.g., $H_{\mu_{Y}}=H_{\mu_{M}}=200$), whereas the BART priors assigned to $\tau_{Y}$ and $\tau_{M}$ representing the heterogeneous effects consist of a small number of trees that regularizes towards the homogeneous effects (e.g., $H_{\tau_{y}}=H_{\tau_{m}}=50$).

To summarize the entire procedure, for the $r$-th MCMC iteration, the BART priors assigned to the $\mu_{M}$ function (prognostic function) and the $\tau_{M}$ function in the $M$ model are sequentially updated. Thereafter, for each $i=1,...,n$, $M_{i}^{\text{mis}}$ is generated, which is an unobservable $M_{i}(1-A_{i})$ value. Through this, the value of the principal stratum membership of the $r$-th MCMC iteration $S_{i}^{(r)}=(M_{i}(1),M_{i}(0))$ is obtained for each $i$-th observation. Similarly, the BART priors corresponding to the $Y$ model are updated. The two potential intermediates, $M_{i}(A_{i})$ and $M_{i}(1-A_{i})$ (forming the principal stratum membership variable, $S_{i}^{(r)}$), along with $\bm{X}_{i}$, serve as independent variables for splitting the tree during the growth of the $\tau_{Y}$ function. After updating all of the BART parameters, the $r$-th iteration is completed by updating the variance parameters $(\sigma_{m}^{2},\sigma_{y}^{2})$ with the Gibbs sampler. It is worth noting that the sampling for potential intermediate variables incorporated during the MCMC takes into account the densities for both the $Y$ and $M$ models. Thus, when sampling $M_{i}^{\text{mis}}$, it is done conditionally based on the structure of the $Y$ model, which includes the pair of potential intermediate variables. The R code for posterior computation is available at https://github.com/lit777/BPCF.

With $R$ of posterior samples for each parameter and unobserved potential outcome, the posterior means of the principal causal effects for a specific principal stratum $\{S_{i}:\epsilon_{\text{lower}}<M_{i}(1)-M_{i}(0)<\epsilon_{\text{upper}}\}$ can be estimated as follows:

where $(M_{i}^{(r)}(0),M_{i}^{(r)}(1))$ and $(Y_{i}^{(r)}(0),Y_{i}^{(r)}(1))$ are the $r$-th MCMC samples of potential intermediate and outcome variables, respectively.

Multiple values of $\bm{\epsilon}$ can be specified to investigate how the effect of $A$ on $Y$ varies as a function of $A$ on $M$ by analyzing the outcome effects within groups categorized based on the strength of the exposure-intermediate relationship. For example, a sequence of $\Delta_{\text{PCE}}(\bm{\epsilon})$ estimate to increase in magnitude as $|\epsilon_{\text{upper}}-\epsilon_{\text{lower}}|$ grows would indicate that larger effects on the intermediate are associated with larger impacts on the outcome. We will motivate and explore a series of $\Delta_{\text{PCE}}(\bm{\epsilon})$ in the context of the analysis of power plant data.

For the first simulation scenario, we consider the following data generating process:

where $h_{1}(x)=(-1)^{I(x\geq 0)}$ and $h_{2}(x)=(-1)^{I(x<0)}$. This simulation scenario is developed to evaluate the accuracy of PCE estimation when the $Y$ and $M$ models are given a complex structure.

In this simulation study, target PCEs were defined based on five different intervals. These intervals were determined using quintiles of observed intermediate values obtained through the data generating process. Specifically, the first interval was (2.00, 2.26), the second interval was (2.26, 2.53), the third interval was (2.53, 2.84), and the fourth interval was (2.84, 3.29), and the last interval was (3.29, 5.09). The results are shown in Table 1. Relative bias and MSE were compared, and the proposed BPCF performed relatively better across all PCEs. In particular, it demonstrated very high accuracy in estimating the causal effect for intermediate variables $[M_{i}(1)-M_{i}(0)]$, and it can be argued that it showed high performance in estimating PCEs by estimating accurate principal stratum membership based on precise $[M_{i}(1)-M_{i}(0)]$ estimates.

A causal effect predictiveness surface (Gilbert and Hudgens, 2008) was created to help interpret the results of each model more precisely. Figure 1 depicts the predictiveness surfaces created using four different methods, including true surface plots. The $x-y$ plane has two axes representing $M(1)$ and $M(0)$ variables, and the $z$-axis represents the corresponding $Y(1)-Y(0)$ value. In terms of the overall response slope, the DP method differed significantly from the actual shape. The dots on the $x-y$ plane represent $M(1)$ and $M(0)$ sample values from one MCMC iteration. Based on this, it is reasonable to interpret that the BART and BPCF methods form more accurate principal strata than does DP. In particular, it can be observed that the DP method failed to accurately estimate $S_{i}=(M_{i}(1),M_{i}(0))$. The reason for this phenomenon occurring is that the data generating model for $M$ has a complex nonlinear structure that includes the interaction of $A$ and $\bm{X}$. In contrast, the DP method (based on Schwartz et al. (2011)) is designed to vary only the intercept parameters of $M(1)$ and $M(0)$. Extensions of Schwartz et al. (2011) to extend the DP model to additional model parameters might improve performance in a scenario where interactions between $A$ and $\bm{X}$ dominate the data generating model for $M$. On the other hand, the BPCF${}_{\text{m only}}$ method relatively accurately estimated the values of $S_{i}$, but the response surface showed deviations from the true values. This can be interpreted as the tree structure shrinking towards a more homogeneous effect when the modifier function only includes $S_{i}$ and not all covariates.

We collected data on 385 coal fuel power plants located in the eastern United States and operating during 2014, incorporating information from the EPA Air Markets Program Data and the Energy Information Administration. The data encompassed plant characteristics, installed emissions control technologies (if applicable), and emissions of various pollutants including SO₂. Additionally, predictions of annual PM_2.5  concentrations at the zip code level were obtained (Wei et al., 2022). Our analysis specifically centered around the year 2014. In terms of confounders, weather variables are major factors that can affect the operation of power plants and the concentration of ambient PM_2.5. We characterized meteorological conditions at the zip code level by measuring relative humidity, temperature, and precipitation at the zip code’s centroid (Tec et al., 2022). Non-local meteorological information from surrounding areas is also a crucial factor that affects local air quality. Recent work in Tec et al. (2022), proposed the use of self-supervised learned latent representations of non-local meteorological information surrounding each zip code as important confounding variables in studies of air pollution. We use the learned representations of non-local meteorological variables within a 100km radius of each zip code as additional confounders which can be extracted from https://huggingface.co/spaces/mauriciogtec/w2vec-app. Furthermore, the 2010 US Census data provided critical demographic information (i.e., total population) for a potential confounder surrounding each power plant.

Using the centroids of the zip codes, we established connections between each power plant and all zip code locations within a specified radius. This approach enabled us to link the average ambient PM_2.5  concentration and the average value of local confounders to each power plant. We used a 50km radius in the main manuscript, with an additional radius (60km) being examined in the Appendix. As of 2014, a power plant was categorized as treated if over half of its EGUs had implemented any control techniques for sulfur dioxide (SO₂). Additionally, only power plants with a total heat input greater than 0 were considered, and power plants with missing values for important confounding factors (especially power plant characteristics) were excluded. Table 2 shows the summary statistics for each variable in the two treatment groups. Before presenting the analytical results, MCMC convergence was assessed with plots of the MCMC log-marginal likelihood (Figure S2 in the Appendix).

A simple average comparison of the two groups showed that the concentration of ambient PM_2.5  was approximately 0.25 ($\mu\text{m}/\text{m}^{3}$) higher within 50 km of the power plants where SO₂  control technologies were installed. The average SO₂  pollution emissions from the power plants with and without any SO₂  control techniques were 713.95 tons and 1063.70 tons, respectively, indicating that the group of power plants without SO₂  control technologies emitted more (roughly 349.76 tons more) SO₂  pollution. Propensity score estimates were obtained using logistic regression with the confounding variables mentioned in the previous section.

Subsequent estimates are reported as posterior means along with 95% credible intervals. The estimated causal effect of SO₂  control strategies on SO₂  emissions was -480.09 (-702.55, -378.54) tons, indicating a significant reduction in SO₂  emissions. Similarly, the estimated causal effect of control strategies on ambient PM_2.5  concentrations was -0.60 (-0.83, -0.37), showing a significant reduction in PM_2.5  concentrations within a 50-kilometer radius. These results coincide with the expectation that introducing SO₂  reduction technology in power plants would consistently result in reduced SO₂  emissions and subsequently lower ambient PM_2.5  levels. The principal stratification approach will more thoroughly examine the influence of power plants’ varying response in SO₂  emissions on the downstream impacts on ambient PM_2.5 .

Notably, Figure 2 indicates that, in both the intermediate and outcome models, there is a positive relationship between propensity score estimates and prognostic function estimates. Particularly for the outcome variable (PM_2.5), an increasing relationship between $\mu$ and $\pi$ indicates that targeted selection bias may exist. This serves as one key motivation for the BCF-based approach.

Principal strata chosen for summarizing the CEP surface were constructed based on the estimated variability among potential intermediate variables. The posterior predictive mean values for each $i$-th $M_{i}(1)-M_{i}(0)$ was estimated, and multiplying the standard deviation of these values by 0.2 and 0.5 resulted in values of 320 and 798, respectively. Using this information, six principal strata were constructed as follows.: interval 1 is ($<$ -798), interval 2 is [-798, -320), interval 3 is [-320, 0), interval 4 is [0, 320), interval 5 is [320, 798), and interval 6 is ($\geq$ 798). The strata defined by the first three intervals demonstrate how the PM_2.5  concentrations vary within units where SO₂  emissions decrease due to SO₂  control installation. The last three intervals explain the opposite impact of SO₂  control techniques on PM_2.5  concentrations in situations where SO₂  emissions actually increase.

In the case of the principal causal effect for interval 1 ($\Delta_{\text{PCE}1}$), 52.3 power plants were assigned to the subgroup on average, and the estimated value of the causal effect was -0.59 (95% CI: -0.84, -0.31). In the case of the principal causal effect for interval 2 ($\Delta_{\text{PCE}2}$), 34.9 power plants were assigned to the subgroup on average, and the estimated value of the causal effect was -0.69 (95% CI: -1.05, -0.32). For interval 3($\Delta_{\text{PCE}3}$), 34.2 power plants were assigned to the subgroup on average, and the estimated value of the causal effect was -0.67 (95% CI: -1.13, -0.22). On the other hand, in the case of the principal causal effects for intervals 4-6 ($\Delta_{\text{PCE}4}$,$\Delta_{\text{PCE}5}$,$\Delta_{\text{PCE}6}$), the groups with increased SO₂  emissions had averages of 19.5, 8.8, 11.4 power plants and estimated effects of -0.61 (95% CI: -1.10, -0.11), -0.24 (95% CI: -1.47, 0.51), -0.45 (95% CI: -1.36, 0.57), respectively. The first three PCEs displayed significant effect estimates, with more power plants assigned to those strata. This suggests that the implementation of SO₂  control techniques led to a substantial reduction in SO₂  emissions and consequently lowered PM_2.5  concentrations. Specifically, based on the fact that the highest number of power plants was assigned to the first principal stratum, it can be inferred that a significant number of power plants achieved a reduction in SO₂  emissions of 798 tons or more through the installation of SO₂  control techniques. In contrast, the last three principal strata, indicating a causal increase in PM_2.5  concentration due to the installation of SO₂  control, did not show significant PCEs except for interval 4. Moreover, the number of power plants within these strata was observed to be relatively small. This suggests that the installation of SO₂  controls tends to not increase SO₂  emissions and, furthermore, exhibits little evidence of an impact on PM_2.5  concentrations.

These results are well illustrated in Figure 3 and Figure 4. Figure 3 displays the principal causal effects corresponding to intervals 1-6 from left to right. The center point of each circle represents the posterior mean of the principal causal effect, and the vertical line represents the 95% credible interval. The size of each circle represents the average value of subjects within the respective stratum, and the numbers below indicate the posterior mean, posterior standard deviation (SD), and average number of subjects in sequential order. The $x$-axis represents the causal effect on SO₂  emissions within each interval.

Figure 4 presents the average response surface. The points on the $x-y$ axis represent the posterior predictive mean values for $S_{i}=(M_{i}(1),M_{i}(0))$. It can be observed that the slope of the response surface generally decreases in regions where $M(1)$ values are lower than $M(0)$ values. This indicates that the reduction in SO₂  emissions corresponds to a decrease in PM_2.5  concentrations, with steeper impacts in areas of the largest SO₂  reductions.