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

[Mika Manninen]

James,

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
https://www.jepusto.com/publication/testing-for-funnel-plot-asymmetry-of-smds/

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)
setting_res

#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
case?
#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 =
"ML")
length_res

### 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
unreported?

#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.R
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)
cd
plot(cd, type="o", pch=19)

#Cooks distance clustered by study

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

#Dfbetas (Change in sds)

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

#Hatvalues

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


### 4. About calculating a modified precision estimate for Egger's type
test
#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
recommendations?

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

mprec <- sqrt(weight)
View(mprec)

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



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

> 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]]