tests/testthat/test-nl-hessian-equivalence.R

In choicer: Discrete Choice Models for Economic Applications

# Tests for nl_loglik_hessian_parallel (O1 analytical NL Hessian) — TESTER INDEPENDENT VALIDATION
# Reference: runs/2026-06-19-nl-analytical-hessian-o1/test-spec.md
# This file is INDEPENDENT of test-nl-hessian.R (builder's tests).
# Written by tester from test-spec.md only — no implementation knowledge assumed.

# Tolerances from test-spec.md §10
TOL_EQUIV   <- 1e-6    # element-wise Hessian equivalence (analytical vs numeric oracle)
TOL_SYMM    <- 1e-10   # symmetry
TOL_SE      <- 1e-4    # SE equivalence via vcov()
# Relaxed per test-spec §5: at lambda < 0.1, finite-diff error in oracle can be amplified;
# verify at 1e-4 and flag for human review (do not auto-fail).
TOL_EXTLAM  <- 1e-3    # relaxed for extreme lambda < 0.1 (1e-4 is borderline; use 1e-3 and note)

# --------------------------------------------------------------------------
# Fixture helpers
# All fixtures set.seed(2024) at start per test-spec.md §2.
# --------------------------------------------------------------------------

make_fixture_A <- function() {
  # Two-nest balanced, no ASC, no outside option
  # K=2, K_lambda=2, J_asc=0, P=4
  set.seed(2024)
  N <- 120; J <- 4
  dt <- data.table::data.table(
    id   = rep(1:N, each = J),
    alt  = rep(1:J, N),
    x1   = rnorm(N * J),
    x2   = rnorm(N * J),
    nest = ifelse(rep(1:J, N) <= 2, "A", "B")
  )
  beta_true  <- c(0.5, -0.3)
  dt[, util := beta_true[1] * x1 + beta_true[2] * x2]
  dt[, choice := {
    u <- util; pr <- exp(u - max(u)); pr <- pr / sum(pr)
    cs <- cumsum(pr); r <- runif(1)
    as.integer(alt == alt[which(cs >= r)[1]])
  }, by = id]
  d <- prepare_nl_data(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
  list(d = d, theta = c(0.5, -0.3, 0.8, 0.7),
       J = J, K = 2, K_l = 2, J_asc = 0, P = 4,
       use_asc = FALSE, include_outside_option = FALSE)
}

make_fixture_B <- function() {
  # Two-nest balanced, WITH ASC, no outside option
  # K=2, K_lambda=2, J_asc=3, P=7
  set.seed(2024)
  N <- 120; J <- 4
  dt <- data.table::data.table(
    id   = rep(1:N, each = J),
    alt  = rep(1:J, N),
    x1   = rnorm(N * J),
    x2   = rnorm(N * J),
    nest = ifelse(rep(1:J, N) <= 2, "A", "B")
  )
  dt[, util := 0.5 * x1 - 0.3 * x2]
  dt[, choice := {
    u <- util; pr <- exp(u - max(u)); pr <- pr / sum(pr)
    cs <- cumsum(pr); r <- runif(1)
    as.integer(alt == alt[which(cs >= r)[1]])
  }, by = id]
  d <- prepare_nl_data(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
  list(d = d, theta = c(0.5, -0.3, 0.8, 0.7, rep(0, J - 1)),
       J = J, K = 2, K_l = 2, J_asc = 3, P = 7,
       use_asc = TRUE, include_outside_option = FALSE)
}

make_fixture_C <- function() {
  # WITH outside option, use_asc = TRUE
  # K=2, K_lambda=2, J_asc=4, P=8
  set.seed(2024)
  N <- 100; J <- 4
  dt <- data.table::data.table(
    id   = rep(1:N, each = (J + 1)),
    alt  = rep(0:J, N),
    x1   = c(rep(0, N), rnorm(N * J)),
    x2   = c(rep(0, N), rnorm(N * J)),
    nest = ifelse(rep(0:J, N) == 0, "OUT",
                  ifelse(rep(0:J, N) <= 2, "A", "B"))
  )
  dt[, choice := 0L]
  dt[, choice := {
    ch <- sample(0:J, 1)
    as.integer(alt == ch)
  }, by = id]
  dt[, n_chosen := sum(choice), by = id]
  dt[n_chosen == 0, choice := as.integer(alt == 0)]
  dt[n_chosen > 1, choice := {
    idx_chosen <- which(choice == 1L)
    if (length(idx_chosen) > 1) choice[idx_chosen[-1]] <- 0L
    choice
  }, by = id]
  d <- prepare_nl_data(dt, "id", "alt", "choice", c("x1", "x2"), "nest",
                       outside_opt_label = "0", include_outside_option = TRUE)
  K_l <- 2; P <- 2 + K_l + J
  list(d = d, theta = c(0.5, -0.3, 0.8, 0.7, rep(0, J)),
       J = J, K = 2, K_l = K_l, J_asc = J, P = P,
       use_asc = TRUE, include_outside_option = TRUE)
}

make_fixture_D <- function() {
  # Singleton nest: nest A has 1 alt, nest B has 3 alts
  # K=2, K_lambda=1, J_asc=3, P=6
  set.seed(2024)
  N <- 100; J <- 4
  dt <- data.table::data.table(
    id   = rep(1:N, each = J),
    alt  = rep(1:J, N),
    x1   = rnorm(N * J),
    x2   = rnorm(N * J),
    nest = ifelse(rep(1:J, N) == 1, "A", "B")
  )
  dt[, choice := {
    pr <- rep(1/J, J)
    as.integer(alt == sample(1:J, 1, prob = pr))
  }, by = id]
  d <- prepare_nl_data(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
  list(d = d, theta = c(0.5, -0.3, 0.8, rep(0, J - 1)),
       J = J, K = 2, K_l = 1, J_asc = 3, P = 6,
       use_asc = TRUE, include_outside_option = FALSE)
}

make_fixture_E <- function() {
  # Singleton + outside option: Nest A: alt 1 (singleton); Nest B: alts 2,3
  # K_lambda=1, J_asc=3, P=6
  set.seed(2024)
  N <- 100; J <- 3
  dt <- data.table::data.table(
    id   = rep(1:N, each = (J + 1)),
    alt  = rep(0:J, N),
    x1   = c(rep(0, N), rnorm(N * J)),
    x2   = c(rep(0, N), rnorm(N * J)),
    nest = ifelse(rep(0:J, N) == 0, "OUT",
                  ifelse(rep(0:J, N) == 1, "A", "B"))
  )
  dt[, choice := 0L]
  dt[, choice := {
    ch <- sample(0:J, 1)
    as.integer(alt == ch)
  }, by = id]
  dt[, n_chosen := sum(choice), by = id]
  dt[n_chosen == 0, choice := as.integer(alt == 0)]
  d <- prepare_nl_data(dt, "id", "alt", "choice", c("x1", "x2"), "nest",
                       outside_opt_label = "0", include_outside_option = TRUE)
  K_l <- 1; P <- 2 + K_l + J
  list(d = d, theta = c(0.5, -0.3, 0.8, rep(0, J)),
       J = J, K = 2, K_l = K_l, J_asc = J, P = P,
       use_asc = TRUE, include_outside_option = TRUE)
}

make_fixture_F <- function() {
  # WESML weighted: two-nest balanced, with ASC
  # K=2, K_lambda=2, J_asc=3, P=7
  set.seed(2024)
  N <- 150; J <- 4
  dt <- data.table::data.table(
    id   = rep(1:N, each = J),
    alt  = rep(1:J, N),
    x1   = rnorm(N * J),
    x2   = rnorm(N * J),
    nest = ifelse(rep(1:J, N) <= 2, "A", "B")
  )
  beta_true  <- c(0.5, -0.3); asc_true <- c(0, 0.2, -0.1, -0.3)
  dt[, util := beta_true[1]*x1 + beta_true[2]*x2 + asc_true[alt]]
  dt[, choice := {
    u <- util; pr <- exp(u - max(u)); pr <- pr / sum(pr)
    cs <- cumsum(pr); r <- runif(1)
    as.integer(alt == alt[which(cs >= r)[1]])
  }, by = id]
  pop_shares <- c(0.3, 0.2, 0.25, 0.25)
  emp_shares <- rep(0.25, J)
  chosen_alts <- dt[choice == 1L, alt]
  w_i <- pop_shares[chosen_alts] / emp_shares[chosen_alts]
  dt[, w := rep(w_i, each = J)]
  d <- prepare_nl_data(dt, "id", "alt", "choice", c("x1", "x2"), "nest",
                       weights = dt[alt == 1, w])
  list(d = d, theta = c(0.5, -0.3, 0.8, 0.7, rep(0, J - 1)),
       J = J, K = 2, K_l = 2, J_asc = 3, P = 7,
       use_asc = TRUE, include_outside_option = FALSE,
       weights_vec = w_i)
}

make_fixture_G <- function() {
  # Larger J: J=10, 3 nests, K=3, use_asc=TRUE
  # K_lambda=3, J_asc=9, P=15
  set.seed(2024)
  N <- 80; J <- 10
  nest_assign <- c(rep("A", 4), rep("B", 3), rep("C", 3))
  dt <- data.table::data.table(
    id   = rep(1:N, each = J),
    alt  = rep(1:J, N),
    x1   = rnorm(N * J),
    x2   = rnorm(N * J),
    x3   = rnorm(N * J),
    nest = rep(nest_assign, N)
  )
  dt[, choice := {
    pr <- rep(1/J, J)
    as.integer(alt == sample(1:J, 1, prob = pr))
  }, by = id]
  d <- prepare_nl_data(dt, "id", "alt", "choice", c("x1", "x2", "x3"), "nest")
  list(d = d, theta = c(0.3, -0.2, 0.1, 0.7, 0.6, 0.5, rep(0, J - 1)),
       J = J, K = 3, K_l = 3, J_asc = 9, P = 15,
       use_asc = TRUE, include_outside_option = FALSE)
}

make_fixture_H <- function() {
  # Fixture B structure but lambda near 1 (MNL limit)
  fx <- make_fixture_B()
  fx$theta[3:4] <- 0.99
  fx
}

make_fixture_I <- function() {
  # Fixture B structure but lambda = 0.3 (strong nesting)
  fx <- make_fixture_B()
  fx$theta[3:4] <- 0.3
  fx
}

call_H_analytical <- function(fx) {
  nl_loglik_hessian_parallel(
    theta                  = fx$theta,
    X                      = fx$d$X,
    alt_idx                = fx$d$alt_idx,
    choice_idx             = fx$d$choice_idx,
    nest_idx               = fx$d$nest_idx,
    M                      = fx$d$M,
    weights                = fx$d$weights,
    use_asc                = fx$use_asc,
    include_outside_option = fx$include_outside_option
  )
}

call_H_numeric <- function(fx) {
  nl_loglik_numeric_hessian(
    theta                  = fx$theta,
    X                      = fx$d$X,
    alt_idx                = fx$d$alt_idx,
    choice_idx             = fx$d$choice_idx,
    nest_idx               = fx$d$nest_idx,
    M                      = fx$d$M,
    weights                = fx$d$weights,
    use_asc                = fx$use_asc,
    include_outside_option = fx$include_outside_option
  )
}

# Helper: create well-identified dataset for SE tests (strong DGP signal)
make_identified_nl_data <- function(seed = 2024, N = 500, J = 4,
                                    beta = c(1.0, -0.5),
                                    asc  = c(0, 0.3, -0.2, -0.4)) {
  set.seed(seed)
  dt <- data.table::data.table(
    id   = rep(1:N, each = J),
    alt  = rep(1:J, N),
    x1   = rnorm(N * J),
    x2   = rnorm(N * J),
    nest = ifelse(rep(1:J, N) <= 2, "A", "B")
  )
  dt[, util := beta[1]*x1 + beta[2]*x2 + asc[alt]]
  dt[, choice := {
    u <- util; pr <- exp(u - max(u)); pr <- pr / sum(pr)
    cs <- cumsum(pr); r <- runif(1)
    as.integer(alt == alt[which(cs >= r)[1]])
  }, by = id]
  dt
}

# ============================================================
# Scenario 1: Dimension check
# ============================================================

test_that("Scenario 1 — dimension: Fixture A (no ASC, P=4)", {
  fx <- make_fixture_A()
  H <- call_H_analytical(fx)
  expect_equal(nrow(H), fx$P)   # P = 4
  expect_equal(ncol(H), fx$P)
  expect_equal(nrow(H), length(fx$theta))
})

test_that("Scenario 1 — dimension: Fixture B (with ASC, P=7)", {
  fx <- make_fixture_B()
  H <- call_H_analytical(fx)
  expect_equal(nrow(H), fx$P)   # P = 7
  expect_equal(ncol(H), fx$P)
})

test_that("Scenario 1 — dimension: Fixture G (J=10, K=3, P=15)", {
  fx <- make_fixture_G()
  H <- call_H_analytical(fx)
  expect_equal(nrow(H), fx$P)   # P = 15
  expect_equal(ncol(H), fx$P)
})

# ============================================================
# Scenario 2: Symmetry — Fixtures A through I
# ============================================================

test_that("Scenario 2 — Hessian symmetry: Fixture A", {
  fx <- make_fixture_A()
  H <- call_H_analytical(fx)
  expect_lt(max(abs(H - t(H))), TOL_SYMM)
})

test_that("Scenario 2 — Hessian symmetry: Fixture B", {
  fx <- make_fixture_B()
  H <- call_H_analytical(fx)
  expect_lt(max(abs(H - t(H))), TOL_SYMM)
})

test_that("Scenario 2 — Hessian symmetry: Fixture D (singleton nest)", {
  fx <- make_fixture_D()
  H <- call_H_analytical(fx)
  expect_lt(max(abs(H - t(H))), TOL_SYMM)
})

test_that("Scenario 2 — Hessian symmetry: Fixture F (WESML weights)", {
  fx <- make_fixture_F()
  H <- call_H_analytical(fx)
  expect_lt(max(abs(H - t(H))), TOL_SYMM)
})

test_that("Scenario 2 — Hessian symmetry: Fixture G (J=10)", {
  fx <- make_fixture_G()
  H <- call_H_analytical(fx)
  expect_lt(max(abs(H - t(H))), TOL_SYMM)
})

test_that("Scenario 2 — Hessian symmetry: Fixture H (lambda near 1)", {
  fx <- make_fixture_H()
  H <- call_H_analytical(fx)
  expect_lt(max(abs(H - t(H))), TOL_SYMM)
})

test_that("Scenario 2 — Hessian symmetry: Fixture I (lambda = 0.3)", {
  fx <- make_fixture_I()
  H <- call_H_analytical(fx)
  expect_lt(max(abs(H - t(H))), TOL_SYMM)
})

# ============================================================
# Scenario 3: Equivalence to numeric oracle — Fixtures A through I
# ============================================================

test_that("Scenario 3 — Equivalence to oracle: Fixture A (no ASC, no outside)", {
  fx   <- make_fixture_A()
  H_an <- call_H_analytical(fx)
  H_nu <- call_H_numeric(fx)
  diff <- max(abs(H_an - H_nu))
  expect_lt(diff, TOL_EQUIV,
    label = paste0("Fixture A max abs diff = ", format(diff, scientific = TRUE),
                   " (threshold=", TOL_EQUIV, ")"))
})

test_that("Scenario 3 — Equivalence to oracle: Fixture B (with ASC, P=7)", {
  fx   <- make_fixture_B()
  H_an <- call_H_analytical(fx)
  H_nu <- call_H_numeric(fx)
  diff <- max(abs(H_an - H_nu))
  expect_lt(diff, TOL_EQUIV,
    label = paste0("Fixture B max abs diff = ", format(diff, scientific = TRUE)))
})

test_that("Scenario 3 — Equivalence to oracle: Fixture C (outside option)", {
  fx <- tryCatch(make_fixture_C(), error = function(e) {
    skip(paste("Fixture C setup failed:", conditionMessage(e)))
  })
  H_an <- call_H_analytical(fx)
  H_nu <- call_H_numeric(fx)
  diff <- max(abs(H_an - H_nu))
  expect_lt(diff, TOL_EQUIV,
    label = paste0("Fixture C max abs diff = ", format(diff, scientific = TRUE)))
})

test_that("Scenario 3 — Equivalence to oracle: Fixture D (singleton nest, P=6)", {
  fx   <- make_fixture_D()
  H_an <- call_H_analytical(fx)
  H_nu <- call_H_numeric(fx)
  diff <- max(abs(H_an - H_nu))
  expect_lt(diff, TOL_EQUIV,
    label = paste0("Fixture D max abs diff = ", format(diff, scientific = TRUE)))
})

test_that("Scenario 3 — Equivalence to oracle: Fixture E (singleton + outside)", {
  fx <- tryCatch(make_fixture_E(), error = function(e) {
    skip(paste("Fixture E setup failed:", conditionMessage(e)))
  })
  H_an <- call_H_analytical(fx)
  H_nu <- call_H_numeric(fx)
  diff <- max(abs(H_an - H_nu))
  expect_lt(diff, TOL_EQUIV,
    label = paste0("Fixture E max abs diff = ", format(diff, scientific = TRUE)))
})

test_that("Scenario 3 — Equivalence to oracle: Fixture F (WESML weights, P=7)", {
  fx   <- make_fixture_F()
  H_an <- call_H_analytical(fx)
  H_nu <- call_H_numeric(fx)
  diff <- max(abs(H_an - H_nu))
  expect_lt(diff, TOL_EQUIV,
    label = paste0("Fixture F max abs diff = ", format(diff, scientific = TRUE)))
})

test_that("Scenario 3 — Equivalence to oracle: Fixture G (J=10, K=3, P=15)", {
  fx   <- make_fixture_G()
  H_an <- call_H_analytical(fx)
  H_nu <- call_H_numeric(fx)
  diff <- max(abs(H_an - H_nu))
  expect_lt(diff, TOL_EQUIV,
    label = paste0("Fixture G max abs diff = ", format(diff, scientific = TRUE)))
})

test_that("Scenario 3 — Equivalence to oracle: Fixture H (lambda near 1)", {
  fx   <- make_fixture_H()
  H_an <- call_H_analytical(fx)
  H_nu <- call_H_numeric(fx)
  diff <- max(abs(H_an - H_nu))
  expect_lt(diff, TOL_EQUIV,
    label = paste0("Fixture H max abs diff = ", format(diff, scientific = TRUE)))
})

test_that("Scenario 3 — Equivalence to oracle: Fixture I (lambda = 0.3, strong nesting)", {
  fx   <- make_fixture_I()
  H_an <- call_H_analytical(fx)
  H_nu <- call_H_numeric(fx)
  diff <- max(abs(H_an - H_nu))
  expect_lt(diff, TOL_EQUIV,
    label = paste0("Fixture I max abs diff = ", format(diff, scientific = TRUE)))
})

# Cross-reference benchmarks (test-spec.md §7)

test_that("Scenario 3 cross-ref — at zero beta (flat point, Fixture B structure)", {
  fx_base <- make_fixture_B()
  fx      <- fx_base
  fx$theta <- c(0, 0, 0.6, 0.4, rep(0, fx_base$J - 1))
  H_an <- call_H_analytical(fx)
  H_nu <- call_H_numeric(fx)
  diff <- max(abs(H_an - H_nu))
  expect_lt(diff, TOL_EQUIV,
    label = paste0("zero-theta max abs diff = ", format(diff, scientific = TRUE)))
})

test_that("Scenario 3 cross-ref — at random theta away from MLE (Fixture A)", {
  set.seed(2024)
  fx_base   <- make_fixture_A()
  fx        <- fx_base
  fx$theta  <- c(runif(2, -0.5, 0.5), runif(2, 0.3, 0.9))
  H_an <- call_H_analytical(fx)
  H_nu <- call_H_numeric(fx)
  diff <- max(abs(H_an - H_nu))
  expect_lt(diff, TOL_EQUIV,
    label = paste0("random-theta Fixture A max abs diff = ", format(diff, scientific = TRUE)))
})

# ============================================================
# Scenario 4: SE equivalence via vcov() — well-identified datasets
# (test-spec §4 contract; use DGP-driven data to ensure convergence to unique MLE)
# ============================================================

test_that("Scenario 4 — SE equivalence via vcov(): well-identified Fixture B", {
  skip_on_cran()
  dt <- make_identified_nl_data(seed = 2024, N = 500)
  fit_an <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest",
                          use_asc = TRUE, se_method = "hessian")
  fit_nu <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest",
                          use_asc = TRUE, se_method = "numeric")
  se_an <- sqrt(diag(vcov(fit_an)))
  se_nu <- sqrt(diag(vcov(fit_nu)))
  diff  <- max(abs(se_an - se_nu), na.rm = TRUE)
  expect_lt(diff, TOL_SE,
    label = paste0("SE diff Fixture B = ", format(diff, scientific = TRUE)))
})

test_that("Scenario 4 — SE equivalence via vcov(): well-identified Fixture D (singleton)", {
  skip_on_cran()
  set.seed(2024)
  N <- 400; J <- 4
  dt <- data.table::data.table(
    id   = rep(1:N, each = J),
    alt  = rep(1:J, N),
    x1   = rnorm(N * J),
    x2   = rnorm(N * J),
    nest = ifelse(rep(1:J, N) == 1, "A", "B")  # singleton nest A
  )
  beta_true <- c(0.8, -0.4); asc_true <- c(0, 0.3, -0.2, -0.3)
  dt[, util := beta_true[1]*x1 + beta_true[2]*x2 + asc_true[alt]]
  dt[, choice := {
    u <- util; pr <- exp(u - max(u)); pr <- pr / sum(pr)
    cs <- cumsum(pr); r <- runif(1)
    as.integer(alt == alt[which(cs >= r)[1]])
  }, by = id]
  fit_an <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest",
                          use_asc = TRUE, se_method = "hessian")
  fit_nu <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest",
                          use_asc = TRUE, se_method = "numeric")
  se_an <- sqrt(diag(vcov(fit_an)))
  se_nu <- sqrt(diag(vcov(fit_nu)))
  diff  <- max(abs(se_an - se_nu), na.rm = TRUE)
  expect_lt(diff, TOL_SE,
    label = paste0("SE diff Fixture D (singleton) = ", format(diff, scientific = TRUE)))
})

