Aggregating Average Treatment Effects with baggr

Data inputs: reported effects or full individual-level data sets

For average effects aggregation, baggr allows 3 types of data inputs. The user may supply, within a data frame environment, any of the following:

1. A set of estimated treatment effects \(\{\hat{\tau_k}\}_{k=1}^{K}\) and their standard errors \(\{\hat{se_k}\}_{k=1}^{K}\) from each study. This should be formatted as two column vectors of length \(K\) within the data frame, where \(\hat{\tau_k}\) is the \(k\)-th entry of the treatment effect vector and \(\hat{se_k}\) is the \(k\)-th entry of the standard errors vector. Columns should be named "tau" and "se". Model will be "rubin" (see below).

2. A set of control group means and estimated treatment effects \(\{\hat{\mu_k},\hat{\tau_k}\}_{k=1}^{K}\), as well as the standard errors for both \(\{\hat{se}_{\mu k}, \hat{se}_{\tau k}\}_{k=1}^{K}\), for each study site This should be formatted as four vectors of length \(K\) within the data frame, analogous to the above. Columns should be names "mu", "tau", "se.mu", "se.tau". Model will be "mutau" (see below).

3. The full data sets from all the original studies \(\{Y_k, T_k\}_{k=1}^{K}\). This should be formatted as three vectors of length \(\sum_{k=1}^K N_{k}\), which we recommend naming "outcome", "treatment", "group" (for site indicators), but the names can also be specified when calling baggr() function. Model will be "rubin_full" (see below).

As an example of an individual-level data set we include in data frames microcredit and microcredit_simplified. The former contains all the microcredit outcome data used in Meager (2019), standardized to USD PPP in 2009 terms per two weeks (a time period is necessary as these are flow variables). It therefore contains NAs and other undesirable features, to allow the user to see how baggr handles these common data issues. The data set microcredit_simplified has these issues cleaned up and contains only one outcome of interest, consumer durables spending.

baggr also has a function that automatically produces summary data from full data sets, in case one wishes to run the comparatively faster summary-data models. The prepare_ma() function applied to a dataframe with columns named "group", "outcome", and "treatment" automatically estimates the control group means, treatment effects, and associated standard errors for each group using an Ordinary Least Squares regression. The resulting output is already formatted as a valid input to the baggr() command itself:

prepare_ma(microcredit_simplified, outcome = "consumption")
#>   group       mu       tau     se.mu    se.tau n.mu n.tau
#> 1     1 303.6065  5.510755  2.559897  4.101140 8298  8262
#> 2     2 280.0887 50.449393 11.141075 22.156317  260   701
#> 3     3 196.4215 -5.171338 14.432604 19.266339  444   551
#> 4     4 276.2791  4.641604  3.730907  5.451094 3248  3579
#> 5     5 327.5246 -2.935731  4.027768  6.022955 2771  2716

ATE aggregation models in baggr

baggr currently contains two different models suitable for aggregating sets of average treatment effects. Consider first the evidence aggregation model from Rubin (1981), discussed extensively in Chapter 5 of Gelman et al. (2013); the model consists of a hierarchical likelihood as follows:

\[\begin{equation} \begin{aligned} \hat{\tau_k} &\sim N(\tau_k, \hat{se_k}^2) \; \forall \; k \\ \tau_k &\sim N(\tau, \sigma_{\tau}^2) \; \forall \; k . \end{aligned} \label{rubin model} \end{equation}\]

The motivation for this model structure is discussed in detail in the sources above and in Meager (2019). To complete the Bayesian model, we now need priors. baggr has a set of default priors for each model (adjusted to data), as well as allowing the user to specify her own priors if desired. In the Rubin model, baggr’s default priors on the hyper-parameters are as follows:

• For \(\tau\), the prior is \(\mathcal{N}(0, 100)\). This is a very weak prior which does little regularization as a default, centered at zero. It assumes that causal effects should not be thought of as large unless data contains evidence to the contrary (a “hand-wavey” form of Occam’s Razor).
• For \(\sigma_{\tau}\) the prior is \(\textrm{Uniform}(0,10\tilde{\sigma})\), where \(\tilde{\sigma}\) is a naive standard deviance estimator for \(\{\hat{\tau_k}\}_{k=1}^{K}\) (generated by the R command sd()).

In case you also have data on the control groups’ mean outcomes and the uncertainty on those,you can augment the Rubin (1981) model to incorporate that information. Following Meager (2019), if one has access to the estimated control means \(\{\hat{\mu_k}\}^K_{k=1}\) and their standard errors \(\{\hat{se}_{\mu k}\}^K_{k=1}\), one can fit a joint Gaussian model on the pairs \(\{\hat{\mu_k},\hat{\tau_k}\}_{k=1}^{K}\):

\[\begin{equation} \begin{aligned} \hat{\tau_k} &\sim N(\tau_k, \hat{se_{\tau k}}^2) \; \forall \; k \\ \hat{\mu_k} &\sim N(\mu_k, \hat{se_{\mu k}}^2) \; \forall \; k \\ \left( \begin{array}{c} \mu_{k}\\ \tau_{k} \end{array} \right) &\sim N\left( \left( \begin{array}{c} \mu\\ \tau \end{array} \right), V \right) \; \text{where} \;V = \left[ \begin{array}{cc} \sigma^2_{\mu} & \sigma_{\tau\mu} \\ \sigma_{\tau\mu} & \sigma_{\tau}^2 \end{array} \right]\forall \; k. \\ \end{aligned} \label{full data model} \end{equation}\]

In baggr this model is referred to as "mutau".

If you have only few groups, the priors on \(V\) will need to be relatively strong to avoid overfitting. See Meager (2019) for more discussion of this issue in particular, or see the Stan Manual on hierarchical priors. The default priors are as follows:

• \(V\) is decomposed into a correlation matrix \(\Omega\), which receives an LKJCorr(3) prior regularizing it towards independence, and a scale parameter \(\theta\) which receives a \(\textrm{Cauchy}(0,s)\) prior, with \(s\) set to either 10 times the empirical SD of \(\mu\) of of \(\tau\), whichever is higher.
• The default prior for \((\mu, \tau)\) is bivariate Gaussian, with zero means and a diagonal covariance matrix with the diagonal (variances) set to 100 times the empirical maximum of \(\mu\) and \(\tau\).

Models, their inputs, likelihood and priors: a summary table

Where inputs are individual-level data, i.e. outcome, treatment, group, you can specify column names in as arguments to baggr() function.

Prior choice in baggr

In the “Rubin” and “\(\mu\) and \(\tau\)” models, the user can specify custom priors beyond the defaults using the prior arguments. These prior arguments are subdivided into 3 categories:

• Priors for the hypermean, the average effect across groups, such as \(\tau\) or the vector \((\mu, \tau)\). This is denoted using the prior_hypermean argument in baggr.

• Priors for the hyper-standard-deviations \(\sigma_{\tau}\), the standard deviation of effects across groups. This also refers to the prior on \(\theta\), the scale parameter of the hyper-variance-covariance matrix \(V\) in the”\(\mu\) and \(\tau\)” model. This type of parameter is denoted prior_hypersd in baggr.

• Priors on hypercorrelations between parameters that may co-vary across the groups. When working with multi-dimensional parameters at the upper level of the model (aka performing multivariate shrinkage), this is the correlation in the distribution of the different parameters across the groups. For example, in the “\(\mu\) and \(\tau\)” model this parameter is the correlation matrix \(\Omega\) that governs the hyper-variance-covariance matrix \(V\). This type of parameter is denoted prior_hypercor in baggr.

The possible prior distributions we allow for in the current version are:

