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

[Yuhang Hu]

Thanks, Wolfgang.

Just curious, if not useful for variance and correlation estimation, then I
wonder when and how BLUPs are needed in practice?

Thanks,
Yuhang

On Wed, Jun 21, 2023 at 8:57 PM Viechtbauer, Wolfgang (NP) <
wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:

> Hi Yuhang,
>
> 1) In the model below, 'outcome' turns into a dummy variable that is 1 for
> PD and 0 for AL. So the intercept random effect represents variation in the
> treatment effects for AL, while the outcomePD random effect represents
> variation in how much the treatment effects differ for PD compared to AL.
> The correlation is then the correlation between these two random effects.
>
> 2) That depends on what you are interested in. But if you want to know how
> much variance there is in a random effect, I would take the variance
> estimate (and not compute the variance of the BLUPs) and similarly if you
> are interested in the correlation between two random effects.
>
> Best,
> Wolfgang
>
> >-----Original Message-----
> >From: Yuhang Hu [mailto:yh342 using nau.edu]
> >Sent: Friday, 16 June, 2023 17:06
> >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
> >
> >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
>


-- 
Yuhang Hu (She/Her/Hers)
Ph.D. Student in Applied Linguistics
Department of English
Northern Arizona University

	[[alternative HTML version deleted]]