[R-meta] Multilevel, multivariate meta analysis on studies with multiple outcomes and dependent effect sizes

[Viechtbauer, Wolfgang (SP)]

Hi Neda,

I don't think you ever received an answer.

1) Since the levels of 'combination' are not always unique within 'studyID' (i.e., the same combination can appear multiple times within the same study), you should add one more level, namely the estimate level as another random effect. So:

res1 <- rma.mv(yi, vi, data=mydata, mods=~factor(A)+factor(B)+C,
                            random = ~1|studyID/combination/effectsizeID,
                            method="REML", tdist=TRUE))

Not trying to make this more complicated than necessary, but you could also consider using crossed random effects here for 'treatment' and 'control':

mydata$treatment.in.study <- paste(mydata$studyID, mydata$treatment)
mydata$control.in.study   <- paste(mydata$studyID, mydata$ control)

res2 <- rma.mv(yi, vi, data=mydata, mods=~factor(A)+factor(B)+C,
                            random = list(~ 1 | studyID, 
                                          ~ 1 | treatment.in.study, 
                                          ~ 1 | control.in.study, 
                                          ~ 1 | effectsizeID),
                            method="REML", tdist=TRUE))

You could compare the fit of these models then:

fitstats(res1, res2)

I don't think these models are nested, so a likelihood ratio test would not be appropriate, but based on the AIC(c)/BIC values, you can still compare their fit.

2) One can think of a network meta-analysis as a hierarchical model, so your second question isn't really answerable.

3) I would always use the outermost level as the clustering variable when using cluster robust inference methods, so yes, using 'studyID' seems appropriate. And the robust() function gives you the same output for the fixed effects as the output for a model, so the interpretation is the same.

Best,
Wolfgang

-----Original Message-----
From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces using r-project.org] On Behalf Of Nasiri Moghadam, Neda
Sent: Wednesday, 29 May, 2019 10:38
To: r-sig-meta-analysis using r-project.org
Subject: [R-meta] Multilevel, multivariate meta analysis on studies with multiple outcomes and dependent effect sizes

Hi everyone,

I’m doing a meta analysis to find the connection between environmental fluctuation and fitness performance of organisms. I have collected data (844 Hedge’s g mean difference effect sizes ) from 43 studies (multiple effect sizes per study).

Dependency between effect sizes within each study arises from one of the following conditions: 

1. several treatments (treatment.type) have been compared with a common control group, 
2. several control conditions (control.type) have been compared with one treatment, 
3. several control and treatment conditions have been compared together. 

The treatment and control conditions are study-specific therefore there is no dependency between studies. 
In each study effect of treatment has been assessed on different traits (growth, yield, survival and so on). The traits and their numbers are different between studies. The applied organism for the experiments is different between studies (insect, bacteria, virus and ...).

Here is a data sample to clear my explanation:
 
Due to the dependency between effect sizes within each study (I don't know the covariances between the different outcomes) I combined information of the treatment.type and control.type in a new column (combination) using the paste function in R. I nested this column in study ID and included in the model as a random factor. 

res <- rma.mv(yi,vi,data=mydata,mods=~factor(A)+factor(B)+C,
>                           random = ~1|studyID/combination,
>                           method="REML",tdist=TRUE))

I've also made the phylogeny tree and correlation matrix between species. I'm going to include the species as a random factor in the model. I checked the likelihood profile plot for sigma^2.1 and sigma^2.2. In both cases the plot peaks at the estimate. 

The distribution of the heterogeneity is as follows:

level1: 59%,
level2: 7%,
level3: 32% 

I obtained cluster robust standard error using the following command: 

robust(res, cluster=mydata$studyID)

The funnel plot is not symmetric and based on the cook's distance plot I have influential data points.

Questions:

1) Is the model correct based on the hierarchical structure in my data set (study > combination > trait > effect size)
2) Isn't it better to use a network meta analysis instead of a three level model for my data set?
3) Is it correct to use studyID as the cluster in robust() function and how should I interpret the output of robust() function?

I’m happy to provide any additional information to facilitate feedback.
Thank you for your help in advance.

Bests,

Neda N. Moghadam 
Postdoctoral Researcher
Department of Biological and Environmental Science
P.O. Box 35, FI-40014 
University of Jyväskylä, Finland
Mobile: +4550393628
Email: neda.n.nasiri-moghadam using jyu.fi
Linkedin: linkedin.com/in/nedanasirimoghadam