Arm-based models in Stan

Arm-based pairwise MA models are a precursor to coding arm-based network meta-analysis, and can also be extended to deal with some problems of the evidence base, such as unreported statistics. This post expands on Chapter 8 of our book.

Let’s give each study a unique number j and do the same over all the arms with k. Although we work with arm-level study statistics, we can’t just mix up all the control arms in one batch and all the intervention arms in another, and compare them. The arms are paired with one another and we need to preserve that pairing. To do this, we estimate mu[j], the control mean in study j, and we compare study j’s intervention mean to that (not to the mean of all studies’ control arms).

There is therefore an array mu of length m (number of studies). These should have a diffuse or flat prior, not a heterogeneity distribution. We don’t want them pulled toward one another or “sharing information”.

Common effect, arm-based

Let’s start with studies reporting means and mean differences first. We ultimately are interested in theta, the mean difference. For m studies, there will be K=2m rows in your dataset. Each reports a mean in “arm_mean”, a standard error of the mean in “arm_se”, and number of participants/subjects in “arm_n”. They also have a variable called “active”, which is 1.0 if that arm had the intervention of interest, 0.0 if it is the control (note that “active” and “arm_n” are declared as real, which is not necessary but is a good habit in Stan if you intend to do any arithmetic with data, just as we use .0 decimal places to be clear to the underlying C++ that these are floating point numbers). Finally, each arm belongs to a study, and these consecutive numbers 1, 2, 3… are in a variable called study_id.

data {
  int K;
  int m;
  array[K] real arm_mean;
  array[K] real active;
  array[K] real arm_se;
  array[K] real arm_n;
  array[K] int study_id;
}
transformed data {
  array[K] real df;
  for(k in 1:K) {
    df[k] = arm_n[k] - 1.0;
  }
}
parameters {
  array[m] real mu;
  real theta;
}
model {
  array[K] real predicted;
  theta ~ ADD_YOUR_PRIOR_HERE
  for(j in 1:m) {
    mu ~ normal(0, 1000); // very diffuse
  }
  for(k in 1:K) {
    predicted[k] = mu[study_id[k]] + theta*active[k];
    arm_mean[k] ~ student_t(predicted[k], arm_se[k], df[k]);
  }
}

Notice that we use double indexing in mu[study_id[k]]. Each row in the data is an arm and it belong to the study given in study_id[k]. That study has an estimated control arm mean stored in mu[study_id[k]].

Warning! Using lower case k and upper case K in the same program increases the risk of a typo (Robert has done it, and been baffled for an hour or so because it is hard to spot). This matches the mathematical notation in the book, but practically, it would be better to use a name like n_arms instead of K.

Because arms are smaller in sample size than whole studies, we have switched to the Student t distribution in the likelihood above, although that is not essential.

In the case of binomial likelihoods, studies report numerators (arm_d) and denominators (arm_n) (for example, the ACE inhibitor versus ARB dataset). Here, mu is the log odds in the control arms, and theta is the log odds ratio:

data {
  int K;
  int m;
  array[K] int arm_n;
  array[K] int arm_d;
  array[K] int study_id;
}
parameters {
  array[m] real mu;
  real theta;
}
model {
  array[K] real pred_logodds;
  for(j in 1:m) {
    mu ~ normal(0, 1000); // very diffuse
  }
  theta ~ ADD_YOUR_PRIOR_HERE
  for(k in 1:K) {
    pred_logodds[k] = mu[study_id[k]] +
                      theta*active[k];
    arm_d[k] ~ binomial_logit(arm_n[k],
                              pred_logodds[k]);
  }
}

Note that we can use the binomial_logit() function to compute a binomial likelihood given the denominator and log odds. This avoid converting the log odds to the corresponding risk (probability), but if you need to do this inverse logistic transformation manually, it is:

prob_risk[k] = 1.0 / (1.0 + exp(-prob_logodds[k]));

You might prefer the equivalent alternative code below. Instead of making the dataset “long”, and having one linear predictor formula for all, regression-style, we can keep it “wide” and have two lines of code that contribute to the likelihood, one for each arm. Here, d0 and n0 are the counts for the control arm, and d1 and n1 are for the intervention arm. You should use whichever code template you feel most comfortable with, as your task is to understand deeply what it does and be able to explain it.

data {
  int m;
  array[m] int d0;
  array[m] int n0;
  array[m] int d1;
  array[m] int n1;
}
parameters {
  array[m] real mu;
  real theta;
}
model {
  mu ~ normal(0, 100);
  theta ~ ADD_YOUR_PRIOR_HERE;
  for(j in 1:m) {
  	d0[j] ~ binomial_logit(n0[j], mu[j]);
  	d1[j] ~ binomial_logit(n1[j], mu[j]+theta);
  }
}

Connection to fixed effects

Now you have seen code for the mu vector with a flat or diffuse prior. This is in fact the same coding approach that we would apply to a vector of theta[j] if we needed to fit a “fixed effects” meta-analysis, where each study has its own theta[j], but they are minimally “shrunken”, “pulled together”, or “sharing information”. In a Bayesian context, all unknowns have priors (even if don’t declare one, there is a flat prior), so the prior on the mu vector might look like a heterogeneity distribution when you read this code for the first time, but is fixed and very diffuse.

Random effects, arm-based

We have seen that, in order to respect the pairing of arms within studies, we have to introduce one more unknown to our model for each study. Further, these can be expected to correlate negatively in the posterior distributions with the intervention effect: if you under-estimate all the studies’ control arm means or log odds, then you have to over-estimate the mean difference or log odds ratio in order to estimate the intervention arms. This all makes the task more challenging for Bayesian algorithms like Metropolis-Hastings and the Gibbs sampler, which we find in BUGS, JAGS, Stata, and parts of JASP. In contrast, the Hamiltonian Monte Carlo algorithm in Stan (hence brms and multiNMA) and PyMC will deal with a higher-dimensional, correlated posterior.