• For prior_hypermean we allow "normal", "uniform", "cauchy", "lognormal", and (generalised) "student_t" with any parameter values which are logically possible given support constraints (e.g. user cannot specify a negative variance on a normal distribution or a negative scale on a Cauchy) When the model has a vector hypermean, “\(\mu\) and \(\tau\)” model, baggr applies the given prior to all the elements of the vector independently. If the user wishes to specify a prior dependence between the components, one can supply prior_hypermean with a multinormal argument (see ?priors for details).
• For prior_hypersd we allow "normal" and "uniform" with any parameter values which are logically possible given support constraints (as above).
• For prior_hypercor, which is only necessary when there is a multivariate shrinkage operation (as in the “\(\mu\) and \(\tau\)” model) we allow an LKJ correlation prior ("lkj") with any parameter values which are logically possible given support constraints. For further details of this strategy see Meager (2019).

Notation for priors is “plain-text”, in that you can write the distributions as normal(), uniform() etc. See ?priors for details, or continue to the example below with the Rubin model.

Running the Rubin Model in baggr

To demonstrate the Rubin model in baggr, consider the 8 schools example from Rubin (1981). In this dataset, 8 schools in the United States performed similar randomized experiments to estimate the causal effect of an SAT tutoring program on learning outcomes. Reported treatment effects and their standard errors are included in baggr as a data frame:

schools
#>      group tau se
#> 1 School A  28 15
#> 2 School B   8 10
#> 3 School C  -3 16
#> 4 School D   7 11
#> 5 School E  -1  9
#> 6 School F   1 11
#> 7 School G  18 10
#> 8 School H  12 18

To fit the model in baggr (having followed the installation instructions and loaded the package):

baggr_schools <- baggr(schools, model = "rubin", pooling = "partial")

This creates a baggr class object, and you can access the underlying stanfit object by calling baggr_schools$fit. If you don’t change the default priors, then baggr will print a message informing you of the priors it has chosen.

Printing baggr_schools returns a summary of the posterior inference. First baggr records the model type and the pooling regime chosen by the user or implemented by default. Second, baggr returns inference on the aggregate treatment effect \(\tau\) by reporting its posterior mean and 95% uncertainty interval, and similar inference on the hyper-SD \(\sigma_{\tau}\). Lastly, it prints the “updated” inference on each of the groups’ treatment effects, displaying their new posterior means, standard deviations and pooling factors (see below).

print(baggr_schools)
#> Model type: Rubin model with aggregate data 
#> Pooling of effects: partial 
#> 
#> Aggregate treatment effect (on mean), 8 groups:
#> Hypermean (tau) =  8.3 with 95% interval -1.8 to 22.0 
#> Hyper-SD (sigma_tau) = 6.93 with 95% interval 0.26 to 23.09 
#> Posterior predictive effect = 8.4 with 95% interval -12.0 to 32.0 
#> Total pooling (1 - I^2) = 0.75 with 95% interval 0.20 to 1.00 
#> 
#> Group-specific treatment effects:
#>          mean  sd   2.5%  50% 97.5% pooling
#> School A 11.7 8.8  -2.47 10.3    33    0.81
#> School B  7.9 6.3  -5.10  8.0    20    0.71
#> School C  6.2 8.0 -11.66  6.8    21    0.83
#> School D  7.6 6.6  -6.07  7.6    21    0.74
#> School E  4.9 6.4  -9.28  5.5    16    0.68
#> School F  5.9 6.5  -8.18  6.3    18    0.74
#> School G 10.9 7.0  -0.63 10.0    27    0.71
#> School H  8.6 7.8  -6.52  8.2    26    0.85

baggr(schools, "rubin", 
      prior_hypermean = normal(-5, 10),
      prior_hypersd   = uniform(0, 5))

custom_priors <- list( hypermean = cauchy(0,25), hypersd = normal(0,30))
baggr(schools, "rubin", pooling = "partial", prior = custom_priors)

Understanding and criticising baggr model performance

