[R-meta] Problem with sigma-square in 3-level meta-analysis

[Viechtbauer, Wolfgang (NP)]

Dear Sebastian,

Just a general note when posting to this list: Please do not reply to a message and then change the subject to something new. There is information stored in the email header that is used to create threads in email clients and the mailing list archives and when doing this, it messes things up. See https://stat.ethz.ch/pipermail/r-sig-meta-analysis/2023-March/thread.html and note how your message is considered a reply to the 'margins command' thread, which is obviously not what you intended. (I am making this reply in an entirely new post, so this should hopefully show up in a proper new thread, although of course your post won't show up at the beginning of this thread, so let me just put the archive URL here: https://stat.ethz.ch/pipermail/r-sig-meta-analysis/2023-March/004447.html).

Moving on to your actual questions:

(a) The observed variances are a combination of the variance in the underlying true effects (apparently r-to-z transformed correlations based on 'zr') and the sampling variances and hence cannot be directly used to examine whether something has gone wrong here.

I would examine profile(overall_mod_cc) to see if the likelihood profile plots are properly peaked at the estimates (which could very well also happen around 0).

(b) I wouldn't put it this way, but yes, in general (and this is not specific to rma.mv()) it is difficult to estimate variance components when the number of levels is small. The REML estimates should still be approximately unbiased, but will then tend to high variance themselves.

(c) Again, I wouldn't automatically conclude here that something is off. It might very well be that the variance in the underlying true effects is due to difference between rows within samples and not between samples.

Speaking of that -- are different rows for the same level of 'Sample_ID' based on the same subjects? In this case, those (r-to-z transformed) correlations are not independent and should not be treated as such. rcalc() is a function that allows to construction of the appropriate var-cov matrix to account for this (but requires information about correlations that may not be available). Alternatively, vcalc() could be used to construct a rough approximation.

And finally -- if you still think something is off here, a possible way forward would be to use a Bayesian model where you can use priors on the variance components to push their posterior distributions in the desired direction (more or less, depending on how much information there is in the data).

Best,
Wolfgang

>-----Original Message-----
>From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces using r-project.org] On
>Behalf Of Röhl, Sebastian via R-sig-meta-analysis
>Sent: Friday, 10 March, 2023 8:40
>To: R Special Interest Group for Meta-Analysis
>Cc: Röhl, Sebastian
>Subject: [R-meta] Problem with sigma-square in 3-level meta-analysis
>
>Dear all,
>
>We are conducting meta-analyses of correlational effects which are clustered in
>studies (or independent samples). Therefore we are using 3-level hierarchical
>random effects models (Sampling error within effect sizes within independent
>samples).
>For one of our sub-analyses which 50 effects of 8 independent samples, the
>results of the rma.mv pointed to nearly no (5.710374e-12) variance component
>between the samples:
>
>> overall_mod_cc <- rma.mv(zr, V=var, random = ~ 1| Sample_ID / nummer, data =
>data_cc, method = "REML")
>> overall_mod_cc
>Multivariate Meta-Analysis Model (k = 50; method: REML)
>Variance Components:
>            estim    sqrt  nlvls  fixed            factor
>sigma^2.1  0.0000  0.0000      8     no         Sample_ID
>sigma^2.2  0.0222  0.1489     50     no  Sample_ID/nummer
>
>However, when I look directly at the manifest variances for the group means
>(0.225) and all effects (0.392), I conclude that there should actually be a much
>larger variance component between the independent samples.
>(a) Am I making a thinking error here?
>(b) Or could this be because the estimation of variance components using rma.mv()
>does not work reliably for the small number of studies (8)?
>(c) And if (b) is correct - what model or function could I use instead for this
>partial analysis?
>
>Thanks a lot for your help!
>
>Best,
>Sebastian