context("testing R2 functions") # data set.seed(12345) p1 <- 10; nsample <- 20; n <- p1 * nsample d <- data.frame(x1 = rnorm(n = n), x2 = rnorm(n = n), u1 = rep(1:p1, each = nsample), u2 = rep(1:p1, times = nsample)) d$u1 <- as.factor(d$u1); d$u2 <- as.factor(d$u2) # LMM: y with random intercept b1 <- 1; b2 <- -1; sd1 <- 1.5 d$y_re_intercept <- b1 * d$x1 + b2 * d$x2 + rep(rnorm(n = p1, sd = sd1), each = nsample) + # random intercept u1 rep(rnorm(n = p1, sd = sd1), times = nsample) + # random intercept u2 rnorm(n = n) # LMM: y with random slope b1 <- 0; sd1 <- 1; sd.x1 <- 2 d$y_re_slope <- b1 * d$x1 + rep(rnorm(n = p1, sd = sd1), each = nsample) + # random intercept u1 d$x1 * rep(rnorm(n = p1, sd = sd.x1), times = nsample) + # random slope u1 rnorm(n = n) # GLMM b1 <- 1; sd1 <- 1.5 prob <- rr2::inv.logit(b1 * d$x1 + rep(rnorm(n = p1, sd = sd1), each = nsample)) # random intercept u1 d$y_binary <- rbinom(n = n, size = 1, prob = prob) # y for PGLS b1 <- 1.5; signal <- 0.7 phy <- ape::compute.brlen(ape::rtree(n = n), method = "Grafen", power = 1) phy.x <- ape::compute.brlen(phy, method = "Grafen", power = .0001) x_trait <- ape::rTraitCont(phy.x, model = "BM", sigma = 1) e <- signal^0.5 * ape::rTraitCont(phy, model = "BM", sigma = 1) + (1-signal)^0.5 * rnorm(n=n) d$x_trait <- x_trait[match(names(e), names(x_trait))] d$y_pgls <- b1 * x_trait + e rownames(d) <- phy$tip.label # y for Phylogenetic logistic regression b1 <- 1.5; signal <- 2 e <- signal * ape::rTraitCont(phy, model = "BM", sigma = 1) e <- e[match(phy$tip.label, names(e))] d$y_phy_binary <- rbinom(n = n, size = 1, prob = rr2::inv.logit(b1 * d$x1 + e)) # data for communityPGLMM # library(dplyr) # comm = phyr::comm_a # comm$site = row.names(comm) # dat = tidyr::gather(comm, key = "sp", value = "freq", -site) %>% # dplyr::left_join(phyr::envi, by = "site") %>% # dplyr::left_join(phyr::traits, by = "sp") # dat$pa = as.numeric(dat$freq > 0) test_that("when missing mod.r, the functions will automatically create one", { testthat::skip_on_cran() # don't run on CRAN since these are time consuming # LM z.f <- lm(y_re_intercept ~ x1 + x2, data = d) z.0 <- lm(y_re_intercept ~ 1, data = d) expect_equal(R2(z.f, z.0), R2(z.f)) # GLM z.f.glm <- glm(y_binary ~ x1 + x2, data = d, family = "binomial") z.v.glm <- glm(y_binary ~ 1, data = d, family = "binomial") expect_equal(R2(mod = z.f.glm, mod.r = z.v.glm), R2(mod = z.f.glm)) # LMM z.f <- lme4::lmer(y_re_intercept ~ x1 + x2 + (1 | u1) + (1 | u2), data = d, REML = F) z.0 <- lm(y_re_intercept ~ 1, data = d) expect_equal(R2(z.f, z.0), R2(z.f)) # GLMM z.f.glmm <- lme4::glmer(y_binary ~ x1 + (1 | u1), data = d, family = "binomial") z.v.glmm <- glm(y_binary ~ 1, data = d, family = "binomial") expect_equal(R2(mod = z.f.glmm, mod.r = z.v.glmm), R2(mod = z.f.glmm)) }) test_that("when missing mod.r, the functions will automatically create one, long-run analyses", { testthat::skip_on_cran() # don't run on CRAN since these are time consuming # PGLS z.x.pgls <- phylolm::phylolm(y_pgls ~ 1, phy = phy, data = d, model = "lambda") lam.x.pgls <- round(z.x.pgls$optpar, digits = 4) z.f.pgls <- phylolm::phylolm(y_pgls ~ x_trait, phy = phy, data = d, model = "lambda") z.v.pgls <- lm(y_pgls ~ 1, data = d) expect_equal(R2(mod = z.f.pgls, mod.r = z.v.pgls, phy = phy), R2(mod = z.f.pgls, phy = phy)) for(md in c("OUrandomRoot", "OUfixedRoot", "BM", "kappa", "delta", "EB")){ z.f.pgls2 <- phylolm::phylolm(y_pgls ~ x_trait, phy = phy, data = d, model = md) z.v.pgls2 <- lm(y_pgls ~ 1, data = d) expect_equal(R2(mod = z.f.pgls2, mod.r = z.v.pgls2, phy = phy), R2(mod = z.f.pgls2, phy = phy)) } # GLS z.f.gls <- nlme::gls(y_pgls ~ x_trait, data = d, correlation = ape::corPagel(1, phy), method = "ML") expect_equal(R2(mod = z.f.gls, mod.r = z.v.pgls), R2(mod = z.f.gls)) # binaryPGLMM z.f.plog <- rr2::binaryPGLMM(y_phy_binary ~ x1, data = d, phy = phy) z.v.plog <- glm(y_phy_binary ~ 1, data = d, family = "binomial") # R2.lik can't be used with binaryPGLMM because it is not a ML method expect_message(t1 <- R2(mod = z.f.plog, mod.r = z.v.plog), "Models of class binaryPGLMM do not have a R2.lik method") expect_message(t2 <- R2(mod = z.f.plog, mod.r = z.v.plog), "Models of class binaryPGLMM do not have a R2.lik method") expect_equal(t1, t2) z.f.plog2 <- phylolm::phyloglm(y_phy_binary ~ x1, data = d, start.alpha = 1, phy = phy) z.v.plog2 <- glm(y_phy_binary ~ 1, data = d, family = "binomial") expect_message(t11 <- R2(mod = z.f.plog2, mod.r = z.v.plog2), "Models of class phyloglm only have a R2.lik method") expect_message(t12 <- R2(z.f.plog2), "Models of class phyloglm only have a R2.lik method") expect_equal(t11, t12) # # communityPGLMM # mod <- phyr::communityPGLMM(freq ~ 1 + shade + (1|sp__) + (1|site) + (1|sp__@site), # dat, tree = phyr::phylotree, REML = F) # mod.r <- lm(freq ~ 1, dat) # expect_message(R2(mod, mod.r)) # expect_equal(R2(mod, mod.r), R2(mod)) }) test_that("when lmer models were fitted with REML = T, R2.lr (but not R2.ls and R2.ce) will change it to FALSE", { z.f2 <- lme4::lmer(y_re_intercept ~ x1 + x2 + (1 | u1) + (1 | u2), data = d, REML = T) expect_warning(R2(z.f2), "mod updated with REML = F") }) test_that("alphaWarn for phylogenetic logistic regression with phyloglm", { testthat::skip_on_cran() # don't run on CRAN since these are time consuming set.seed(12345) n <- 100 b1 <- 1.5 signal <- 1 phy <- ape::compute.brlen(ape::rtree(n = n), method = 'Grafen', power = 1) phy.x <- ape::compute.brlen(phy, method = 'Grafen', power = .0001) # Generate random data x <- rnorm(n) d <- data.frame(x = x, y = 0) e <- signal * ape::rTraitCont(phy, model = 'BM', sigma = 1) e <- e[match(phy$tip.label, names(e))] d$y <- rbinom(n = n, size = 1, prob = rr2::inv.logit(b1 * d$x + e)) rownames(d) <- phy$tip.label # Use the function phyloglm() from the phylolm package. z.f.phyloglm <- phylolm::phyloglm(y ~ x, data = d, start.alpha = 1, phy = phy) z.x.phyloglm <- phylolm::phyloglm(y ~ 1, data = d, phy = phy, start.alpha = min(20, z.f.phyloglm$alpha)) z.v <- glm(y ~ x, data = d, family = 'binomial') expect_warning(R2(mod = z.f.phyloglm, mod.r = z.x.phyloglm), "In mod, alphaWarn = 2, so model refit with glm()") expect_warning(R2(mod = z.f.phyloglm, mod.r = z.v), "In mod, alphaWarn = 2, so model refit with glm()") })