This is essentially an addendum to the previous post where I simulated data from multiple RCTs to explore an analytic method to pool data across different studies. In that post, I used the nlme package to conduct a meta-analysis based on individual level data of 12 studies. Here, I am presenting an alternative hierarchical modeling approach that uses the Bayesian package rstan.
Create the data set
We’ll use the exact same data generating process as described in some detail in the previous post.
library(simstudy)
library(rstan)
library(data.table)defS <- defData(varname = "a.k", formula = 3, variance = 2, id = "study")
defS <- defData(defS, varname = "d.0", formula = 3, dist = "nonrandom")
defS <- defData(defS, varname = "v.k", formula = 0, variance = 6, dist= "normal")
defS <- defData(defS, varname = "s2.k", formula = 16, variance = .2, dist = "gamma")
defS <- defData(defS, varname = "size.study", formula = ".3;.5;.2", dist = "categorical")
defS <- defData(defS, varname = "n.study",
formula = "(size.study==1) * 20 + (size.study==2) * 40 + (size.study==3) * 60",
dist = "poisson")
defI <- defDataAdd(varname = "y", formula = "a.k + x * (d.0 + v.k)", variance = "s2.k")
RNGkind(kind = "L'Ecuyer-CMRG")
set.seed(12764)
ds <- genData(12, defS)
dc <- genCluster(ds, "study", "n.study", "id", )
dc <- trtAssign(dc, strata = "study", grpName = "x")
dc <- addColumns(defI, dc)
d.obs <- dc[, .(study, id, x, y)]Build the Stan model
There are multiple ways to estimate a Stan model in R, but I choose to build the Stan code directly rather than using the brms or rstanarm packages. In the Stan code, we need to define the data structure, specify the parameters, specify any transformed parameters (which are just a function of the parameters), and then build the model - which includes laying out the prior distributions as well as the likelihood.
In this case, the model is slightly different from what was presented in the context of a mixed effects model. This is the mixed effects model:
\[ y_{ik} = \alpha_k + \delta_k x_{ik} + e_{ik} \\ \\ \delta_k = \delta_0 + v_k \\ e_{ik} \sim N(0, \sigma_k^2), v_k \sim N(0,\tau^2) \] In this Bayesian model, things are pretty much the same: \[ y_{ik} \sim N(\alpha_k + \delta_k x_{ik}, \sigma_k^2) \\ \\ \delta_k \sim N(\Delta, \tau^2) \]
The key difference is that there are prior distributions on \(\Delta\) and \(\tau\), introducing an additional level of uncertainty into the estimate. I would expect that the estimate of the overall treatment effect \(\Delta\) will have a wider 95% CI (credible interval in this context) than the 95% CI (confidence interval) for \(\delta_0\) in the mixed effects model. This added measure of uncertainty is a strength of the Bayesian approach.
data {
int<lower=0> N; // number of observations
int<lower=1> K; // number of studies
real y[N]; // vector of continuous outcomes
int<lower=1,upper=K> kk[N]; // study for individual
int<lower=0,upper=1> x[N]; // treatment arm for individual
}
parameters {
vector[K] beta; // study-specific intercept
vector[K] delta; // study effects
real<lower=0> sigma[K]; // sd of outcome dist - study specific
real Delta; // average treatment effect
real <lower=0> tau; // variation of treatment effects
}
transformed parameters{
vector[N] yhat;
for (i in 1:N)
yhat[i] = beta[kk[i]] + x[i] * delta[kk[i]];
}
model {
// priors
sigma ~ normal(0, 2.5);
beta ~ normal(0, 10);
tau ~ normal(0, 2.5);
Delta ~ normal(0, 10);
delta ~ normal(Delta, tau);
// outcome model
for (i in 1:N)
y[i] ~ normal(yhat[i], sigma[kk[i]]);
}Generate the posterior distributions
With the model in place, we transform the data into a list so that Stan can make sense of it:
N <- nrow(d.obs) ## number of observations
K <- dc[, length(unique(study))] ## number of studies
y <- d.obs$y ## vector of continuous outcomes
kk <- d.obs$study ## study for individual
x <- d.obs$x ## treatment arm for individual
ddata <- list(N = N, K = K, y = y, kk = kk, x = x)And then we compile the Stan code:
rt <- stanc("model.stan")
sm <- stan_model(stanc_ret = rt, verbose=FALSE)Finally, we can sample data from the posterior distribution:
fit <- sampling(sm, data=ddata, seed = 3327, iter = 10000, warmup = 2500,
control=list(adapt_delta=0.9))Check the diagonstic plots
Before looking at any of the output, it is imperative to convince ourselves that the MCMC process was a stable one. The trace plot is the most basic way to assess this. Here, I am only showing these plots for \(\Delta\) and \(\tau\), but the plots for the other parameters looked similar, which is to say everything looks good:
pname <- c("Delta", "tau")
stan_trace(object = fit, pars = pname)
Look at the results
It is possible to look inspect the distribution of any or all parameters. In this case, I am particularly interested in the treatment effects at the study level, and overall. That is, the focus here is on \(\Delta\), \(\delta_k\), and \(\tau\).
pname <- c("delta", "Delta","tau")
print(fit, pars=pname, probs = c(0.05, 0.5, 0.95))## Inference for Stan model: model.
## 4 chains, each with iter=10000; warmup=2500; thin=1;
## post-warmup draws per chain=7500, total post-warmup draws=30000.
##
## mean se_mean sd 5% 50% 95% n_eff Rhat
## delta[1] 6.39 0.01 1.13 4.51 6.41 8.22 29562 1
## delta[2] -0.78 0.01 1.62 -3.45 -0.78 1.85 28188 1
## delta[3] -0.14 0.01 1.39 -2.37 -0.16 2.18 28909 1
## delta[4] 3.08 0.00 0.59 2.09 3.08 4.05 34277 1
## delta[5] -0.16 0.01 1.01 -1.77 -0.18 1.52 27491 1
## delta[6] 3.87 0.00 0.86 2.47 3.87 5.27 35079 1
## delta[7] 4.04 0.01 1.11 2.21 4.03 5.87 32913 1
## delta[8] 5.23 0.01 1.29 3.12 5.23 7.36 33503 1
## delta[9] 1.79 0.01 1.25 -0.27 1.78 3.82 30709 1
## delta[10] 1.38 0.01 1.12 -0.46 1.38 3.21 30522 1
## delta[11] 4.47 0.01 1.25 2.43 4.47 6.54 34573 1
## delta[12] 0.79 0.01 1.45 -1.60 0.80 3.16 33422 1
## Delta 2.48 0.00 0.89 1.01 2.50 3.89 31970 1
## tau 2.72 0.00 0.71 1.72 2.64 4.01 24118 1
##
## Samples were drawn using NUTS(diag_e) at Sat Jun 27 15:47:15 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).The forest plot is quite similar to the one based on the mixed effects model, though as predicted, the 95% CI is considerably wider:
As a comparison, here is the plot from the mixed effects model estimated using the nlme package in the previous post. The bootstrapped estimates of uncertainty at the study level are quite close to the Bayesian measure of uncertainty; the difference really lies in the uncertainty around the global estimate.
