[R-meta] Queries re: rma.mv

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 at bristol.ac.uk]
Sent: Friday, 13 May, 2022 11:38
To: Viechtbauer, Wolfgang (NP); r-sig-meta-analysis at 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 at maastrichtuniversity.nl>

Date: Wednesday, 20 April 2022 at 19:42
To: Adelina Artenie <adelina.artenie at bristol.ac.uk>,

r-sig-meta-analysis at r-

project.org <r-sig-meta-analysis at 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 at r-project.org] On

Behalf Of Adelina Artenie
Sent: Wednesday, 20 April, 2022 19:15
To: r-sig-meta-analysis at 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)