Adapting inlaNMA to BRMS

[timdisher]

Hello everyone,

I am trying to adapt code from nmaINLA to allow for an easy interface for running and expanding these models from stan. Network meta-analysis is an extension of pairwise meta-analysis for more than two treatments. The code for these models is primarily based off code from NICE TSD 2 with some extensions to allow for inconsistency between direct (A vs B) and indirect (A vs B through C vs A and C vs B).

Stepping through the nma_inla function from their example code

library(nmaINLA)
library(INLA)
SmokdatINLA ← create_INLA_dat(dat = Smokdat, armVars = c(‘treatment’ = ‘t’, ‘responders’ = ‘r’
,‘sampleSize’ = ‘n’), nArmsVar = ‘na’)
if(requireNamespace(‘INLA’, quietly = TRUE)){
require(‘INLA’, quietly = TRUE)
fit.Smok.cons.INLA ← nma_inla(SmokdatINLA, likelihood = ‘binomial’, type = ‘consistency’,
tau.prior = ‘uniform’, tau.par = c(0, 5))
}

They create a long format data frame with study indicators for non-baseline arms (stored in het below) and indicators for multi-arm trials (g) and the basic parameters (d12 …).

study treatment responders sampleSize na baseline mu d12 d13 d14 g het inc
1 1 9 140 3 1 1 0 0 0 NA NA NA
1 3 23 140 3 1 1 0 1 0 1 1 1
1 4 10 138 3 1 1 0 0 1 2 1 1
2 2 11 78 3 2 2 0 0 0 NA NA NA
2 3 12 85 3 2 2 -1 1 0 1 2 2
2 4 29 170 3 2 2 -1 0 1 2 2 2
3 1 75 731 2 1 3 0 0 0 NA NA NA
3 3 363 714 2 1 3 0 1 0 1 3 3

The multi-arm bit is used to re-build equation 14 from NICE TSD2 (p.35). Ultimately, the model call to INLA is:

INLA::inla(stats::as.formula(inla.form), data = datINLA,
family = “binomial”, verbose = verbose, control.fixed = list(expand.factor.strategy = “inla”,
mean = fixed.par[1], prec = 1/fixed.par[2]), Ntrials = Ntrials,
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE,
mlik = TRUE, config = TRUE), control.inla = list(strategy = inla.strategy,
correct = correct, correct.factor = correct.factor))

Where fixed.par is a vague prior on treatment effects, Ntrials is sample size in the binomial case and inla.form is

“Y ~ -1 + mu + d12+d13+d14 + f(het, model="iid", hyper = list(theta = list(prior = unif.prior.table)), group = g, control.group = list(model = "exchangeable", hyper = list(rho = list(fixed = TRUE, initial = cor.inla.init))))”

Where Y is the number of events in arm i of study j and mu is the outcome in the baseline treatment of that arm. Where I’m stuck is how to translate the f() portion of the call to something for BRMS. I had originally assumed 1 | studyid would give me the single between-study variance parameter I was looking for but the multi-arm piece has me stumped.

[paul.buerkner]

I am not used to reading INLA code. Could you write down the statistical model itself? Then I can better tell you whether such a model is possible in brms.

[anon75146577]

I’m on a phone but it’s a random effects model here the civariance matrix is the kroneker product of an iid random effect and an exchangeable random effect (sigma^2R) where R is a correlation matrix where all the off diagonal elements are the same. The fixed=TRUE part means that the correlation parameter is fixed rather than inferred.

I’m travelling today but let me know if you need more.

[paul.buerkner]

Thanks Dan! I don’t yet fully understand all the details, but based on what Dan writes,
you might be able to solve this via argument cov_ranef of brm, which allows you to specify fixed correlation matrices between random effects of the same level. For usage of cov_ranef you may also look at vignette("brms_phylogenetics") although the application presented there is very different from your use case.

[timdisher]

dat.csv (2.1 KB)

Have the time (and courage, thanks to a nudge from @avehtari) to re-visit this and thought it best to work from a toy example up to something more realistic.

Using the attached simulated data I have fit a fixed an random effect meta-analysis that returns the correct value/should be equivalent to my target model:

