# [R-meta] Calculating conditional and marginal r-squared values for rma.mv models in metafor

Michael Dewey lists at dewey.myzen.co.uk
Tue Oct 3 12:46:54 CEST 2017

```Dear Danielle

I think your interpretation of the values of sigma i correct, most of
the variation is between levels of ACC rather than between id within ACC

On 02/10/2017 20:37, Danielle Claar wrote:
> Dear Wolfgang,
>
> Thank you very much for your reply! I apologize for my delay in responding, I have been ill this past week.
>
> Thanks for pointing out the common mistake in three-level models. Adding in random=~1|ACC/id definitely makes sense, and I have gone ahead and done this for all of my models. Can you confirm that when I include this term in my models that sigma^2.1 (factor=ACC) demonstrates between-study variance, and that sigma^2.2 (factor=ACC/id) demonstrates within-study variance? So if I get this result:
> Variance Components:
>              estim    sqrt  nlvls  fixed  factor
> sigma^2.1  2.2596  1.5032     29     no     ACC
> sigma^2.2  0.1508  0.3883    192     no  ACC/id
> That means that more variance is explained by differences among papers (between-study variance) than within individual papers (within-study variance)? This result would make sense, I just want to make sure I am understanding/reporting the results correctly.
>
> The calculation of a pseudo-R^2 seems reasonable as well. If I'm understanding correctly, it doesn't give you an 'absolute' quantification of model fit, but it does tell you how much better your full model is, once the chosen moderators are included (i.e. how much variance your moderators explain).
>
> Thanks again, and all the best,
> Danielle
>
>
> -----Original Message-----
> From: Viechtbauer Wolfgang (SP) [mailto:wolfgang.viechtbauer at maastrichtuniversity.nl]
> Sent: Wednesday, September 27, 2017 1:54 PM
> To: Danielle Claar <danielle.claar at gmail.com>; r-sig-meta-analysis at r-project.org
> Subject: RE: [R-meta] Calculating conditional and marginal r-squared values for rma.mv models in metafor
>
> 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)
>
>
> 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 Ni o 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
>
> _______________________________________________
> R-sig-meta-analysis mailing list
> R-sig-meta-analysis at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis
>
> ---
> This email has been checked for viruses by AVG.
> http://www.avg.com
>
>

--
Michael
http://www.dewey.myzen.co.uk/home.html

```