Solved – Main effects and interaction in multivariate meta-analysis (network meta-analysis) in R

meta-analysismultivariate analysisnetwork-meta-analysisr

I would like to conduct a multivariate meta-analysis (or network meta-analysis) to investigate the effect of two factors (A & B, 1 & 2). For every combination of the two factors I have a number of studies that report results that describe the difference in effect size as well as the variance thereof (A1 vs B1, A1 vs B2, A2 vs B1, A2 vs B2). Every effect size comes from an separate study. I would like to investigate the main-effects as well as the interaction of my two factors. Here I provide a dummy example of the data in R that I create initially for trying out the gemtc package:

set.seed(100)
treat1 <- c(rep("A1", 10),rep("A1", 8),rep("A2", 15),rep("A2", 5))
treat2 <- c(rep("B1", 10),rep("B2", 8),rep("B1", 15),rep("B2", 5))
TE <- c(rnorm(10,5,1),rnorm(8,2,1),rnorm(15,3,1),rnorm(5,1,1))
seTE <- c(rnorm(10,1,1),rnorm(8,2,1),rnorm(15,0.5,1),rnorm(5,0.2,1))
studlab <- paste(rep("study_", 38), as.character(c(1:38)), sep="")
my_data <- data.frame(treat1, treat2, TE, seTE, studlab)
my_data

I would be thankful for a pointer to an appropriate R package (gemtc, netmeta, mvmeta ?) that could handle such data and how to model not only main effects but also interactions of my factors. Some hints how to get started would be fantastic…

Best Answer

Given that there is a single effect size estimate from each study (see comments above), the analysis can be carried out with regular meta-regression methods. You can carry out such an analysis with the metafor package. The "trick" is to code variables that indicate what treatments have been compared within a particular study:

library(metafor)

my_data$A1 <- ifelse(treat1 == "A1", 1, 0)
my_data$A2 <- ifelse(treat1 == "A2", 1, 0)
my_data$B1 <- ifelse(treat2 == "B1", -1, 0)
my_data$B2 <- ifelse(treat2 == "B2", -1, 0)

res <- rma(TE, sei=seTE, mods = ~ A1 + A2 + B1 - 1, data=my_data)
res

yields:

Mixed-Effects Model (k = 38; tau^2 estimator: REML)

tau^2 (estimated amount of residual heterogeneity):     0.2898 (SE = 0.1578)
tau (square root of estimated tau^2 value):             0.5384
I^2 (residual heterogeneity / unaccounted variability): 59.02%
H^2 (unaccounted variability / sampling variability):   2.44

Test for Residual Heterogeneity: 
QE(df = 35) = 93.5215, p-val < .0001

Test of Moderators (coefficient(s) 1,2,3): 
QM(df = 3) = 435.5223, p-val < .0001

Model Results:

     estimate      se     zval    pval    ci.lb    ci.ub     
A1    2.2446  0.2837   7.9123  <.0001   1.6886   2.8006  ***
A2    0.9060  0.3387   2.6751  0.0075   0.2422   1.5699   **
B1   -2.2983  0.3467  -6.6294  <.0001  -2.9778  -1.6188  ***

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Since variable B2 has been left out, this becomes the "reference" treatment. So, the coefficient for A1 is the estimated average effect when comparing treatment A1 against B2. The coefficient for A2 is the estimated average effect when comparing treatment A2 against B2. And the coefficient for B1 is the estimated average effect when comparing treatment B1 against B2.

The network that is analyzed here looks like this:

    A1   A2
    |\   /|
    | \ / |
    |  X  |
    | / \ |
    |/   \|
    B1   B2

So, the comparison between B1 and B2 is based purely on indirect evidence.

There are 3 more comparisons that can be obtained here besides the ones above (i.e., A1 vs A2, A1 vs B1, and A2 vs B1). You can obtain those by changing the "reference" treatment.

An assumption made here is that the amount of heterogeneity is the same regardless of the comparison. This may or may not be true.

An article that describes this type of analysis is:

Salanti et al. (2008). Evaluation of networks of randomized trials. Statistical Methods in Medical Research, 17, 279-301.

Edit: To test whether the effect of the first factor (A & B) depends on the second factor (1 & 2), that is, whether (A1 vs B1) = (A2 vs B2) or not, first note that:

(A1-B2) - (A2-B2) - (B1-B2) = (A1-A2) - (B1-B2) 
                            = (A1-B1) - (A2-B2).

So, you just have to test whether b1 - b2 - b3 = 0. You can do this with:

predict(res, newmods=c(1,-1,-1))

or install/load the multcomp package and use:

summary(glht(res1, linfct=rbind(c(1,-1,-1))), test=adjusted("none"))

which yields:

         Simultaneous Tests for General Linear Hypotheses

Linear Hypotheses:
       Estimate Std. Error z value Pr(>|z|)    
1 == 0   3.6369     0.5779   6.293 3.12e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- none method)
Related Question