https://github.com/cran/bayestestR
Tip revision: 01482dc32c49cc56111762e704c854b5f287966a authored by Dominique Makowski on 20 April 2020, 05:10:28 UTC
version 0.6.0
version 0.6.0
Tip revision: 01482dc
test-emmGrid.R
if (require("rstanarm") && require("emmeans")) {
context("emmGrid_*")
library(rstanarm)
library(emmeans)
set.seed(300)
model <- stan_glm(extra ~ group, data = sleep, refresh = 0)
em_ <- emmeans(model, ~group)
c_ <- pairs(em_)
all_ <- rbind(em_, c_)
all_summ <- summary(all_)
test_that("emmGrid ci", {
testthat::skip_on_travis()
xci <- ci(all_, ci = 0.9)
testthat::expect_equal(xci$CI_low, c(-0.236749774206338, 1.23103419307697, -2.99025704276072), tolerance = 0.2)
testthat::expect_equal(xci$CI_high, c(1.83, 3.35, -0.02), tolerance = 0.2)
})
test_that("emmGrid equivalence_test", {
testthat::skip_on_travis()
xeqtest <- equivalence_test(all_, ci = 0.9, range = c(-0.1, 0.1))
testthat::expect_equal(xeqtest$ROPE_Percentage, c(0.0553, 0, 0.0183), tolerance = 0.2)
testthat::expect_equal(xeqtest$ROPE_Equivalence, c("Undecided", "Rejected", "Undecided"))
})
test_that("emmGrid estimate_density", {
testthat::skip_on_travis()
xestden <- estimate_density(c_, method = "logspline", precision = 5)
testthat::expect_equal(xestden$x, c(-4.67, -2.91, -1.16, 0.60, 2.35), tolerance = 0.2)
testthat::expect_equal(log(xestden$y), c(-6.18, -2.12, -0.86, -3.62, -7.90), tolerance = 0.2)
})
test_that("emmGrid hdi", {
testthat::skip_on_travis()
xhdi <- hdi(all_, ci = 0.95)
testthat::expect_equal(xhdi$CI_low, c(-0.41, 0.99, -3.23), tolerance = 0.2)
testthat::expect_equal(xhdi$CI_high, c(2.06, 3.56, 0.28), tolerance = 0.2)
testthat::expect_equal(xhdi$CI_low, all_summ$lower.HPD, tolerance = 0.2)
testthat::expect_equal(xhdi$CI_high, all_summ$upper.HPD, tolerance = 0.2)
})
test_that("emmGrid p_direction", {
testthat::skip_on_travis()
xpd <- p_direction(all_, method = "direct")
testthat::expect_equal(xpd$pd, c(0.9025, 0.999, 0.952), tolerance = 0.01)
})
test_that("emmGrid p_map", {
testthat::skip_on_travis()
xpmap <- p_map(all_, precision = 2^9)
testthat::expect_equal(xpmap$p_MAP, c(0.42, 0, 0.26), tolerance = 0.1)
})
test_that("emmGrid mhdior", {
testthat::skip_on_travis()
xprope <- mhdior(all_, range = c(-0.1, 0.1), precision = 0.5)
testthat::expect_equal(xprope$mhdior, c(0.695, 1, 0.87), tolerance = 0.1)
})
test_that("emmGrid point_estimate", {
testthat::skip_on_travis()
xpest <- point_estimate(all_, centrality = "median", dispersion = TRUE)
testthat::expect_equal(xpest$Median, c(0.78, 2.29, -1.52), tolerance = 0.1)
testthat::expect_equal(xpest$MAD, c(0.60, 0.61, 0.88), tolerance = 0.1)
testthat::expect_equal(xpest$Median, all_summ$emmean, tolerance = 0.1)
})
test_that("emmGrid rope", {
testthat::skip_on_travis()
xrope <- rope(all_, range = "default", ci = .9)
testthat::expect_equal(xrope$ROPE_Percentage, c(0.0553, 0, 0.0183), tolerance = 0.1)
})
test_that("emmGrid bayesfactor_parameters", {
testthat::skip_on_travis()
testthat::skip_on_cran()
set.seed(4)
xsdbf <- bayesfactor_parameters(all_, prior = model)
testthat::expect_equal(log(xsdbf$BF), c(-2.5756125848835, 1.69713280431204, -0.212277519930343), tolerance = .1)
testthat::expect_warning(bayesfactor_parameters(all_))
# error - cannot deal with regrid / transform
testthat::expect_error(bayesfactor_parameters(regrid(all_), prior = model))
})
test_that("emmGrid bayesfactor_restricted", {
testthat::skip_on_travis()
testthat::skip_on_cran()
set.seed(4)
hyps <- c("`1` < `2`", "`1` < 0")
xrbf <- bayesfactor_restricted(em_, prior = model, hypothesis = hyps)
testthat::expect_equal(log(xrbf$BF), c(0.667225563308528, -1.62521030757486), tolerance = .1)
testthat::expect_equal(xrbf$Prior_prob, c(0.49775, 0.504), tolerance = 1e-2)
testthat::expect_equal(xrbf$Posterior_prob, c(0.952, 0.0975), tolerance = 1e-2)
testthat::expect_warning(bayesfactor_restricted(em_, hypothesis = hyps))
})
test_that("emmGrid si", {
testthat::skip_on_travis()
testthat::skip_on_cran()
set.seed(4)
xrsi <- si(em_, prior = model)
testthat::expect_equal(xrsi$CI_low, c(-0.8479125, 0.5738828), tolerance = .1)
testthat::expect_equal(xrsi$CI_high, c(2.387275, 4.004303), tolerance = .1)
})
test_that("emmGrid describe_posterior", {
testthat::skip_on_travis()
testthat::skip_on_cran()
set.seed(4)
xpost <- describe_posterior(
all_,
centrality = "median", dispersion = TRUE,
ci = 0.95, ci_method = "hdi",
test = c("pd", "rope", "bf"),
rope_range = "default", rope_ci = 0.89,
bf_prior = model
)
testthat::expect_equal(log(xpost$BF), c(-2.58, 2.00, -0.25), tolerance = 0.1)
testthat::expect_warning(describe_posterior(all_, test = "bf"))
})
## For non linear models
set.seed(333)
df <- data.frame(
G = rep(letters[1:3], each = 2),
Y = rexp(6)
)
fit_bayes <- stan_glm(Y ~ G,
data = df,
family = Gamma(link = "identity"),
refresh = 0
)
fit_bayes_prior <- update(fit_bayes, prior_PD = TRUE)
bayes_sum <- emmeans(fit_bayes, ~G)
bayes_sum_prior <- emmeans(fit_bayes_prior, ~G)
test_that("emmGrid bayesfactor_restricted2", {
testthat::skip_on_travis()
testthat::skip_on_cran()
hyps <- c("a < b", "b < c")
xrbf1 <- bayesfactor_restricted(bayes_sum, fit_bayes, hypothesis = hyps)
xrbf2 <- bayesfactor_restricted(bayes_sum, bayes_sum_prior, hypothesis = hyps)
testthat::expect_equal(xrbf1, xrbf2, tolerance = 0.1)
})
test_that("emmGrid bayesfactor_parameters", {
set.seed(333)
testthat::skip_on_travis()
testthat::skip_on_cran()
xsdbf1 <- bayesfactor_parameters(bayes_sum, prior = fit_bayes)
xsdbf2 <- bayesfactor_parameters(bayes_sum, prior = bayes_sum_prior)
# testthat::expect_equal(log(xsdbf1$BF), log(xsdbf2$BF), tolerance = 0.1)
})
}