I am running a random-effects multivariate meta-analysis with the package metafor. I have binary data with one treatment and one control group. I am interested in the overall effect of the treatment. I did not want to average the effect sizes, but use all effect sizes I have, as of by now I just include 15 effect sizes. This post comes close to this one here, however, I am not quite sure if I do it correctly. So I would be more than happy if someone could help. Please see the generic data below. Regarding the dataset: Study 5 has 3 different treatment groups always compared to the no treatment group, therefore these effect sizes are correlated. For study 6 the effect sizes are also correlated, as the same participants are measured 2 weeks later again.
data <- structure(list(es_id = c(1, 2, 3, 4, 5, 6 ,7 ,8, 9, 10, 11, 12, 13, 14, 15),
study = c(1, 1, 2, 2, 3, 4, 5, 5, 5, 6, 6, 6, 6, 7, 7),
obs= c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 3L, 1L, 2L, 1L, 2L, 1L, 2L),
ai = c(22, 10, 16, 17, 19, 124, 15, 9, 9, 27, 14, 14, 13, 70, 79),
n1i = c(32, 20, 26, 27, 29, 134, 25, 19, 22, 37, 24, 24, 23, 80, 90),
ci = c(18, 6, 12, 13, 15, 120, 5, 5, 5, 23, 10, 10, 9, 66, 66),
n2i = c(31, 19, 24, 24, 30, 144, 25, 25, 25, 38, 29, 38, 29, 80, 80)),
.Names = c("Number", "Study", "obs", "ai", "n1i", "ci", "n2i"), row.names = c(NA, -15L), class = "data.frame")
library(metafor)
d1.1 <- escalc(measure = "OR", ai = ai, n1i = n1i, ci = ci, n2i = n2i, data = data)
I am running 2 different options.
The problem with Option #1 are the profile plots, which have a gap for sigma 2:
Option1 <- rma.mv(yi = yi,
V = vi,
random = list(~ 1 | obs,
~ 1 | Number,
~ 1 | Study),
tdist = F,
data = d1.1,
method = "REML")
summary(Option1)
robust(Option1, cluster = Study)
par(mfrow=c(3,1))
profile(Option1, sigma2=1, xlim=c(0,.6))
profile(Option1, sigma2=2, xlim=c(0,.6))
profile(Option1, sigma2=3, xlim=c(0,.6))
par(mfrow=c(1,1))
For the second option, I was thinking of just including Study/obs, which does not show any gaps in the profile plots
Option2 <- rma.mv(yi = yi,
V = vi,
slab = Study,
data = d1.1,
random = ~ 1 | Study/obs)
summary(Option2)
robust(Option2, cluster = Study)
par(mfrow=c(2,1))
profile(Option2, sigma2=1, xlim=c(0,.6))
profile(Option2, sigma2=2, xlim=c(0,.6))
par(mfrow=c(1,1))
I was also thinking of changing it to a long format as I saw in the example section of the GitHub page by Prof. Viechtbauer (here). Option 1 does not show any profile plots however almost no heterogeneity for sigma2 and sigma3.
dat.long <- to.long(measure="OR", ai=ai, bi=n1i, ci=ci, di=n2i, data=data)
dat.long$obs <- c(rep(1,14), 1, 2, 1, 3, 1, 1, 2, 2, 1, 1, 2, 2, 1, 1, 1 , 2) # reassigning number of observations
dat.long[,c("ai","n1i","ci","n2i","study")] <- NULL # deleting variables which are not needed any longer
dat.long <- escalc(measure="PLO", xi=out1, mi=out2, data=dat.long)
res <- rma.mv(yi, vi, mods = ~ group, random = ~ 1 | Study/group/obs, data=dat.long)
par(mfrow=c(3,1))
profile(res, sigma2=1, xlim=c(0,.6))
profile(res, sigma2=2, xlim=c(0,.6))
profile(res, sigma2=3, xlim=c(0,.6))
par(mfrow=c(1,1))
Another option would also be to delete the reference groups that are included multiple times and to average the results of study 6, as this is the only study that includes measures with a time gap in between. Not quite sure, if this is ok, as it is the only study that would be averaged.
I hope someone could help me out.
Update:
I tried to construct the variance-covariance matrix as described on the metafor page here by the following code. However, I am not quite sure if this takes care of the dependent effect sizes (same control group and one time participants got measured two weeks after the first measurement).
data$pti <- with(data, ai / n1i)
data$pci <- with(data, ci / n2i)
data <- escalc(measure = "OR", ai = ai, n1i = n1i, ci = ci, n2i = n2i, data = data)
calc.v <- function(x) {
v <- matrix(1/(x$n2i[1]x$pci[1](1-x$pci[1])), nrow=nrow(x), ncol=nrow(x))
diag(v) <- x$vi
v
}
V <- bldiag(lapply(split(data, data$Study), calc.v))
res <- rma.mv(yi, V, random = ~ 1 | Study/Number, data=data)
summary(res)
robust(res, cluster = Study)
Thanks again in advance for any help.
As I am quite new to meta-analysis, so is the following approach also correct, if I stick to metafor:
data <- structure(list(es_id = c(1, 2, 3, 4, 5, 6 ,7 ,8, 9, 10, 11, 12, 13, 14, 15),
study = c(1, 1, 2, 2, 3, 4, 5, 5, 5, 6, 6, 6, 6, 7, 7),
obs= c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 3L, 1L, 2L, 1L, 2L, 1L, 2L),
ai = c(22, 10, 16, 17, 19, 124, 15, 9, 9, 27, 14, 14, 13, 70, 79),
n1i = c(32, 20, 26, 27, 29, 134, 25, 19, 22, 37, 24, 24, 23, 80, 90),
ci = c(18, 6, 12, 13, 15, 120, 5, 5, 5, 23, 10, 10, 9, 66, 66),
n2i = c(31, 19, 24, 24, 30, 144, 25, 25, 25, 38, 29, 38, 29, 80, 80)),
.Names = c("Number", "Study", "obs", "ai", "n1i", "ci", "n2i"), row.names = c(NA, -15L), class = "data.frame")
library(metafor)
data$pti <- with(data, ai / n1i)
data$pci <- with(data, ci / n2i)
data <- escalc(measure = "OR", ai = ai, n1i = n1i, ci = ci, n2i = n2i, data = data)
calc.v <- function(x) {
v <- matrix(1/(x$n2i[1]x$pci[1](1-x$pci[1])), nrow=nrow(x), ncol=nrow(x))
diag(v) <- x$vi
v
}
V <- bldiag(lapply(split(data, data$Study), calc.v))
res <- rma.mv(yi, V2, random = ~ Number | Study, struct = "UN", data=data)
summary(res)
robust(res, cluster = Study)
I also saw the vcalc function by metafor. I get the same final results as well as the same results on the diagonal, however, when constructing the covariance matrix I get the following warning for all my clusters: 1: The var-cov matrix appears to be not positive definite in cluster 1.
V2 <- vcalc(vi, cluster = Study, data = data)
Can I use this and does it take care of my problem that I sometimes have the same control group (study 5 and 7) and for one study, the same participants are measured two weeks after the first one (study 6)?
