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.