skip_if_not_installed("sandwich")
skip_on_cran()

# standard errors -------------------------------------
test_that("robust-se glm warn with profile-CI", {
  mglm <- glm(mpg ~ wt, data = mtcars)
  expect_message(
    ci(mglm, vcov = "HC3"),
    regex = "available"
  )
  expect_message(
    model_parameters(mglm, vcov = "HC3", ci_method = "profile"),
    regex = "modifies"
  )
})


# standard errors -------------------------------------
test_that("robust-se lm", {
  m <- lm(Petal.Length ~ Sepal.Length * Species, data = iris)
  se1 <- standard_error(m, vcov = "HC")
  se2 <- sqrt(diag(sandwich::vcovHC(m)))
  expect_equal(se1$SE, se2, tolerance = 1e-4, ignore_attr = TRUE)
})

test_that("robust-se polr", {
  skip_if_not_installed("MASS")
  data(housing, package = "MASS")
  m <- MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
  se1 <- standard_error(m, vcov = "vcovCL")
  se2 <- sqrt(diag(sandwich::vcovCL(m)))
  expect_equal(se1$SE, se2, tolerance = 1e-4, ignore_attr = TRUE)

  se1 <- standard_error(m, vcov = "vcovOPG")
  se2 <- sqrt(diag(sandwich::vcovOPG(m)))
  expect_equal(se1$SE, se2, tolerance = 1e-4, ignore_attr = TRUE)
})


test_that("robust-se zeroinfl", {
  skip_if_not_installed("pscl")
  skip_if_not_installed("clubSandwich")
  data("bioChemists", package = "pscl")
  m <- pscl::zeroinfl(art ~ fem + mar + kid5 + ment | kid5 + phd, data = bioChemists)
  se1 <- standard_error(m, vcov = "vcovCL")
  se2 <- sqrt(diag(sandwich::vcovCL(m)))
  expect_equal(se1$SE, se2, tolerance = 1e-4, ignore_attr = TRUE)

  se1 <- standard_error(m, vcov = "vcovOPG")
  se2 <- sqrt(diag(sandwich::vcovOPG(m)))
  expect_equal(se1$SE, se2, tolerance = 1e-4, ignore_attr = TRUE)
})

test_that("robust-se ivreg", {
  skip_if_not_installed("AER")
  skip_if_not_installed("clubSandwich")
  data(CigarettesSW, package = "AER")
  CigarettesSW$rprice <- with(CigarettesSW, price / cpi)
  CigarettesSW$rincome <- with(CigarettesSW, income / population / cpi)
  CigarettesSW$tdiff <- with(CigarettesSW, (taxs - tax) / cpi)

  m <- AER::ivreg(
    log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax / cpi),
    data = CigarettesSW,
    subset = year == "1995"
  )

  se1 <- standard_error(m, vcov = "vcovCL")
  se2 <- sqrt(diag(sandwich::vcovCL(m)))
  expect_equal(se1$SE, se2, tolerance = 1e-4, ignore_attr = TRUE)

  se1 <- standard_error(m, vcov = "vcovOPG")
  se2 <- sqrt(diag(sandwich::vcovOPG(m)))
  expect_equal(se1$SE, se2, tolerance = 1e-4, ignore_attr = TRUE)
})


test_that("robust-se survival", {
  skip_if_not_installed("survival")
  set.seed(123)
  m <- survival::survreg(
    formula = survival::Surv(futime, fustat) ~ ecog.ps + rx,
    data = survival::ovarian,
    dist = "logistic"
  )

  se1 <- standard_error(m, vcov = "vcovOPG")
  se2 <- sqrt(diag(sandwich::vcovOPG(m)))
  expect_equal(se1$SE, se2, tolerance = 1e-4, ignore_attr = TRUE)
})




# p-values -------------------------------------

test_that("robust-p lm", {
  m <- lm(Petal.Length ~ Sepal.Length * Species, data = iris)
  p1 <- p_value(m, vcov = "HC")
  # robust p manually
  se <- sqrt(diag(sandwich::vcovHC(m)))
  dof <- degrees_of_freedom(m, method = "wald", verbose = FALSE)
  stat <- coef(m) / se
  p2 <- 2 * pt(abs(stat), df = dof, lower.tail = FALSE)
  expect_equal(p1$p, p2, tolerance = 1e-4, ignore_attr = TRUE)
})