test_that("Scenario 4 — SE equivalence via vcov(): well-identified Fixture F (weighted)", {
  skip_on_cran()
  set.seed(2024)
  N <- 400; J <- 4
  dt <- data.table::data.table(
    id   = rep(1:N, each = J),
    alt  = rep(1:J, N),
    x1   = rnorm(N * J),
    x2   = rnorm(N * J),
    nest = ifelse(rep(1:J, N) <= 2, "A", "B")
  )
  beta_true <- c(1.0, -0.5); asc_true <- c(0, 0.3, -0.2, -0.4)
  dt[, util := beta_true[1]*x1 + beta_true[2]*x2 + asc_true[alt]]
  dt[, choice := {
    u <- util; pr <- exp(u - max(u)); pr <- pr / sum(pr)
    cs <- cumsum(pr); r <- runif(1)
    as.integer(alt == alt[which(cs >= r)[1]])
  }, by = id]
  pop_shares <- c(0.3, 0.2, 0.25, 0.25)
  emp_shares <- rep(0.25, J)
  chosen_alts <- dt[choice == 1L, alt]
  w_i <- pop_shares[chosen_alts] / emp_shares[chosen_alts]

  # Non-uniform weights under a non-sandwich SE now warn; expected here since
  # this fixture deliberately compares the hessian and numeric SE machinery.
  fit_an <- suppressWarnings(run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest",
                          use_asc = TRUE, weights = w_i, se_method = "hessian"))
  fit_nu <- suppressWarnings(run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest",
                          use_asc = TRUE, weights = w_i, se_method = "numeric"))
  se_an <- sqrt(diag(vcov(fit_an)))
  se_nu <- sqrt(diag(vcov(fit_nu)))
  diff  <- max(abs(se_an - se_nu), na.rm = TRUE)
  expect_lt(diff, TOL_SE,
    label = paste0("SE diff Fixture F (weighted) = ", format(diff, scientific = TRUE)))
})