We start again with arms reporting means. For small sample sizes, you can use the student_t() likelihood with arm_n[k]-1 degress of freedom, where k is an index of all the arms, and arm k belongs to study study_id[k].

data {
  int m;
  int K;
  array[K] real arm_mean;
  array[K] real arm_se;
  array[K] int arm_n;
  array[K] int study_id;
}
transformed data {
  array[K] real df;
  for(k in 1:K) {
    df[k] = arm_n[k] - 1.0;
  }
}
parameters {
  real theta;
  real<lower=0> tau;
  array[m] u;
  array[m] mu;
}
model {
  array[K] real predicted;

  theta ~ ADD_YOUR_PRIOR_HERE;
  tau ~ ADD_YOUR_PRIOR_HERE;
 
  for(j in 1:m){
    mu[j] ~ normal(0, 1000); # very diffuse
    u[j] ~ normal(0, tau); # heterogeneity
  }

  for(k in 1:K){
    predicted[k] = mu[study_id[k]] + (theta+u[study_id[k]])*active[k];
    arm_mean[k] ~ student_t(predicted[k], arm_se[k], df[k])
  }
}

The changes necessary for binary outcomes are only those previously shown:

data {
  int m;
  int K;
  array[K] int arm_d;
  array[K] int arm_n;
  array[K] int study_id;
}
parameters {
  real theta;
  real<lower=0> tau;
  array[m] u;
  array[m] mu;
}
model {
  array[K] real pred_logodds;

  theta ~ ADD_YOUR_PRIOR_HERE;
  tau ~ ADD_YOUR_PRIOR_HERE;
 
  for(j in 1:m){
    mu[j] ~ normal(0, 1000); # very diffuse
    u[j] ~ normal(0, tau); # heterogeneity
  }

  for(k in 1:K){
    pred_logodds[k] = mu[study_id[k]] + (theta+u[study_id[k]])*active[k];
    arm_d[k] ~ binomial_logit(arm_n[k], pred_logodds[k])
  }
}

And this is the code for a “wide” data format, as shown in common effects above:

data {
  int m;
  array[m] int d0;
  array[m] int n0;
  array[m] int d1;
  array[m] int n1;
}
parameters {
  array[m] real mu;
  real theta;
  array[m] real u;
  real<lower=0> tau;
}
transformed parameters {
	real or;
	or = exp(theta);
}
model {
  mu ~ normal(0, 100);
  theta ~ normal(0, 3);
  tau ~ gamma(1, 1);
  for(j in 1:m) {
  	u[j] ~ normal(0, tau);
  	d0[j] ~ binomial_logit(n1[j], mu[j]);
  	d1[j] ~ binomial_logit(n2[j], mu[j]+theta+u[j]);
  }
}

Once we have mu[j] and u[j], theta and tau, we have 2m+2 unknowns but 2m arms, so we are slicing the study information very thinly indeed, and trusting in the sharing of information in the model to get us through. When the studies are all the same size, none is too small, and the events not too rare or too common, which we can achieve in simulation (see the “all pairwise models” posts), this model will work. However, real-life evidence bases are not so obliging. The ACEI vs ARB dataset has one very large study and some of the others are too small to contribute much information at all. Further, the event is fairly rare, so the absolute numbers get very small. If you try these models with that dataset, you may see some problems in MCMC convergence.

Fully Bayesian models

The “fully Bayesian” model has a likelihood contribution from the arm mean and another from the arm standard deviation (see Section 8.4 in the book for explanation of the sampling distribution of standard deviations and variances).

There are important modelling assumptions to be chosen and justified regarding the population standard deviation of the outcome variable. In this random effects example, we have heterogeneity on the study means but we assume a common population SD shared by all control arms (sigma_ctl), and another for all intervention arms (sigma_int). This can be changed in various ways, as long as the assumptions can be justified. For example, there could be one sigma common to all arms in all studies, or some form of heterogeneity distribution could be applied across studies.

data {
  int K; // number of arms in the data
  int m; // number of studies
  array[K] int study_id;  // in (1:m)
  array[K] real mean;
  array[K] real<lower=0> sd;
  array[K] real treat;  // 0.0 or 1.0
  array[K] real<lower=0> n;
}
parameters {
  real theta;
  real<lower=0> tau;
  array[m] real mu;
  array[m] real u;
  real<lower=0> sigma;
}
model {
  array[K] real<lower=0> varx2;
  array[K] real<lower=0> arm_sigma;

  # priors:
  theta ~ ADD_YOUR_PRIOR_HERE;
  tau ~ ADD_YOUR_PRIOR_HERE;
  sigma_int ~ ADD_YOUR_PRIOR_HERE;
  sigma_ctl ~ ADD_YOUR_PRIOR_HERE;

  for(j in 1:m) {
    # very diffuse:
    mu[j] ~ normal(0, 1000);  
    # heterogeneity:
    u[j] ~ normal(0, tau);
  }

  # likelihood:
  for(k in 1:K) {
    # likelihood of mu, u, theta, given sigmas
    mean[k] ~ normal(mu[study_id[k]] +
                      treat[k]*theta,
                     arm_sigma[k]);
    # likelihood of sigmas
    arm_sigma[k] = treat*sigma_int +
                      (1.0-treat)*sigma_ctl;
    varx2[k] = sd[k]*sd[k]*(n[k]-1.0) /
                (arm_sigma[k]*arm_sigma[k]);
    varx2[k] ~ chi_square(n[k]-1.0);
  }
}

tags: #stan, #stan-repo, #repo