In this vignette, we show how to use BCF to estimate treatment effects of a simulated intervention.
library(bcf)
library(ggplot2)
library(latex2exp)
library(rpart)
library(rpart.plot)
library(partykit)First, we simulate some data for testing. This data set has an outcome variable \(y\), a treatment indicator \(z\), and three covariates \(x_1\), \(x_2\) and \(x_3\). Of the covariates, two affect the outcome \(y\) at both levels of treatment, while the third is an effect moderator.
We draw three random \(x\)s for each unit and generate each unit’s expected outcome without treatment, \(\mu\), as a function of \(x_1\) and \(x_2\). Each unit’s probability of joining the intervention, \(\pi\), is also a function of \(\mu\), so that units with larger responses under control are more likely to participate in the intervention. We then assign units to treatment (\(z = 1\)) or comparison (\(z = 0\)) status as a function of \(\pi\).
Then we generate the true treatment effect for each unit, \(\tau\). As noted above, \(\tau\) is a function of \(x_3\). The observed outcome, \(y\), is a function of \(\mu\), \(\tau\), a random error term with variance \(\sigma^{2}\), and weights \(w\) if applicable.
set.seed(1)
p <- 3 # two control variables and one effect moderator
n <- 1000 # number of observations
n_burn <- 2000
n_sim <- 1000
x <- matrix(rnorm(n*(p-1)), nrow=n)
x <- cbind(x, x[,2] + rnorm(n))
weights <- abs(rnorm(n))
# create targeted selection, whereby a practice's likelihood of joining the intervention (pi)
# is related to their expected outcome (mu)
mu <- -1*(x[,1]>(x[,2])) + 1*(x[,1]<(x[,2])) - 0.1
# generate treatment variable
pi <- pnorm(mu)
z <- rbinom(n,1,pi)
# tau is the true treatment effect. It varies across practices as a function of
# X3 and X2
tau <- 1/(1 + exp(-x[,3])) + x[,2]/10
# generate the expected response using mu, tau and z
y_noiseless <- mu + tau*z
# set the noise level relative to the expected mean function of Y
sigma <- diff(range(mu + tau*pi))/8
# draw the response variable with additive error
y <- y_noiseless + sigma*rnorm(n)/sqrt(weights)In this data set we have observed \(y\), \(x\), and \(\pi\) values to which we can fit our BCF model. With BCF, we can distinguish between control variables – which affect the outcome at both levels of treatment – and moderator variables – which affect the estimated treatment effect.
Note that we are using the n_chains argument to
bcf(), which allows us to run several MCMC chains in
parallel and assess whether they have converged to the posterior
distribution.
After fitting the BCF model, we can compare the \(\hat{\tau}\) estimates from BCF to the true \(\tau\) from the data-generating process.
We use the summary.bcf function to obtain posterior
summary statistics and MCMC diagnostics. We use those diagnostics to
assess convergence of our run.
summary(bcf_out)
#> Summary statistics for each Markov Chain Monte Carlo run
#>
#> Iterations = 1:1000
#> Thinning interval = 1
#> Number of chains = 4
#> Sample size per chain = 1000
#>
#> 1. Empirical mean and standard deviation for each variable,
#> plus standard error of the mean:
#>
#> Mean SD Naive SE Time-series SE
#> sigma 0.37556 0.008577 0.0001356 0.0001689
#> tau_bar 0.49308 0.040953 0.0006475 0.0012522
#> mu_bar -0.08981 0.023568 0.0003726 0.0005772
#> yhat_bar 0.19119 0.012922 0.0002043 0.0002069
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> sigma 0.3592 0.3698 0.37533 0.38120 0.39302
#> tau_bar 0.4131 0.4656 0.49293 0.52112 0.57317
#> mu_bar -0.1362 -0.1058 -0.08954 -0.07331 -0.04391
#> yhat_bar 0.1657 0.1825 0.19128 0.20000 0.21660
#>
#>
#> ----
#> Effective sample size for summary parameters
#> Reverting to coda's default ESS calculation. See ?summary.bcf for details.
#>
#> sigma tau_bar mu_bar yhat_bar
#> 2676.224 1077.885 1696.041 3911.178
#>
#> ----
#> Gelman and Rubin's convergence diagnostic for summary parameters
#> Potential scale reduction factors:
#>
#> Point est. Upper C.I.
#> sigma 1 1.00
#> tau_bar 1 1.01
#> mu_bar 1 1.01
#> yhat_bar 1 1.00
#>
#> Multivariate psrf
#>
#> 1
#>
#> ----Since our “\(\hat{R}s\)” (Gelman and Rubin’s convergence factor or Potential Scale Reduction Factor) are between 0.9 and 1.1, we don’t see any obvious mixing issues. Note that mixing on other parameters (like subgroup average treatment effects defined below, or evaluations of \(\tau(x)\) or \(\mu(x)\) at particular values of \(x\)) might be worse than these aggregate summaries. But in our experience good mixing on \(\sigma\) tends to indicate good mixing on other quantities of interest.
Now that we’ve successfully fit our model, let’s explore the output. First, since this is a simulation, let’s compare the unit-specific treatment effect estimates \(\hat{\tau}\) to the true unit-specific treatment effects \(\tau\).
We plot the true and estimated \(\tau\) versus \(x_3\), since \(x_3\) is one of our primary effect moderators (see the data creation step above).
tau_ests <- data.frame(Mean = colMeans(bcf_out$tau),
Low95 = apply(bcf_out$tau, 2, function(x) quantile(x, 0.025)),
Up95 = apply(bcf_out$tau, 2, function(x) quantile(x, 0.975)))
ggplot(NULL, aes(x = x[,3])) +
geom_pointrange(aes(y = tau_ests$Mean, ymin = tau_ests$Low95, ymax = tau_ests$Up95), color = "forestgreen", alpha = 0.5) +
geom_smooth(aes(y = tau), se = FALSE) +
xlab(TeX("$x_3$")) +
ylab(TeX("$\\hat{\\tau}$")) +
xlim(-4, 6) +
geom_segment(aes(x = 3, xend = 4, y = 0.2, yend = 0.2), color = "blue", alpha = 0.9) +
geom_text(aes(x = 4.5, y = 0.2, label = "Truth"), color = "black") +
geom_segment(aes(x = 3, xend = 4, y = 0.1, yend = 0.1), color = "forestgreen", alpha = 0.7) +
geom_text(aes(x = 5.2, y = 0.1, label = "Estimate (95% CI)"), color = "black")
#> `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
#> Warning: Removed 6 rows containing non-finite values (`stat_smooth()`).
#> Warning: Removed 6 rows containing missing values (`geom_pointrange()`).