[R-meta] Plotting the correlation among true/random effects across categories

Viechtbauer, Wolfgang (NP) wo||g@ng@v|echtb@uer @end|ng |rom m@@@tr|chtun|ver@|ty@n|
Wed May 3 09:51:01 CEST 2023


You can extract the BLUPs of the random effects and create a scatterplot based on them:

re <- ranef(model)
re

paired <- do.call(rbind, split(re[[1]]$intrcpt, dat.berkey1998$trial))
paired

plot(paired, pch=21, bg="gray", cex=1.5, lwd=1.2)

And before somebody asks why cor(paired) does not yield 0.7752 (or why the values in var(paired) do not match up with the variances as estimated from the model), see for example:

https://stats.stackexchange.com/q/69882/1934

Or let me give a more technical explanation based on the standard RE model:

dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
res <- rma(yi, vi, data=dat)
res$tau2
var(ranef(res)$pred)

You will notice that the latter is smaller than tau^2. By the law of total variance:

tau^2 = var(u_i) = E(var(u_i|y_i)) + var(E(u_i|y_i)).

The conditional means of the random effects (which is what ranef() provides estimates of) are E(u_i|y_i) and hence their variance is only part of the total variance. Therefore, the estimate of tau^2 and the estimated variance of the BLUPs of the random effects will not match up.

In more complex models, this then also affects things like the correlation between the BLUPs of the random effects.

Best,
Wolfgang

>-----Original Message-----
>From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces using r-project.org] On
>Behalf Of Yuhang Hu via R-sig-meta-analysis
>Sent: Wednesday, 03 May, 2023 1:21
>To: R meta
>Cc: Yuhang Hu
>Subject: [R-meta] Plotting the correlation among true/random effects across
>categories
>
>Hello Colleagues,
>
>I was wondering if there is a way to scatterplot the correlation between
>the categories of variable "outcome" (AL and PD) which is estimated to be
>rho = .7752 in my model below?
>
>model <- rma.mv(yi~ outcome, vi, data =  dat.berkey1998,
>                random = ~ outcome | trial, struct = "UN")
>
>    rho.AL  rho.PD    AL  PD
>AL       1             -   5
>PD  0.7752       1    no   -
>
>Thanks,
>Yuhang



More information about the R-sig-meta-analysis mailing list