2
$\begingroup$

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 :

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 function by . 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)?

$\endgroup$

1 Answer 1

1
$\begingroup$

You either need to construct the covariance matrix, as stated, which you could quite easily do by a) analyzing the data from each study e.g. using logistic regression and outputting the estimates along with their covariance matrix, or b) based on theoretical asymptotic considerations, which should give you a very similar results with a reasonable number of events.

Alternatively, you can just switch to logistic regression instead of trying to figure out the correlations, for which the input format would be one record per study arm (containing a studyid, the number of cases and the denominator). E.g. a (exact, if you worry about small sample size effects with very few events in some arms of some studies) logistic regression with treatment main (fixed) effects and study main (fixed) effects is a way of doing a fixed effects meta-analysis for binomial data. If you want a random effects meta-analysis, you would additionally specify a study random effect on the treatment effect.

$\endgroup$
1
  • $\begingroup$ Thanks a lot for your reply. I used the answer function due to character restrictions. $\endgroup$
    – Joshua
    Jun 8, 2022 at 12:47

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.