Skip to content

Commit

Permalink
update tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@strengejacke
strengejacke committed Jun 13, 2024
1 parent 0fd39ed commit bad32c1
Show file tree
Hide file tree
Showing 12 changed files with 918 additions and 821 deletions.
2 changes: 2 additions & 0 deletions R/get_variances.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,8 @@
#' @param approximation Character string, indicating the approximation method
#' for the distribution-specific (residual) variance. Can be `"lognormal"`
#' (default), `"delta"` or `"trigamma"`. See _Nakagawa et al. 2017_ for details.
#' @param model_component For models that can have a zero-inflation component,
#' specify for which component variances should be returned.
#' @param ... Currently not used.
#'
#' @return A list with following elements:
Expand Down
3 changes: 3 additions & 0 deletions man/get_variance.Rd
View file

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

821 changes: 0 additions & 821 deletions tests/testthat/test-r2_nakagawa_MuMIn.R
View file

This file was deleted.

197 changes: 197 additions & 0 deletions tests/testthat/test-r2_nakagawa_bernoulli.R
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,197 @@
skip_on_cran()

skip_if_not_installed("glmmTMB")
skip_if_not_installed("MuMIn")
skip_if_not_installed("lme4")
skip_if_not_installed("performance")
skip_if_not_installed("datawizard")


# ==============================================================================
# Bernoulli mixed models, glmmTMB ----
# ==============================================================================

test_that("glmmTMB, bernoulli", {
# dataset ---------------------------------
set.seed(123)
dat <- data.frame(
outcome = rbinom(n = 500, size = 1, prob = 0.3),
var_binom = as.factor(rbinom(n = 500, size = 1, prob = 0.3)),
var_cont = rnorm(n = 500, mean = 10, sd = 7)
)
dat$var_cont <- datawizard::standardize(dat$var_cont)
dat$group <- NA
dat$group[dat$outcome == 1] <- sample(
letters[1:5],
size = sum(dat$outcome == 1),
replace = TRUE,
prob = c(0.1, 0.2, 0.3, 0.1, 0.3)
)
dat$group[dat$outcome == 0] <- sample(
letters[1:5],
size = sum(dat$outcome == 0),
replace = TRUE,
prob = c(0.3, 0.1, 0.1, 0.4, 0.1)
)

# glmmTMB, no random slope -------------------------------------------------
m <- glmmTMB::glmmTMB(
outcome ~ var_binom + var_cont + (1 | group),
data = dat,
family = binomial(link = "logit")
)
out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
out2 <- performance::r2_nakagawa(m)
# matches theoretical values
expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)

# glmmTMB, probit, no random slope -----------------------------------------
m <- glmmTMB::glmmTMB(
outcome ~ var_binom + var_cont + (1 | group),
data = dat,
family = binomial(link = "probit")
)
out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
out2 <- performance::r2_nakagawa(m)
# matches theoretical values
expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)

# glmmTMB, cloglog, no random slope -----------------------------------------
m <- glmmTMB::glmmTMB(
outcome ~ var_binom + var_cont + (1 | group),
data = dat,
family = binomial(link = "cloglog")
)
out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
out2 <- performance::r2_nakagawa(m)
# matches theoretical values
expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)

# glmmTMB, probit, random slope -------------------------------------------------
m <- suppressWarnings(glmmTMB::glmmTMB(
outcome ~ var_binom + var_cont + (1 + var_cont | group),
data = dat,
family = binomial(link = "probit")
))
out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
out2 <- performance::r2_nakagawa(m, tolerance = 1e-8)
# matches theoretical values
expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)

# glmmTMB, random slope -------------------------------------------------
m <- glmmTMB::glmmTMB(
outcome ~ var_binom + var_cont + (1 + var_cont | group),
data = dat,
family = binomial(link = "logit")
)
out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
out2 <- performance::r2_nakagawa(m)
# matches theoretical values
expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)

# glmmTMB, cloglog, random slope -------------------------------------------------
m <- glmmTMB::glmmTMB(
outcome ~ var_binom + var_cont + (1 + var_cont | group),
data = dat,
family = binomial(link = "cloglog")
)
out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
out2 <- performance::r2_nakagawa(m)
# matches theoretical values
expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
})


# ==============================================================================
# Bernoulli mixed models, lme4 ----
# ==============================================================================

