[R-meta] Metafor modeling question

Tue Mar 28 10:36:49 CEST 2023

I would start by running:

par(mfrow=c(3,2))
profile(mm)

which will indicate that none of the var-cor components are identifiable.

The only study with any variability in SES is study 18 (all other studies have SES = 1), which makes analyzing the association between the effect sizes and SES a difficult endeavor.

Not sure how I would approach this. Possibly, I might consider just looking at study 18 by itself to examine this relationship:

plot(yi ~ SES, data=dat, subset=study==18, pch=21, bg="gray", cex=0.5/sqrt(vi))

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: Tuesday, 28 March, 2023 6:35
>To: R meta
>Cc: Yuhang Hu
>Subject: Re: [R-meta] Metafor modeling question
>
>Dear Wolfgang,
>
>Thanks you so much for your response. Regarding the fitted vs residuals
>plot, if we see a violation (e.g., heteroskedasticity and/or non-random
>pattern in residuals' distribution around the regression line) as in the
>case below, what should we do?
>
>I would imagine that if I was working with an rma.uni() model, I could
>model the heteroskedasticity using the 'scale=' argument?
>
>Thanks,
>Yuhang
>
>mm <- rma.mv(yi ~ SES, vi, random = list(~SES|study, ~SES|outcome),
>             struct = c("GEN","GEN"), data = dat)
>
>plot(fitted(mm), resid(mm))
>
>dat <- structure(list(study = c(1L, 1L, 2L, 2L, 3L, 3L, 4L, 5L, 6L,
>7L, 7L, 8L, 8L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 11L, 12L,
>13L, 13L, 14L, 14L, 15L, 15L, 16L, 16L, 17L, 18L, 18L, 18L, 18L,
>18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L,
>18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L,
>18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 19L, 19L, 20L,
>20L), group = c(1, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1,
>2, 2, 1, 1, 2, 2, 1, 1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1,
>1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2,
>2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1), outcome = c("A",
>"A", "B", "B", "B", "B", "B", "A", "A", "A", "B", "A", "A", "A",
>"B", "A", "B", "A", "B", "A", "B", "B", "A", "A", "B", "B", "B",
>"A", "A", "B", "B", "A", "B", "B", "B", "B", "B", "B", "B", "B",
>"B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B",
>"B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B",
>"B", "B", "B", "B", "B", "B", "A", "B", "A", "B"), SES = c(1,
>1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
>1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 7, 13, 10, 8, 16, 3, 2, 18, 1,
>5, 4, 20, 11, 9, 14, 17, 15, 6, 12, 19, 7, 13, 10, 8, 16, 3,
>2, 18, 1, 5, 4, 20, 11, 9, 14, 17, 15, 6, 12, 19, 1, 1, 1, 1),
>    yi = c(0.661693, 0.594366, 1.199301, 1.626664, 1.216272,
>    -0.751534, -1.170659, 0.959518, 0.353214, 1.595024, 1.689419,
>    0.462459, 1.324527, 0.61639, 0.615806, 0.616544, 0.617199,
>    1.420133, 1.542125, 1.338082, 1.523729, 0.569801, 1.041618,
>    6.796397, 3.241257, 0.97541, 1.369298, 1.706266, 1.964219,
>    1.554362, 1.693888, 0.939487, 0.349941, 0.364232, 0.41826,
>    0.37037, 0.378151, 0.265958, 0.250613, 0.32379, 0.318541,
>    0.331738, 0.373367, 0.242215, 0.420524, 0.153225, 0.287845,
>    0.43818, 0.315411, 0.264576, 0.33054, 0.35668, 0.393478,
>    0.468176, 0.342097, 0.39018, 0.332325, 0.414179, 0.303209,
>    0.290351, 0.398231, 0.283192, 0.290464, 0.352266, 0.323642,
>    0.357534, 0.18504, 0.349599, 0.322425, 0.330017, 0.273626,
>    0.484511, 3.087989, -0.427958, 1.133235, 1.073414), vi = c(0.030703,
>    0.355305, 0.190678, 0.525746, 0.257342, 0.371465, 0.65578,
>    0.970636, 0.605049, 0.584372, 0.185532, 0.140877, 0.618397,
>    0.880055, 0.29737, 0.186943, 0.376905, 0.327622, 0.232547,
>    0.288472, 0.88999, 0.203553, 0.202477, 0.94609, 0.694122,
>    0.79155, 0.260697, 0.689051, 0.332882, 0.883863, 0.265871,
>    0.443643, 0.71795, 0.325523, 0.720786, 0.905436, 0.777423,
>    0.424495, 0.449885, 0.436236, 0.551392, 0.412329, 0.841006,
>    0.03359, 0.592244, 0.159267, 0.430026, 0.044465, 0.341106,
>    0.163984, 0.037114, 0.9678, 0.215364, 0.111536, 0.292115,
>    0.548628, 0.25953, 0.217442, 0.615006, 0.108288, 0.80669,
>    0.789254, 0.312538, 0.996466, 0.277414, 0.975438, 0.710994,
>    0.028242, 0.488088, 0.921909, 0.622881, 0.145062, 0.240907,
>    0.461285, 0.881521, 0.974101)), row.names = c(NA, -76L), class =
>"data.frame")
>
>On Mon, Mar 27, 2023 at 12:49 AM Viechtbauer, Wolfgang (NP) <
>wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
>
>> Dear Yuhang,
>>
>> Please see below for my answers.
>>
>> 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: Sunday, 26 March, 2023 6:28
>> >To: R meta
>> >Cc: Yuhang Hu
>> >Subject: [R-meta] Metafor modeling question
>> >
>> >Dear All,
>> >
>> >I have two questions about the following model:
>> >
>> >mm <- rma.mv(y ~ SES, V, random = list(~SES|study, ~SES|outcome), struct
>> =
>> >c("GEN","GEN"))
>> >
>> >fit to a dataset structurally similar to:
>> >
>> >study outcome SES
>> >1     A       1.5
>> >1     A       1.7
>> >1     B       1.5
>> >1     B       1.7
>> >2     A       2.1
>> >3     A       1.1
>> >3     A       2.3
>> >3     B       1.1
>> >3     B       2.3
>> >
>> >Question 1: I know, if I used: "random = list(~1 | study, ~1| outcome)",
>> >then 'study' and 'outcome' would be crossed, random effects. But are
>> >'study' and 'outcome' in model 'mm' still crossed,
>> >random effects? (In other words, can crossed random-effects have varying
>> >coefficients in them)
>>
>> Yes, they are crossed in this example.
>>
>> >Question 2: Is it possible (or useful) to plot residuals vs. fitted values
>> >for model 'mm'? If yes, how can we do that in metafor?
>>
>> fitted(mm) gives you the fitted values and resid(mm) gives you the
>> residuals, so you can just plot them against each other with
>> plot(fitted(mm), resid(mm)) or you can use the standardized residuals on
>> the y-axis with plot(fitted(mm), rstandard(mm)$z), since the residuals
>> themselves are heteroscedastic, but the standardized residuals are ...
>> well, standardized!
>>
>> Whether such a plot is useful is difficult to say in general, but I
>> suppose it can be used in the same manner as it is used with primary data.
>>
>> >Many thanks for your help,
>> >Yuhang