R2_resid function

Calculate R2_resid

Calculate partial and total R2s for LMM, GLMM, PGLS, and PGLMM using R2_resid, an extension of ordinary least-squares (OLS) R2s. For LMMs and GLMMs, R2_resid is related to the method proposed by Nakagawa and Schielzeth (2013).

R2_resid(
  mod = NULL,
  mod.r = NULL,
  sigma2_d = c("s2w", "NS", "rNS"),
  phy = NULL
)

Arguments

• mod: A regression model with one of the following classes: 'lm', 'glm', 'lmerMod', 'glmerMod', 'phylolm', 'gls', 'pglmm_compare' or 'binaryPGLMM'. For 'glmerMod', only family = c('binomial', 'poisson') are supported.
• mod.r: A reduced model; if not provided, the total R2 will be given by setting 'mod.r' to the model corresponding to 'mod' with the intercept as the only predictor.
• sigma2_d: Distribution-specific variance $sigma2d$ (see Details). For binomial GLMs, GLMMs and PGLMMs with logit link functions, options are c('s2w', 'NS', 'rNS'). For binomial GLMs, GLMMs and PGLMMs with probit link functions, options are c('s2w', 'NS'). Other families use 's2w'.
• phy: The phylogeny for phylogenetic models (as a 'phylo' object), which must be specified for models of class phylolm.

Returns

R2_resid value.

Details

R2_resid works with classes 'lm', 'glm', 'lmerMod', 'glmerMod', 'phylolm', 'pglmm_compare', and 'binaryPGLMM'.

LMM (lmerMod):

$$partial R^2 = 1 - \sigma^2_{e.f}/\sigma^2_{e.r}partial R2 = 1 - σ2e.f/σ2e.r$$ $$total R^2 = 1 - \sigma^2_{e.f}/var(y)total R2 = 1 - σ2e.f/var(y)$$

where $σ2e.f$ and $σ2e.r$ are the estimated residual variances from the full and reduced LMM, and var(y) is the total variance of the response (dependent) variable.

GLMM (glmerMod):

$$total R^2 = 1 - \sigma^2_d/(\sigma^2_x + \sigma^2_b + \sigma^2_d)total R2 = 1 - σ2d/(σ2x + σ2b + σ2d)$$

where $σ2x$ and $σ2b$ are the estimated variances associated with the fixed and random effects. $σ2d$ is a term that scales the implied 'residual variance' of the GLMM (see Ives 2018, Appendix 1). The default used for $σ2d$ is $σ2w$ which is computed from the iterative weights of the GLMM. Specifically,

$$\sigma_{w}^{2}=var(g'(\mu)*(y-\mu))σ2w = var(g'(μ) * (y – μ))$$

where g'() is the derivative of the link function, and $(y – μ)$ is the difference between the data y and their predicted values $μ$. This is the default option specified by sigma2_d = 's2w'. For binomial models with a logit link function, sigma2_d = 'NS' gives the scaling $σ2d = π^2/3$ from Nakagawa and Schielzeth (2013), and sigma2_d = 'rNS' gives $σ2d = 0.8768809 * π^2/3$ which is a "corrected" version of 'NS' (see Ives 2018, Appendix 1). For binomial models with a probit link function, sigma2_d = 'NS' gives the scaling $σ2d = 1$. In general, option sigma2_d = 's2w' will give values lower than sigma2_d = 'NS' and 'rNS', but the values will be closer to R2_lik() and R2_pred(). For other forms of sigma2_d from Nakagawa and Schielzeth (2013) and Nakagawa et al. (2017), see the MuMIn package.

Partial R2s are given by the standard formula

$$partial R^2 = 1 - (1 - R^2_{.f})/(1 - R^2_{.r})partial R2 = 1 - (1 - R2.f)/(1 - R2.r)$$

where R2.f and R2.r are the total R2s for full and reduced models, respectively.

PGLS (phylolm, pglmm_compare):

$$partial R^2 = 1 - c.f * \sigma^2_{.f}/(c.r * \sigma^2_{.r})partial R2 = 1 - c.f * σ2.f/(c.r * σ2.r)$$

