Edit: changed code to include test = "t"
I'm hoping to better understand how the updated dev version of Metafor 2.5-101 will help me to adjust my degrees of freedom in a multi-level model to provide some protection against type 1 error.
My understanding of this has come from the Nakagawa preprint "Methods for testing publication bias in ecological and evolutionary meta-analyses" https://ecoevorxiv.org/k7pmz/ and their "Supplemental_Impleentation_Example.Rmd" file, following along with their line 133-142:
Before moving on to some useful corrections, users should be aware that the most up-to-date version of
metafor
(version 2.5-101) does now provide users with some protection against Type I errors. Instead of using the number of effect sizes in the calculation of the degrees of freedom we can instead make use of the total numbers of papers instead. We show in our simulations that a "papers-1" degrees of freedom can be fairly good. This can be implemented as follows after installing the development version ofmetafor
(see "R Packages Required" above):
mod_multilevel_pdf = rma.mv(yi = yi, V = vi, mods = ~1,
random=list(~1|study_id,~1|obs),
data=data, test="t", dfs = "contain")
summary(mod_multilevel_pdf)
We can see here that the df for the model has changes from 149 to 29, and the p-value has been adjusted accordingly.
So my understanding here is that the model now shows df as 29 (the original no. of papers (30) -1, instead of the no. of papers x no. of effects (30 papers with 5 effects each (150) -1)
Adapting this to my code, where I have n=18 papers and total of n=24 effects, I would expect using the above code would adjust my df to 17 (the original no. of papers (18) -1), however I still have df as 23 (total no. of effects (24) -1).
The output using the df code:
mod_multilevel_pdf = rma.mv(yi = yi, V = vi, mods = ~1,
random=list(~1|study_id,~1|es_id),
data=dat, test="t", dfs = "contain")
summary(mod_multilevel_pdf)
Is:
Multivariate Meta-Analysis Model (k = 24; method: REML)
logLik Deviance AIC BIC AICc
-30.2270 60.4540 66.4540 69.8604 67.7171
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.6783 0.8236 18 no study_id
sigma^2.2 0.1416 0.3763 24 no es_id
Test for Heterogeneity:
Q(df = 23) = 167.2145, p-val < .0001
Model Results:
estimate se tval df pval ci.lb ci.ub
-0.3508 0.2219 -1.5809 17 0.1323 -0.8190 0.1174
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Quite stumped on this one! Any help would be majorly appreciated.
You neither have df=17 nor df=23, since you did not specify that you want a t-test. With test="t", dfs = "contain"
, you will get the expected t-test with df=17.