test_that("lme4, bernoulli", {
# dataset ---------------------------------
set.seed(123)
dat <- data.frame(
outcome = rbinom(n = 500, size = 1, prob = 0.3),
var_binom = as.factor(rbinom(n = 500, size = 1, prob = 0.3)),
var_cont = rnorm(n = 500, mean = 10, sd = 7)
)
dat$var_cont <- datawizard::standardize(dat$var_cont)
dat$group <- NA
dat$group[dat$outcome == 1] <- sample(
letters[1:5],
size = sum(dat$outcome == 1),
replace = TRUE,
prob = c(0.1, 0.2, 0.3, 0.1, 0.3)
)
dat$group[dat$outcome == 0] <- sample(
letters[1:5],
size = sum(dat$outcome == 0),
replace = TRUE,
prob = c(0.3, 0.1, 0.1, 0.4, 0.1)
)

# lme4, no random slope ----------------------------------------------------
m <- lme4::glmer(
outcome ~ var_binom + var_cont + (1 | group),
data = dat,
family = binomial(link = "logit")
)
out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
out2 <- performance::r2_nakagawa(m)
# matches theoretical values
expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)

# lme4, probit, no random slope ---------------------------------------------
m <- lme4::glmer(
outcome ~ var_binom + var_cont + (1 | group),
data = dat,
family = binomial(link = "probit")
)
out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
out2 <- performance::r2_nakagawa(m)
# matches theoretical values
expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)

# lme4, cloglog, no random slope ---------------------------------------------
m <- lme4::glmer(
outcome ~ var_binom + var_cont + (1 | group),
data = dat,
family = binomial(link = "cloglog")
)
out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
out2 <- performance::r2_nakagawa(m)
# matches theoretical values
expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)

# lme4, random slope -------------------------------------------------
m <- lme4::glmer(
outcome ~ var_binom + var_cont + (1 + var_cont | group),
data = dat,
family = binomial(link = "logit")
)
out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
out2 <- performance::r2_nakagawa(m)
# matches theoretical values
expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)

# lme4, cloglog, random slope -------------------------------------------------
m <- lme4::glmer(
outcome ~ var_binom + var_cont + (1 + var_cont | group),
data = dat,
family = binomial(link = "cloglog")
)
out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
out2 <- performance::r2_nakagawa(m)
# matches theoretical values
expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
})
27 changes: 27 additions & 0 deletions tests/testthat/test-r2_nakagawa_beta.R
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,27 @@
skip_on_cran()

skip_if_not_installed("glmmTMB")
skip_if_not_installed("MuMIn")
skip_if_not_installed("performance")


# ==============================================================================
# beta mixed models, glmmTMB
# ==============================================================================

skip_if_not_installed("betareg")

test_that("glmmTMB, beta_family", {
# dataset ---------------------------------
data(FoodExpenditure, package = "betareg")
m <- glmmTMB::glmmTMB(
I(food / income) ~ income + (1 | persons),
data = FoodExpenditure,
family = glmmTMB::beta_family()
)
out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
out2 <- suppressWarnings(performance::r2_nakagawa(m, verbose = FALSE))
# matches theoretical values
expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
})
50 changes: 50 additions & 0 deletions tests/testthat/test-r2_nakagawa_binomial.R
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,50 @@
skip_on_cran()

skip_if_not_installed("glmmTMB")
skip_if_not_installed("MuMIn")
skip_if_not_installed("lme4")
skip_if_not_installed("performance")


# ==============================================================================
# Binomial mixed models, lme4 ----
# ==============================================================================

test_that("lme4, binomial", {
# dataset
data(cbpp, package = "lme4")

# lme4, no random slope ----------------------------------------------------
m <- lme4::glmer(
cbind(incidence, size - incidence) ~ period + (1 | herd),
data = cbpp,
family = binomial()
)
out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
out2 <- performance::r2_nakagawa(m)
# matches theoretical values
expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-3)
expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-3)
})


# ==============================================================================
# Binomial mixed models, glmmTMB ----
# ==============================================================================

test_that("glmmTMB, binomial", {
# dataset
data(cbpp, package = "lme4")

# lme4, no random slope ----------------------------------------------------
m <- glmmTMB::glmmTMB(
cbind(incidence, size - incidence) ~ period + (1 | herd),
data = cbpp,
family = binomial()
)
out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
out2 <- performance::r2_nakagawa(m)
# matches theoretical values
expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-3)
expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-3)
})
79 changes: 79 additions & 0 deletions tests/testthat/test-r2_nakagawa_gamma.R
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,79 @@
skip_on_cran()

skip_if_not_installed("MuMIn")
skip_if_not_installed("performance")

# ==============================================================================
# Gamma mixed models, glmmTMB ----
# ==============================================================================

