I'm working on a Bayesian Network Meta-Analysis using the gemtc package on a dataset similar to the following:
df <- data.frame(study = c("A", "A", "B", "B", "C", "C", "D", "D", "E", "E", "F", "F",
"G", "G", "H", "H", "I", "I", "J", "J", "K", "K"),
treatment = c("A", "B", "B", "C", "B", "C", "A", "B", "B",
"C", "B", "C", "A", "B", "B", "C", "B", "C",
"A", "C", "B", "C"),
responders = c(1, 5, 0, 0, 3, 1, 0, 2, 0, 2, 0, 2, 0,
0, 1, 2, 0, 0, 2, 9, 1, 1),
sampleSize = c(32, 33, 30, 30, 18, 20, 15, 15, 20,
20, 30, 30, 36, 32, 15, 15, 23, 22, 24, 23, 18, 16))
While I have been able to set up the network model and run the analysis just fine, I have been struggling with specifying the order in which I would like the treatments to be compared in the node-splitting consistency analysis. For example, I want the odds ratios and 95% credible intervals to be calculated using the "B" treatment as the reference group when comparing "B" with "A" and "C" as the reference group when comparing "A" with "C" and "B" with "C". Below is the code I have tried:
library(gemtc)
library(rjags)
# Create mtc.network element to be used in modeling ------
network <- mtc.network(data.ab = df)
# Compile model ------
network.mod <- mtc.model(network,
linearModel = "random", # random effects model
n.chain = 4) # 4 Markov chains
# Assess network consistency using nodesplit method ------
nodesplit <- mtc.nodesplit(network.mod,
linearModel = "random", # random effects model
n.adapt = 5000, # burn-in iterations
n.iter = 100000, # actual simulation iterations
thin = 10) # extract values of every 10th iteration
summary(nodesplit) # High p-values indicate consistent results
plot(summary(nodesplit))
My results provide ORs (95% CrIs) for:
- "A" vs. "C"
- "B" vs. "C"
- "B" vs. "A"
I have created a separate data frame specifying that I want "A" vs. "B" comparisons via:
# Specify desired comparisons ------
comparisons = data.frame(t1 = "A", t2 = "B")
# Assess network consistency using nodesplit method, adding comparisons argument ------
nodesplit <- mtc.nodesplit(network.mod,
comparisons = comparisons,
linearModel = "random", # random effects model
n.adapt = 5000, # burn-in iterations
n.iter = 100000, # actual simulation iterations
thin = 10) # extract values of every 10th iteration
summary(nodesplit) # High p-values indicate consistent results
But I still get "B" vs. "A" results. I have also tried specify t1="B", t2="A", and I get the same results. Any assistance with this would be greatly appreciated. Thanks in advance.