where $σ2.f$ and $σ2.r$ are the variances estimated for the PGLS full and reduced models, and c.f and c.r are the scaling values for full and reduce models that equal the total sum of phylogenetic branch length estimates. Note that the phylogeny needs to be specified in R2_resid.

phylolm() can have difficulties in finding solutions when there is no phylogenetic signal; when the estimate indicates no phylogenetic signal, you should refit the model with the corresponding LM.

PGLMM (pglmm_compare, binaryPGLMM):

The R2_resid for PGLMMs is computed in the same way as the GLMM (glmer), with options sigma_d = c('s2w', 'NS', 'rNS'). The estimated variance of the random effect associated with the phylogeny, $σ2b$, is multiplied by the diagonal elements of the phylogenetic covariance matrix. For binary models, this covariance matrix should be standardized so that all diagonal elements are the same (a contemporaneous or ultrametric phylogenetic tree) (Ives and Garland 2014). In case this is not done, however, the code takes the geometric average of the diagonal elements.

Note that the version of binaryPGLMM() in the package ape is replaced by a version contained within rr2 that outputs all of the required information for the calculation of R2_resid()

LM (lm) and GLM (glm):

For compatibility and generating reduced models, rr2 will compute R2_resid() for LM and GLM that correspond to LMM/PGLS and GLMM/PGLMM.

Examples

library(ape)
library(phylolm)
library(lme4)
library(nlme)
library(phyr)

set.seed(12345)
p1 <- 10
nsample <- 10
n <- p1 * nsample

d <- data.frame(x1 = 0, x2 = 0, u1 = rep(1:p1, each = nsample),
                u2 = rep(1:p1, times = nsample))
d$u1 <- as.factor(d$u1)
d$u2 <- as.factor(d$u2)

b1 <- 1
b2 <- -1
sd1 <- 1.5

d$x1 <- rnorm(n = n)
d$x2 <- rnorm(n = n)
d$y.lmm <- b1 * d$x1 + b2 * d$x2 + 
  rep(rnorm(n = p1, sd = sd1), each = nsample) +
  rep(rnorm(n = p1, sd = sd1), times = nsample) + 
  rnorm(n = n)

prob <- inv.logit(b1 * d$x1 + rep(rnorm(n = p1, sd = sd1), each = nsample))
d$y.glmm <- rbinom(n = n, size = 1, prob = prob)

# LMM with two fixed and two random effects ----
z.f <- lmer(y.lmm ~ x1 + x2 + (1 | u1) + (1 | u2), data = d, REML = FALSE)
z.x <- lmer(y.lmm ~ x1 + (1 | u1) + (1 | u2), data = d, REML = FALSE)
z.v <- lmer(y.lmm ~ 1 + (1 | u2), data = d, REML = FALSE)
z.0 <- lm(y.lmm ~ 1, data = d)

R2_resid(z.f, z.x)
R2_resid(z.f, z.v)
R2_resid(z.f)

# GLMM with one fixed and one random effect ----
z.f <- glmer(y.glmm ~ x1 + (1 | u1), data = d, family = 'binomial')
z.x <- glmer(y.glmm ~ 1 + (1 | u1), data = d, family = 'binomial')
z.v <- glm(y.glmm ~ x1, data = d, family = 'binomial')

R2_resid(z.f, z.x)
R2_resid(z.f, z.v)
R2_resid(z.f)

# PGLS with a single fixed effect ----
n <- 100
d <- data.frame(x = array(0, dim = n), y = 0)

b1 <- 1.5
signal <- 0.7

phy <- compute.brlen(rtree(n = n), method = 'Grafen', power = 1)
phy.x <- compute.brlen(phy, method = 'Grafen', power = .0001)

# Generate random data
x <- rTraitCont(phy.x, model = 'BM', sigma = 1)
e <- signal ^ 0.5 * rTraitCont(phy, model = 'BM', sigma = 1) +
  (1 - signal) ^ 0.5 * rnorm(n = n)
d$x <- x[match(names(e), names(x))]
d$y <- b1 * x + e
rownames(d) <- phy$tip.label
d$sp <- phy$tip.label

z.x <- pglmm_compare(y ~ 1, phy = phy, data = d, REML=FALSE)
z.f <- pglmm_compare(y ~ x, phy = phy, data = d, REML=FALSE)
z.v <- lm(y ~ x, data = d)

