[R-meta] Calculating conditional and marginal r-squared values for rma.mv models in metafor
Viechtbauer Wolfgang (SP)
wolfgang.viechtbauer at maastrichtuniversity.nl
Wed Sep 27 23:14:51 CEST 2017
Dear Danielle,
First of all, I think you should use:
cover.hedges$id <- 1:nrow(cover.hedges)
cover.rma <- rma.mv(yi, vi, random=~1|ACC/id, method="REML", data=cover.hedges)
See also: http://www.metafor-project.org/doku.php/analyses:konstantopoulos2011#a_common_mistake_in_the_three-level_model
I also don't understand what you mean by "We report the fsn for the model without coefficients, as well as the AICc for the top models with coefficients (within 2AICc units of the best model)." What do you mean by "fsn" here?
As for a pseudo-R^2, one can easily extend the way of computing R^2 for 'standard' mixed-effects meta-regression models to this type of model. That is, let 'res0' be the model without coefficients and 'res1' the model with (both fitted as shown above). Then:
(res0$sigma2 - res1$sigma2) / res0$sigma2
gives the proportional reduction in each variance component, so R^2 for the study level and the estimate level. A 'total' R^2 would be:
(sum(res0$sigma2) - sum(res1$sigma2)) / sum(res0$sigma2)
Best,
Wolfgang
-----Original Message-----
From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces at r-project.org] On Behalf Of Danielle Claar
Sent: Sunday, 24 September, 2017 20:28
To: r-sig-meta-analysis at r-project.org
Subject: [R-meta] Calculating conditional and marginal r-squared values for rma.mv models in metafor
Hello,
Thanks to Wolfgang for pointing me in the direction of this list-serve!
I am currently revising a manuscript that uses metafor to investigate how
sea surface temperature stress affects coral bleaching and mortality during
El Nino events. We are using linear mixed effects models to take into
account within-study variance, with study (ACC) as a random factor in our
model. To do this, we are using rma.mv. In our simplest model, this looks
like this:
> cover.rma <- rma.mv(yi, vi, random=~1|ACC, method="REML",
data=cover.hedges)
> cover.rma
Multivariate Meta-Analysis Model (k = 192; method: REML)
Variance Components:
estim sqrt nlvls fixed factor
sigma^2 2.3650 1.5378 29 no ACC
Test for Heterogeneity:
Q(df = 191) = 1355.9701, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
-1.5225 0.2946 -5.1682 <.0001 -2.1000 -0.9451 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
We have also used glmulti (as explained here
http://www.metafor-project.org/doku.php/tips:model_selection_with_glmulti),
to conduct coefficient selection, which looks like this (output truncated
for brevity).
> cov <- glmulti(yi~DHWtodate+Timelagtodate+tmean+tvar,
random=list("~1|ACC"), # First try the full model with all 4 potential
moderators
+ data=cover.hedges, fitfunction = rma.glmulti, level=2,
crit="aicc") # Level 2 looks at the interactions between all moderators.
AICc allows for correction of a small N
Initialization...
TASK: Exhaustive screening of candidate set.
Fitting...
After 50 models:
Best model: yi~1+tmean:DHWtodate+tvar:tmean Crit= 780.334181647103 Mean
crit= 806.821778458654 .
After 1050 models:
Best model: yi~1+DHWtodate+tmean+Timelagtodate:DHWtodate+tmean:DHWtodate
Crit= 773.219300469166
Mean crit= 776.66700421609
Completed.
> print(cov) # View summary of cov model
glmulti.analysis
Method: h / Fitting: rma.glmulti / IC used: aicc
Level: 2 / Marginality: FALSE
From 100 models:
Best IC: 773.219300469166
Best model:
[1] "yi ~ 1 + DHWtodate + tmean + Timelagtodate:DHWtodate + tmean:DHWtodate"
Evidence weight: 0.0488927941262297
Worst IC: 777.81840821166
8 models within 2 IC units.
90 models to reach 95% of evidence weight.
> tmp <- weightable(cov) # Save tmp weightable for cov2 tmp <-
> tmp[tmp$aicc <= min(tmp$aicc) + 10000,] tmp # print AICc and weights for
all models
model aicc weights
1 yi ~ 1 + DHWtodate + tmean + Timelagtodate:DHWtodate + tmean:DHWtodate
773.2193 0.048892794
2 yi ~ 1 + DHWtodate + tmean + tmean:DHWtodate + tvar:DHWtodate +
tvar:Timelagtodate 773.9826 0.033380831
3 yi ~ 1 + DHWtodate + tvar + tmean:DHWtodate + tvar:tmean
774.2527 0.029163932
4 yi ~ 1 + DHWtodate + tvar + tmean:DHWtodate
774.8099 0.022072363
5 yi ~ 1 + DHWtodate + tmean + tvar + Timelagtodate:DHWtodate +
tmean:DHWtodate 774.9541 0.020537024
6 yi ~ 1 + DHWtodate + tmean + Timelagtodate:DHWtodate + tmean:DHWtodate +
tvar:tmean 775.0963 0.019127989
7 yi ~ 1 + DHWtodate + tmean + Timelagtodate:DHWtodate + tmean:DHWtodate +
tvar:DHWtodate 775.1307 0.018801511
8 yi ~ 1 + DHWtodate + tmean + Timelagtodate:DHWtodate + tmean:DHWtodate +
tvar:Timelagtodate 775.1505 0.018616442
> cov.summ <- summary(cov at objects[[1]]) # create summary object for data
> cov.summ
Multivariate Meta-Analysis Model (k = 192; method: ML)
logLik Deviance AIC BIC AICc
-380.3826 588.5923 772.7652 792.3102 773.2193
Variance Components:
estim sqrt nlvls fixed factor
sigma^2 2.2116 1.4872 29 no ACC
Test for Residual Heterogeneity:
QE(df = 187) = 1117.8900, p-val < .0001
Test of Moderators (coefficient(s) 2,3,4,5):
QM(df = 4) = 82.9952, p-val < .0001
Model Results:
estimate se intrcpt 14.4108
3.0184 DHWtodate -0.6322 0.2191 tmean
-0.5737 0.1103 DHWtodate:Timelagtodate -0.0001 0.0000 DHWtodate:tmean
0.0228 0.0086
zval pval intrcpt 4.7743
<.0001 DHWtodate -2.8858 0.0039 tmean
-5.2030 <.0001 DHWtodate:Timelagtodate -2.1840 0.0290 DHWtodate:tmean
2.6613 0.0078
ci.lb ci.ub intrcpt 8.4948
20.3268 *** DHWtodate -1.0617 -0.2028 ** tmean
-0.7898 -0.3576 *** DHWtodate:Timelagtodate -0.0001 -0.0000 *
DHWtodate:tmean 0.0060 0.0395 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
We report the fsn for the model without coefficients, as well as the AICc
for the top models with coefficients (within 2AICc units of the best model).
Our reviewer has asked for a measure of model performance:
"I don't see a clear measure of model performance beyond AICc. This is
concerning as most readers will be looking for an R2, deviance explained or
something similar, particularly from linear models. Were the fits poor?
Either way you need to report these values, otherwise it's hard for the
reader to judge how much we should be interpreting the models."
I understand that it is not possible to calculate standard R2 values from
linear mixed effects models, and that pseudo-r2s are potentially
questionable (see Robert Long's paragraph here
https://stats.stackexchange.com/questions/111150/calculating-r2-in-mixed-mod
els-using-nakagawa-schielzeths-2013-r2glmm-me), but that conditional-r2 and
marginal-r2 (Nakagawa & Schielzeth 2012) are potential candidates for
demonstrating model fit. I see that these metrics are implemented in MuMIn
(https://ecologyforacrowdedplanet.wordpress.com/2013/08/27/r-squared-in-mixe
d-models-the-easy-way/, r.squaredGLMM(x)
https://cran.r-project.org/web/packages/MuMIn/MuMIn.pdf) and in piecewiseSEM
(sem.model.fits
https://jonlefcheck.net/2013/03/13/r2-for-linear-mixed-effects-models/),
however neither of these packages/functions can utilize a rma (or rma.mv)
object. I have also looked at the original equations from Nakagawa and
Schielzeth, but am having trouble understanding how they and the output of
our models fit together, given that only one metric of variance is supplied
in the output (sigma^2).
Can anyone point me in the right direction of either a) why it isn't
possible to calculate conditional & marginal-r2 for an rma.mv object, b) how
to calculate conditional & marginal-r2 for our model, or c) if there is a
better metric to demonstrate model fit in this situation?
Please let me know if I have misunderstood anything, or if I can provide any
additional information that would assist with discussing my question.
Thanks in advance!
All the best,
Danielle
--
Danielle Claar
PhD Candidate - University of Victoria
PO Box 1700 STN CSC
Victoria, BC, V8W 2Y2 Canada
dclaar at uvic.ca
danielleclaar.weebly.com
More information about the R-sig-meta-analysis
mailing list