Skip to content

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

3 messages · James Pustejovsky, Mika Manninen

Mika Manninen
# 
Hello,

I am doing a meta-analysis looking at the effect of a teaching intervention
(versus control) on six types of motivation/behavioural regulation.
Theoretically and empirically these constructs form a continuum in which
the continuum neighbours are most strongly positively correlated and the
furthest from one another most negatively correlated.

I received help from Professors Wolfgang Viechtbauer and James Pustejovsky
already in July with the variance-covariance matrix computation. Thank you
again!

Now with all data collected, I have a few additional questions.regarding.

1. Warning message "positive definite"
2. Getting Q and I squared values from a multivariate meta-analysis output
3. and 4.Interpreting moderator analysis with categorical moderators
(Egger-type test, too for publication bias)
5. Suggestion for reporting forest plot and funnel plot for this type of
analysis (137 effects and 37 studies in the original datafile).

A subsample of data and relevant analyses can be done with the code
provided below. The aforementioned questions are blended into the text to
make answering them easier.
__________

#Install metafor package before proceeding with meta-analysis

library("metafor", lib.loc="~/R/R-4.0.0/library")

#Create the dataframe (without covariance matrix)

meta <- structure(list(Author = c("Abula et al., 2018", "Abula et al.,
2018",
                          "Abula et al., 2018", "Abula et al., 2018",
"Abula et al., 2018",
                          "Barkoukis et al., 2020", "Barkoukis et al.,
2020", "Barkoukis et al., 2020",
                          "Barkoukis et al., 2020", "Behzadnia et al.,
2017", "Behzadnia et al., 2017",
                          "Behzadnia et al., 2017", "Behzadnia et al.,
2017", "Chang et al., 2016",
                          "Chang et al., 2016", "Chang et al., 2016",
"Chang et al., 2016",
                          "Chang et al., 2016", "Cheon et al., 2010",
"Cheon et al., 2010",
                          "Cheon et al., 2010", "Cheon et al., 2010",
"Cheon et al., 2010",
                          "Cheon et al., 2012", "Cheon et al., 2015",
"Cheon et al., 2016",
                          "Edmunds et al., 2008", "Edmunds et al., 2008",
"Edmunds et al., 2008",
                          "Edmunds et al., 2008", "Edmunds et al., 2008",
"Edmunds et al., 2008",
                          "Fin et al., 2019", "Fin et al., 2019", "Fin et
al., 2019", "Fin et al., 2019",
                          "Fin et al., 2019", "Franco et al., 2017",
"Fransen etal., 2018",
                          "Gillison et al., 2013", "Gillison et al., 2013",
"Gillison et al., 2013",
                          "Gillison et al., 2013", "Gillison et al., 2013",
"Gonz?lez-Cutre et al., 2018",
                          "Gonz?lez-Cutre et al., 2018", "Gonz?lez-Cutre et
al., 2018",
                          "Gonz?lez-Cutre et al., 2018", "Gonz?lez-Cutre et
al., 2018",
                          "Gonz?lez-Cutre et al., 2018", "Ha et al., 2019",
"Ha et al., 2019"),
study = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4,
                         4, 4, 5, 5, 5, 5, 5, 6, 7, 8, 9, 9, 9, 9, 9, 9,
10, 10, 10, 10,
                         10, 11, 12, 13, 13, 13, 13, 13, 14, 14, 14, 14,
14, 14, 15, 15),
outcome = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,
                              15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25,
26, 27, 28, 29, 30,
                              31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41,
42, 43, 44, 45, 46,
                              47, 48, 49, 50, 51, 52),
motivation = c(1, 3, 4, 5, 6, 1, 3, 4, 5, 1, 3, 4, 5, 1, 3, 4, 5, 6, 1, 3,
4, 5, 6, 6, 6, 6, 1, 2,
                              3, 4, 5, 6, 1, 3, 4, 5, 6, 1, 1, 1, 3, 4, 5,
6, 1, 2, 3, 4, 5,
                              6, 1, 3),
