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 at gmail.com)
kirjoitti:
Hi Mika,
Comments below.
James
On Thu, Sep 24, 2020 at 12:45 PM Mika Manninen <mixu89 at 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.