[R-meta] Multivariate meta-analysis - moderator analysis and tau squared

Mika Manninen m|xu89 @end|ng |rom gm@||@com
Tue Sep 29 20:57:06 CEST 2020


Thank you very much for your detailed answer!

I took some time to learn (first time doing the multivariate analyses) and
I would have 4 quick follow-up questions.

1. Follow-up on moderators (especially categorical moderators with more
than two levels).
2. About calculating the I squared and the need for Qs
3. About influential study diagnostics (reestimate = FALSE)
4. About calculating the modified precision estimate for the Egger's type
test based on

The data and main analysis can be generated with the code from my original
message (sorry about the correlations being all over the place).

The follow-up questions and related code are presented below

###1. About moderators

#The moderators were selected a priori and the purpose was to all of the
moderators for each level of motivation (increased risk noted!)

#In my understanding here the latter half of the output is the difference
to the corresponding motivation level on the first half (categorical
moderator level 0)
setting_res <- rma.mv(g, V, mods = ~ factor(motivation)*I(setting) -1,
random = ~ factor(motivation) | study, struct="UN", data=meta)

#Is there a simple way to get the mean effect and its CIs for all the
motivations levels and setting levels? Usually the -1 addition
#provides this, but is it that with the multivariate model, the level 1
moderator effects need to be subtracted/added from the level 0 to have
#mean effects for all the levels of motivation with different moderators?

#Also, there are a few moderators with more than two levels (example with
three). Would the following be the correct way to run the analysis in this
#Again, would the level 1 and level 2 estimates be the differences from the
level 0 estimates? So, the mean effects sizes for all the levels can be
gained by subtraction/addition to and from the level 0

length_res <- rma.mv(g, V, mods = ~ factor(motivation)*factor(length)-1,
random = ~ factor(motivation) | study, struct="UN", data=meta, method =

### 2. About calculating the I squared
#Which of these three methods would you recommend using? Also, if I am
reporting taus and I squared, would not this suffice and Qs could be left

#Constructing a block diagonal of the variance-covariance matrix (V from
the original message).

c <- bldiag(V)

##I2 computations from metafor-project.org

#1st option (solving W, X, and P)

W <- solve(c)
X <- model.matrix(res)
P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W
100 * res$tau2 / (res$tau2 + (res$k-res$p)/sum(diag(P)))

#2nd option. P comes from the code above (not sure if I computed the W
correct by assigning the bldiag function - constructs a block diagonal
matrix) )

c(100 * res$tau2[1] / (res$tau2[1] + (sum(meta$motivation ==
1)-1)/sum(diag(P)[meta$motivation == 1])),
  100 * res$tau2[2] / (res$tau2[2] + (sum(meta$motivation ==
2)-1)/sum(diag(P)[meta$motivation == 2])),
  100 * res$tau2[3] / (res$tau2[3] + (sum(meta$motivation ==
3)-1)/sum(diag(P)[meta$motivation == 3])),
  100 * res$tau2[4] / (res$tau2[4] + (sum(meta$motivation ==
4)-1)/sum(diag(P)[meta$motivation == 4])),
  100 * res$tau2[5] / (res$tau2[5] + (sum(meta$motivation ==
5)-1)/sum(diag(P)[meta$motivation == 5])),
  100 * res$tau2[6] / (res$tau2[6] + (sum(meta$motivation ==
6)-1)/sum(diag(P)[meta$motivation == 6])))

#3rd option (Jackson 2012) ( %s lower than with the previous two approaches)

res.R <- rma.mv(g, V, mods = ~ factor(motivation) - 1, random = ~
factor(motivation) | study, struct="UN", data=meta)
res.F <- rma.mv(g, V, mods = ~ factor(motivation) - 1, data=meta)

c(100 * (vcov(res.R)[1,1] - vcov(res.F)[1,1]) / vcov(res.R)[1,1],
  100 * (vcov(res.R)[2,2] - vcov(res.F)[2,2]) / vcov(res.R)[2,2],
  100 * (vcov(res.R)[3,3] - vcov(res.F)[3,3]) / vcov(res.R)[3,3],
  100 * (vcov(res.R)[4,4] - vcov(res.F)[4,4]) / vcov(res.R)[4,4],
  100 * (vcov(res.R)[5,5] - vcov(res.F)[5,5]) / vcov(res.R)[5,5],
  100 * (vcov(res.R)[6,6] - vcov(res.F)[6,6]) / vcov(res.R)[6,6])

### 3. About influential studies
#Is setting the reestimate to FALSE on cooks.distance and dfbetas okay?
Otherwise, the diagnostics resulted in NAs.

#Cooks distance on each outcome (Mahalanobis distance)

cd <- cooks.distance(res, reestimate = FALSE)
plot(cd, type="o", pch=19)

#Cooks distance clustered by study