test_that("robust-p polr", {
  skip_if_not_installed("MASS")
  data(housing, package = "MASS")
  m <- MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
  p1 <- p_value(m, vcov = "vcovCL")
  # robust p manually
  se <- sqrt(diag(sandwich::vcovCL(m)))
  dof <- degrees_of_freedom(m, method = "wald", verbose = FALSE)
  stat <- c(m$coefficients, m$zeta) / se
  p2 <- 2 * pt(abs(stat), df = dof, lower.tail = FALSE)
  expect_equal(p1$p, p2, tolerance = 1e-4, ignore_attr = TRUE)

  p1 <- p_value(m, vcov = "vcovOPG")
  # robust p manually
  se <- sqrt(diag(sandwich::vcovOPG(m)))
  dof <- degrees_of_freedom(m, method = "wald", verbose = FALSE)
  stat <- c(m$coefficients, m$zeta) / se
  p2 <- 2 * pt(abs(stat), df = dof, lower.tail = FALSE)
  expect_equal(p1$p, p2, tolerance = 1e-4, ignore_attr = TRUE)
})


test_that("robust-p ivreg", {
  skip_if_not_installed("AER")
  data(CigarettesSW, package = "AER")
  CigarettesSW$rprice <- with(CigarettesSW, price / cpi)
  CigarettesSW$rincome <- with(CigarettesSW, income / population / cpi)
  CigarettesSW$tdiff <- with(CigarettesSW, (taxs - tax) / cpi)

  m <- AER::ivreg(
    log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax / cpi),
    data = CigarettesSW,
    subset = year == "1995"
  )
  p1 <- p_value(m, vcov = "vcovCL")
  # robust p manually
  se <- sqrt(diag(sandwich::vcovCL(m)))
  dof <- degrees_of_freedom(m, method = "wald", verbose = FALSE)
  stat <- coef(m) / se
  p2 <- 2 * pt(abs(stat), df = dof, lower.tail = FALSE)
  expect_equal(p1$p, p2, tolerance = 1e-4, ignore_attr = TRUE)

  p1 <- p_value(m, vcov = "vcovOPG")
  # robust p manually
  se <- sqrt(diag(sandwich::vcovOPG(m)))
  dof <- degrees_of_freedom(m, method = "wald", verbose = FALSE)
  stat <- coef(m) / se
  p2 <- 2 * pt(abs(stat), df = dof, lower.tail = FALSE)
  expect_equal(p1$p, p2, tolerance = 1e-4, ignore_attr = TRUE)
})


test_that("robust-p zeroinfl", {
  skip_if_not_installed("pscl")

  data("bioChemists", package = "pscl")
  m <- pscl::zeroinfl(art ~ fem + mar + kid5 + ment | kid5 + phd, data = bioChemists)

  p1 <- p_value(m, vcov = "vcovCL")
  # robust p manually
  se <- sqrt(diag(sandwich::vcovCL(m)))
  dof <- degrees_of_freedom(m, method = "wald", verbose = FALSE)
  stat <- coef(m) / se
  p2 <- 2 * pt(abs(stat), df = dof, lower.tail = FALSE)
  expect_equal(p1$p, p2, tolerance = 1e-4, ignore_attr = TRUE)

  p1 <- p_value(m, vcov = "vcovOPG")
  # robust p manually
  se <- sqrt(diag(sandwich::vcovOPG(m)))
  dof <- degrees_of_freedom(m, method = "wald", verbose = FALSE)
  stat <- coef(m) / se
  p2 <- 2 * pt(abs(stat), df = dof, lower.tail = FALSE)
  expect_equal(p1$p, p2, tolerance = 1e-4, ignore_attr = TRUE)
})


test_that("robust-p survival", {
  skip_if_not_installed("survival")

  set.seed(123)
  m <- survival::survreg(
    formula = survival::Surv(futime, fustat) ~ ecog.ps + rx,
    data = survival::ovarian,
    dist = "logistic"
  )
  p1 <- p_value(m, vcov = "vcovOPG")
  # robust p manually
  se <- sqrt(diag(sandwich::vcovOPG(m)))
  dof <- degrees_of_freedom(m, method = "wald", verbose = FALSE)
  stat <- insight::get_parameters(m)$Estimate / se
  p2 <- 2 * pt(abs(stat), df = dof, lower.tail = FALSE)
  expect_equal(p1$p, p2, tolerance = 1e-4, ignore_attr = TRUE)
})




# CI -------------------------------------


test_that("robust-ci lm", {
  data(iris)
  m <- lm(Petal.Length ~ Sepal.Length * Species, data = iris)
  ci1 <- ci(m, vcov = "HC")
  # robust CI manually
  params <- insight::get_parameters(m)
  se <- sqrt(diag(sandwich::vcovHC(m)))
  dof <- degrees_of_freedom(m, method = "wald", verbose = FALSE)
  fac <- suppressWarnings(stats::qt(0.975, df = dof))
  ci2 <- as.data.frame(cbind(
    CI_low = params$Estimate - se * fac,
    CI_high = params$Estimate + se * fac
  ))
  expect_equal(ci1$CI_low, ci2$CI_low, tolerance = 1e-4, ignore_attr = TRUE)
  expect_equal(ci1$CI_high, ci2$CI_high, tolerance = 1e-4, ignore_attr = TRUE)
})


