## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(hypothesize)

## -----------------------------------------------------------------------------
# Setup: testing whether a Poisson rate equals 5
# Observed: 60 events in 10 time units
# MLE: lambda_hat = 6, se = sqrt(6/10)

# Wald: evaluate at the MLE
w <- wald_test(estimate = 6, se = sqrt(6/10), null_value = 5)

# Score: evaluate at the null
# Score function at lambda0=5: n*xbar/lambda0 - n = 60/5 - 10 = 2
# Fisher info at lambda0=5: n/lambda0 = 10/5 = 2
s <- score_test(score = 2, fisher_info = 2)

# LRT: evaluate at both
# loglik(lambda) = sum(x)*log(lambda) - n*lambda
ll_null <- 60 * log(5) - 10 * 5
ll_alt  <- 60 * log(6) - 10 * 6
l <- lrt(null_loglik = ll_null, alt_loglik = ll_alt, dof = 1)

# All three give similar (not identical) answers:
data.frame(
  test = c("Wald", "Score", "LRT"),
  statistic = c(test_stat(w), test_stat(s), test_stat(l)),
  p_value = round(c(pval(w), pval(s), pval(l)), 4)
)

## -----------------------------------------------------------------------------
# A significant Wald test
w <- wald_test(estimate = 3.0, se = 1.0)
pval(w)                    # small — rejects H0

# Its complement
cw <- complement_test(w)
pval(cw)                   # large — fails to reject

## -----------------------------------------------------------------------------
pval(complement_test(complement_test(w)))
pval(w)

## -----------------------------------------------------------------------------
class(cw)

## -----------------------------------------------------------------------------
# Both must be significant
intersection_test(0.01, 0.03)          # both < 0.05
is_significant_at(intersection_test(0.01, 0.03), 0.05)

# One non-significant component blocks rejection
intersection_test(0.01, 0.80)
is_significant_at(intersection_test(0.01, 0.80), 0.05)

## -----------------------------------------------------------------------------
# Drug effect estimate: 1.05, with SE = 0.08
# Must show: effect > 0.8 AND effect < 1.25
t_lower <- wald_test(estimate = 1.05, se = 0.08, null_value = 0.8)
t_upper <- wald_test(estimate = 1.05, se = 0.08, null_value = 1.25)

bioequiv <- intersection_test(t_lower, t_upper)
is_significant_at(bioequiv, 0.05)

## ----eval=FALSE---------------------------------------------------------------
# # This is (essentially) the actual source code of union_test:
# union_test <- function(...) {
#   complement_test(
#     do.call(intersection_test, lapply(tests, complement_test))
#   )
# }

## -----------------------------------------------------------------------------
# Three p-values
p1 <- 0.03; p2 <- 0.15; p3 <- 0.07

# Direct union
ut <- union_test(p1, p2, p3)

# Manual De Morgan construction
tests <- list(
  hypothesis_test(stat = 0, p.value = p1, dof = 1),
  hypothesis_test(stat = 0, p.value = p2, dof = 1),
  hypothesis_test(stat = 0, p.value = p3, dof = 1)
)
dm <- complement_test(
  do.call(intersection_test, lapply(tests, complement_test))
)

# Identical
c(union = pval(ut), de_morgan = pval(dm))

## -----------------------------------------------------------------------------
# Compose: take the union, then combine with another test
ut <- union_test(0.01, 0.05)
combined <- fisher_combine(ut, wald_test(estimate = 2, se = 1))
pval(combined)

## -----------------------------------------------------------------------------
# Invert a Wald test into a confidence interval
cs <- invert_test(
  test_fn = function(theta) wald_test(estimate = 2.5, se = 0.8, null_value = theta),
  grid = seq(0, 5, by = 0.01)
)
cs

## -----------------------------------------------------------------------------
# Analytical CI from the Wald test
confint(wald_test(estimate = 2.5, se = 0.8))

# Numerical CI from test inversion
c(lower = lower(cs), upper = upper(cs))

## -----------------------------------------------------------------------------
# A custom test: reject if |observation - theta| > 2
# (This has no analytical CI, but invert_test handles it)
my_test <- function(theta) {
  obs <- 5.0
  hypothesis_test(
    stat = (obs - theta)^2,
    p.value = if (abs(obs - theta) > 2) 0.01 else 0.5,
    dof = 1
  )
}

cs <- invert_test(my_test, grid = seq(0, 10, by = 0.1))
c(lower = lower(cs), upper = upper(cs))

## -----------------------------------------------------------------------------
# Test two parameters jointly
est <- c(2.0, 3.0)
V <- matrix(c(1.0, 0.3, 0.3, 1.0), 2, 2)  # correlated

w_joint <- wald_test(estimate = est, vcov = V, null_value = c(0, 0))
w_joint

## -----------------------------------------------------------------------------
# Independent parameters: diagonal vcov
se1 <- 0.8; se2 <- 1.2
V_diag <- diag(c(se1^2, se2^2))

w_multi <- wald_test(estimate = c(2.0, 3.0), vcov = V_diag)
w1 <- wald_test(estimate = 2.0, se = se1)
w2 <- wald_test(estimate = 3.0, se = se2)

# The joint statistic equals the sum of individual statistics
c(joint = test_stat(w_multi), sum = test_stat(w1) + test_stat(w2))