# ==============================================================================
# Validate against Nakagawa et al. 2017 paper!
test_that("glmmTMB, Gamma", {
# example data from Nakagawa et al. 2017
Population <- gl(12, 80, 960)
Container <- gl(120, 8, 960)
Sex <- factor(rep(rep(c("Female", "Male"), each = 8), 60))
Habitat <- factor(rep(rep(c("Dry", "Wet"), each = 4), 120))
Treatment <- factor(rep(c("Cont", "Exp"), 480))
Data <- data.frame(
Population = Population, Container = Container, Sex = Sex,
Habitat = Habitat, Treatment = Treatment
)
DataFemale <- Data[Data$Sex == "Female", ]
set.seed(777)
PopulationE <- rnorm(12, 0, sqrt(0.4))
ContainerE <- rnorm(120, 0, sqrt(0.05))
EggL <- with(DataFemale, 1.1 + 0.5 * (as.numeric(Treatment) - 1) + 0.1 * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container] + rnorm(480, 0, sqrt(0.1)))
DataFemale$Egg <- rpois(length(EggL), exp(EggL))
DataAll <- Data
PopulationE <- rnorm(12, 0, sqrt(0.5))
ContainerE <- rnorm(120, 0, sqrt(0.8))
ParasiteL <- with(DataAll, 1.8 + 2 * (-1) * (as.numeric(Sex) - 1) + 0.8 * (-1) * (as.numeric(Treatment) - 1) + 0.7 * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container])
DataAll$Parasite <- rnbinom(length(ParasiteL), size = 5, mu = exp(ParasiteL))
PopulationE <- rnorm(12, 0, sqrt(1.3))
ContainerE <- rnorm(120, 0, sqrt(0.3))
DataAll$BodyL <- 15 + 3 * (-1) * (as.numeric(Sex) - 1) + 0.4 * (as.numeric(Treatment) - 1) + 0.15 * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container] + rnorm(960, 0, sqrt(1.2))
PopulationE <- rnorm(12, 0, sqrt(0.2))
ContainerE <- rnorm(120, 0, sqrt(0.2))
ExplorationL <- with(DataAll, 4 + 1 * (-1) * (as.numeric(Sex) - 1) + 2 * (as.numeric(Treatment) - 1) + 0.5 * (-1) * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container])
DataAll$Exploration <- rgamma(length(ExplorationL), shape = exp(ExplorationL) * 0.3, rate = 0.3)

sizemodeGLMERr <- lme4::glmer(
BodyL ~ 1 + (1 | Population) + (1 | Container),
family = Gamma(link = log),
data = DataAll
)
# Fit alternative model including fixed and all random effects
sizemodeGLMERf <- lme4::glmer(
BodyL ~ Sex + Treatment + Habitat + (1 | Population) + (1 | Container),
family = Gamma(link = log), data = DataAll
)

VarF <- var(as.vector(model.matrix(sizemodeGLMERf) %*% lme4::fixef(sizemodeGLMERf)))
nuF <- 1 / attr(lme4::VarCorr(sizemodeGLMERf), "sc")^2
VarOdF <- 1 / nuF # the delta method
VarOlF <- log(1 + 1 / nuF) # log-normal approximation
VarOtF <- trigamma(nuF) # trigamma function

# lognormal
R2glmmM <- VarF / (VarF + sum(as.numeric(lme4::VarCorr(sizemodeGLMERf))) + VarOlF)
R2glmmC <- (VarF + sum(as.numeric(lme4::VarCorr(sizemodeGLMERf)))) / (VarF + sum(as.numeric(lme4::VarCorr(sizemodeGLMERf))) + VarOlF)
out <- performance::r2_nakagawa(sizemodeGLMERf, null_model = sizemodeGLMERr)
expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)

# delta
R2glmmM <- VarF / (VarF + sum(as.numeric(lme4::VarCorr(sizemodeGLMERf))) + VarOdF)
R2glmmC <- (VarF + sum(as.numeric(lme4::VarCorr(sizemodeGLMERf)))) / (VarF + sum(as.numeric(lme4::VarCorr(sizemodeGLMERf))) + VarOdF)
out <- performance::r2_nakagawa(sizemodeGLMERf, null_model = sizemodeGLMERr, approximation = "delta")
expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)

# trigamma
R2glmmM <- VarF / (VarF + sum(as.numeric(lme4::VarCorr(sizemodeGLMERf))) + VarOtF)
R2glmmC <- (VarF + sum(as.numeric(lme4::VarCorr(sizemodeGLMERf)))) / (VarF + sum(as.numeric(lme4::VarCorr(sizemodeGLMERf))) + VarOtF)
out <- performance::r2_nakagawa(sizemodeGLMERf, null_model = sizemodeGLMERr, approximation = "trigamma")
expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
})
Loading

0 comments on commit bad32c1

Please sign in to comment.