Baggr models are run in Stan, and the fit and results need to be checked, understood and criticised as you would any stan model or indeed any MCMC model fitting exercise. You must pay attention to printed warnings about the Rhat criterion: if you see a warning that Rhat statistic exceeds 1.05 for any parameter, you MUST NOT use the results for inference. This warning means the MCMC chains have not converged, and it is exceedingly unlikely that the “posterior inference” printed out corresponds to anything close to the true posterior implied by your model and data. If you use results from which the Rhat statistic exceeds 1.05 YOUR INFERENCE WILL BE WRONG.

If you see this warning, try re-running the model with the option iter set to a large number such as 10,000, as below. It is also good practice to run many chains, such as 8 rather than the default 4, to have a greater chance to detect pathological MCMC behaviour. You do this by passing baggr the stan arguments iter = 10000 and chains = 8, like so:

baggr_schools <- baggr(schools, model = "rubin", pooling = "partial", 
                       iter = 10000, chains = 8)

Other warnings you may see involve “divergent transitions”. While not as serious as high Rhat, this can signal problems with the model. As the stan message that you will see suggests, try adjusting adapt_delta argument above 0.8, e.g. to 0.99. You cannot pass this parameter directly to stan and thus you cannot pass it directly to baggr, so instead you must pass the argument control = list(adapt_delta = 0.99).

Measuring “pooling”

It is often useful to measure the extent to which the hierarchical model is “pooling” or sharing information across the groups in the process of aggregation. Baggr automatically computes and prints such a metric, as seen above. You can access more details by writing pooling(baggr_schools) or heterogeneity(baggr_schools).

pooling(baggr_schools)
#> , , 1
#> 
#>            [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
#> 2.5%  0.2967497 0.1579246 0.3243730 0.1849549 0.1318760 0.1849549 0.1579246
#> mean  0.8139490 0.7089556 0.8282654 0.7357306 0.6781800 0.7357306 0.7089556
#> 97.5% 0.9996918 0.9993067 0.9997291 0.9994270 0.9991442 0.9994270 0.9993067
#>            [,8]
#> 2.5%  0.3779680
#> mean  0.8525197
#> 97.5% 0.9997859

Plotting and model comparison in baggr

A fundamental step to understanding the model is to plot the posterior distributions. baggr has several automatic plot functions which you can access by calling baggr_plot() or, equivalently, using the default plot function; these visuals are based on bayesplot package. Plotting functions always take baggr class object as their first argument. By default, means and 95% posterior intervals of the effects in each group are shown. Extra options are available, such as whether to order the results by effect size. For the 8 schools Rubin model we have

plot(baggr_schools, order = FALSE)

Similarly important is the plot of potential treatment effects (posterior predictive distribution), obtained by repeatedly drawing hypermean and hyper-SD values from baggr’s Markov Chain model, then drawing new value of treatment effect conditional on these:

effect_plot(baggr_schools)

The values underlying the plot can be obtained by effect_draw, which can be used e.g. to draw 1 new value:

effect_draw(baggr_schools, draws = 1)
#> [1] -0.7846135

You can also summarise over the entire distribution to obtain mean, SD, and quantiles (effect_draw(baggr_schools, summary = TRUE)). For meta-regression models (models with covariates), you can also use this function to make predictions for new data, at set values of covariates.

However, from the viewpoint of model building, the most important plots are the ones that directly compare many possible models.

my_baggr_comparison <- baggr_compare(schools)

plot(my_baggr_comparison) + 
  ggtitle("8 schools: model comparison")

baggr_schools_v2 <- baggr(schools, prior_hypermean = normal(10, 2.5))

effect_plot("Default model" = baggr_schools, "normal(10, 2.5)" = baggr_schools_v2) +
  coord_cartesian(xlim = c(-10, 30)) + theme(legend.position = "top")

