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.
- effect-size as response
- environmental variables as fixed effects,
- sites nested in studies as random effects and
- 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)
- Is the global approach valid?
- Which approach would you suggest?
- Is my R code correct? - with special attention to the weights argument varIdent(~vi) in the lme function.