#### How to implement the design pattern using a Bayesian estimator: We have now modeled the same conceptual model on the same data using [[Statistical modeling of RD|local polynomial regression]] and [[(Double) Machine learning using RD|DoubleML with flexible covariate adjustment]]. Both gave us a treatment effect at the cutoff of approximately \$1,567. What if we model the same problem using Bayesian estimation? As we discussed in the [[IV the Bayesian way|IV Bayesian chapter]], we should expect *consistent* results, not *identical* results. The output of a frequentist analysis is a point estimate with a confidence interval; the output of a Bayesian analysis is a posterior distribution. They answer different questions about the same parameter. The posterior distribution is what makes the Bayesian approach distinctive here. With it, we can compute probability statements like "what is the probability that the true treatment effect exceeds \$1,000?" We can also visualize the distribution of the jump at the cutoff and explore heterogeneity by interacting the treatment with covariates of interest. As a reminder, the frequentist result from [[Statistical modeling of RD|Chapter 2.2.2]] using `rdrobust` is as follows: ```r RD Effect 1567.178 36.225 0.000 [1496.639 , 1667.855] ``` A \$1,567 jump in `future_spending` at the \$10,000 cutoff, with a robust 95% confidence interval of [\$1,497, \$1,668]. **Bayesian RD model** We will use the [[brms]] package in R as we did in the [[IV the Bayesian way|IV chapter]]. As discussed there, we will not replicate the Bayesian models in Python due to the computational cost; `brms` is the more mature tool for this purpose. A simple Bayesian analog of the local linear estimator is a linear regression fit to the observations within a bandwidth around the cutoff. We reuse the [[MSE-optimal bandwidth]] from the frequentist chapter (\$1,849) so the comparison is fair. The model centers the running variable at the cutoff and adds an interaction between treatment and the centered running variable, which lets the slope differ on each side of the cutoff: ```r df_local <- df %>% filter(abs(past_spending - 10000) <= 1849) %>% mutate(x_c = past_spending - 10000) model_rd <- brm( future_spending ~ status_platinum + x_c + status_platinum:x_c, data = df_local, family = gaussian(), iter = 4000, cores = 4, seed = 2026 ) summary(model_rd) ``` The coefficient on `status_platinum` is the [[Local Average Treatment Effect|treatment effect at the cutoff]]: it is the jump in fitted values at `x_c = 0`. We rely on the default priors in `brms` here; the follow-up note on [[Oh my! Priors in the Bayesian RD model|priors in the Bayesian RD model]] explores how explicit prior choices can affect the result. > [!NOTE] > A faithful Bayesian analog of `rdrobust` would use a triangular kernel that down-weights observations farther from the cutoff. We use uniform weights within the bandwidth here for simplicity and because the goal of this chapter is comparison rather than precision. Adding triangular kernel weights to the model would bring the two estimators closer if stricter alignment is needed. ***Bayesian results*** ```r Family: gaussian Links: mu = identity Formula: future_spending ~ status_platinum + x_c + status_platinum:x_c Data: df_local (Number of observations: 13535) Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1; total post-warmup draws = 8000 Regression Coefficients: Estimate Est.Error l-95% CI u-95% CI Rhat Intercept 8581.54 25.25 8530.84 8630.92 1.00 status_platinumTRUE 1574.20 35.27 1505.51 1643.45 1.00 x_c 0.03 0.02 -0.01 0.08 1.00 status_platinumTRUE:x_c -0.02 0.03 -0.09 0.05 1.00 Further Distributional Parameters: Estimate Est.Error l-95% CI u-95% CI Rhat sigma 1046.44 6.45 1034.01 1058.84 1.00 ``` The mean of the posterior for `status_platinumTRUE` is \$1,574, close to the frequentist point estimate of \$1,567. The 95% credible interval [\$1,506, \$1,643] is narrower than the frequentist robust confidence interval [\$1,497, \$1,668] but overlaps it substantially. The slope coefficients on `x_c` and `status_platinumTRUE:x_c` both have credible intervals that comfortably include zero, which simply reflects that within the narrow MSE-optimal window of \$1,849 around the cutoff, the local linear trend is essentially flat. Now to the part where Bayesian estimation pays off. ***Probability statements*** Once we have posterior samples for the treatment effect, we can ask probability questions directly: ```r posterior <- as_draws_df(model_rd) %>% rename(tau = b_status_platinumTRUE) mean(posterior$tau > 1000) mean(posterior$tau > 1500) ``` ```r [1] 1 [1] 0.983 ``` The first quantity is the posterior probability that the true treatment effect exceeds \$1,000 — a business-relevant threshold. With probability essentially equal to 1, the posterior places virtually no mass below \$1,000. The second is the posterior probability that the effect exceeds the ground-truth effect of \$1,500 from the [[Synthetic data for the RDD|data generating process]]: 0.983, or 98.3%. ***Posterior density plot*** ```r ggplot(posterior, aes(x = tau)) + stat_halfeye(fill = "skyblue", alpha = .5, point_interval = mode_hdi, .width = .95) + geom_vline(xintercept = 1500, linetype = "dashed", color = "grey40") + xlab("Treatment effect at the cutoff (USD)") + ylab("Density") + theme_classic() + theme(text = element_text(size = 20)) + labs(title = "Posterior density of the treatment effect at the cutoff") ``` ![[posterior-density-rd.png]] The dashed line marks the ground-truth effect of \$1,500. The dark line on the x-axis is the highest posterior density interval at the 95% probability level, with the mode in the middle. Unlike the confidence interval (a statement about the long-run frequency of intervals containing the true value), the HDI is a statement about where the parameter most plausibly lies given the data and the model. ***Heterogeneity in the treatment effect*** A natural follow-up question: does Platinum Status affect future spending equally for all customer segments? In the [[IV the Bayesian way|IV chapter]], we explored heterogeneity across U.S. states using random effects. Here we run a simple exploratory heterogeneity model by interacting the treatment with `primary_region` (North America, Europe, Asia-Pacific) and refit: ```r df_local <- df_local %>% mutate(region = factor(primary_region, levels = c("North America", "Europe", "Asia-Pacific"), labels = c("NorthAmerica", "Europe", "AsiaPacific"))) model_het <- brm( future_spending ~ status_platinum * region + x_c + status_platinum:x_c, data = df_local, family = gaussian(), iter = 4000, cores = 4, seed = 2026 ) ``` We then derive the region-specific treatment effects from the posterior draws by combining the main treatment coefficient with the treatment-by-region interactions (the reference region is North America, so its $\tau$ is just the main effect): ```r het_draws <- as_draws_df(model_het) %>% mutate( NorthAmerica = b_status_platinumTRUE, Europe = b_status_platinumTRUE + `b_status_platinumTRUE:regionEurope`, AsiaPacific = b_status_platinumTRUE + `b_status_platinumTRUE:regionAsiaPacific` ) %>% select(NorthAmerica, Europe, AsiaPacific) %>% tidyr::pivot_longer(everything(), names_to = "region", values_to = "tau") %>% mutate(region = factor(region, levels = c("NorthAmerica", "Europe", "AsiaPacific"), labels = c("North America", "Europe", "Asia-Pacific"))) ggplot(het_draws, aes(y = region, x = tau, fill = region)) + stat_halfeye(point_interval = mode_hdi, .width = .95) + geom_vline(xintercept = 1500, linetype = "dashed", color = "grey40") + labs(y = "Region", x = "Treatment effect at the cutoff (USD)") + theme_classic() + theme(legend.position = "none", text = element_text(size = 20)) ``` ![[posterior-heterogeneity-rd.png]] The three posteriors form a visible ordering: **North America** at \$1,545 [HDI \$1,475, \$1,629], **Europe** at \$1,587 [HDI \$1,497, \$1,673], and **Asia-Pacific** at \$1,626 [HDI \$1,514, \$1,723]. Asia-Pacific's mode sits roughly \$80 above North America's, and its posterior is the widest because it has the smallest sample share (~15% of customers). Two of the three HDIs cover the ground-truth \$1,500 from the [[Synthetic data for the RDD|data generating process]]; Asia-Pacific's lower bound at \$1,514 sits just above it, suggesting a slight upward shift in the regional posterior relative to the ground truth. The coefficient for the interaction between Platinum Status and Asia-Pacific is +\$71 with a 95% credible interval of [-\$27, \$169], leaning positive but covering zero, so even though Asia-Pacific's marginal HDI excludes \$1,500, the *joint uncertainty* in the comparison with North America does not let us conclude that the two regions respond differently. The point of the heterogeneity exercise is the same as in the [[IV the Bayesian way|IV chapter]]: visualize differences in the causal effect across segments using full posterior distributions instead of competing point estimates. Whether the segments differ "significantly" in the frequentist sense matters less than whether the posterior distributions tell a coherent business story. Finally, for more on how priors influence the Bayesian RD result, see [[Oh my! Priors in the Bayesian RD model]]! > [!info]- Last updated: May 14, 2026