[R-meta] Network Meta-analysis metafor

[Viechtbauer Wolfgang (SP)]

Dear Caroline,

The example in help(dat.senn2013) uses mean differences to illustrate a network MA with metafor but is easy to adjust to standardized mean differences. For SMDs, you have to use a contrast-based model. The dat.senn2013 dataset is given in an arm-based format, so it first needs to be restructured. Note that Willms (1999) is the only multiarm study in this dataset, but the code should work regardless of how many multiarm studies there are. Just make sure that the 'reference' treatment for multiarm studies always goes under 'trt1'.

library(metafor)

############################################################################

### load data
dat <- get(data(dat.senn2013, package="metafor"))

### restructure dataset to a contrast-based format
dat <- dat.senn2013[,c(1,4,3,2,5,6)]
dat.c <- lapply(split(dat, dat$study),
                function(x) cbind(x[rep(1,nrow(x)-1),], x[-1,c(3:6)]))
dat.c <- do.call(rbind, dat.c)
names(dat.c)[3:10] <- c("trt1", "n1i", "m1i", "sd1i", "trt2", "n2i", "m2i", "sd2i")
rownames(dat.c) <- 1:nrow(dat.c)
dat.c$id <- 1:nrow(dat.c)
dat <- dat.c
rm(dat.c)

############################################################################

### make sure dataset is sorted by 'study'
dat <- dat[order(dat$study),]

### compute the total sample sizes of the studies
dat$Ni <- unlist(lapply(split(dat, dat$study), function(x) rep(x$n1i[1] + sum(x$n2i), each=nrow(x))))

### compute the pooled SDs of the studies
dat$sdpi <- unlist(lapply(split(dat, dat$study), function(x) sqrt(((x$n1i[1]-1)*x$sd1i^2 + sum((x$n2i-1)*x$sd2i^2)) / ((x$n1i[1] - 1) + sum(x$n2i - 1)))))

### calculate SMDs and corresponding sampling variances
dat$yi <- with(dat, (m1i-m2i)/sdpi)
dat$vi <- with(dat, 1/n1i + 1/n2i + yi^2/(2*Ni))

### calculate the variance-covariance matrix of the SMDs for multitreatment studies
calc.v <- function(x) {
   v <- matrix(1/x$n1i[1] + outer(x$yi, x$yi, "*")/(2*x$Ni[1]), nrow=nrow(x), ncol=nrow(x))
   diag(v) <- x$vi
   v
}
V <- bldiag(lapply(split(dat, dat$study), calc.v))

### convert trt1 and trt2 variables to factors and set levels
lvls <- levels(factor(c(dat$trt1, dat$trt2)))
dat$trt1 <- factor(dat$trt1, levels=lvls)
dat$trt2 <- factor(dat$trt2, levels=lvls)

### create variables to indicate the contrasts examined
dat <- cbind(dat, model.matrix(~ dat$trt1 - 1) - model.matrix(~ dat$trt2 - 1))
names(dat)[(ncol(dat)-9):ncol(dat)] <- lvls

############################################################################

### network meta-analysis using a contrast-based random-effects model
### by setting rho=1/2, tau^2 reflects the amount of heterogeneity for all treatment comparisons
### the treatment left out (placebo) becomes the reference level for the treatment comparisons
res <- rma.mv(yi, V, mods = ~ acarbose + benfluorex + metformin + miglitol + pioglitazone +
                              rosiglitazone + sitagliptin + sulfonylurea + vildagliptin - 1,
              random = ~ factor(id) | study, rho=1/2, data=dat)
res

### estimate/test all pairwise differences between treatments
### need 'multcomp' package for contrMat() function
library(multcomp)
contr <- contrMat(setNames(rep(1,res$p), colnames(res$X)), type="Tukey")
sav <- predict(res, newmods=contr)
sav[["slab"]] <- rownames(contr)
sav$pval <- anova(res, L=contr)$pval
sav

############################################################################

Two points:

1) Note that I am using the pooled SD of all arms within each trial to compute the SMDs. That way, results within trials are internally consistent (i.e., A vs B = (A vs C) - (B vs C) for a three-arm study). If you just use the pooled SD of each pairwise comparison, then this won't be the case.

2) I see that you have the reference treatment listed under treat2. If you want to use the code above, just rename your variables to: study trt2 trt1 n2i m2i sd2i n1i m1i sd1i.

Best,
Wolfgang

-- 
Wolfgang Viechtbauer, Ph.D., Statistician | Department of Psychiatry and 
Neuropsychology | Maastricht University | P.O. Box 616 (VIJV1) | 6200 MD 
Maastricht, The Netherlands | +31 (43) 388-4170 | http://www.wvbauer.com 

>-----Original Message-----
>From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces at r-
>project.org] On Behalf Of Caroline Woehl
>Sent: Wednesday, 14 February, 2018 8:45
>To: r-sig-meta-analysis at r-project.org
>Subject: [R-meta] Network Meta-analysis metafor
>
>Dear all,
>
>I've only calucated network meta-analyses with the package netmeta yet.
>Due to the possibility to calculate subgroup analyses, I'd like to have
>a try at a network meta-analysis with metafor.
>
>The analysis is based on 20 randomized controlled trials. They compare
>physical activity or cognitive activity with usual care. One of them is
>a multi-arm trial comparing physical and cognitive activity with usual
>care. I could extract sample size, mean und standard deviation
>(posttreatment) from each study arm. Thus, the file contains the study,
>treat1, treat2, n1i, m1i, sd1i, n2i, m2i and sd2i.
>
>Study                                                treat1       treat2
>        n1i          m1i      sd1i        n2i          m2i       sd2i
>       1      physical activity    usual care     37   17,16    4,33
>36   15,17    4,48
>       2                    physical activity    usual care     39
>19,40    8,68    35   18,48    8,40
>       2                    cognitive activity   usual care    36
>19,98    8,16    35   18,48    8,40
>       3      cognitive activity usual care    35    42,34   9,30    33
>  44,49   9,53
>
>So far, I’ve calculated the standardized mean differences (yi) and
>variances (vi) for each pairwise comparison.
>nma <- escalc(measure="SMD", m1i=m1i, sd1i=sd1i, n1i=n1i, m2i=m2i,
>sd2i=sd2i, n2i=n2i, data=nma)
>
>How can I define usual care as the reference group in the next step?
>
>Many thanks for your help.
>
>Best regards,
>Caroline