nma.ab generates a summary file for effect sizes by conducting the arm-based approach, proposed by Zhang et al (2013), in network meta-analysis. The generated summary file contains odds ratio (OR) and patient-centered paramenters, such as risk ratio (RR), risk difference (RD), and absolute risk (AR). Also, it can provide DIC statistics for checking the goodness of fit, and give trace plot and posterior density plot for the population-averaged event rates.nma.ab(s.id, t.id, event.n, total.n, o.path = getwd(), f.name = "",
model = "het2", sigma2.a = 0.001, sigma2.b = 0.001, mu.prec = 0.001,
R, param = c("probt","RR","RD","OR","rk","best"), ity = "estimate",
n.iter = 200000, n.burnin = floor(n.iter/2), n.chains = 3,
n.thin = max(1, floor((n.iter - n.burnin)/100000)), dic = TRUE,
trace = FALSE, postdens = FALSE)dic = TRUE), trace plot (if trace = TRUE) and posterior density plot (if p"hom" (homogeneous variance, denoted as model HOM), "het1" (model HET1), or "het2" (mosigma2.a and sigma2.b are set as the default, the prior is weatN by tN covariance matrix for Wishart prior in the model HET2 (model = "het2"), where tN is the number of treatments. The default is a matrix with diagonal elements being 1 and off-diagonmu in the JAGS model. It can be set as "estimate" (the default) and "same". See "Details".n.iter/2.TRUE (the default), the function would generate a file containing the DIC statistics, and a node named "deviance" would be contained in the summary result file; otherwise, the DIC statistics would not be calculated.FALSE (the default), the function would not give the trace plot; otherwise, trace plots for population-averaged event rate ("probt") would be given, and they would be saved in the working directory where o.path specifFALSE (the default), the function would not give the posterior density plot; otherwise, posterior density plots for population-averaged event rate ("probt") would be given, and they would be saved in the working directory where nma.ab generates a summary statistics file using the arm-based method. Furthermore, this function would give a DIC statistics file if dic = TRUE, a trace plot file if trace = TRUE, a posterior density file if postdens = TRUE.
In the summary file, each row contains statistics for corresponding OR, RD, RR, population-averaged event rate ("probt"), rank of treatment ("rk"), probability of being the best treatment ("best"), etc. Note that RR[i, j], RD[i, j] or OR[i, j] means that treatment i is compared with treatment j, e.g., RD[i,j] = probt[i] - probt[j]. The columns show the statistics of these nodes, including mean, standard deviance, 2.5% percentile, median, and 97.5% percentile.
The DIC file contains the value of pD and DIC; the trace plot and posterior density plot file contain the corresponding plots for the node "probt" (population-averaged event rate).ity is specified as "estimate", the initial value for fixed effect would be estimated (in probit scale) as the proportion of sum of event numbers for a treatment over sum of total numbers in corresponding study; otherwise, the initial value of mu in the JAGS model would be set as 0. Setting ity as "estimate" will be helpful for the Markov chain converging more quickly.
The homogeneous model (model = "hom") considers a common variance for the random effects, and it assumes that the random effects for different treatments are the same in each study. The JAGS model is given as follows:
model{
for(i in 1:sN){
p[i] <- phi(mu[t[i]] + sigma*vi[s[i]])
r[i] ~ dbin(p[i], totaln[i])
}
for(j in 1:tS){
vi[j] ~ dnorm(0, 1)
}
sigma <- 1/sqrt(tau)
tau ~ dgamma(sigma2.a, sigma2.b)
for(j in 1:tN){
mu[j] ~ dnorm(0, mu.prec)
probt[j] <- phi(mu[j]/sqrt(1 + 1/tau))
}
for(j in 1:tN){
for(k in 1:tN){
RR[j, k] <- probt[j]/probt[k]
RD[j, k] <- probt[j] - probt[k]
OR[j, k] <- probt[j]/(1 - probt[j])/probt[k]*(1 - probt[k])
}
}
rk[1:tN] <- tN + 1 - rank(probt[])
best[1:tN] <- equals(rk[], 1)
}
The first heterogeneous model (model = "het1") accounts for the heterogeneity for the variances of random effects, but it still assumes that the random effects for different treatments are the same in each study. The following shows the corresponding JAGS model:
model{
for(i in 1:sN){
p[i] <- phi(mu[t[i]] + sigma[t[i]]*vi[s[i]])
r[i] ~ dbin(p[i], totaln[i])
}
for(j in 1:tS){
vi[j] ~ dnorm(0, 1)
}
for(j in 1:tN){
mu[j] ~ dnorm(0, mu.prec)
tau[j] ~ dgamma(sigma2.a, sigma2.b)
sigma[j] <- 1/sqrt(tau[j])
probt[j] <- phi(mu[j]/sqrt(1 + 1/tau[j]))
}
for(j in 1:tN){
for(k in 1:tN){
RR[j, k] <- probt[j]/probt[k]
RD[j, k] <- probt[j] - probt[k]
OR[j, k] <- probt[j]/(1 - probt[j])/probt[k]*(1 - probt[k])
}
}
rk[1:tN] <- tN + 1 - rank(probt[])
best[1:tN] <- equals(rk[], 1)
}
The second heterogeneous model (model = "het2") covers the most general cases, and it employs a Wishart prior for the inverse of covariance matrix for random effects. The JAGS model is defined as follows:
model{
for(i in 1:sN){
p[i] <- phi(mu[t[i]] + vi[s[i], t[i]])
r[i] ~ dbin(p[i], totaln[i])
}
for(j in 1:tS){
vi[j, 1:tN] ~ dmnorm(mn[1:tN], T[1:tN, 1:tN])
}
invT[1:tN, 1:tN] <- inverse(T[,])
for(j in 1:tN){
mu[j] ~ dnorm(0, mu.prec)
sigma[j] <- sqrt(invT[j, j])
probt[j] <- phi(mu[j]/sqrt(1 + invT[j, j]))
}
T[1:tN, 1:tN] ~ dwish(R[1:tN, 1:tN], tN)
for(j in 1:tN){
for(k in 1:tN){
RR[j, k] <- probt[j]/probt[k]
RD[j, k] <- probt[j] - probt[k]
OR[j, k] <- probt[j]/(1 - probt[j])/probt[k]*(1 - probt[k])
}
}
rk[1:tN] <- tN + 1 - rank(probt[])
best[1:tN] <- equals(rk[], 1)
}data(Ara09)
attach(Ara09)
set.seed(12345)
## nma.ab(s.id = Study.ID, t.id = Treatment, event.n = r, total.n = n,
## model = "hom", f.name = "Ara09_hom_", n.iter = 500)
## nma.ab(s.id = Study.ID, t.id = Treatment, event.n = r, total.n = n,
## model = "het1", f.name = "Ara09_het1_", n.iter = 500)
nma.ab(s.id = Study.ID, t.id = Treatment, event.n = r, total.n = n,
model = "het2", f.name = "Ara09_het2_", n.iter = 500)
detach(Ara09)
data(Lam07)
attach(Lam07$data)
set.seed(12345)
nma.ab(s.id = Study.ID, t.id = Treatment, event.n = r, total.n = n, model = "hom",
f.name = "Lam07_", ity = "same", n.iter = 500, n.chains = 2,
param = c("probt", "RR", "RD", "OR", "rk", "best", "mu", "sigma", "vi"),
trace = TRUE, postdens = TRUE)
detach(Lam07$data)Run the code above in your browser using DataLab