[R-meta] Queries re: rma.mv
James Pustejovsky
jepu@to @end|ng |rom gm@||@com
Fri May 13 16:57:37 CEST 2022
Hi Adelina,
I agree with Wolfgang that looking at the between-study and within-study
variation in the predictor would probably be a helpful step for diagnosing
what's going on here. To do so, try calculating group means of the
predictor and a group-centered version of the predictor, for instance as
follows:
trend_dat_Inci <- escalc(measure="IRLN", xi=cases, ti=prs_yrs_100,
data=Inci_trends)
trend_dat_Inci_centered <-
trend_dat_Inci %>%
group_by(slab) %>%
mutate(
midpoint_b = mean(midpoint),
midpoint_w = midpoint - midpoint_b
)
Then include both predictors in the meta-regression model:
rma.mv(yi = yi,
V = vi,
slab = author,
data = trend_dat_Inci,
random = ~ 1 |slab/ID,
test = "z",
method = "REML",
mods = ~ midpoint_b + midpoint_w)
In your original meta-regression model, the coefficient on midpoint is
actually a weighted average of the coefficients on midpoint_b and
midpoint_w. So, if there is a big difference (especially a sign difference)
in the coefficients on midpoint_b and midpoint_w, then that would help to
explain why the visual appearance of the bubble plot doesn't match up with
the meta-regression result.
James
On Fri, May 13, 2022 at 5:33 AM Viechtbauer, Wolfgang (NP) <
wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
> 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)
>
> _______________________________________________
> R-sig-meta-analysis mailing list @ R-sig-meta-analysis using r-project.org
> To manage your subscription to this mailing list, go to:
> https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis
>
[[alternative HTML version deleted]]
More information about the R-sig-meta-analysis
mailing list