brm_fe <- brms::brm(r | trials(n) ~ -1 + factor(study_id) + t, family = binomial(), data = long_dat)` 

brm_re <- brms::brm(r | trials(n) ~ -1 + factor(study_id) + t + (-1 + t | study_id), 
                    control = list(adapt_delta = 0.99), cores = 4, iter = 5000, family = binomial(), data = long_dat)

I needed to change adapt delta from 0.8 to 0.9 to 0.99 in the re model because of divergent transitions.

Now the first issue would be when extending the random effect model to more than two treatments, it needs to be specified so that:

delta_t1 <- 0
delta_t2 ~ dnorm(t2, sd_shared)
delta_t3 ~ dnorm(t3, sd_shared)

So the group level effects on the treatments are independent in the sense that they don’t come from a common distribution (i.e. I don’t want to shrink estimates for t2, t3 together) but they have a common sd that’s estimated.

My understanding is that the current re specification I have won’t accomplish that. I will get 1 sd but my estimates for t2/t3 will be shrunken together. Is that true? Is there an easy workaround if so?

[paul.buerkner]

Fist, in brm_re you should remove factor(study_id) as you now have the study_id grouping variable include via (... | study_id).

Is t the treatment variable? you may model what you have in mind by writing the data in long format to have (you already do I assume) to have a single varying effect for the treatments coded as 0, 1, 1 for t1, t2, and t3, respectively. That is use (1 + td | study_id) for this dummy variable. That would ensure the equal SD on t2 and t3. Then, in the fixed effects part, just go with 1 + t to have varying overall effects of t2 and t3.

[timdisher]

@paul.buerkner

Fist, in brm_re you should remove factor(study_id) as you now have the study_id grouping variable include via (... | study_id) .

Will specifying it this way keep a study specific fixed effect baseline? This model is based on current industry-standard guidelines

@paul.buerkner

Is t the treatment variable? you may model what you have in mind by writing the data in long format to have (you already do I assume) to have a single varying effect for the treatments coded as 0, 1, 1 for t1, t2, and t3, respectively.

Yes that’s correct. So the idea is I would have a dataframe that looks like the one below, correct?

study  t  td
1      0  0
1      1  1
1      2  1

That is use (1 + td | study_id) for this dummy variable.

This would introduce a random baseline (mu[i] in notation above) correct? I think to replicate the target model I would need to avoid this as per:

"An important feature of all the meta-analytic models presented here is that no model is assumed
for the trial-specific baselines µi. They are regarded as nuisance parameters which are
estimated in the model. An alternative is to place a second hierarchical model on the trial
baselines, or to put a bivariate normal model on both.32,33 However, unless this model is correct,
the estimated relative treatment effects will be biased"

The final piece of this would be to incorporate correlation in random effects for multi-arm trials.

Is my assumption that translation of (td | study) will accomplish this?

[paul.buerkner]

If we just had two groups, (3) with a hierarchial prior would be expressed as

... ~ 1 + t + (1 + t | study_id)

@timdisher

“An important feature of all the meta-analytic models presented here is that no model is assumed for the trial-specific baselines µi. They are regarded as nuisance parameters which are estimated in the model. An alternative is to place a second hierarchical model on the trial baselines, or to put a bivariate normal model on both.32,33 However, unless this model is correct, the estimated relative treatment effects will be biased”

"An important feature of all the meta-analytic models presented here is that no model is assumed
for the trial-specific baselines µi. They are regarded as nuisance parameters which are
estimated in the model. An alternative is to place a second hierarchical model on the trial
baselines, or to put a bivariate normal model on both.32,33 However, unless this model is correct,
the estimated relative treatment effects will be biased"

This is a weird way off expressing the situation and I disagree with the description, but if this is what you want to follow, you need to use the 0 + factor(study_id) approach (which I argue against) but then you cannot reasonably also include (1 | study_id).

I am not sure if (td | study_id) would achieve what you want since you don’t distinquish between the “delta” groups and hence do not estimate a correlation (other than the ones with the baseline which you don’t want to include apparently).

[timdisher]

@paul.buerkner

I am not sure if (td | study_id) would achieve what you want since you don’t distinquish between the “delta” groups and hence do not estimate a correlation (other than the ones with the baseline which you don’t want to include apparently).

So I think you are perhaps more in line with the “arm-based” approach that has been popping up and is (as far as I can tell) more firmly connected with standard approaches to hierarchical models.

The model for the arm-based approach is:

where i and k index studies and treatments. Would this just be

brms::brm(r | trials(n) ~ t + (1 + t | study_id),
control = list(adapt_delta = 0.99), cores = 4, iter = 5000, family = binomial(), data = long_dat)

Example code for target model in JAGS would be

for(ss in 1:nStudies){
alpha[ss] ~ dnorm(0,1.0E-4)
for (tt in 1:nTx){
re[ss,tt] ~dnorm(0,reTau)
}
}

 #fit data 
 for(ii in 1:nObs ){
  x[ii] <- alpha[study[ii]] + beta[tx[ii]]+re[study[ii],tx[ii]]
 
  #logit link for probability of response in control and treatment arms
  logit(prob[ii]) <- x[ii]

  #binomial likelihood
  r[ii] ~ dbin(prob[ii], n[ii])

[paul.buerkner]

Yes indeed the brms model looks reasonable to me.

As a side note, your Jags code looks slightly incorrect as you don’t seem to model the correlations there. Or perhaps I am overlooking something.

[timdisher]

@paul.buerkner

As a side note, your Jags code looks slightly incorrect as you don’t seem to model the correlations there. Or perhaps I am overlooking something.

No, I think you’re right typically these are modelled with mvnorm for random effects but it’s just dnorm in that code for some reason, thanks for point it out.

One thing I often see written of these models is that the first set I showed you are referred to as contrast models and only rely on relative effects being constant across treatments while the second set of models are said to assume absolute effects. I am always slightly confused by this because if, for example, I were analyzing a multi-site trial where I thought treatment effect would vary by site I would do something like:

y ~ t + (1 + t | site)

This seems analogous to the arm-based models but I wasn’t under the impression that I was modeling anything other than relative effects.

[paul.buerkner]

Me neither. Not sure what they mean here exactly. Of course you could do

y ~ 0 + t + (0 + t | site)

to model absolute effects.

[timdisher]

@paul.buerkner

to model absolute effects.

Ah yes, that seems to align with the model specification above. I wonder why the specification you suggested hasn’t been used.

One last question:
With:
y ~ t + (1 + t | study)

Is it still possible to get a single SD for each t while estimating correlation within multi-arm studies? perhaps by specifying them as exchangeable?

[paul.buerkner]

With this approach you get a single SD per treatment.

[timdisher]

Thanks Paul! This worked, returned the correct parameter from my simulated dataset and resulted in loo actually working.

[paul.buerkner]

Glad I could help!