4

I have a collection of continuous data from the literature, including the mean, the standard deviation and the number of observations for both experimental and control groups, as well some environmental variables. A meta-analysis coupled with a meta-regression could be done. However, some studies had several experimental sites, i.e. each line is not a single study, but a single site taken from a study: a random site effect could thus be nested in studies. The following meta mixed-model can be described.

  1. effect-size as response
  2. environmental variables as fixed effects,
  3. sites nested in studies as random effects and
  4. each observation is weighted on the sampling variance.

I tried two approaches with R, that returned somewhat different results.

As first step, let's create a dummy data frame.

set.seed(40)
Study = factor(c(1,1,1,1,1,2,2,3,3,3,4,5,5,5,5))
Site = factor(c(1,2,3,4,5,1,2,1,2,3,1,1,2,3,4))
n_case = length(Study)
experimental_mean = rnorm(n_case,5,2)
experimental_sd = abs(rnorm(n_case,1,0.1))
experimental_n = round(runif(n_case, 2, 10))
control_mean = experimental_mean+runif(n_case,0,5)
control_sd = abs(rnorm(n_case,1,0.1))
control_n = round(runif(n_case, 2, 10))
A = experimental_mean/control_mean * runif(n_case, 0.5, 0.8)
B = rnorm(n_case,-1,1)
C = factor(letters[round(runif(n_case, 1, 3))])

data_table = data.frame(Study, Site, 
                        experimental_mean, experimental_sd, experimental_n,
                        control_mean, control_sd, control_n,
                        A, B, C)

Then, compute the effect size and the sampling variance.

library(metafor)
meta_table = escalc(measure='ROM', data=data_table,
                  m1i=experimental_mean, m2i=control_mean,
                  sd1i=experimental_sd, sd2i=control_sd,
                  n1i=experimental_n, n2i=control_n)

The rma.mv function from the metafor package could be used to run the meta mixed-model.

res = rma.mv(yi, vi, data=meta_table, method="REML", level=95,
              mods = ~ A+B+C, random=~1|Study/Site)
summary(res)

Another option is to run a mixed-model on the effect-size.

library(nlme)
mixed_meta_model = lme(data = meta_table,
                       fixed = yi ~ A+B+C, # yi is the effect size
                       random = ~1|Study/Site,
                       weights = varIdent(~vi)) # vi is the sampling variance
summary(mixed_meta_model)

Three questions:

  1. Is the global approach valid?
  2. Which approach would you suggest?
  3. Is my R code correct? - with special attention to the weights argument varIdent(~vi) in the lme function.
amoeba
  • 93,463
  • 28
  • 275
  • 317
essicolo
  • 185
  • 2
  • 7

1 Answers1

6

Indeed,

res <- rma.mv(yi, vi, data=meta_table, method="REML", level=95,
              mods = ~ A+B+C, random=~1|Study/Site)

will fit a three-level meta-analytic model. For a pretty extensive illustration of this model (based on Konstantopoulos, 2011), take a look at http://www.metafor-project.org/doku.php/analyses:konstantopoulos2011 (the districts in that example corresponds to the studies in your case and the studies in that example correspond to the sites in your case).

On the metafor package website, I have also written up a rather extensive discussion as to what lme() does when you try to fit meta-analytic models with it (see: http://www.metafor-project.org/doku.php/tips:rma_vs_lm_and_lme). The varIdent(~vi) part is in principle correct, but the lme() function will assume that the sampling variances are known up to a proportionality constant (what is estimated as the residual variance by the function). That is different than what we typically want to do in meta-analyses (where the sampling variances are pretty much known).

References

Konstantopoulos, S. (2011). Fixed effects and variance components estimation in three-level meta-analysis. Research Synthesis Methods, 2(1), 61–76.

Wolfgang
  • 15,542
  • 1
  • 47
  • 74