[R-meta] Clarifications about workflow for complex meta-analyses with dependent effect sizes

Tue, Nov 7, 2023 1:16 AM

Dear Alberto,

Please see below for some responses.

Best,
Wolfgang

You can get prediction intervals for the 4 timepoints with:

predict(res, newmods=diag(4), tau2.levels=1:4)

To add these to the plot, you could do:

pred <- predict(res, newmods=diag(4), tau2.levels=1:4)
arrows(1:4, pred$pi.lb, 1:4, pred$pi.ub, lwd=3, code=3, angle=90, length=0.1)

The 0.8 in the presentation was chosen a bit arbitrarily. In the example code here:

https://wviechtb.github.io/metadat/reference/dat.ishak2007.html

V is calculated using phi=0.97.

Technically, there is a different value of the autocorrelation coefficient for each study. If one had the raw data for a study, one could easily compute this. However, in practice, this is often not available. Often, we just use a 'guestimate' (which may be based on one or two studies for which the raw data are actually available) and assume that this applies to all studies.

In the actual analyses of these data, the authors (Ishak et al., 2007) actually estimated the value of phi from the studies at hand. In essence, we just profile the likelihood over different values of phi and check for which value of phi the likelihood is maximized:

phis <- seq(0.8, 0.99, by=0.01)
lls <- sapply(phis, function(phi) {
V <- vcalc(vi, cluster=study, time1=time, data=dat, phi=phi)
res <- rma.mv(yi, V, mods = ~ factor(time) - 1, random = ~ time | study,
              struct="HAR", data=dat, digits=2)
logLik(res)
})

plot(phis, lls, type="o", pch=21, bg="lightgray")
abline(v=phis[which.max(lls)])
phis[which.max(lls)]

This leads to the estimate of 0.97. In my experience, trying to do this with other datasets often doesn't work so well and the estimate of the correlation drifts towards -1 or +1. In fact, the estimated value above seems awfully high to me. So this is why I just (somewhat naively) used 0.8 in the presentation. However, the assumed value of phi has only a minor impact on the results, so it doesn't matter (at least in this example) too much what we use in the first place.

In general, whenever we construct a V matrix that is really just an approximation, I would recommend to use cluster-robust inference methods as the last step:

robust(res, cluster=study, clubSandwich=TRUE)

Here, the SEs are barely changed and the degrees of freedom are quite high, so this provides evidence that the assumed model is an adequate approximation.

I don't know this paper, but if they have repeated Log.RR values based on the same sample of subjects, then this model seems to ignore this (since it doesn't use a proper V matrix and just the var.log values for the sampling variances, in essence assuming that V is diagonal). Also, 'random = ~ 1 | Study.ID/effect.number' implies that the correlation between multiple effect sizes within studies is constant, but this is often not so sensible for longitudinal data (where the correlation tends to become weaker for effect sizes further apart in time).

[R-meta] Clarifications about workflow for complex meta-analyses with dependent effect sizes

Thread (3 messages)