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

[Viechtbauer, Wolfgang (NP)]

Please see below for responses.

Best,
Wolfgang

>-----Original Message-----
>From: Yuhang Hu [mailto:yh342 using nau.edu]
>Sent: Wednesday, 03 May, 2023 22:48
>To: Viechtbauer, Wolfgang (NP)
>Cc: R Special Interest Group for Meta-Analysis
>Subject: Re: [R-meta] Plotting the correlation among true/random effects across
>categories
>
>Thank you very much, Wolfgang. Two quick follow-ups:
>
>1) To convert these estimated random deviations to true effects, I should add the
>fixed effect estimates (assuming I use 'outcome + 0' in model formula) for each
>category of outcome to its relevant column, right?
>
>AL = paired[,1] + model$b[1,1]
>PD = paired[,2] + model$b[2,1]
>plot(PD~AL, pch=21, bg="gray", cex=1.5, lwd=1.2)

Correct.

>2) When using categorical variables (with "UN") to the left of |, I think we drop
>the intercept in the random-effects design matrix, so what is actually allowed to
>vary across the trials given that eac trial has only one instance of AL and PD in
>it?

Not entirely sure what you mean by this question. The help file explains what the different structures do:

https://wviechtb.github.io/metafor/reference/rma.mv.html#specifying-random-effects

>Thank you for your time.
>Yuhang
>
>On Wed, May 3, 2023 at 12:51 AM Viechtbauer, Wolfgang (NP)
><wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
>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