[R-meta] Queries re: rma.mv

[Adelina Artenie]

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)

	[[alternative HTML version deleted]]