test_that("robust-ci polr", {
  skip_if_not_installed("MASS")
  data(housing, package = "MASS")
  m <- MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
  ci1 <- ci(m, vcov = "vcovCL")
  # robust CI manually
  se <- sqrt(diag(sandwich::vcovCL(m)))
  dof <- degrees_of_freedom(m, method = "wald", verbose = FALSE)
  fac <- suppressWarnings(stats::qt(0.975, df = dof))
  ci2 <- as.data.frame(cbind(
    CI_low = c(m$coefficients, m$zeta) - se * fac,
    CI_high = c(m$coefficients, m$zeta) + se * fac
  ))
  expect_equal(ci1$CI_low, ci2$CI_low, tolerance = 1e-4, ignore_attr = TRUE)
  expect_equal(ci1$CI_high, ci2$CI_high, tolerance = 1e-4, ignore_attr = TRUE)

  ci1 <- ci(m, vcov = "vcovOPG")
  # robust CI manually
  se <- sqrt(diag(sandwich::vcovOPG(m)))
  dof <- degrees_of_freedom(m, method = "wald", verbose = FALSE)
  fac <- suppressWarnings(stats::qt(0.975, df = dof))
  ci2 <- as.data.frame(cbind(
    CI_low = c(m$coefficients, m$zeta) - se * fac,
    CI_high = c(m$coefficients, m$zeta) + se * fac
  ))
  expect_equal(ci1$CI_low, ci2$CI_low, tolerance = 1e-4, ignore_attr = TRUE)
  expect_equal(ci1$CI_high, ci2$CI_high, tolerance = 1e-4, ignore_attr = TRUE)
})


test_that("robust-ci ivreg", {
  skip_if_not_installed("AER")
  data(CigarettesSW, package = "AER")
  CigarettesSW$rprice <- with(CigarettesSW, price / cpi)
  CigarettesSW$rincome <- with(CigarettesSW, income / population / cpi)
  CigarettesSW$tdiff <- with(CigarettesSW, (taxs - tax) / cpi)

  m <- AER::ivreg(
    log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax / cpi),
    data = CigarettesSW,
    subset = year == "1995"
  )
  ci1 <- ci(m, vcov = "vcovCL")
  # robust CI manually
  se <- sqrt(diag(sandwich::vcovCL(m)))
  dof <- degrees_of_freedom(m, method = "wald", verbose = FALSE)
  fac <- suppressWarnings(stats::qt(0.975, df = dof))
  ci2 <- as.data.frame(cbind(
    CI_low = coef(m) - se * fac,
    CI_high = coef(m) + se * fac
  ))
  expect_equal(ci1$CI_low, ci2$CI_low, tolerance = 1e-4, ignore_attr = TRUE)
  expect_equal(ci1$CI_high, ci2$CI_high, tolerance = 1e-4, ignore_attr = TRUE)

  ci1 <- ci(m, vcov = "vcovOPG")
  # robust CI manually
  se <- sqrt(diag(sandwich::vcovOPG(m)))
  dof <- degrees_of_freedom(m, method = "wald", verbose = FALSE)
  fac <- suppressWarnings(stats::qt(0.975, df = dof))
  ci2 <- as.data.frame(cbind(
    CI_low = coef(m) - se * fac,
    CI_high = coef(m) + se * fac
  ))
  expect_equal(ci1$CI_low, ci2$CI_low, tolerance = 1e-4, ignore_attr = TRUE)
  expect_equal(ci1$CI_high, ci2$CI_high, tolerance = 1e-4, ignore_attr = TRUE)
})


test_that("robust-ci zeroinfl", {
  skip_if_not_installed("pscl")
  data("bioChemists", package = "pscl")
  m <- pscl::zeroinfl(art ~ fem + mar + kid5 + ment | kid5 + phd, data = bioChemists)
  ci1 <- ci(m, vcov = "vcovCL")
  # robust CI manually
  se <- sqrt(diag(sandwich::vcovCL(m)))
  dof <- degrees_of_freedom(m, method = "wald", verbose = FALSE)
  fac <- suppressWarnings(stats::qt(0.975, df = dof))
  ci2 <- as.data.frame(cbind(
    CI_low = coef(m) - se * fac,
    CI_high = coef(m) + se * fac
  ))
  expect_equal(ci1$CI_low, ci2$CI_low, tolerance = 1e-4, ignore_attr = TRUE)
  expect_equal(ci1$CI_high, ci2$CI_high, tolerance = 1e-4, ignore_attr = TRUE)

  ci1 <- ci(m, vcov = "vcovOPG")
  # robust CI manually
  se <- sqrt(diag(sandwich::vcovOPG(m)))
  dof <- degrees_of_freedom(m, method = "wald", verbose = FALSE)
  fac <- suppressWarnings(stats::qt(0.975, df = dof))
  ci2 <- as.data.frame(cbind(
    CI_low = coef(m) - se * fac,
    CI_high = coef(m) + se * fac
  ))
  expect_equal(ci1$CI_low, ci2$CI_low, tolerance = 1e-4, ignore_attr = TRUE)
  expect_equal(ci1$CI_high, ci2$CI_high, tolerance = 1e-4, ignore_attr = TRUE)
})


