Skip to content

[R-meta] Help specifying a three-level model

3 messages · Frank Bosco, Wolfgang Viechtbauer, _m@i@@ccou@t m@iii@g oii ir@@kbosco@com

Original
Frank Bosco
# 
Greetings to all,

I seek help specifying a three-level model. I suspect I need just syntax 
help, but it might be more complicated. I?ve examined some previous 
examples, but none seem to address my specific need. Hoping there?s a 
three-level guru out there willing to help. Any help would be greatly 
appreciated.

Here?s the situation: I have 94 effect sizes pertaining to a particular 
bivariate relation. The 94 effect sizes are nested within 56 independent 
samples. In turn, the 56 independent samples are nested within 9 
countries. The data file is at the URL in the script below.

For each country, I have values for three Hofstede cultural dimensions 
(these are continuous variables): Individualism-collectivism, power 
distance, and masculinity/femininity. For each independent sample, I 
have a control variable: publication year.

I would like to ascertain the incremental variance (which, I presume, 
would be a pseudo R^2) for the addition of the three cultural factors 
(i.e., all three cultural factors entered simultaneously).

At any rate, as shown below, I?m somewhat lost at the moment (I'm trying 
to address a reviewer comment). I paste below the variable names and my 
script skeletons.

Thanks in advance for any help provided.

Frank

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Variable names
yi
vi
SampleID
Country
PubYear
HOF_powerDist
HOF_indivCol
HOF_masFem

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

#-- Script attempt (Fit0 contains only the control variable [PubYear]; 
Fit1 contains the control variable and the three cultural factors [HOF_*])

library(metafor)
dat <- read.csv("http://www.frankbosco.com/test.csv")

Fit0 <- rma.mv(yi, vi,
 ????????????? mods = ~ PubYear + ??? -1,???? #-- (not sure whether I 
need ?-1? here)
 ????????????? random = ???,? #-- I've tried several nesting approaches 
- some cause an error
 ????????????? data = dat,
 ????????????? struct = ??????? #-- (Not sure what I need here)
)

Fit1 <- rma.mv(yi, vi,
 ????????????? mods = ~ PubYear + HOF_powerDist + HOF_indivCol + 
HOF_masFem -1,????? #-- (not sure whether I need ?-1? here)
 ????????????? random = ???,? #-- I've tried several nesting approaches 
- some cause an error
 ????????????? data = dat,
 ????????????? struct = ???????? #-- (Not sure what I need here)
)

anova(Fit0, Fit1)

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-


Frank Bosco, Ph.D.
Director, metaBUS.org
Associate Professor
School of Business
Department of Management

Virginia Commonwealth University
Snead Hall
301 West Main Street, Room B4151
Richmond, Virginia 23284
804 828-3197 ? Fax: 804 828-1602
business.vcu.edu
7 days later
7 days later
_m@i@@ccou@t m@iii@g oii ir@@kbosco@com
# 
Thanks, Wolfgang!

This is just what we needed.

Others looking for existing examples of 3-, 4-, and 5-level models: You might consider the following review of multi-level metas:

https://doi.org/10.3758/s13428-020-01373-9

(In their supplemental materials, you'll find some good sources, including a reference to a 4-level model reported in Psychological Bulletin [http://dx.doi.org/10.1037/bul0000142] whose supplemental materials, in turn, contain R script using metafor.)

Cheers,

Frank

-----Original Message-----
From: Viechtbauer, Wolfgang (SP) <wolfgang.viechtbauer at maastrichtuniversity.nl> 
Sent: Thursday, February 18, 2021 4:00 AM
To: meta at frankbosco.com; r-sig-meta-analysis at r-project.org
Subject: RE: [R-meta] Help specifying a three-level model

Dear Frank,
random = ~ 1 | Country / SampleID / EstimateID

would be a starting point, where 'EstimateID' is just 

dat$EstimateID <- 1:nrow(dat)

Argument 'struct' is irrelevant here.

Don't use -1 in 'mods'. See here:

https://www.metafor-project.org/doku.php/tips:models_with_or_without_intercept

With anova(Fit0, Fit1) you can do a LRT comparing the two models. For comparing models with different fixed effects (as you are doing), you would have to use ML estimation (i.e., method="ML") instead of REML estimation (which is the default). But you can also just do a Wald-type test instead with:

anova(Fit1, btt=3:5)

which should give you the omnibus test of the coefficients for HOF_powerDist, HOF_indivCol, HOF_masFem.

To get pseudo R^2 values, you can do:

(Fit0$sigma2 - Fit1$sigma2) / Fit0$sigma2

to estimate the proportional reduction in each variance component as the additional fixed effects are added to the model or:

(sum(Fit0$sigma2) - sum(Fit1$sigma2)) / sum(Fit0$sigma2)

to get the proportional reduction of the total amount of heterogeneity.

Given that the moderators are measured at the country level, the largest reductions (if any) should be observable for the country-level variance component, although in practice things don't divide that neatly.

Best,
Wolfgang