baggr_compare("Default model" = baggr_schools, "normal(10, 2.5)" = baggr_schools_v2)
#> 
#> Mean treatment effects:
#>                     2.5%    mean   97.5%  median      sd
#> Default model   -1.77837 8.29902 22.0418 7.91293 5.93786
#> normal(10, 2.5)  5.18573 9.51974 13.9187 9.48642 2.21891
#> 
#> SD for treatment effects:
#>                     2.5%    mean   97.5%  median      sd
#> Default model   0.263394 6.92524 23.0915 5.54017 5.93912
#> normal(10, 2.5) 0.219007 5.84283 17.6026 4.71528 4.67683
#> 
#> Posterior predictive effects:
#>                      2.5%    mean   97.5%  median       sd
#> Default model   -13.28690 8.28680 31.9200 7.76553 11.24630
#> normal(10, 2.5)  -7.91404 9.45364 27.1755 9.39610  8.08842

forest_plot(baggr_schools)

forest_plot(baggr_schools, show = "both")

Cross-validation in baggr

Baggr has built-in, automated leave-one-out cross validation for its models. The values returned by loocv() can be used to understand how any one group affects the overall result, as well as how well the model predicts the omitted group.

This function automatically runs \(K\) baggr models, leaving out one group at a time, and then calculates expected log predictive density (ELPD) for that group (see Gelman, Hwang, and Vehtari (2014)). The main output is the cross-validation information criterion, or -2 times the ELPD averaged over \(K\) models. This is related to, and often approximated by, the Watanabe-Akaike Information Criterion. A value closer to zero (i.e. a smaller number in magnitude) means a better fit. For more information on cross-validation see this overview article.

Therefore, if you have \(K\) groups, using the loocv() will run your model of choice \(K\) times. Be aware that this may take a while even for simple models. (You should see a progress bar in your terminal window.) The loocv() function takes in the same arguments as baggr(), plus an option (return_models) for whether to return all the models or just the summary statistics. For the 8 schools example we can do

loocv_res <- loocv(schools, return_models = FALSE, 
                   iter = 1000, #just to make it a bit faster -- don't try it at home!
                   model = "rubin", pooling = "partial")

loocv_res
#> LOO estimate based on 8-fold cross-validation
#> 
#>       Estimate Standard Error
#> elpd     -31.7          0.954
#> looic     63.5          1.910

The loocv() output contains more than just the matrix it prints, and this additional information can be accessed via attributes(), e.g. the mean treatment effects, their variability and elpd for each model that are stored in the attribute df:

names(attributes(loocv_res))
#> [1] "names" "class"
attr(loocv_res, "df")
#> NULL

This data frame can then be used to examine or compute the variation in the inference on \(\tau\) in the absence of each group. If the user is interested in manually checking the consequences of excluding a particular group or set of groups, this is also possible in baggr using subsetting. For example, suppose that we want to run the Rubin model on school groups 1-7 and predict the effect in the 8th school. The code below shows how you can specify a subset of the dataframe as your “data” argument, and then designate another subset as the “testing” holdout set by assigning it to the argument “test_data” in the baggr command. Here we have done it for both partial and full pooling:

We can compare the performance of the two models using the mean log predictive density. This itself is a density as the name suggests so we will compute the log expected value of this density: here, as before, a number closer to zero is better. In this case the full pooling model actually does slightly better:

fit1$mean_lpd
#> [1] 7.945999
fit2$mean_lpd
#> [1] 7.728722

If we run the full loocv for both models, we can estimate the difference between the ELPD of those models as well as the standard error of that difference. As an example, let’s consider the two models above.

loocv_full <- loocv(data = schools, 
                    model = "rubin", 
                    pooling = "full")

We can compare those fits with loo_compare, which gives use the differences in ELPD for the various models we pass it. Generally it is important to pay attention to both the difference in the predictive density and also the standard error of that difference.

loo_compare(loocv_res, loocv_full)
#> Comparison of cross-validation
#> 
#>                     ELPD ELPD SE
#> Model 1 - Model 2 -0.697   0.355
#> 
#> Positive ELPD indicates the reference group is preferred.