R2_resid(z.f, z.x)
R2_resid(z.f, z.v)
R2_resid(z.f)

z.x <- phylolm(y ~ 1, phy = phy, data = d, model = 'lambda')
z.f <- phylolm(y ~ x, phy = phy, data = d, model = 'lambda')
z.v <- lm(y ~ x, data = d)

R2_resid(z.f, z.x, phy = phy)
R2_resid(z.f, z.v, phy = phy)
R2_resid(z.f, phy = phy)

# This also works for models fit with gls() in {nlme}
z.f <- gls(y ~ x, data = d, correlation = corPagel(1, phy, form = ~sp), method = "ML")
z.x <- gls(y ~ 1, data = d, correlation = corPagel(1, phy, form = ~sp), method = "ML")
z.v <- lm(y ~ x, data = d)

R2_resid(z.f, z.x)
R2_resid(z.f, z.v)
R2_resid(z.f)

# But note that you need to define weights for gls() with non-ultrametric trees;
# if not, you will get a error "Matrix is not block-diagonal"

phy.nu <- rtree(n = n)

# Generate random data
e <- signal ^ 0.5 * rTraitCont(phy.nu, model = 'BM', sigma = 1) +
  (1 - signal) ^ 0.5 * rnorm(n = n)
d$x <- x[match(names(e), names(x))]
d$y <- b1 * x + e
rownames(d) <- phy.nu$tip.label
d$sp <- phy.nu$tip.label

weights <- diag(vcv.phylo(phy.nu))
z.x <- gls(y ~ 1,data = d,
         correlation = corPagel(1, phy.nu, form = ~sp),
         weights=varFixed(~weights), method = "ML")
z.f <- gls(y ~ x,data = d,
         correlation = corPagel(1, phy.nu, form = ~sp),
         weights=varFixed(~weights), method = "ML")
z.v <- lm(y ~ x, weights = 1/weights, data = d)

R2_resid(z.f, z.x)
R2_resid(z.f, z.v)
R2_resid(z.f)

# PGLMM with one fixed effect ----

n <- 100
b1 <- 1.5
signal <- 2

phy <- compute.brlen(rtree(n = n), method = 'Grafen', power = 1)
phy.x <- compute.brlen(phy, method = 'Grafen', power = .0001)

# Generate random data
x <- rnorm(n)
d <- data.frame(x = x, y = 0)

e <- signal * rTraitCont(phy, model = 'BM', sigma = 1)
e <- e[match(phy$tip.label, names(e))]

d$y <- rbinom(n = n, size = 1, prob = inv.logit(b1 * d$x + e))
rownames(d) <- phy$tip.label

# Use the function pglmm_compare in {phyr}.
z.f <- pglmm_compare(y ~ x, data = d, family = "binomial", phy = phy)
z.x <- pglmm_compare(y ~ 1, data = d, family = "binomial", phy = phy)
z.v <- glm(y ~ x, data = d, family = 'binomial')

R2_resid(z.f, z.x)
R2_resid(z.f, z.v)
R2_resid(z.f)

References

Ives A.R. and Li D. 2018. rr2: An R package to calculate R2s for regression models. Journal of Open Source Software. DOI:10.21105/joss.01028

Ives A.R. 2018. R2s for Correlated Data: Phylogenetic Models, LMMs, and GLMMs. Systematic Biology, Volume 68, Issue 2, March 2019, Pages 234-251. DOI:10.1093/sysbio/syy060

Ives A. R., Garland T., Jr. 2014. Phylogenetic regression for binary dependent variables. In: Garamszegi LZ editor. Modern Phylogenetic Comparative Methods and Their Application in Evolutionary Biology. Berlin Heidelberg, Springer-Verlag, p. 231-261.

Nakagawa S., Schielzeth H. 2013. A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods in Ecology and Evolution, 4:133-142.

Nakagawa S., Johnson P. C. D., Schielzeth H. 2017. The coefficient of determination R2 and intra-class correlation coefficient from generalized linear mixed-effects models revisited and expanded. Journal of the Royal Society Interface, 14.

See Also

MuMIn, lme4, ape, phylolm, pez

Author(s)

Anthony R. Ives