g = c(0.528464772324917, 0.0589694070444635, 0.282008465683339,
      -0.0920405520247558, 0.0404267845152442, 0.419921031000001,
0.288692999072491,
      -0.129114905405604, 0.0153234677481859, 1.692130349854,
1.5655112070775,
      0.797525006036895, 0.202220243923423, 0.560898338638158,
0.120018930689234,
      0.320629566196145, -0.0255557448041566, -0.356052957510261,
0.836181226696161,
      0.945619432746331, -1.14695724875324, -0.316929747502631,
-0.444939499881025,
      -0.103827982870712, -0.124551198480087, -0.270481778614604,
0.193279463967818,
      -0.220165204655694, 0.328515459564001, -0.661405951944594,
-0.0588677781991213,
      0.381443142567869, 0.405279105271576, -0.0689560292055143,
-0.406763178319009,
      -1.28981574597883, -0.12834586771615, 0.493626759881918,
0.527629708371335,
      0.14237493801129, 0.181922479878477, -0.012193145480936,
-0.476411871637917,
      -0.239292055322439, 0.135428894231014, -0.140303735152328,
-0.125317079792252,
      -0.199539217658822, 0.347519633564979, -0.0569303505457916,
0.0331393420421724,
      0.127165941804818),
v = c(0.0126374422978132, 0.0136501896827314,
      0.0172092970523336, 0.0158304677686747, 0.015196984404496,
0.0117016875107911,
      0.0129905670084342, 0.0160782230295901, 0.0130349120410749,
0.166244292906422,
      0.172788985120758, 0.169998384821774, 0.147236242254395,
0.026050596807523,
      0.0280225511437046, 0.0353688036495172, 0.0324162425207384,
0.0316517208270417,
      0.125246076172345, 0.143299200337245, 0.182870602701095,
0.148364422975138,
      0.144468441612776, 0.00274503202263367, 0.00657178221924943,
      0.00378625733110702, 0.0567043412519983, 0.0698128695544358,
      0.064578154151857, 0.083501070217046, 0.0737350489050993,
0.0721490273295586,
      0.0536371594302375, 0.0589692454871346, 0.0750760396309753,
0.0825343404161793,
      0.0657659083731156, 0.0622812720119753, 0.0317147286821705,
0.0165661561004418,
      0.0186223247379937, 0.023286686815142, 0.0220852914411238,
0.0207926243477567,
      0.0402229479672577, 0.0494884645771515, 0.045350892784268,
0.0568053388291199,
      0.0531589046798887, 0.0504220144875667, 0.00506476243264833,
      0.00572614282481596),
setting = c(0, 0, 0, 0, 0, 0, 0, 0, 0,
            0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1,
            1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0),
scope = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 2, 2, 2, 2,
2, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 0, 1, 1, 1, 1, 1, 3, 3, 3,
3, 3, 3, 3, 3),
typei = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 0, 0),
typec = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0,
0, 0, 0, 0, 0),
length = c(2, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 1, 1,
1, 1, 1, 1, 1, 1),
age = c(1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1,
1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0),
oq = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2,
1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 0, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1,
1, 1, 1, 1),
mqi = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 2, 2, 2, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1,
1, 1, 1, 1, 1),
tq = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1,
1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1,
1, 1, 1, 1)),
row.names = c(NA,-52L), class = c("tbl_df", "tbl", "data.frame"))

View(meta)

corlist <- list(structure(c(1, NA, 0.96, 0.34, -0.33, -0.66, NA, NA, NA,
                 NA, NA, NA, 0.96, NA, 1, 0.55, -0.12, -0.5, 0.34, NA,
0.55, 1,
                 0.53, 0.05, -0.33, NA, -0.12, 0.53, 1, 0.72, -0.66, NA,
-0.5,
                 0.05, 0.72, 1), .Dim = c(6L, 6L)), structure(c(1, NA,
0.374,
                                                                0.034,
0.108, NA, NA, NA, NA, NA, NA, NA, 0.374, NA, 1, 0.261,
                                                                0.22, NA,
0.034, NA, 0.261, 1, 0.381, NA, 0.108, NA, 0.22, 0.381,
                                                                1, NA, NA,
NA, NA, NA, NA, NA), .Dim = c(6L, 6L)), structure(c(1,

                                                   NA, 0.96, 0.34, -0.33,
NA, NA, NA, NA, NA, NA, NA, 0.96, NA,

                                                   1, 0.55, -0.12, NA,
0.34, NA, 0.55, 1, 0.53, NA, -0.33, NA, -0.12,

                                                   0.53, 1, NA, NA, NA, NA,
NA, NA, NA), .Dim = c(6L, 6L)), structure(c(1,


                                            NA, 0.96, 0.34, -0.33, -0.66,
NA, NA, NA, NA, NA, NA, 0.96, NA,


                                            1, 0.55, -0.12, -0.5, 0.34, NA,
0.55, 1, 0.53, 0.05, -0.33, NA,


                                            -0.12, 0.53, 1, 0.72, -0.66,
NA, -0.5, 0.05, 0.72, 1), .Dim = c(6L,



                                6L)), structure(c(1, NA, 0.85, 0.35, 0.06,
-0.44, NA, NA, NA,



                                                  NA, NA, NA, 0.85, NA, 1,
0.59, 0.46, -0.14, 0.35, NA, 0.59, 1,



                                                  0.71, 0.27, 0.06, NA,
0.46, 0.71, 1, 0.52, -0.44, NA, -0.14,



                                                  0.27, 0.52, 1), .Dim =
c(6L, 6L)), structure(c(NA, NA, NA, NA,




                     NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA,




                     NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, 1



                                                  ), .Dim = c(6L, 6L)),
structure(c(NA, NA, NA, NA, NA, NA, NA,




        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,




        NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1), .Dim = c(6L,




                                                                     6L)),
structure(c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,





           NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,





           NA, NA, NA, NA, NA, NA, NA, NA, 1), .Dim = c(6L, 6L)),