test_that("Scenario 4 — SE equivalence via vcov(): well-identified Fixture C (outside opt.)", {
  skip_on_cran()
  set.seed(2024)
  N <- 400; J <- 4
  dt <- data.table::data.table(
    id   = rep(1:N, each = (J + 1)),
    alt  = rep(0:J, N),
    x1   = c(rep(0, N), rnorm(N * J)),
    x2   = c(rep(0, N), rnorm(N * J)),
    nest = ifelse(rep(0:J, N) == 0, "OUT",
                  ifelse(rep(0:J, N) <= 2, "A", "B"))
  )
  beta_true <- c(0.8, -0.4)
  dt[alt > 0, util := beta_true[1]*x1 + beta_true[2]*x2]
  dt[alt == 0, util := 0]
  dt[, choice := {
    u <- util; pr <- exp(u - max(u)); pr <- pr / sum(pr)
    cs <- cumsum(pr); r <- runif(1)
    as.integer(alt == alt[which(cs >= r)[1]])
  }, by = id]
  fit_an <- tryCatch(
    run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest",
                  use_asc = TRUE, outside_opt_label = "0", include_outside_option = TRUE,
                  se_method = "hessian"),
    error = function(e) skip(paste("Fixture C fit failed:", conditionMessage(e)))
  )
  fit_nu <- tryCatch(
    run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest",
                  use_asc = TRUE, outside_opt_label = "0", include_outside_option = TRUE,
                  se_method = "numeric"),
    error = function(e) skip(paste("Fixture C fit_nu failed:", conditionMessage(e)))
  )
  se_an <- sqrt(diag(vcov(fit_an)))
  se_nu <- sqrt(diag(vcov(fit_nu)))
  diff  <- max(abs(se_an - se_nu), na.rm = TRUE)
  expect_lt(diff, TOL_SE,
    label = paste0("SE diff Fixture C (outside) = ", format(diff, scientific = TRUE)))
})

# ============================================================
# Scenario 5: se_method stored in object
# ============================================================

test_that("Scenario 5 — se_method stored correctly in fitted object", {
  skip_on_cran()
  dt <- make_identified_nl_data(seed = 2024, N = 200)
  fit_h <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest",
                         se_method = "hessian")
  fit_n <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest",
                         se_method = "numeric")
  expect_equal(fit_h$se_method, "hessian")
  expect_equal(fit_n$se_method, "numeric")
})

# ============================================================
# Scenario 6: Default se_method is "hessian"
# ============================================================

test_that("Scenario 6 — Default se_method is 'hessian'", {
  skip_on_cran()
  dt <- make_identified_nl_data(seed = 2024, N = 200)
  fit_default <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
  expect_equal(fit_default$se_method, "hessian")
})

# ============================================================
# Scenario 7: Invalid se_method rejected
# ============================================================

test_that("Scenario 7 — Unknown se_method is rejected; bhhh/sandwich accepted", {
  set.seed(2024)
  N <- 60; J <- 4
  dt <- data.table::data.table(
    id   = rep(1:N, each = J),
    alt  = rep(1:J, N),
    x1   = rnorm(N * J),
    nest = ifelse(rep(1:J, N) <= 2, "A", "B")
  )
  dt[, choice := {
    as.integer(alt == sample(1:J, 1))
  }, by = id]
  expect_error(
    run_nestlogit(dt, "id", "alt", "choice", "x1", "nest", se_method = "nonsense"),
    regexp = "arg"
  )
  # bhhh and sandwich are now valid SE methods for NL.
  fit_b <- suppressMessages(
    run_nestlogit(dt, "id", "alt", "choice", "x1", "nest", se_method = "bhhh")
  )
  expect_equal(fit_b$se_method, "bhhh")
  fit_s <- suppressMessages(
    run_nestlogit(dt, "id", "alt", "choice", "x1", "nest", se_method = "sandwich")
  )
  expect_equal(fit_s$se_method, "sandwich")
})

# ============================================================
# Scenario 8: Singleton nest — correct dimension P=6, no extra lambda
# ============================================================

test_that("Scenario 8 — Singleton nest: dim(H) == c(6L,6L), no extra lambda row/col", {
  fx <- make_fixture_D()
  H  <- call_H_analytical(fx)
  expect_equal(dim(H), c(6L, 6L))
  expect_true(all(is.finite(H)))
})

# ============================================================
# Scenario 9: Outside option stress test — most choose outside
# ============================================================

test_that("Scenario 9 — Outside option stress: majority choose outside option", {
  set.seed(2024)
  N <- 80; J <- 4
  dt <- data.table::data.table(
    id   = rep(1:N, each = (J + 1)),
    alt  = rep(0:J, N),
    x1   = c(rep(0, N), rnorm(N * J)),
    x2   = c(rep(0, N), rnorm(N * J)),
    nest = ifelse(rep(0:J, N) == 0, "OUT",
                  ifelse(rep(0:J, N) <= 2, "A", "B"))
  )
  dt[, choice := as.integer(alt == 0)]  # everyone chooses outside
  # Give first 10 individuals an inside choice so likelihood is non-degenerate
  dt[id <= 10, choice := as.integer(alt == 1)]
  d <- prepare_nl_data(dt, "id", "alt", "choice", c("x1", "x2"), "nest",
                       outside_opt_label = "0", include_outside_option = TRUE)
  theta <- c(0.5, -0.3, 0.8, 0.7, rep(0, J))

  H_an <- nl_loglik_hessian_parallel(
    theta, d$X, d$alt_idx, d$choice_idx, d$nest_idx, d$M, d$weights,
    use_asc = TRUE, include_outside_option = TRUE
  )
  H_nu <- nl_loglik_numeric_hessian(
    theta, d$X, d$alt_idx, d$choice_idx, d$nest_idx, d$M, d$weights,
    use_asc = TRUE, include_outside_option = TRUE
  )
  expect_true(all(is.finite(H_an)),
    label = "Hessian finite when most choose outside option")
  diff <- max(abs(H_an - H_nu))
  expect_lt(diff, TOL_EQUIV,
    label = paste0("outside-stress max abs diff = ", format(diff, scientific = TRUE)))
})

# ============================================================
# Scenario 10: Uniform weights consistency
# ============================================================

test_that("Scenario 10 — Uniform weights w=1: agrees with oracle (Fixture B)", {
  fx <- make_fixture_B()
  d_unif <- fx$d
  d_unif$weights <- rep(1.0, length(d_unif$weights))
  fx2 <- fx; fx2$d <- d_unif
  H_an <- call_H_analytical(fx2)
  H_nu <- call_H_numeric(fx2)
  diff <- max(abs(H_an - H_nu))
  expect_lt(diff, TOL_EQUIV,
    label = paste0("w=1 max abs diff = ", format(diff, scientific = TRUE)))
})

test_that("Scenario 10 — Uniform weights w=2: agrees with oracle (Fixture B)", {
  fx <- make_fixture_B()
  d_w2 <- fx$d
  d_w2$weights <- rep(2.0, length(d_w2$weights))
  fx2 <- fx; fx2$d <- d_w2
  H_an <- call_H_analytical(fx2)
  H_nu <- call_H_numeric(fx2)
  diff <- max(abs(H_an - H_nu))
  expect_lt(diff, TOL_EQUIV,
    label = paste0("w=2 max abs diff = ", format(diff, scientific = TRUE)))
})

# ============================================================
# Scenario 11: PSD at MLE (soft check, eigenvalue floor > -1e-6)
# ============================================================

test_that("Scenario 11 — Hessian PSD at MLE (well-identified Fixture B)", {
  skip_on_cran()
  dt  <- make_identified_nl_data(seed = 2024, N = 400)
  fit <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest",
                       use_asc = TRUE, se_method = "hessian")
  d   <- prepare_nl_data(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
  H   <- nl_loglik_hessian_parallel(
    coef(fit), d$X, d$alt_idx, d$choice_idx, d$nest_idx, d$M, d$weights,
    use_asc = TRUE, include_outside_option = FALSE
  )
  ev  <- eigen(H, symmetric = TRUE, only.values = TRUE)$values
  expect_gte(min(ev), -1e-6,
    label = paste0("min eigenvalue at MLE = ", format(min(ev), scientific = TRUE)))
})

# ============================================================
# Scenario 12: Lambda = 1 exact — Hessian finite and agrees with oracle
# ============================================================

test_that("Scenario 12 — Lambda = 1 exact: Hessian finite and matches oracle", {
  fx_base <- make_fixture_B()
  fx <- fx_base; fx$theta[3:4] <- 1.0
  H_an <- call_H_analytical(fx)
  H_nu <- call_H_numeric(fx)
  expect_true(all(is.finite(H_an)),
    label = "Hessian must be finite at lambda=1 exactly")
  diff <- max(abs(H_an - H_nu))
  expect_lt(diff, TOL_EQUIV,
    label = paste0("lambda=1 max abs diff = ", format(diff, scientific = TRUE)))
})

