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

[Yuhang Hu]

Hi Wolfgang,

Thank you very much for your response. I wanted to ask two
questions regarding your responses.

1-- If instead of struct="UN", I use struct="GEN" (below), then what does
the correlation reported between AL and PD represent?

My understanding is that it represents the correlation between the
'differences' (not between the means of AL and PD) between AL and PD across
the trials, right?

model <- rma.mv(yi~ outcome, vi, data =  dat.berkey1998,
                random = ~ outcome | trial, struct = "GEN")

2-- You mentioned that the estimated variation by the model (say tau2)
accounts for an additional piece (i.e., E(var(u_i|y_i))) in random effects
(ui) that is not present in the BLUPs.

My question is then, is the variation, or correlation obtained by the model
preferred over (and more reliable than) that obtained by hand-calculating
them from the BLUPs?

Many thanks,
Yuhang

On Thu, Jun 1, 2023 at 8:18 PM Viechtbauer, Wolfgang (NP) <
wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:

> You get the estimated variation in the underlying true means of each
> category (around the fixed effects, which correspond to the expected value
> of the underlying true means of each category) and the correlation of that
> variation between the two categories. I have now added the model to the
> write-up of this example here:
>
> https://www.metafor-project.org/doku.php/analyses:berkey1998
>
> so you can see exactly what the code (given on that page) corresponds to
> model-wise.
>
> Best,
> Wolfgang
>
> >-----Original Message-----
> >From: Yuhang Hu [mailto:yh342 using nau.edu]
> >Sent: Thursday, 04 May, 2023 23:39
> >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. Regarding my second (unclear) follow-up, I
> might
> >have a conceptual misunderstanding that I hope to clear up.
> >
> >When we use a "categorical" moderator (with struct="UN") to the left of |
> (e.g.,
> >~outcome | trial_id), conceptually it means that in each unique trial_id,
> we fit
> >"effect_size ~ outcome+0", and obtain the Mean for each outcome category
> (say two
> >categories). Then in the rma.mv() output, we get the variation in means
> of each
> >category as well as the correlation between the means of the two
> categories
> >across all unique trial_ids.
> >
> >Is my conceptual understanding correct?
> >
> >Thanks,
> >Yuhang
> >
> >On Wed, May 3, 2023 at 1:59 PM Viechtbauer, Wolfgang (NP)
> ><wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
> >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
>

	[[alternative HTML version deleted]]