R2_resid: Calculate R2_resid
In arives/rr2: R2s for Regression Models
R2_resid R Documentation
Calculate R2_resid
Description
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).
Usage
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 \sigma^2_d (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'.
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}
total R^2 = 1 - \sigma^2_{e.f}/var(y)
where \sigma^2_{e.f} and \sigma^2_{e.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)
where \sigma^2_x and \sigma^2_b are the estimated variances associated with the fixed and random effects. \sigma^2_d is a term that scales the implied 'residual variance' of the GLMM (see Ives 2018, Appendix 1). The default used for \sigma^2_d is \sigma^2_w which is computed from the iterative weights of the GLMM. Specifically,
\sigma_{w}^{2}=var(g'(\mu)*(y-\mu))
where g'() is the derivative of the link function, and (y-\mu) is the difference between the data y and their predicted values \mu. This is the default option specified by sigma2_d = 's2w'. For binomial models with a logit link function, sigma2_d = 'NS' gives the scaling \sigma^2_d = \pi^2/3 from Nakagawa and Schielzeth (2013), and sigma2_d = 'rNS' gives \sigma^2_d = 0.8768809 * \pi^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 \sigma^2_d = 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})
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})
where \sigma^2_{.f} and \sigma^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, \sigma^2_b, 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.
Value
R2_resid value.
Author(s)
Anthony R. Ives
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
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)
arives/rr2 documentation built on Aug. 20, 2023, 3:24 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| R2_resid | R Documentation |
Calculate R2_resid
Description
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).
Usage
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 |
phy |
The phylogeny for phylogenetic models (as a 'phylo' object), which must be specified for models of class 'phylolm'. |
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}
total R^2 = 1 - \sigma^2_{e.f}/var(y)
where \sigma^2_{e.f} and \sigma^2_{e.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)
where \sigma^2_x and \sigma^2_b are the estimated variances associated with the fixed and random effects. \sigma^2_d is a term that scales the implied 'residual variance' of the GLMM (see Ives 2018, Appendix 1). The default used for \sigma^2_d is \sigma^2_w which is computed from the iterative weights of the GLMM. Specifically,
\sigma_{w}^{2}=var(g'(\mu)*(y-\mu))
where g'() is the derivative of the link function, and (y-\mu) is the difference between the data y and their predicted values \mu. This is the default option specified by sigma2_d = 's2w'. For binomial models with a logit link function, sigma2_d = 'NS' gives the scaling \sigma^2_d = \pi^2/3 from Nakagawa and Schielzeth (2013), and sigma2_d = 'rNS' gives \sigma^2_d = 0.8768809 * \pi^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 \sigma^2_d = 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})
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})
where \sigma^2_{.f} and \sigma^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, \sigma^2_b, 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.
Value
R2_resid value.
Author(s)
Anthony R. Ives
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
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)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.