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()
    testthat::skip_on_ci()
    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()
    testthat::skip_on_ci()
    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()
    testthat::skip_on_ci()
    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()
    testthat::skip_on_ci()
    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()
    testthat::skip_on_ci()

    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()
    testthat::skip_on_ci()

    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)
  })

  test_that("emmGrid bayesfactor_parameters / describe w/ nonlinear models", {
    testthat::skip_on_travis()
    testthat::skip_on_cran()
    testthat::skip_on_ci()
    set.seed(333)

    model <- stan_glm(vs ~ mpg,
                      data = mtcars,
                      family = "binomial",
                      refresh = 0)

    probs <- emmeans(model, "mpg", type = "resp")
    link <- emmeans(model, "mpg")

    bfp1 <- bayesfactor_parameters(probs, prior = model, null = 0.5)
    bfp2 <- bayesfactor_parameters(link, prior = model, null = 0)

    testthat::expect_equal(bfp1$BF, 0.083, tolerance = 0.01)
    testthat::expect_equal(bfp2$BF, 0.049, tolerance = 0.01)
    testthat::expect_error(bayesfactor_parameters(regrid(link), prior = model))

    hdip1 <- hdi(probs, ci = 0.9)
    hdip2 <- hdi(link, ci = 0.9)
    testthat::expect_equal(hdip1$CI_low, 0.265, tolerance = 0.01)
    testthat::expect_equal(hdip1$CI_high, 0.646, tolerance = 0.01)
    testthat::expect_equal(hdip2$CI_low, -1.016, tolerance = 0.01)
    testthat::expect_equal(hdip2$CI_high, 0.607, tolerance = 0.01)
  })
}