# Tests for proportion confidence interval functions

# --- clopper_pearson ---------------------------------------------------------

test_that("clopper_pearson returns a named numeric vector of length 3", {
  result <- clopper_pearson(20, 40)
  expect_type(result, "double")
  expect_length(result, 3)
  expect_named(result, c("low", "observed", "high"))
})

test_that("clopper_pearson observed matches num/den", {
  result <- clopper_pearson(20, 40)
  expect_equal(unname(result["observed"]), 20 / 40)
})

test_that("clopper_pearson matches binom.test in base R", {
  # The implementation is documented as producing the same results as
  # binom.test(), which uses the Clopper-Pearson method.
  for (num in c(0, 1, 5, 20, 39, 40)) {
    result <- clopper_pearson(num, 40)
    bt <- binom.test(num, 40)$conf.int
    expect_equal(unname(result["low"]),     bt[1], tolerance = 1e-10)
    expect_equal(unname(result["high"]),    bt[2], tolerance = 1e-10)
  }
})

test_that("clopper_pearson interval is valid (low <= observed <= high)", {
  result <- clopper_pearson(5, 10)
  expect_lte(result["low"],      result["observed"])
  expect_gte(result["high"],     result["observed"])
})

test_that("clopper_pearson handles edge cases", {
  # num = 0: lower bound should be exactly 0
  zero_result <- clopper_pearson(0, 10)
  expect_equal(unname(zero_result["low"]),    0, tolerance = 1e-15)
  expect_gt(unname(zero_result["high"]),      0)
  # Upper bound should match binom.test exactly
  bt_zero <- binom.test(0, 10)$conf.int[2]
  expect_equal(unname(zero_result["high"]),   unname(bt_zero), tolerance = 1e-10)

  # num = den: upper bound should be exactly 1 (observed is at the boundary)
  full_result <- clopper_pearson(10, 10)
  expect_equal(unname(full_result["high"]),   1, tolerance = 1e-15)
  # Lower bound for all-successes case is well below observed (asymmetric interval)
  bt_full <- binom.test(10, 10)$conf.int
  expect_equal(unname(full_result["low"]),    unname(bt_full[1]), tolerance = 1e-10)
})

test_that("clopper_pearson respects conf.level", {
  wide <- clopper_pearson(20, 40, conf.level = 0.99)
  narrow <- clopper_pearson(20, 40, conf.level = 0.80)

  # Higher confidence level produces a wider interval
  expect_gt(unname(wide["high"]) - unname(wide["low"]),
            unname(narrow["high"]) - unname(narrow["low"]))
})


# --- z_univariate ------------------------------------------------------------

test_that("z_univariate returns a single numeric value", {
  result <- z_univariate(0.13, 0.11, 2500)
  expect_type(result, "double")
  expect_length(result, 1)
})

test_that("z_univariate equals the formula by hand calculation", {
  # z = (p_hat - p_0) / sqrt(p_0 * (1 - p_0) / n)
  unit_prop <- 0.13
  global_prop <- 0.11
  unit_denom <- 2500

  expected <- (unit_prop - global_prop) /
    sqrt((global_prop * (1 - global_prop)) / unit_denom)

  result <- z_univariate(unit_prop, global_prop, unit_denom)
  expect_equal(result, expected, tolerance = 1e-12)
})

test_that("z_univariate is zero when proportions are equal", {
  expect_equal(z_univariate(0.5, 0.5, 100), 0, tolerance = 1e-15)
})

test_that("z_univariate sign follows the direction of deviation", {
  # When unit_prop > global_prop, z should be positive
  expect_gt(z_univariate(0.2, 0.1, 100), 0)

  # When unit_prop < global_prop, z should be negative
  expect_lt(z_univariate(0.1, 0.2, 100), 0)
})


# --- waldInterval ------------------------------------------------------------

test_that("waldInterval returns a named numeric vector of length 2", {
  result <- waldInterval(x = 20, n = 40)
  expect_type(result, "double")
  expect_length(result, 2)
  expect_named(result, c("lwr", "upr"))
})

test_that("waldInterval matches documented example values", {
  # The roxygen @examples comment says: waldInterval(x = 20, n = 40)
  # returns approximately 0.345 and 0.655
  result <- waldInterval(20, 40)

  p_hat <- 20 / 40                         # 0.5
  se <- sqrt(p_hat * (1 - p_hat) / 40)      # ~0.0791
  z_crit <- qnorm(0.975)                    # ~1.96

  expect_equal(unname(result["lwr"]), p_hat - z_crit * se, tolerance = 1e-12)
  expect_equal(unname(result["upr"]), p_hat + z_crit * se, tolerance = 1e-12)
})

test_that("waldInterval interval is centered on the sample proportion", {
  result <- waldInterval(30, 50)
  midpoint <- (unname(result["lwr"]) + unname(result["upr"])) / 2
  expect_equal(midpoint, 30 / 50, tolerance = 1e-12)
})

test_that("waldInterval respects conf.level", {
  wide <- waldInterval(20, 40, conf.level = 0.99)
  narrow <- waldInterval(20, 40, conf.level = 0.80)

  expect_gt(unname(wide["upr"]) - unname(wide["lwr"]),
            unname(narrow["upr"]) - unname(narrow["lwr"]))
})


# --- agresti_coull_interval --------------------------------------------------

test_that("agresti_coull_interval returns a named numeric vector of length 3", {
  result <- agresti_coull_interval(20, 40)
  expect_type(result, "double")
  expect_length(result, 3)
  expect_named(result, c("low", "observed", "high"))
})

test_that("agresti_coull_interval observed matches num/den", {
  result <- agresti_coull_interval(20, 40)
  expect_equal(unname(result["observed"]), 20 / 40)
})

test_that("agresti_coull_interval interval is valid (low <= observed <= high)", {
  result <- agresti_coull_interval(15, 30)
  expect_lte(result["low"],      result["observed"])
  expect_gte(result["high"],     result["observed"])
})

test_that("agresti_coull_interval matches manual formula calculation", {
  num <- 20
  den <- 40
  conf_level <- 0.95

  z <- qnorm(1 - (1 - conf_level) / 2)
  n_tilde <- den + z^2
  p_tilde <- (num + z^2 / 2) / n_tilde
  margin <- z * sqrt(p_tilde * (1 - p_tilde) / n_tilde)

  result <- agresti_coull_interval(num, den, conf.level = conf_level)
  expect_equal(unname(result["low"]),    p_tilde - margin, tolerance = 1e-12)
  expect_equal(unname(result["high"]),   p_tilde + margin, tolerance = 1e-12)
})

test_that("agresti_coull_interval respects conf.level", {
  wide <- agresti_coull_interval(20, 40, conf.level = 0.99)
  narrow <- agresti_coull_interval(20, 40, conf.level = 0.80)

  expect_gt(unname(wide["high"]) - unname(wide["low"]),
            unname(narrow["high"]) - unname(narrow["low"]))
})