# R2 for Bayesian models v. Nakagawa delta method

**URL:** <https://discourse.mc-stan.org/t/r2-for-bayesian-models-v-nakagawa-delta-method/29112>\
**Category:** General\
**Tags:** performance, loo, techniques\
**Created:** [October 12, 2022, 9:59am UTC](https://discourse.mc-stan.org/t/r2-for-bayesian-models-v-nakagawa-delta-method/29112 "2022-10-12T09:59:20Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![RobertoCerina](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/robertocerina/32/20357_2.png) [@RobertoCerina](https://discourse.mc-stan.org/u/RobertoCerina)\
**Post date:** [October 12, 2022, 9:59am UTC](https://discourse.mc-stan.org/t/r2-for-bayesian-models-v-nakagawa-delta-method/29112/1 "2022-10-12T09:59:21Z")

</div>

I am looking to sensibly partition the variance from a mixed-effects logit model fit via stan. I understand typical Bayesian R^2 ([Bayesian R2 and LOO-R2](https://avehtari.github.io/bayes_R2/bayes_R2.html)) is:

R^2 = \frac{\sigma^2\_{\mu}}{\sigma^2\_{\mu} + \sigma^2\_{\epsilon}}.

I am however interested in obtaining estimates of the share of variance explained by each model component - so suppose a set of fixed effects with variance \sigma^2\_f and a random effect with variance \sigma^2\_a, then I’d like to calculate:

R\_{f}^2 = \frac{\sigma^2\_{f}}{\sigma^2\_{f} + \sigma^2\_{a} + \sigma^2\_{\epsilon}}, \mbox{ and } R\_{a}^2 = \frac{\sigma^2\_{a}}{\sigma^2\_{f} + \sigma^2\_{a} + \sigma^2\_{\epsilon}}. 

Classic Bayesian approaches for logistic regression calculate \sigma^2\_{\epsilon} at the predicted-probability level, as opposed to the latent logit level. I try to stick to this and calculate the component-specific variance with an approach which uses the predicted probabilities at each level (similarly to what is suggested in the old [Gelman (2006)](http://www.stat.columbia.edu/~gelman/research/published/rsquared.pdf) paper in the context of linear models). For every posterior simulation, I calculate:

\widehat{\sigma^2}\_a = Var \mbox{ ilogit}( \theta\_a),

But this has the very undesirable property to produce R^2 s that do not sum to 1 ( and at times are larger than 1 on their own), regardless of which method for calculating the residual variance is chosen.

The [Nakagawa et al. (2017)](https://royalsocietypublishing.org/doi/full/10.1098/rsif.2017.0213) paper offers an interesting solution using the corrected delta method to estimate the residual variance at the latent (logit) level. Based on other posts in this forum (e.g. [here](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461/19)) I have seen people have a familiarity with the paper. However it is formulated from a frequentist prospective, and i’m having trouble applying this to a Bayesian model.

The formula they derive for the latent-scale variance is:

\sigma^2\_\epsilon \approx \frac{1}{\hat{p}(1-\hat{p})};\\ \hat{p} \approx \mbox{ilogit}(\beta\_0) + 0.5\sigma^2\_\tau \frac{exp(\beta\_0)(1-exp(\beta\_0)}{(1+exp(\beta\_0))^3} ;

where \beta\_0 is the model intercept; \sigma^2\_\tau is the sum of all random-effect variance components in the model, and \hat{p} is the best guess for the expected value of our bernoulli-distributed outcome variable y.

If I calculate the R^2 replacing \sigma^2\_\tau with the variance of my linear predictor Var{(\mu)} , I get the following:

 ![R2_calculation](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/3/314fad8c7694fb2df62d7131e6dcb829df6646d3.png)

Which would seem to suggest this is a sensible approach which both gives R2 which is consistent with other R2 measures, and allows me to partition the variance as I want to.

I have the following questions which I hope folks on the forum may help with:

1. The Nakagawa formula is very insistent on this distinction between random and fixed effects - namely \sigma^2\_\tau is meant to be the sum of the variances of random effect components (not including the fixed-effects variance) - but from a Bayesian prospective I do not think there is any actual difference – everything is a ‘random’ effect (or better a ‘modeled’ effect) – so I thought it would be sensible to also include the variance of the fixed effect components in \sigma^2\_\tau – hence my choice to replace it with Var{(\mu)} .  
**Is this a sensible approach in your estimation ?**

2. Relatedly, the paper suggests it is possible to replace \beta\_0 with the fixed effects – namely \beta\_0 + \sum^K\_{k=1}\beta\_k x\_{i,k}. Here I do not understand why I could not / should not include the random effects in this predictor — again I’m not sure about how the conceptual difference between frequentist and bayesian random effects interacts with the Nakagawa formula;

3. Is it somehow possible to obtain the variance partitioning I want - namely a set of R^2 for each model component of interest which also sum to 1 - within the loo framework ?

Thanks in advance for your help and time.

---

<div class="post-metadata">

**Author:** ![prototaxites](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/prototaxites/32/16846_2.png) [@prototaxites](https://discourse.mc-stan.org/u/prototaxites)\
**Post date:** [October 12, 2022, 11:37am UTC](https://discourse.mc-stan.org/t/r2-for-bayesian-models-v-nakagawa-delta-method/29112/2 "2022-10-12T11:37:00Z")

</div>

I don’t have answers to your specific questions, but if you haven’t seen it, [this paper](https://www.biorxiv.org/content/10.1101/2021.10.17.464682v2.full) proposes a full framework for variance partitioning in a Bayesian context that you may find useful, and which addresses corner cases found in a lot of more complex models.

Doors and corners of variance partitioning in statistical ecology  
Torsti Schulz, Marjo Saastamoinen, Jarno Vanhatalo  
bioRxiv 2021.10.17.464682; doi: [Doors and corners of variance partitioning in statistical ecology | bioRxiv](https://doi.org/10.1101/2021.10.17.464682)

---

<div class="post-metadata">

**Author:** ![RobertoCerina](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/robertocerina/32/20357_2.png) [@RobertoCerina](https://discourse.mc-stan.org/u/RobertoCerina)\
**Post date:** [October 12, 2022, 1:29pm UTC](https://discourse.mc-stan.org/t/r2-for-bayesian-models-v-nakagawa-delta-method/29112/3 "2022-10-12T13:29:16Z")

</div>

Thanks for sharing that paper - this is very interesting and relevant - especially the methods to account for covariance amongst linear terms. On the terms of this paper, my example above would be a strictly ‘diagonal variance partition’.

Still, if I understand the paper correctly, the authors do not deal with a full variance-partitioning exercise, but rather stick to partitioning the linear predictor. This is cool and very useful, but we still cant calculate R^2 if the residuals are non-lineariseable.

The authors say:

> The approach is not limited to partitioning the variation in the linear predictor, but can also be used for the more traditional partitioning of explained variation, R^2, when working with a linear model or a model with linearizable residuals.

I guess one could combine the Nakagawa’s formula to linearise the variance of a Bernoulli-distributed variable, and the linear-predictor partitioning proposed in this paper to obtain a fully partitionable measure of R^2, but my questions about how exactly to do that in a Bayesian context remain, unfortunately, unanswered.

---

<div class="post-metadata">

**Author:** ![RobertoCerina](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/robertocerina/32/20357_2.png) [@RobertoCerina](https://discourse.mc-stan.org/u/RobertoCerina)\
**Post date:** [April 23, 2025, 7:58am UTC](https://discourse.mc-stan.org/t/r2-for-bayesian-models-v-nakagawa-delta-method/29112/4 "2025-04-23T07:58:04Z")

</div>

Hi folks,

I’ve come back to this a few years later and I think things are a bit clearer for me – I’ll put down the reasoning in case it’s helpful for anyone.

* * *

In a binomial logit GLMM, I have

x\_i \mid \beta\_{j[i]} \sim Binomial(n\_i,\pi\_i);  
\mu\_i = \log\left(\frac{\pi\_i}{1-\pi\_i}\right) = \alpha + \beta\_{j[i]};  
\beta\_j \sim N(0,{\sigma\_\beta}^2).

E[x\_i \mid \beta\_{j[i]} ] = n\_i\pi\_i;  
Var[x\_i \mid \beta\_{j[i]} ] = n\_i\pi\_i(1-\pi\_i).

We need to estimate the residual variance to calculate R^2. By the law of total variance, we can decompose total variation in the binomial counts in observation-level (residual) variance + random-effect variance:

Var[x\_i] = E\_\beta[Var(x\_i\mid \beta\_{j[i]})] + Var\_\beta[E[x\_i \mid \beta\_{j[i]}]]

Where E\_\beta[Var(x\_i\mid \beta\_{j[i]})] describes how much the counts fluctuate after we fix the random (and fixed) effects to their realised values; whilst Var\_\beta[E[x\_i \mid \beta\_{j[i]}]] describes the average difference in expected counts between the j groups.

The residual variance across counts is computable, though not immediately useful – we wish our R^2 metric to be computed on the scale of the linear predictor, where different sources of variance can be compared on the same scale.

Given the number of Bernoulli trials n\_i is constant, we can focus exclusively on variation arising from conditional probabilities \pi\_i. We do so on the logit scale for the reasons mentioned above:

Var[\mu\_i] = Var[logit(\pi\_i)] =\\ E\_\beta[Var(logit(\pi\_i)\mid \beta\_{j[i]})] + Var\_\beta[E[logit(\pi\_i)\mid \beta\_{j[i]}]] =\\ = \sigma^2\_\epsilon + \sigma^2\_\beta

Hence the quantity we want to estimate is the residual variance on the logit scale: \sigma^2\_\epsilon = E\_\beta[Var(logit(\pi\_i)\mid \beta\_{j[i]})].

One issue we face is that we never observe \pi\_i – we do not sample this quantity directly, so it’s challenging to describe its sampling variability. To overcome this, define the following sample estimator for the conditional probability \pi\_i:

\hat{\pi}\_i = \frac{x\_i}{n};  
E[\hat{\pi}\_i \mid \beta\_{j[i]}] = \pi\_i;  
Var[\hat{\pi}\_i\mid \beta\_{j[i]}] =\frac{\pi\_i(1-\pi\_i)}{n\_i}.

Now that we have this, we can get Var(logit(\hat{\pi}\_i)\mid \beta\_{j[i]}) via delta method:

g(x) \approx g(E[x]) + g’(E[x])(x - E[x]);  
Var\left[g(x) \right] \approx ( g’(E[x])) ^2 Var[x].

Var[g(\hat{\pi}\_i) \mid \beta\_j] = Var\left[log\left( \frac{\hat{\pi}\_i}{1-\hat{\pi}\_i} \right) \right] \approx \left( \frac{1}{\pi\_i(1-\pi\_i)} \right)^2 Var[\hat{\pi}\_i] = \frac{1}{n\_i\pi\_i(1-\pi\_i)}.

This is the conditional residual variance.

* * *

To get the marginal residual variance we need to compute the expected value of Var[g(\hat{\pi}\_i) \mid \beta\_j] with respect to the distribution of betas.

This is non-trivial. As best as I can tell, Nakagawa et al. propose multiple approaches:

1. approximate \pi\_i with the sample estimator \hat{\pi}\_i, or better with its predicted value conditional on simulated \beta\_{j[I]}, and take the average across observations:  
\hat{\sigma}^2\_{\epsilon} = \frac{1}{N} \sum\_i \frac{1}{n\_i\hat{\pi}\_i(1-\hat{\pi}\_i)}

2. replace \pi\_i with \bar{\pi}, the expectation of \pi\_i with respect to \beta:  
\bar{\pi} = E\_\beta[\pi\_i \mid \beta\_{j[i]}]

With this second simplification, the objective then becomes to get an unbiased estimate of \bar{\pi}.  
– One option is the sample average, though Nakagawa et al. caution against this because taking average of the observed conditional probabilities is not equivalent marginalising over the distribution of \beta.  
– Another intuitive option is logit^{-1}(\alpha), given \alpha is the global mean on the logit scale. But unfortunately this is equivalent to logit^{-1}(E[\mu\_i])], not E[logit^{-1}(\mu\_i)] , and due to Jensen’s inequality, those two expectations are not the same ^1.

Nakagawa et al. propose a second order Taylor expansion of logit^{-1}(\mu), around the expected value of \mu that can be used to correct for this bias:

g(x) \approx g(E[x]) + g’(E[x])(x-E[x]) + \frac{1}{2} g’’(E[x])(x-E[x])^2  
E[g(x)] \approx g(E[x]) + \frac{1}{2} Var[x] g’’(E[x])

Taking here the expectation with respect to \beta, we have:

E\_{\beta}[\pi\_i] = E\_{\beta}[logit^{-1}(\mu\_i)] \approx \\ \approx \frac{exp(E[\mu\_i])}{1+exp(E[\mu\_i])} + \frac{1}{2} Var[\mu\_i] \frac{exp(E[\mu\_i])\left(1-exp(E[\mu\_i])\right)}{(1+exp(E[\mu\_i]))^3} = \\ = \frac{exp(\alpha)}{1+exp(\alpha)} + \frac{1}{2} \sigma\_\beta^2 \frac{exp(\alpha)\left(1-exp(\alpha)\right)}{(1+exp(\alpha))^3}

Here we can then plug-in the estimates of \alpha and \sigma\_\beta to get the approximately unbiased estimate of \pi\_i, and in turn plug that into the variance formula obtained from the 1st delta method to get the observation-level, logit-scale, residual variance.

* * *

So what should a Bayesian do ? The 2nd delta method formula seems fairly inconsequential from a Bayesian prospective: we have the luxury of sampling from the posterior of \pi\_i, conditional on distribution of \beta, so the 1st delta method formula is where I would stop, and simply simulate the posterior of \sigma^2\_{\epsilon} = \frac{1}{N} \sum\_i s^2\_{\epsilon,i}, where s^2\_{\epsilon,i} is the conditional variance described above.

Even better, in a MrP setting one could get the weighted posterior predictive distribution of \sigma^2\_{\epsilon} by weighted average according to a stratification frame.

Apologies for any typos / mistakes etc. , if I got anything wrong would love to hear.

* * *

^1 It does seem strange to worry about Jensen’s inequality only here, when E\_{\beta}[Var(logit{-1}(\pi\_i) \mid \beta\_{j[i]})] = E\_{\beta}\left[\frac{1}{n\_i\pi\_i(1-\pi\_i)}\right] suffers from the similar pathology if we just sub E\_\beta[p\_i] in for \pi\_i.