test_that("robust-ci survival", {
  skip_if_not_installed("survival")
  set.seed(123)
  m <- survival::survreg(
    formula = survival::Surv(futime, fustat) ~ ecog.ps + rx,
    data = survival::ovarian,
    dist = "logistic"
  )

  ci1 <- ci(m, vcov = "vcovOPG")
  # robust CI manually
  se <- sqrt(diag(sandwich::vcovOPG(m)))
  dof <- degrees_of_freedom(m, method = "wald", verbose = FALSE)
  fac <- suppressWarnings(stats::qt(0.975, df = dof))
  ci2 <- as.data.frame(cbind(
    CI_low = insight::get_parameters(m)$Estimate - se * fac,
    CI_high = insight::get_parameters(m)$Estimate + se * fac
  ))
  expect_equal(ci1$CI_low, ci2$CI_low, tolerance = 1e-4, ignore_attr = TRUE)
  expect_equal(ci1$CI_high, ci2$CI_high, tolerance = 1e-4, ignore_attr = TRUE)
})




# mixed models ----------------------

skip_if_not_installed("clubSandwich")
skip_if_not_installed("lme4")

test_that("robust-se lmer", {
  data(iris)
  set.seed(1234)
  iris$grp <- as.factor(sample(1:3, nrow(iris), replace = TRUE))

  m <- lme4::lmer(
    Sepal.Length ~ Species * Sepal.Width + Petal.Length + (1 | grp),
    data = iris
  )
  se1 <- standard_error(m, vcov = "vcovCR", vcov_args = list(type = "CR1", cluster = iris$grp))
  se2 <- sqrt(diag(clubSandwich::vcovCR(m, type = "CR1", cluster = iris$grp)))
  expect_equal(se1$SE, se2, tolerance = 1e-4, ignore_attr = TRUE)
})

test_that("robust-p lmer", {
  data(iris)
  set.seed(1234)
  iris$grp <- as.factor(sample(1:3, nrow(iris), replace = TRUE))

  m <- lme4::lmer(
    Sepal.Length ~ Species * Sepal.Width + Petal.Length + (1 | grp),
    data = iris
  )
  p1 <- p_value(m, vcov = "vcovCR", vcov_args = list(type = "CR1", cluster = iris$grp))
  se <- sqrt(diag(clubSandwich::vcovCR(m, type = "CR1", cluster = iris$grp)))
  dof <- degrees_of_freedom(m, method = "wald", verbose = FALSE)
  stat <- lme4::fixef(m) / se
  p2 <- 2 * pt(abs(stat), df = dof, lower.tail = FALSE)
  expect_equal(p1$p, p2, tolerance = 1e-4, ignore_attr = TRUE)
})

test_that("robust-ci lmer", {
  data(iris)
  set.seed(1234)
  iris$grp <- as.factor(sample(1:3, nrow(iris), replace = TRUE))

  m <- lme4::lmer(
    Sepal.Length ~ Species * Sepal.Width + Petal.Length + (1 | grp),
    data = iris
  )
  ci1 <- ci(m, vcov = "vcovCR", vcov_args = list(type = "CR1", cluster = iris$grp))
  # robust CI manually
  params <- insight::get_parameters(m)
  se <- sqrt(diag(clubSandwich::vcovCR(m, type = "CR1", cluster = iris$grp)))
  dof <- degrees_of_freedom(m, method = "wald", verbose = FALSE)
  fac <- suppressWarnings(stats::qt(0.975, df = dof))
  ci2 <- as.data.frame(cbind(
    CI_low = params$Estimate - se * fac,
    CI_high = params$Estimate + se * fac
  ))
  expect_equal(ci1$CI_low, ci2$CI_low, tolerance = 1e-4, ignore_attr = TRUE)
  expect_equal(ci1$CI_high, ci2$CI_high, tolerance = 1e-4, ignore_attr = TRUE)
})