[R-meta] Queries re: rma.mv

[Viechtbauer, Wolfgang (NP)]

Hi Adelina,

In multilevel models, one has to distinguish the between- and the within-study relationship between some predictor and the outcome. I don't have access to your data, but possibly those trends differ in your data. So, it might appear visually that the trend is increasing, but not when analyzed using a multilevel model.

The weight structure is more complex in models with moderators, so it's not just a matter of computing row-sum weights to get the same results. In general, there is an entire weight matrix. On the page that I linked to, I only discuss the special case of an "intercept-only model". James Pustejovsky has a blog entry where he goes further into discussing weights also in models with moderators:

https://www.jepusto.com/weighting-in-multivariate-meta-analysis/

I am not aware of a general way of reducing the weight matrix down to a vector of weights, such that a standard meta-analysis model will give the same results as the more complex multilevel/multivariate model from which the weight matrix is derived. So, since it is not generally possible (as far as I know) to compute 'estimate-specific weights', forest and bubble plots just use the diagonal of the weight matrix as a rough indication of how much weight is assigned to individual points (but you are of course free to specify other point sizes).

Best,
Wolfgang

>-----Original Message-----
>From: Adelina Artenie [mailto:adelina.artenie using bristol.ac.uk]
>Sent: Friday, 13 May, 2022 11:38
>To: Viechtbauer, Wolfgang (NP); r-sig-meta-analysis using r-project.org
>Subject: Re: Queries re: rma.mv
>
>Hi Wolfgang,
>
>Thank you for your detailed response. It has been very helpful. I am still
>pondering one thing though about the weights, if I could have your advice please:
>
>I have re-calculated the weights as you suggested (shown in my code below) and if
>I run a bubble plot using regplot and use these weights as the size of the
>bubbles, it would appear that, indeed, the correlated data are down-weighted.
>
>The issue is the following (note: my moderator is a measure of time): looking at
>the bubble plot, one can tell visually that the trend is clearly increasing over
>time. However, the coefficient from rma.mv and the best fit line in regplot show
>a decreasing trend. The only way this would occur is if the weights of the
>correlated data used in the model are not, in fact, down-weighted but up-weighted
>(the correlated data have a decreasing trend over time). In other words, it seems
>as if rma.mv is not using the row-sum weights but the weights that are based on
>the diagonal of the matrix (or some other weights..).
>
>Relatedly, I also wonder why, if the default weights shown in the forest plot or
>the bubble/regplot, are not the ones used by the rma.mv model, why are they the
>default weights plotted? This too makes me think that these are the ones used
>internally in the model estimation by rma.mv
>
>Apologies for the long e-mail. This package has been incredibly useful so far.
>I've just hit a point where I do not fully understand how rma.mv works and how to
>make sense of my results.
>
>Many thanks for your help
>Adelina
>
>My code:
>
>Inci_trends$prs_yrs_100 = Inci_trends$prs_yrs/100
>trend_dat_Inci <- escalc(measure="IRLN", xi=cases, ti=prs_yrs_100,
>data=Inci_trends)
>summary.escalc(trend_dat_Inci)
>
>multi_level_m_Inci <- rma.mv(yi = yi,
>                            V = vi,
>                            slab = author,
>                            data = trend_dat_Inci,
>                            random = ~ 1 |slab/ID,
>                            test = "z",
>                            method = "REML",
>                            mods = ~ midpoint )
>summary(multi_level_m_Inci)
>round(exp(coef(summary(multi_level_m_Inci))[-1,c("estimate", "ci.lb", "ci.ub")]),
>2)
>
>w_reg_Inci <- weights(multi_level_m_Inci)
>w_m_Inci <- weights(multi_level_m_Inci, type = "matrix")
>
>w_m_Inci_adjw <- rowSums(w_m_Inci) / sum(w_m_Inci) * 100
>
>regplot(multi_level_m_Inci, label = F, transf=exp,
>        ylab = "Rate",
>        xlab = "midpoint",
>        col="grey",
>        bg="red",
>        labsize = 0.5,
>        legend = F,
>        psize=c(w_m_Inci_adjw)
>)
>
>From: Viechtbauer, Wolfgang (SP) <wolfgang.viechtbauer using maastrichtuniversity.nl>
>Date: Wednesday, 20 April 2022 at 19:42
>To: Adelina Artenie <adelina.artenie using bristol.ac.uk>, r-sig-meta-analysis using r-
>project.org <r-sig-meta-analysis using r-project.org>
>Subject: RE: Queries re: rma.mv
>Dear Adelina,
>
>See below for my responses.
>
>Best,
>Wolfgang
>
>>-----Original Message-----
>>From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces using r-project.org] On
>>Behalf Of Adelina Artenie
>>Sent: Wednesday, 20 April, 2022 19:15
>>To: r-sig-meta-analysis using r-project.org
>>Subject: [R-meta] Queries re: rma.mv
>>
>>Hello
>>
>>I am trying to do a meta-regression to explore whether incidence rates are
>>changing over time. Some of the estimates are correlated so I am using the
>rma.mv
>>function. The moderator is a measure of time.
>>
>>I have two questions please:
>>
>>  1.  I am trying to understand why is it that the correlated estimates are
>>systematically up-weighted. I would have thought that any correlated data would
>>be down-weighted or, at the very least, have a similar weight as the
>uncorrelated
>>data. I appreciate that the weight calculation in this case is complex
>>(https://www.metafor-project.org/doku.php/tips:weights_in_rma.mv_models ). I
>>wonder if, conceptually, this would be expected? In my mock example, the
>>correlated estimates have weights of 7-10% whereas the uncorrelated ones have
>>weights of 4%. I see a similar trend in my own data, and I find the difference
>in
>>the weights to be quite considerable.
>
>I assume you are referring to the weights shown in the forest plot. These are
>just based on the diagonal of the weight matrix. As you note, the weighting in
>such models is more complex, so looking at these weights doesn't really allow you
>to draw such a conclusion. In fact, generally, the effect will be as you expect,
>namely that correlated estimates will be downweighted. Consider the
>dat.konstantopoulos2011 dataset and let's make all sampling variances equal to
>each other to see this more clearly. Also, we will compute the row-sum weights,
>since these reflect more directly the weight each estimate is getting in this
>model.
>
>dat <- dat.konstantopoulos2011
>
>dat$vi <- median(dat$vi)
>res <- rma.mv(yi, vi, random = ~ 1 | district/school, data=dat)
>res
>
>wi <- weights(res, type="rowsum")
>
>data.frame(k = c(table(dat$district)),
>           weight = tapply(wi, dat$district, mean))
>
>So, as you can see, in clusters with more estimates, each individual estimate
>gets less weight.
>
>Things will be more complicated when sampling variances differ and also when
>including moderators in the model.
>
>>  2.  Could I confirm with you that my interpretation of the coefficient for the
>>moderator is correct? In this case, the estimate is 0.0846, which represents the
>>mean/median log IRR (incidence rate ratio) for a one-unit increase in the
>>moderator. Therefore, the IRR is 1.088: for each unit increase in time, the
>>incidence rate increases on average by 8.8% (ignoring the confidence intervals
>>for now).
>
>Correct.
>
>>Many thanks in advance,
>>Adelina
>>
>>Adelina Artenie, MSc, PhD
>>CIHR Postdoctoral Research Fellow
>>Population Health Sciences
>>Bristol Medical School
>>University of Bristol
>>
>>Code
>>***
>>
>>ID <- c(1:18)
>>author <- c("AA", "AA", "BB", "CC", "DD", "EE", "EE", "EE", "FF", "GG", "HH",
>>"II", "JJ",
>>            "KK", "LL", "MM", "NN", "OO")
>>cases <- c(8:25)
>>prs_yrs <-
>>c(200,150,3000,400,500,100,400,300,1200,456,177,296,664,123,123,432,67,3045)
>>moderator<-c(1:18)
>>
>>mydata<-data.frame(ID,author,cases,prs_yrs,moderator)
>>print(mydata)
>>
>>my_example <- escalc(measure="IRLN", xi=cases, ti=prs_yrs, data=mydata)
>>my_example
>>my_meta_reg <- rma.mv(yi = yi,
>>                            V = vi,
>>                            slab = author,
>>                            data = my_example,
>>                            random = ~ 1 | author/ID,
>>                            test = "z",
>>                            method = "REML",
>>                            mods = ~ moderator)
>>
>>summary(my_meta_reg)
>>regplot(my_meta_reg, label = T, transf=exp,
>>        ylab = "Incidence rate",
>>        xlab = "Moderator",
>>        col="black",
>>        plim = c(1,NA),
>>        labsize = 0.5,
>>        legend = F)
>>forest(my_meta_reg, showweights=TRUE)