structure(c(1,






  0.79, 0.84, 0.24, -0.2, -0.5, 0.79, 1, 0.93, 0.47, 0, -0.31,






  0.84, 0.93, 1, 0.57, -0.07, -0.52, 0.24, 0.47, 0.57, 1, 0.43,






  0.01, -0.2, 0, -0.07, 0.43, 1, 0.55, -0.5, -0.31, -0.52, 0.01,






  0.55, 1), .Dim = c(6L, 6L)), structure(c(1, NA, 0.96, 0.34, -0.33,






                                           -0.66, NA, NA, NA, NA, NA, NA,
0.96, NA, 1, 0.55, -0.12, -0.5,






                                           0.34, NA, 0.55, 1, 0.53, 0.05,
-0.33, NA, -0.12, 0.53, 1, 0.72,






                                           -0.66, NA, -0.5, 0.05, 0.72, 1),
.Dim = c(6L, 6L)), structure(c(1,







                               NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA,







                               NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA,







                               NA, NA, NA), .Dim = c(6L, 6L)),
structure(c(1, NA, NA, NA, NA,








 NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,








 NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA), .Dim = c(6L,








                                                                     6L)),
structure(c(1, NA, 0.96, 0.34, -0.33, -0.66, NA, NA, NA,









           NA, NA, NA, 0.96, NA, 1, 0.55, -0.12, -0.5, 0.34, NA, 0.55, 1,









           0.53, 0.05, -0.33, NA, -0.12, 0.53, 1, 0.72, -0.66, NA, -0.5,









           0.05, 0.72, 1), .Dim = c(6L, 6L)), structure(c(1, 0.76, 0.96,









                                                          0.34, -0.33,
-0.66, 0.76, 1, 0.62, 0.25, -0.09, -0.41, 0.96,









                                                          0.62, 1, 0.55,
-0.12, -0.5, 0.34, 0.25, 0.55, 1, 0.53, 0.05,









                                                          -0.33, -0.09,
-0.12, 0.53, 1, 0.72, -0.66, -0.41, -0.5, 0.05,









                                                          0.72, 1), .Dim =
c(6L, 6L)), structure(c(1, NA, 0.79, NA, NA,










                       NA, NA, NA, NA, NA, NA, NA, 0.79, NA, 1, NA, NA, NA,
NA, NA,










                       NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA









                                                          ), .Dim = c(6L,
6L)))
# A correlation matrix for the levels of motivation actually observed -
omitting NAs (by W.Viechtbauer, July2020)
mots <- split(meta$motivation, meta$study)
corlist <- mapply(function(rmat, mots) rmat[mots,mots], corlist, mots)
corlist

#Creating the covariance matrix - V-matrix function (w. Viechtbauer, July
2020)

Vfun <- function(vi,rmat) {
  if (length(vi) == 1L) {
    matrix(vi)
  } else {
    S <- diag(sqrt(vi))
    S %*% rmat %*% S
  }
}

V <- mapply(Vfun, split(meta$v, meta$study), corlist)
V

#Here I sometimes get the following warning message: "V appears to be not
positive definite". 1. Should this be ignored?

#The multivariate model (all clear so far, except for positive definite
warning):

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

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

#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"? )

#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

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

###Thank you in advance!###

Mika
James Pustejovsky
# 
Hi Mika,

Comments below.

James
On Thu, Sep 24, 2020 at 12:45 PM Mika Manninen <mixu89 at gmail.com> wrote:

            

      
      
      
    

        
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

            

      
      
      
    

        
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?

            

      
      
      
    

        
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.

            

      
      
      
    

        
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.

            

      
      
      
    

        
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.

  
  


  


      

    

    
5 days later

    

      
      

        


  

        

      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 at gmail.com)
kirjoitti: