Extending R 2 and intra-class correlation coefficient from generalized linear mixed-effects models: capturing and characterizing biological variation

The coefficient of determination R2 quantifies the proportion of variance explained by a statistical model and is an important summary statistic of biological interest. However, estimating R2 for (generalized) linear mixed models (GLMMs) remains challenging. We have previously introduced a version of R2 that we called R2GLMM for Poisson and binomial GLMMs using biological examples, but not for other distributional families. Similarly, we earlier discussed how to estimate intra-class correlation coefficients ICC using only Poisson and binomial GLMMs. In this article, we expand our methods to all the other non-Gaussian distributions such as negative binomial and gamma distributions, which are common in biological data. While expanding our approach, we highlight two useful concepts for biologists, Jensen9s inequality and the delta method, both of which help us in understanding the properties of GLMMs. Jensen9s inequality has important implications for biologically more meaningful interpretation of GLMMs, while the delta method allows a general derivation of distribution-specific variances. We also discuss some special considerations for binomial GLMMs with binary or proportion data. We illustrate the implementation of our extension by worked examples from the field of ecology and evolution in the R environment although our method can be used regardless of statistical environments.

For overdispersed Poisson, negative binomial and gamma GLMMs with log link, the observation- level variance can be obtained via the variance of the log-normal distribution, as described above (see Appendix S1). There are two more alternative methods to obtain the same target: the 1 4 1 delta method and the trigamma function. The two alternatives have different advantages and will be 1 4 2 discussed in some detail below. The delta method for variance approximation uses a first order Taylor series expansion, which is where x is a random variable (typically represented by observations), f represents a function (e.g. log or square-root), var denotes variance, and d/dx is a (first) derivative with respect to variable x. is also possible (see Appendix S1; cf. Table 1).
The use of the trigamma function is limited to distributions with log link, but it should provide We note that in calculations of heritability (which can be seen as a type of ICC although in a strict Imagine a Poisson GLMM with log link and additive overdispersion fitted as an observation-level other symbols are the same as above. Using the log-normal approximation R 2 GLMM(m) and (adjusted) ICC GLMM can be calculated as: the median value of y ij rather than the mean of y ij , assuming a Poisson distribution. Thus, the use of compared to when using averaged y ij , which is a better estimate of E[y ij ]. Quantitative differences 2 1 0 between the two approaches may often be negligible, but when λ is small, the difference can be 2 1 1 substantial so the choice of the method needs to be reported for reproducibility (Appendix S2). Our  This recommendation for obtaining λ also applies to negative binomial GLMMs (see Table 1).
The comparison between equation (5.8) (exact) and equation (6.2) (approximate) is shown in Appendix S3. The approximation is most useful when the exact formula is not available as in the 2 3 9 case of a binomial GLMM with logit link (model 6): where y ij is the number of 'success' in n ij trials by the ith individual at the jth occasion (for binary 2 4 5 data, n ij is always 1), p ij is the underlying probability of success, and the other symbols are the same 2 4 6 as above.

4 7
To obtain corresponding values between the latent scale and data (observation) scale, we need to 2 4 8 account for Jensen's inequality (note the logit function combines of concave and convex sections).

4 9
For example, the overall intercept, on the latent scale could be transformed not just with the 2 5 0 inverse (anti) logit function ( logit −1 (x) = exp(x) / (1+ exp(x))) but also the bias corrected approximation. For the case of the binomial GLMM, we can use this approximation below given We can replace with any value obtained from the fixed part of the model (i.e. ).

5 5
Another approximation proposed by Zeger and colleagues [28] produces similar (but slightly better) 2 5 6 estimates than equation (6.7). Using our notation, this approximation can be written as: A comparison between equations (6.7) and (6.8) is also shown in Appendix S3. This approximation 2 5 9 uses the exact solution for the inverse probit function, which can be written for a model like model  The observation-level variance σ ε 2 can be thought of as being added to the latent scale on which In the following, we present a worked example by expanding the beetle dataset that was generated areas as determined by soil moisture. Larvae are exposed to two different dietary treatments  We have data on the five phenotypes, two of them sex-limited: (i) the number of eggs laid by each 3 2 7 female after random mating which we had generated previously using Poisson distributions (with Accordingly, adjusted ICC [container] is smaller in the parasite and size models but not in the 3 5 8 exploration model. The last thing to note is that for the morph models (binomial mixed models), 3 5 9 both R 2 and ICC values are larger when using the distribution-specific variance rather than the 3 6 0 observation-level variance, as discussed above (Table 3; also see Appendix S4). writing and editing of the manuscript. We thank Losia Lagisz for help in making Figure 1. This work has been benefited from discussion             Wiley & Sons.