[R-meta] variance explained by fixed & random effects

[Viechtbauer, Wolfgang (SP)]

(I get so excited when I talk about this stuff that I (again) hit 'Send' too soon! Sorry -- continued properly below.)

Yes, this is a good reference.

I will try to address a few of your questions.

1) What does it means if the pseudo R^2 statistic is negative?

R^2 type measures in meta-analysis are typically computed as proportional reductions in variance components (that reflect heterogeneity) when predictors/moderators are added to the model. However, for the models we typically use for a meta-analysis, there is no guarantee that the variance components always decrease (or remain unchanged) when we add predictors. So, it can happen that a variance component increases, in which case the R^2 value is technically negative. An example:

library(metafor)

dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)

res0 <- rma(yi, vi, data=dat)
res1 <- rma(yi, vi, mods = ~ alloc, data=dat)
res1

# proportional 'reduction' in tau^2 (how R^2 is computed here)
(res0$tau2 - res1$tau2) / res0$tau2

If you look at 'res1', R^2 is given as 0%, since negative R^2 values are simply set to 0.

So what does this mean? I would say it simply means that there is no reduction in the estimate of tau^2, so 'alloc' does not account for any heterogeneity.

But also keep in mind that the value of R^2 is an estimate that can be quite inaccurate esp. if the number of studies is small (https://onlinelibrary.wiley.com/doi/abs/10.1111/bmsp.12002). One may need maybe 30 or even more studies to get a fairly decent estimate (but don't quote me on that number).

Another important point here is that we compute R^2 based on how much *heterogeneity* is accounted for, not how much of the *total variance* is accounted for. So, for example:

res1 <- rma(yi, vi, mods = ~ ablat, data=dat)
res1

This suggests that ~76% of the heterogeneity is accounted for by 'ablat' (but again, this estimate could be way off with k=13). But the total variance in this example is tau^2 plus the amount of sampling variability in the estimates (the 'vi' values). Since each study has a different amount of sampling variability, it isn't super clear how we then should think of the total variance, but one commonly used approach is to compute a 'typical' sampling variance and use that to define the total variance. For example, Higgins and Thompson (2002) suggest to compute the 'typical' sampling variance with:

k <- res1$k
wi <- 1/dat$vi
vt <- (k-1) / (sum(wi) - sum(wi^2)/sum(wi)) 
vt

Others (https://academic.oup.com/aje/article/150/2/206/55325) have suggested using the harmonic mean, which would be:

1/mean(wi)

Regardless, we can then use tau^2 + vt as the total variance. Then how much of the *total variance* does 'ablat' account?

((res0$tau2 + vt) - (res1$tau2 + vt)) / (res0$tau2 + vt)

That yields about 70%. Why don't we do that? For one, we like larger R^2 values :) But study-level moderators cannot account for sampling variability, so one could also argue it is a bit unfair to see how much of the total variability is accounted for.

2) How can we think of variance accounted for by random effects?

Another issue in this discussion is the question of how to think of variance accounted for by random effects. It makes no sense to ask how much heterogeneity is accounted for by tau^2 here, since that is always 100%. It makes sense to ask how much of the total variability is due to the random effects. In the context of model 'res0', that would be:

res0$tau2 / (res0$tau2 + vt)

If you look at res0, you will see that this is I^2.

For model 'res1', things are a bit more tricky. One approach would be as follow. We again define res0$tau2 + vt as the total variability (note I am using res0$tau2 here, so the amount of heterogeneity before we include moderators, plus sampling variability). Then:

res1$tau2 / (res0$tau2 + vt)

indicates how much of this total variability is due to the random effects, that is, residual heterogeneity (about 22%). And:

vt / (res0$tau2 + vt)

indicates how much of this total variability is due to sampling variability (about 8%). And by subtraction, we then have:

1 - res1$tau2 / (res0$tau2 + vt) - vt / (res0$tau2 + vt)

which indicates how much of the total variability is due to 'ablat' (about 70%).

This is *not* something that has been done or suggested before, but makes sense and is completely defensible.

I hope this helps.

Best,
Wolfgang
-----Original Message-----
From: Theresa Stratmann [mailto:theresa.stratmann using senckenberg.de] 
Sent: Tuesday, 26 March, 2019 9:58
To: Viechtbauer, Wolfgang (SP); r-sig-meta-analysis using r-project.org
Subject: Re: [R-meta] variance explained by fixed & random effects

Hi all,

A paper from Nakagawa & Santos (2012) seems to answer my question (specifically the section "Heterogeneity - Quantifying Heterogeneity"). 

Nakagawa, S. & Santos, E.S.A. Evol Ecol (2012) 26: 1253. https://doi.org/10.1007/s10682-012-9555-5

Theresa

> On March 14, 2019 at 5:13 PM Theresa Stratmann <theresa.stratmann using senckenberg.de> wrote:
> 
> Dear Wolfgang,
> 
> Thank-you for the fast reply! 
> 
> So the code is easy enough, but I am struggling to understand some of the nuance/ how to properly think about the different numbers. I have added some questions, between the dashed lines, below your answers. If they are within the scope of the mailing list and not to cumbersome, I would be curious for some explanations to help clear up my confusion. 
> 
> Thanks!
> 
> Theresa
> 
> One way of estimating an R^2-type measure for the fixed effects is to fit the reduced model without the predictor and then compute the proportional reduction in the variance components:
> 
> res0 <- rma.mv(yi, vi, random = ~ 1 | district/school, data=dat) 
> pmax(0, (res0$sigma2 - res$sigma2) / res0$sigma2) * 100
> max(0, (sum(res0$sigma2) - sum(res$sigma2)) / sum(res0$sigma2)) * 100
> 
> One can compute this per component or overall. There is no guarantee that the value is >= 0, so I use pmax()/max() to set negative values to 0. One could debate whether it even makes sense here to estimate how much of the school-level heterogeneity is accounted for by 'year'(since this is a district-level predictor), but that's another issue.
> 
> -----------------------------------------------------------------------------------------------------------------------------------------
> ***
> So we have: 
> 
> sigma^2 for res0 .... heterogeneity in district & school:
> 0.06506194 0.03273652
> sigma^2 for res1 .... heterogeneity in district & school, accounting for year:
> 0.11482670 0.03208976
> 
> (res0$sigma2 - res$sigma2) / res0$sigma2
> -0.76488265  0.01975659
> 
> So it looks like accounting for year increases the heterogeneity in district.
> 
> Acknowledging and moving on from the fact that year is not the best fixed effect... What does it mean if the value is negative/zero? 
> There is more variance in the random intercepts when you include the fixed effect? 
> Adding this particular fixed effect does not improve the model? 
> If I was writing this up, could I say: "Accounting for year does not explain any additional variation in the estimated variance of the random intercepts."
> More specifically: 
> "The heterogeneity accounted for by year is 0% for district level variance, and 1.98% for school level variance."
> "The overall heterogeneity accounted for by year is 0%." ???
> 
> Thinking about marginal and conditional R^2 for mixed-effects models, there is a separation of fixed and random effects: fixed effects account for X% of the variation
> and random effects accounting for Y% of the variation and Z% residual variance. But here the fixed effect is somehow included 
> in the sigma^2 (estimated variance of the random intercepts)?  This just confuses me a bit/ I missed the thought process. 
> 
> ***
> -----------------------------------------------------------------------------------------------------------------------------------------
> 
> I am not sure how to think of the question how much variance is explained by a random effect here. One could ask how much of the total heterogeneity 
> is due to district- and school-level variance, which would be:
> 
> res0$sigma2 / sum(res0$sigma2) * 100
> 
> Or one could ask how much of the unaccounted for heterogeneity (so heterogeneity not accounted for by 'year') is due to district- and school-level variance, which would be:
> 
> res$sigma2 / sum(res$sigma2) * 100
> 
> -----------------------------------------------------------------------------------------------------------------------------------------
> ***
> So, to double check, the heterogeneity accounted for by year is, what you showed above: 
> 
> max(0, (sum(res0$sigma2) - sum(res$sigma2)) / sum(res0$sigma2)) * 100
> 
> We also know that: 
> 
> res0$sigma2 / sum(res0$sigma2) * 100  = 66.52655; 33.47345
> 
> max(0, (sum(res0$sigma2) - sum(res$sigma2)) / sum(res0$sigma2)) * 100 = 0
> 
> Since these have the same denominator... is it true that: 
> [max(0, (sum(res0$sigma2) - sum(res$sigma2)) / sum(res0$sigma2)) * 100] + sum(res0$sigma2 / sum(res0$sigma2) * 100) = 100
> (I check on my data and the answer seems to be no...)
> ***
> -----------------------------------------------------------------------------------------------------------------------------------------
> 
> Or one could ask how much of the total variance (so heterogeneity + sampling variance) is due to district- and school-level variance, which would be:
> 
> k <- res$k
> wi <- 1/dat$vi
> s2 <- (k-1) * sum(wi) / (sum(wi)^2 - sum(wi^2))
> 
> res0$sigma2 / (sum(res0$sigma2) + s2) * 100
> 
> (here I am using the Higgins & Thompson, 2002, definition of a 'typical' sampling variance).
> 
> Or one could ask how much of the unaccounted for variance (so unaccounted for heterogeneity [heterogeneity not accounted for by 'year'] + sampling variance)
> is due to district- and school-level variance, which would be:
> 
> res1$sigma2 / (sum(res1$sigma2) + s2) * 100
> 
> -----------------------------------------------------------------------------------------------------------------------------------------
> ***
> Can one ask here how much of the total variance (estimate of variance in true effects + sample variation) is accounted for by year? This is really what I want. 
> ***
> -----------------------------------------------------------------------------------------------------------------------------------------
> 
> > On March 13, 2019 at 6:00 PM "Viechtbauer, Wolfgang (SP)" <wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
> > 
> > 
> > Dear Theresa,
> > 
> > You forgot some parts in your code, so to make this example reproducible, let's use this model:
> > 
> > dat <- dat.konstantopoulos2011 
> > res <- rma.mv(yi, vi, mods = ~ factor(year), random = ~ 1 | district/school, data=dat)
> > 
> > And I would consider schools nested within districts here, but this isn't pertinent to the issue.
> > 
> > One way of estimating an R^2-type measure for the fixed effects is to fit the reduced model without the predictor and then compute the proportional reduction in the variance components:
> > 
> > res0 <- rma.mv(yi, vi, random = ~ 1 | district/school, data=dat) 
> > pmax(0, (res0$sigma2 - res$sigma2) / res0$sigma2) * 100
> > max(0, (sum(res0$sigma2) - sum(res$sigma2)) / sum(res0$sigma2)) * 100
> > 
> > One can compute this per component or overall. There is no guarantee that the value is >= 0, so I use pmax()/max() to set negative values to 0. One could debate whether it even makes sense here to estimate how much of the school-level heterogeneity is accounted for by 'year' (since this is a district-level predictor), but that's another issue.
> > 
> > I am not sure how to think of the question how much variance is explained by a random effect here. One could ask how much of the total heterogeneity is due to district- and schol-level variance, which would be:
> > 
> > res0$sigma2 / sum(res0$sigma2) * 100
> > 
> > Or one could ask how much of the unaccounted for heterogeneity (so heterogeneity not accounted for by 'year') is due to district- and schol-level variance, which would be:
> > 
> > res$sigma2 / sum(res$sigma2) * 100
> > 
> > Or one could ask how much of the total variance (so heterogenity + sampling variance) is due to district- and schol-level variance, which would be:
> > 
> > k <- res$k
> > wi <- 1/dat$vi
> > s2 <- (k-1) * sum(wi) / (sum(wi)^2 - sum(wi^2))
> > res0$sigma2 / (sum(res0$sigma2) + s2) * 100
> > 
> > (here I am using the Higgins & Thompson, 2002, definition of a 'typical' sampling variance).
> > 
> > Or one could ask how much of the unaccounted for variance (so unaccounted for heterogenity + sampling variance) is due to district- and schol-level variance, which would be:
> > 
> > res1$sigma2 / (sum(res1$sigma2) + s2) * 100
> > 
> > These are at least some possibilities.
> > 
> > Best,
> > Wolfgang
> > 
> > -----Original Message-----
> > From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces using r-project.org] On Behalf Of Theresa Stratmann
> > Sent: Wednesday, 13 March, 2019 17:29
> > To: r-sig-meta-analysis using r-project.org
> > Subject: [R-meta] variance explained by fixed & random effects
> > 
> > Dear Meta-Analysis Community,
> > 
> > I am using the metafor package to help me summarize some data on habitat selection. I have a question whose answer I have not found yet on the useful metafor website, and was hoping that someone could help me (or point me to an online resource I missed). 
> > 
> > As part of the description of my meta-analysis I would like to explain how much of the variation in habitat selection is due to my fixed effects (season) versus my random effects (individual & year). 
> > 
> > For a normal mixed-effects model I know how to compute the marginal and conditional R^2 (Nakagawa & Schielzeth 2013), but I am wondering how to do this correctly for a meta-analysis, or maybe, better said, talk about variance explained by the fixed and random effects. 
> > 
> > First I looked into I^2 and the post on this for random effect models. But then I read more from Borenstein et al. (https://www.meta-analysis-workshops.com/download/common-mistakes1.pdf) and was not quiet sure anymore if this is what I want (I^2 = proportion of observed variance that reflects variance in the true effects and not the sampling error ... so this seems more of a summary for the entire model (?), which also needs to be treated with care and presented properly). 
> > 
> > Then I found this question on this mailing list: https://stat.ethz.ch/pipermail/r-sig-meta-analysis/2017-September/000232.html  ...but cannot really get this code to produce anything sensible, because I guess I am not clear on the formula for res0 vs. res1. 
> > 
> > So I will give a different example, which closely approximates what I am trying to do (a model with fixed & random effects): 
> > 
> > dat <- dat.konstantopoulos2011 
> > res <- rma.mv(yi, vi, mods = ~ factor(year), list( ~ 1 | district , ~ 1|school)  data=dat) 
> > res 
> > 
> > How would you calculate / correctly describe variance explained by the fixed effects (year)? By the random effects (district, school)? I am not tied to a particular statistic, just want to make sure that whatever I do is an honest representation of my results and gives readers (& myself) the information they need to interpret the results. 
> > 
> > I am new to meta-analyses and this is just a small part of my project, so I have been able to do some reading, but not an in-depth dive. Therefore I would appreciate any help you can offer.
> > 
> > Many thanks,
> > 
> > Theresa Stratmann 
> > 
> > B.S. Ecology, The University of Georgia 
> > M.S. Wildlife & Fisheries Biology, Clemson University 
> > 
> > PhD Student 
> > Goethe University 
> > and 
> > Senckenberg Biodiversity and Climate Research Center