cd <- cooks.distance(res, ncpus = 1, cluster = meta$study, reestimate =
plot(cd, type="o", pch=19)

#Dfbetas (Change in sds)

dfb <- dfbetas(res, ncpus = 1, reestimate = FALSE)
plot(dfb, type="o", pch=19)


hats <- hatvalues(res, ncpus = 1)
plot(hats, type="o", pch=19)

### 4. About calculating a modified precision estimate for Egger's type
#Could mprc (after running the code) be used in the Egger's test as a
moderator/modified precision estimate to keep up with the articles

weight <- weights(res, type = "diagonal")

mprec <- sqrt(weight)

###Thank you very much in advance!###

to 24. syysk. 2020 klo 19.21 James Pustejovsky (jepusto using gmail.com)

> Hi Mika,
> Comments below.
> James
> On Thu, Sep 24, 2020 at 12:45 PM Mika Manninen <mixu89 using gmail.com> wrote:
>> #Here I sometimes get the following warning message: "V appears to be not
>> positive definite". 1. Should this be ignored?
> This means that one of your input correlation matrices is not positive
> definite, i.e., it is not a valid correlation matrix. It could be due to a
> typo or rounding of the entries. You can find the offending matrix using
> the following
> isPosDef <- function(x) all(eigen(x)$values > 0)
> sapply(corlist, isPosDef)
> Although it's probably a minor issue, I would still suggest correcting
> before proceeding.
> #Each level of motivation has its own variance component in the output
>> (tau). 2. How could I obtain the Q-values (and consequently I squared) to
>> report?
> I would recommend reporting the tau estimates directly. Because they are
> model parameters, they are more meaningful and more directly interpretable
> than Q or I-squared.
> As far as Q, there are a number of different ways to define Q in
> multivariate models (or generally, models with multiple variance
> components), some of which are global rather than being specific to the
> dimension of the outcome. Do you want to report Q as a global description
> of excess heterogeneity, or do you want something specific to the outcome
> dimension?
>> #If I now wanted to do a moderator analysis with each level of motivation
>> using the categorical variable "setting"
>> #to see if there is significant difference between the two settings for
>> motivation 1,2,3 etc. 3. Would the following be correct?
>> res <-  rma.mv(g, V, mods = ~ factor(motivation)+
>> factor(motivation):I(setting) -1, random = ~ factor(motivation) | study,
>> struct="UN", data=meta)
>> res
>> res <- rma.mv(g, V, mods = ~ factor(motivation)*I(setting) -1, random = ~
>> factor(motivation) | study, struct="UN", data=meta)
>> res
>> #I am not exactly sure what I should be interpreting. I would like to know
>> if setting moderates the effect of the intervention on motivation levels
>> 1,2,3,4,5,6 (or can I do that - is the comparison done "overall"? )
>> Is your question about setting moderates for *any* motivation level? Or
> whether it moderates for  *specific* motivation levels?
> If it is the former, then you can do a likelihood ratio test comparing the
> model with moderator to the model without. You would, however, have to
> switch to using ML rather than REML estimation for the variance components.
> Syntax as follows:
> res_ML <- rma.mv(g, V, mods = ~ factor(motivation) - 1,
>                  random = ~ factor(motivation) | study,
>                  struct="UN", data = meta,
>                  method = "ML")
> mod_ML <- update(res_ML, mods = ~ factor(motivation)*I(setting) -1)
> anova(res_ML, mod_ML)
> If it is the latter, then the first set of syntax gives you coefficient
> estimates for the difference in average effect size between setting = 1
> versus setting = 0, for each distinct level of motivation, so you can
> interpret those coefficients (CIs, t-tests) directly. If you're looking at
> the t-tests for all six levels, it would be prudent to use a correction for
> multiple comparisons.
>> #Same goes for the Egger's type test, 4. how should the results of the
>> moderator analyses be interpreted if I were to use the variance estimate
>> or
>> its inverse (v) as the predictor to test for publication bias/funnel plot
>> asymmetry
> If these effect sizes are standardized mean differences, then you'll need
> to use a modified measure of precision rather than the variance or standard
> error of the effect sizes. Details here:
> https://www.jepusto.com/publication/testing-for-funnel-plot-asymmetry-of-smds/
> There are at least two ways to implement Egger's test in this setting. One
> would be to simply add the modified measure of precision as a predictor. A
> significant slope coefficient would be indicative of small-study effects.
> Alternately, you could interact the predictor with the levels of motivation
> and then report the likelihood ratio test, as with the previous question.
> It is hard to say which approach is more powerful generally. Perhaps others
> on the listserv have insights.
>> #5. Any other advice to report the results (e.g., how to display the
>> forest
>> plot and is there any sense to do a funnel plot?)
> You could make a forest plot with the points and whiskers in different
> color shades corresponding to the levels of motivation. Alternately, make
> separate forest plots per level of motivation. Similarly, a funnel plot
> with points in different colors corresponding to the levels of motiviation,
> or make separate funnels per level of motivation.

	[[alternative HTML version deleted]]

More information about the R-sig-meta-analysis mailing list