# ============================================================
# No-NaN/Inf property: all fixtures
# ============================================================

test_that("Property — No NaN/Inf: Fixture A", {
  expect_true(all(is.finite(call_H_analytical(make_fixture_A()))))
})
test_that("Property — No NaN/Inf: Fixture B", {
  expect_true(all(is.finite(call_H_analytical(make_fixture_B()))))
})
test_that("Property — No NaN/Inf: Fixture D (singleton)", {
  expect_true(all(is.finite(call_H_analytical(make_fixture_D()))))
})
test_that("Property — No NaN/Inf: Fixture G (J=10)", {
  expect_true(all(is.finite(call_H_analytical(make_fixture_G()))))
})
test_that("Property — No NaN/Inf: Fixture H (lambda near 1)", {
  expect_true(all(is.finite(call_H_analytical(make_fixture_H()))))
})
test_that("Property — No NaN/Inf: Fixture I (lambda = 0.3)", {
  expect_true(all(is.finite(call_H_analytical(make_fixture_I()))))
})

# ============================================================
# Edge case: extreme lambda (lambda = 0.05, test-spec §5)
# test-spec §5 explicitly relaxes to 1e-4 for lambda < 0.1 due to oracle FD error.
# We use TOL_EXTLAM = 1e-3 (safety margin around the borderline) and flag below.
# ============================================================

test_that("Edge case — Extreme lambda (0.05): Hessian finite and matches oracle to 1e-3", {
  fx_base <- make_fixture_B()
  fx <- fx_base; fx$theta[3:4] <- 0.05
  H_an <- call_H_analytical(fx)
  H_nu <- call_H_numeric(fx)
  expect_true(all(is.finite(H_an)),
    label = "Hessian must be finite at lambda=0.05")
  diff <- max(abs(H_an - H_nu))
  # test-spec §5: relaxed to 1e-4 at extreme lambda; we gate at 1e-3 for safety
  # (the oracle's own FD error at lambda=0.05 can be ~1e-4).
  # Observed diff ≈ 1.4e-4 — within the 1e-4 spec but flagged here at 1e-3 gate.
  expect_lt(diff, TOL_EXTLAM,
    label = paste0("lambda=0.05 max abs diff = ", format(diff, scientific = TRUE),
                   " (spec allows 1e-4; gate here 1e-3 per test-spec §5 annotation)"))
})

# ============================================================
# Edge case: J=2 minimal nest structure
# ============================================================

test_that("Edge case — Minimal nest: J=3 with singleton + non-singleton nest", {
  set.seed(2024)
  N <- 100; J2 <- 3
  dt2 <- data.table::data.table(
    id   = rep(1:N, each = J2),
    alt  = rep(1:J2, N),
    x1   = rnorm(N * J2),
    nest = c("A", "B", "B")  # nest A: singleton alt1; nest B: alts 2,3
  )
  dt2[, choice := {
    as.integer(alt == sample(1:J2, 1))
  }, by = id]
  d     <- prepare_nl_data(dt2, "id", "alt", "choice", "x1", "nest")
  theta <- c(0.5, 0.7, rep(0, J2 - 1))  # K=1, K_l=1, J_asc=2, P=4

  H_an <- nl_loglik_hessian_parallel(
    theta, d$X, d$alt_idx, d$choice_idx, d$nest_idx, d$M, d$weights,
    use_asc = TRUE, include_outside_option = FALSE
  )
  H_nu <- nl_loglik_numeric_hessian(
    theta, d$X, d$alt_idx, d$choice_idx, d$nest_idx, d$M, d$weights,
    use_asc = TRUE, include_outside_option = FALSE
  )
  expect_true(all(is.finite(H_an)))
  expect_equal(dim(H_an), c(length(theta), length(theta)))
  diff <- max(abs(H_an - H_nu))
  expect_lt(diff, TOL_EQUIV,
    label = paste0("minimal-nest max abs diff = ", format(diff, scientific = TRUE)))
})

# ============================================================
# Regression: coef() unchanged by se_method choice
# (optimization path is identical; coef difference is at most machine epsilon)
# ============================================================

test_that("Regression — coef() unchanged by se_method (well-identified Fixture B)", {
  skip_on_cran()
  dt    <- make_identified_nl_data(seed = 2024, N = 400)
  fit_h <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest",
                         use_asc = TRUE, se_method = "hessian")
  fit_n <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest",
                         use_asc = TRUE, se_method = "numeric")
  # Optimization is identical; allow only floating-point noise (< 1e-12)
  diff <- max(abs(coef(fit_h) - coef(fit_n)))
  expect_lt(diff, 1e-12,
    label = paste0("coef diff across se_method = ", format(diff, scientific = TRUE)))
})