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

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

## -----------------------------------------------------------------------------
# A manufacturer claims widgets weigh 100g on average
# We test 30 widgets (population SD is known to be 5g from historical data)
set.seed(42)
weights <- rnorm(30, mean = 102, sd = 5)

# Test H0: mu = 100 vs H1: mu != 100
z <- z_test(weights, mu0 = 100, sigma = 5)
z

## -----------------------------------------------------------------------------
test_stat(z)   # The z-statistic
pval(z)        # The p-value
dof(z)         # Degrees of freedom (Inf for z-test)

## -----------------------------------------------------------------------------
# Suppose we fit a regression and got coefficient = 2.3 with SE = 0.8
# Test H0: beta = 0 (no effect) vs H1: beta != 0
w <- wald_test(estimate = 2.3, se = 0.8, null_value = 0)
w

is_significant_at(w, 0.05)

## -----------------------------------------------------------------------------
w$z            # The z-score
test_stat(w)   # z² (the Wald statistic)

## -----------------------------------------------------------------------------
# From raw log-likelihoods (specify dof manually)
test <- lrt(null_loglik = -250, alt_loglik = -240, dof = 3)
test

# Or from fitted models -- dof is derived automatically
set.seed(42)
x <- 1:50
y <- 2 + 3 * x + rnorm(50, sd = 5)
m0 <- lm(y ~ 1)       # intercept only
m1 <- lm(y ~ x)       # add slope
lrt(logLik(m0), logLik(m1))

## -----------------------------------------------------------------------------
w <- wald_test(estimate = 5.2, se = 1.1)
confint(w)               # 95% CI
confint(w, level = 0.99) # 99% CI

## -----------------------------------------------------------------------------
# Is 0 in the 95% CI?
ci <- confint(w)
ci["lower"] < 0 && 0 < ci["upper"]

# Is the test significant at 5%?
is_significant_at(w, 0.05)

## -----------------------------------------------------------------------------
# Three independent studies test the same hypothesis
# None is individually significant, but together...
p1 <- 0.08
p2 <- 0.12
p3 <- 0.06

combined <- fisher_combine(p1, p2, p3)
combined

is_significant_at(combined, 0.05)

## -----------------------------------------------------------------------------
t1 <- wald_test(estimate = 1.2, se = 0.8)
t2 <- wald_test(estimate = 0.9, se = 0.6)
t3 <- wald_test(estimate = 1.5, se = 1.0)

fisher_combine(t1, t2, t3)

## -----------------------------------------------------------------------------
combined <- fisher_combine(0.05, 0.10, 0.15)
pval(combined)
test_stat(combined)
dof(combined)

## -----------------------------------------------------------------------------
# Perform 5 tests
tests <- list(
  wald_test(estimate = 2.5, se = 1.0),
  wald_test(estimate = 1.8, se = 0.9),
  wald_test(estimate = 1.2, se = 0.7),
  wald_test(estimate = 0.9, se = 0.8),
  wald_test(estimate = 2.1, se = 1.1)
)

# Original p-values
vapply(tests, pval, numeric(1))

# Bonferroni-adjusted p-values
adjusted <- adjust_pval(tests, method = "bonferroni")
vapply(adjusted, pval, numeric(1))

# Less conservative: Benjamini-Hochberg (FDR control)
adjusted_bh <- adjust_pval(tests, method = "BH")
vapply(adjusted_bh, pval, numeric(1))

## -----------------------------------------------------------------------------
adj <- adjusted[[1]]
test_stat(adj)           # Same as original
adj$original_pval        # Can access unadjusted p-value
adj$adjustment_method    # Method used

## -----------------------------------------------------------------------------
# First adjust, then combine
adjusted <- adjust_pval(tests[1:3], method = "bonferroni")
fisher_combine(adjusted[[1]], adjusted[[2]], adjusted[[3]])

## -----------------------------------------------------------------------------
# Example: Create a simple chi-squared goodness-of-fit test wrapper
chisq_gof <- function(observed, expected) {
  stat <- sum((observed - expected)^2 / expected)
  df <- length(observed) - 1
  p.value <- pchisq(stat, df = df, lower.tail = FALSE)

  hypothesis_test(
    stat = stat,
    p.value = p.value,
    dof = df,
    superclasses = "chisq_gof_test",
    observed = observed,
    expected = expected
  )
}

# Use it
result <- chisq_gof(
  observed = c(45, 35, 20),
  expected = c(40, 40, 20)
)
result

# It works with all the standard functions
pval(result)
is_significant_at(result, 0.05)