I would like to use network meta-analysis (or multiple treatment comparison) to investigate the main-effect and interaction of two factors on a continuous outcome variable. As an example the effect of hair-color (dark/light) and gender (male/female) on a blood pressure (continuous). I have mean, variance and sample size from studies that compared dark-males to dark-females (studies 1 & 2), studies that compared dark-males to light-males (study 5), studies that compared dark-females with light-females, studies that compared dark-males with light-males and light-females (study 7). Most importantly I also have studies that compared males to females irrespective of hair-color (studies 3 and 4).

my_data <- data.frame(study = c("Study_1","Study_1", "Study_2", "Study_2", "Study_3", "Study_3", "Study_4", "Study_4", "Study_5", "Study_5", "Study_6", "Study_6", "Study_7", "Study_7", "Study_7"), treatment = c("dark_male", "dark_female","dark_male", "dark_female","male", "female","male", "female", "dark_male", "light_male", "dark_female", "light_female", "dark_male", "light_male", "light_female"), mean=c(10,12,11,14,11,21,14,25,10,12,14,20,12,14,20), std.dev=c(1,1.2,1,1.2,2,2.1,1,1.3,1.3,1,1,1,1.2,1.2,1.3), sampleSize=c(10,10,20,20,15,15,35,35,18,19,20,22,40,41,40))


I am not sure if I can include studies 3 and 4 since they include light- and dark-haired subjects.

Is it possible to investigate the main effects of hair-color, gender as well as their interaction using one of the available r packages for network meta-analysis (gemtc, netmeta, metafor)? What would be a suitable way to calculate post-hoc tests?


1 Answer 1


I think you will have to make some kind of assumption about the proportion of dark/light hair-color individuals in studies 3 and 4. Then you can treat the 'dark' color variable as a proportion that is equal to 0 in samples composed entirely of light hair-color individuals, 1 in samples composed entirely of dark hair-color individuals, and something in between when it is a mix (as in studies 3 and 4). For gender, you just set up a standard dummy variable. Then you can analyze these data with the rma.mv function from the metafor package (which allows fitting network-type meta-analytic models). Here is how this would be done, assuming that 50% of the individuals in studies 3 and 4 were dark hair-color individuals:

# dummy variable for gender and proportion of dark hair-color individuals
my_data$male <- c(1,0,1,0,1,0,1,0,1,1,0,0,1,1,0)
my_data$dark <- c(1,1,1,1,.5,.5,.5,.5,1,0,1,0,1,0,0)

# compute sampling variances of the means
my_data$vi <- my_data$std.dev^2 / my_data$sampleSize

# need a variable to indicate the arm within each study 
my_data$arm <- sequence(table(my_data$study))

# meta-analysis model using correlated random effects within each study
res <- rma.mv(mean, vi, mods = ~ male*dark, random = list(~ arm | study), data=my_data)

The output from this is:

Multivariate Meta-Analysis Model (k = 15; method: REML)

Variance Components: 

outer factor: study (nlvls = 7)
inner factor: arm   (nlvls = 3)

            estim    sqrt  fixed
tau^2      6.8909  2.6250     no
rho        0.3316             no

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

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

Model Results:

           estimate      se     zval    pval     ci.lb    ci.ub     
intrcpt     21.9548  1.5913  13.7965  <.0001   18.8359  25.0738  ***
male        -9.0630  2.0448  -4.4322  <.0001  -13.0708  -5.0552  ***
dark        -6.8660  2.1176  -3.2424  0.0012  -11.0164  -2.7156   **
male:dark    4.9911  2.7210   1.8343  0.0666   -0.3419  10.3241    .

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

Maybe there is an interaction, but strictly speaking, it is not significant at $\alpha = .05$ (two-sided). Actually, let's get the predicted means for the 'female-light', 'female-dark', 'male-light', and 'male-dark' combinations:

predict(res, newmods=rbind(c(0,0,0), c(0,1,0), c(1,0,0), c(1,1,1)))

     pred     se   ci.lb   ci.ub   cr.lb   cr.ub
1 21.9548 1.5913 18.8359 25.0738 15.9383 27.9714
2 15.0888 1.3667 12.4101 17.7676  9.2882 20.8894
3 12.8918 1.5816  9.7919 15.9917  6.8851 18.8985
4 11.0169 1.2026  8.6599 13.3739  5.3577 16.6761

So, for females, going from 'light' to 'dark' results in a pretty substantial decrease, but not so for males.

Again, I made an assumption about the proportion of dark hair-color individuals in studies 3 and 4. My suggestion now would be to carry out a sensitivity analysis, varying those proportions within a sensible range (e.g., .3 to .7), checking whether the conclusions depend on the specific proportions assumed.

  • $\begingroup$ Thanks for the fantastic & comprehensive answer - as always ;) Would the suggested analysis change if the above presented data represented raw means and their variances instead of effect sizes? $\endgroup$
    – jokel
    Jan 5, 2014 at 16:46
  • $\begingroup$ I actually assumed that these are raw means and that the SDs given are just the SDs of the raw scores within the groups (vi, that is, the sampling variances, were computed under that assumption, since the variance of a mean is equal to the variance of the raw scores divided by the sample size). $\endgroup$
    – Wolfgang
    Jan 5, 2014 at 18:29
  • 1
    $\begingroup$ By the way, if you feel my answer adequately addresses your question, feel free to accept it (the same applies to stats.stackexchange.com/q/71404/1934). $\endgroup$
    – Wolfgang
    Jan 5, 2014 at 18:31
  • $\begingroup$ Did so - thanks again! Could you maybe suggest a way to do a post-hoc test to compare different factor levels (e.g. dark-males vs. dark females). Also could you elaborate on the columns specification in the predict() call - is it hair-color, gender, interaction? $\endgroup$
    – jokel
    Jan 9, 2014 at 23:03
  • $\begingroup$ The order of the values specified via the newmods argument corresponds to the order of the variables under the Model Results. So, it is male, dark, and then the interaction male:dark. So, newmods=c(1,1,1) would be male=1, dark=1, and hence male:dark=1. And newmods=c(0,1,0) would be male=0 (i.e., female), dark=1, and hence male:dark=0. $\endgroup$
    – Wolfgang
    Jan 10, 2014 at 21:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.