skip_on_os(c("mac", "solaris"))
skip_if_not_installed("marginaleffects")
skip_if_not_installed("ggplot2")
skip_if_not_installed("lme4")
skip_if_not_installed("nlme")
set.seed(123)
n <- 200
d <- data.frame(
outcome = rnorm(n),
groups = as.factor(sample(c("treatment", "control"), n, TRUE)),
episode = as.factor(sample.int(3, n, TRUE)),
ID = as.factor(rep(1:10, n / 10)),
sex = as.factor(sample(c("female", "male"), n, TRUE, prob = c(0.4, 0.6)))
)
model1 <- lm(outcome ~ groups * episode, data = d)
test_that("test_predictions, error", {
expect_error(
test_predictions(model1, c("groups", "episode"), engine = "ggeffects"),
regex = "Argument `engine` must be"
)
pr <- predict_response(model1, c("groups", "episode"))
expect_silent(test_predictions(pr, engine = "ggeffects"))
})
test_that("test_predictions, categorical, pairwise", {
out <- test_predictions(model1, c("groups", "episode"))
expect_named(out, c("groups", "episode", "Contrast", "conf.low", "conf.high", "p.value"))
expect_equal(
out$Contrast,
c(
0.4183, -0.2036, -0.1482, 0.0709, 0.1211, -0.6219, -0.5666,
-0.3475, -0.2972, 0.0554, 0.2745, 0.3247, 0.2191, 0.2694, 0.0503
),
tolerance = 1e-3,
ignore_attr = FALSE
)
expect_identical(
out$groups,
c(
"control-treatment", "control-control", "control-treatment",
"control-control", "control-treatment", "treatment-control",
"treatment-treatment", "treatment-control", "treatment-treatment",
"control-treatment", "control-control", "control-treatment",
"treatment-control", "treatment-treatment", "control-treatment"
)
)
expect_equal(
attributes(out)$standard_error,
c(
0.23286, 0.21745, 0.23533, 0.22247, 0.21449, 0.23558, 0.25218,
0.24022, 0.23286, 0.23803, 0.22532, 0.21745, 0.24262, 0.23533,
0.22247
),
tolerance = 1e-3
)
})
test_that("test_predictions, categorical, pairwise, p_adjust", {
out1 <- test_predictions(model1, c("groups", "episode"))
out2 <- test_predictions(model1, c("groups", "episode"), p_adjust = "tukey")
expect_equal(
out1$p.value,
c(
0.074, 0.3503, 0.5295, 0.7504, 0.5729, 0.009, 0.0258, 0.1497,
0.2033, 0.8163, 0.2247, 0.137, 0.3676, 0.2538, 0.8215
),
tolerance = 1e-3,
ignore_attr = FALSE
)
expect_equal(
out2$p.value,
c(
0.4704, 0.9366, 0.9887, 0.9996, 0.9931, 0.0927, 0.2215, 0.6985,
0.7976, 0.9999, 0.8276, 0.6689, 0.9453, 0.862, 0.9999
),
tolerance = 1e-3,
ignore_attr = FALSE
)
})
test_that("test_predictions, categorical, NULL", {
out <- test_predictions(model1, c("groups", "episode"), test = NULL)
expect_named(out, c("groups", "episode", "Predicted", "conf.low", "conf.high", "p.value"))
expect_equal(out$Predicted, c(0.028, -0.3903, 0.2316, 0.1763, -0.0428, -0.0931),
tolerance = 1e-3,
ignore_attr = FALSE
)
expect_identical(
out$groups,
structure(c(1L, 2L, 1L, 2L, 1L, 2L), levels = c("control", "treatment"), class = "factor")
)
})
test_that("test_predictions, interaction", {
data(iris)
model2 <- lm(Sepal.Width ~ Sepal.Length * Species, data = iris)
out <- test_predictions(model2, c("Sepal.Length", "Species"))
expect_named(out, c("Sepal.Length", "Species", "Contrast", "conf.low", "conf.high", "p.value"))
expect_equal(out$Contrast, c(0.4788, 0.5666, 0.0878),
tolerance = 1e-3,
ignore_attr = FALSE
)
expect_identical(out$Sepal.Length, c("slope", "slope", "slope"))
})
test_that("test_predictions, by-argument", {
skip_if_not_installed("datawizard")
data(efc, package = "ggeffects")
efc$c161sex <- datawizard::to_factor(efc$c161sex)
efc$c172code <- datawizard::to_factor(efc$c172code)
mfilter <- lm(neg_c_7 ~ c161sex * c172code + e42dep + c12hour, data = efc)
prfilter <- ggpredict(mfilter, "c172code")
out <- test_predictions(prfilter, by = "c161sex")
expect_identical(nrow(out), 6L)
expect_identical(
out$c172code,
c(
"low level of education-intermediate level of education",
"low level of education-high level of education",
"intermediate level of education-high level of education",
"low level of education-intermediate level of education",
"low level of education-high level of education",
"intermediate level of education-high level of education"
)
)
expect_equal(out$p.value, c(0.3962, 0.6512, 0.7424, 0.9491, 0.0721, 0.0288), tolerance = 1e-3)
out <- test_predictions(prfilter, by = "c161sex", p_adjust = "tukey")
expect_equal(out$p.value, c(0.6727, 0.8934, 0.9422, 0.9978, 0.1699, 0.0734), tolerance = 1e-3)
prfilter <- ggpredict(mfilter, c("c172code", "c161sex"))
out <- test_predictions(prfilter, p_adjust = "tukey")
out <- out[out$c161sex %in% c("Male-Male", "Female-Female"), , drop = FALSE]
expect_equal(out$p.value, c(0.9581, 0.9976, 0.9995, 1, 0.4657, 0.2432), tolerance = 1e-3)
expect_error(test_predictions(mfilter, "c161sex", by = "c12hour"), regex = "categorical")
})
model3 <- suppressMessages(lme4::lmer(outcome ~ groups * episode + sex + (1 | ID), data = d))
test_that("test_predictions, categorical, pairwise", {
out <- test_predictions(model3, c("groups", "episode"))
expect_named(out, c("groups", "episode", "Contrast", "conf.low", "conf.high", "p.value"))
expect_equal(
out$Contrast,
c(
0.4199, -0.2051, -0.1528, 0.0666, 0.1187, -0.6251, -0.5727,
-0.3533, -0.3012, 0.0524, 0.2718, 0.3239, 0.2194, 0.2715, 0.0521
),
tolerance = 1e-3,
ignore_attr = FALSE
)
expect_identical(
out$groups,
c(
"control-treatment", "control-control", "control-treatment",
"control-control", "control-treatment", "treatment-control",
"treatment-treatment", "treatment-control", "treatment-treatment",
"control-treatment", "control-control", "control-treatment",
"treatment-control", "treatment-treatment", "control-treatment"
)
)
expect_identical(
out$episode,
c(
"1-1", "1-2", "1-2", "1-3", "1-3", "1-2", "1-2", "1-3", "1-3",
"2-2", "2-3", "2-3", "2-3", "2-3", "3-3"
)
)
})
test_that("test_predictions, categorical, NULL", {
out <- test_predictions(model3, c("groups", "episode"), test = NULL)
out <- out[order(out$groups, out$episode), ]
expect_named(out, c("groups", "episode", "Predicted", "conf.low", "conf.high", "p.value"))
expect_equal(out$Predicted, c(0.0559, 0.2611, -0.0107, -0.364, 0.2087, -0.0628),
tolerance = 1e-3,
ignore_attr = FALSE
)
expect_identical(
out$groups,
structure(c(1L, 1L, 1L, 2L, 2L, 2L), levels = c("control", "treatment"), class = "factor")
)
expect_identical(
out$episode,
structure(c(1L, 2L, 3L, 1L, 2L, 3L), levels = c("1", "2", "3"), class = "factor")
)
})
d <- nlme::Orthodont
m <- lme4::lmer(distance ~ age * Sex + (1 | Subject), data = d)
test_that("test_predictions, numeric, one focal, pairwise", {
out <- test_predictions(m, "age")
expect_named(out, c("age", "Slope", "conf.low", "conf.high", "p.value"))
expect_equal(out$Slope, 0.6602, tolerance = 1e-3, ignore_attr = FALSE)
})
test_that("test_predictions, numeric, one focal, NULL", {
out <- test_predictions(m, "age", test = NULL)
expect_named(out, c("age", "Slope", "conf.low", "conf.high", "p.value"))
expect_equal(out$Slope, 0.6602, tolerance = 1e-3, ignore_attr = FALSE)
})
test_that("test_predictions, categorical, one focal, pairwise", {
out <- test_predictions(m, "Sex")
expect_named(out, c("Sex", "Contrast", "conf.low", "conf.high", "p.value"))
expect_equal(out$Contrast, 2.321023, tolerance = 1e-3, ignore_attr = FALSE)
})
test_that("test_predictions, categorical, one focal, NULL", {
out <- test_predictions(m, "Sex", test = NULL)
expect_named(out, c("Sex", "Predicted", "conf.low", "conf.high", "p.value"))
expect_equal(out$Predicted, c(24.9688, 22.6477), tolerance = 1e-3, ignore_attr = FALSE)
})
test_that("test_predictions, masked chars in levels", {
set.seed(123)
n <- 200
d <- data.frame(
outcome = rnorm(n),
groups = as.factor(sample(c("ta-ca", "tb-cb"), n, TRUE)),
episode = as.factor(sample.int(3, n, TRUE)),
ID = as.factor(rep(1:10, n / 10)),
sex = as.factor(sample(c("1", "2"), n, TRUE, prob = c(0.4, 0.6)))
)
model <- suppressMessages(lme4::lmer(outcome ~ groups * sex + episode + (1 | ID), data = d))
out <- test_predictions(model, c("groups", "sex"))
expect_named(out, c("groups", "sex", "Contrast", "conf.low", "conf.high", "p.value"))
expect_equal(
out$Contrast,
c(-0.1854, -0.4473, -0.2076, -0.2619, -0.0222, 0.2397),
tolerance = 1e-3,
ignore_attr = FALSE
)
expect_identical(
out$groups,
c(
"ta-ca-ta-ca", "ta-ca-tb-cb", "ta-ca-tb-cb", "ta-ca-tb-cb",
"ta-ca-tb-cb", "tb-cb-tb-cb"
)
)
# ggeffects can be passed directly when model is named "model"
pr <- predict_response(model, c("groups", "sex"))
expect_silent(test_predictions(pr))
})
test_that("test_predictions, don't drop single columns", {
data(iris)
iris$Sepal.Width.factor <- factor(as.numeric(iris$Sepal.Width >= 3))
m <- lme4::lmer(Petal.Length ~ Petal.Width * Sepal.Width.factor + (1 | Species), data = iris)
expect_s3_class(
test_predictions(m, c("Sepal.Width.factor", "Petal.Width [0.5]")),
"ggcomparisons"
)
})
test_that("test_predictions, make sure random effects group is categorical", {
data(sleepstudy, package = "lme4")
set.seed(123)
sleepstudy$grp <- as.factor(sample(letters[1:3], nrow(sleepstudy), replace = TRUE))
sleepstudy$ID <- as.numeric(sleepstudy$Subject)
m <- lme4::lmer(Reaction ~ Days + (1 | ID), sleepstudy)
out <- test_predictions(ggpredict(m, "Days"))
expect_equal(out$Slope, 10.467285959584, tolerance = 1e-4)
m <- lme4::lmer(Reaction ~ Days * grp + (1 | ID), sleepstudy)
out <- test_predictions(ggpredict(m, c("Days", "grp")))
expect_equal(out$Contrast, c(-0.0813, -1.26533, -1.18403), tolerance = 1e-4)
})
test_that("test_predictions, works with glmmTMB and w/o vcov", {
skip_if_not_installed("glmmTMB")
skip_if_not_installed("datawizard")
data(efc, package = "ggeffects")
efc <- datawizard::to_factor(efc, select = c("c161sex", "c172code", "c175empl"))
efc <- datawizard::recode_values(
efc,
select = "c160age",
recode = list(`1` = "min:40", `2` = 41:64, `3` = "65:max")
)
efc <- datawizard::data_rename(
efc,
pattern = c("c161sex", "c160age", "quol_5", "c175empl"),
replacement = c("gender", "age", "qol", "employed")
)
efc <- datawizard::data_modify(efc, age = factor(age, labels = c("-40", "41-64", "65+")))
m_null <- glmmTMB::glmmTMB(qol ~ 1 + (1 | gender:employed:age), data = efc)
predictions <- ggpredict(m_null, c("gender", "employed", "age"), type = "random", ci_level = NA)
out1 <- test_predictions(predictions, verbose = FALSE)[1:5, ]
out2 <- test_predictions(predictions, vcov = TRUE, verbose = FALSE)[1:5, ]
expect_equal(out1$conf.low, out2$conf.low, tolerance = 1e-4)
expect_equal(out1$Contrast, out2$Contrast, tolerance = 1e-4)
# validate against raw values
expect_equal(out1$Contrast, c(0.06846, -0.87857, -0.79452, 0.30375, 1.48621), tolerance = 1e-4)
expect_equal(out1$conf.low, c(0.06846, -0.87857, -0.79452, 0.30375, 1.48621), tolerance = 1e-4)
})
test_that("test_predictions, correct order of character vectors", {
skip_if_not_installed("marginaleffects", minimum_version = "0.20.0")
skip_if_not_installed("datawizard")
set.seed(1234)
dat <- data.frame(
outcome = rbinom(n = 100, size = 1, prob = 0.35),
var_binom = as.factor(rbinom(n = 100, size = 1, prob = 0.3)),
var_cont = rnorm(n = 100, mean = 10, sd = 7),
groups = sample(letters[1:2], size = 100, replace = TRUE)
)
m1 <- glm(outcome ~ var_binom * groups + var_cont, data = dat, family = binomial())
pr1 <- predict_response(m1, c("var_binom", "groups"), verbose = FALSE)
out1 <- test_predictions(pr1, engine = "ggeffects")
out2 <- test_predictions(pr1)
out2 <- datawizard::data_arrange(out2, c("var_binom", "groups"))
expect_equal(out1$Contrast, out2$Contrast, tolerance = 1e-4)
expect_equal(
as.data.frame(out1)[c("var_binom", "groups")],
as.data.frame(out2)[c("var_binom", "groups")],
ignore_attr = TRUE
)
})
test_that("test_predictions, zero-inflated models", {
skip_if_not_installed("glmmTMB")
data(Salamanders, package = "glmmTMB")
m1 <- glmmTMB::glmmTMB(count ~ mined + (1 | site),
ziformula = ~mined,
family = poisson, data = Salamanders
)
pr1 <- predict_response(m1, "mined", margin = "empirical")
out1 <- test_predictions(pr1)
out2 <- test_predictions(m1, "mined", scale = "conditional")
expect_equal(out1$Contrast, out2$Contrast, tolerance = 1e-4)
pr1 <- predict_response(m1, "mined", type = "zero_inflated", margin = "empirical")
out1 <- test_predictions(pr1)
out2 <- test_predictions(m1, "mined", scale = "response")
expect_equal(out1$Contrast, out2$Contrast, tolerance = 1e-4)
pr1 <- predict_response(m1, "mined", type = "zi_prob")
out1 <- test_predictions(pr1)
out2 <- test_predictions(m1, "mined", scale = "zprob")
expect_equal(out1$Contrast, out2$Contrast, tolerance = 1e-4)
# validate against emmeans
skip_if_not_installed("emmeans")
emm <- emmeans::emmeans(m1, "mined", regrid = "response", component = "zi")
out3 <- as.data.frame(confint(emmeans::contrast(emm, method = "pairwise")))
expect_equal(attributes(out1)$standard_error, out3$SE, tolerance = 1e-1)
expect_equal(out1$conf.low, out3$asymp.LCL, tolerance = 1e-2)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.