## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
have_ggplot2 <- requireNamespace("ggplot2", quietly = TRUE)
points <- function(x) sprintf("%.0f", 100 * x)

## ----library, message = FALSE-------------------------------------------------
library(foundryR)
library(dplyr)

## ----simulate-----------------------------------------------------------------
set.seed(20260927)
n_cases <- 120
tickets <- tibble(
  case_id = seq_len(n_cases),
  segment = rep(c("billing", "technical", "account"), times = c(50, 40, 30)),
  difficulty = rnorm(n_cases)
)

simulate_run <- function(tickets, skill, technical_bonus = 0) {
  shift <- c(billing = 0.6, technical = -0.5, account = 0.2)[tickets$segment]
  bonus <- ifelse(tickets$segment == "technical", technical_bonus, 0)
  p_pass <- plogis(skill - 1.2 * tickets$difficulty + shift + bonus)
  tickets |>
    mutate(.grader = "label-match", .passed = runif(n()) < p_pass, p_true = p_pass) |>
    select(-difficulty)
}

simulated <- bind_rows(
  "model-a" = simulate_run(tickets, skill = 1.0),
  "model-b" = simulate_run(tickets, skill = 1.3, technical_bonus = 0.5),
  .id = "target"
)
expected <- simulated |>
  group_by(target, segment) |>
  summarise(expected_pass_rate = mean(p_true), .groups = "drop")
results <- select(simulated, -p_true)
results

## ----pass-rates---------------------------------------------------------------
wilson_interval <- function(passed, n, level = 0.95) {
  z <- qnorm(1 - (1 - level) / 2)
  p <- passed / n
  center <- (p + z^2 / (2 * n)) / (1 + z^2 / n)
  half <- z * sqrt(p * (1 - p) / n + z^2 / (4 * n^2)) / (1 + z^2 / n)
  tibble(lower = center - half, upper = center + half)
}

pass_rates <- results |>
  group_by(target) |>
  summarise(passed = sum(.passed), cases = n(), .groups = "drop") |>
  mutate(pass_rate = passed / cases, wilson_interval(passed, cases))

pass_rates

## ----pass-rate-widths, include = FALSE----------------------------------------
independent_half_width <- {
  p <- mean(pass_rates$pass_rate)
  qnorm(0.975) * sqrt(2 * p * (1 - p) / n_cases)
}

## ----paired-------------------------------------------------------------------
paired <- inner_join(
  results |> filter(target == "model-a") |> select(case_id, segment, a = .passed),
  results |> filter(target == "model-b") |> select(case_id, b = .passed),
  by = "case_id"
)
count(paired, a, b)

only_b <- sum(!paired$a & paired$b)
only_a <- sum(paired$a & !paired$b)
mcnemar <- binom.test(only_b, only_b + only_a)
mcnemar$p.value

## ----paired-interval----------------------------------------------------------
difference <- mean(paired$b) - mean(paired$a)
bootstrap <- replicate(2000, {
  rows <- sample.int(nrow(paired), replace = TRUE)
  mean(paired$b[rows]) - mean(paired$a[rows])
})
paired_ci <- quantile(bootstrap, c(0.025, 0.975), names = FALSE)
unpaired_ci <- prop.test(
  c(sum(paired$b), sum(paired$a)),
  c(nrow(paired), nrow(paired)),
  correct = FALSE
)$conf.int
paired_wald_ci <- difference + c(-1, 1) * qnorm(0.975) *
  sqrt((only_a + only_b - (only_b - only_a)^2 / nrow(paired)) / nrow(paired)^2)

tibble(
  method = c(
    "paired bootstrap over tickets",
    "paired Wald interval",
    "unpaired two-sample interval"
  ),
  difference = difference,
  lower = c(paired_ci[[1]], paired_wald_ci[[1]], unpaired_ci[[1]]),
  upper = c(paired_ci[[2]], paired_wald_ci[[2]], unpaired_ci[[2]])
)

## ----true-difference, include = FALSE-----------------------------------------
expected_overall <- simulated |>
  group_by(target) |>
  summarise(expected_pass_rate = mean(p_true), .groups = "drop")
true_difference <- diff(expected_overall$expected_pass_rate[match(c("model-a", "model-b"), expected_overall$target)])

## ----by-segment---------------------------------------------------------------
by_segment <- results |>
  group_by(target, segment) |>
  summarise(passed = sum(.passed), cases = n(), .groups = "drop") |>
  mutate(pass_rate = passed / cases, wilson_interval(passed, cases))

by_segment

## ----segment-reversals, include = FALSE---------------------------------------
segment_gaps <- inner_join(by_segment, expected, by = c("target", "segment")) |>
  group_by(segment) |>
  summarise(
    observed_gap = pass_rate[target == "model-b"] - pass_rate[target == "model-a"],
    expected_gap = expected_pass_rate[target == "model-b"] - expected_pass_rate[target == "model-a"],
    .groups = "drop"
  )
reversed <- segment_gaps[sign(segment_gaps$observed_gap) != sign(segment_gaps$expected_gap), ]
reversal_sentence <- if (nrow(reversed) > 0) {
  sprintf(
    paste(
      "The %s segment shows how easily that happens. By construction `model-b` is",
      "better on %s tickets, by %s points in expectation, yet this sample shows it",
      "%s points lower."
    ),
    reversed$segment[[1]],
    reversed$segment[[1]],
    points(reversed$expected_gap[[1]]),
    points(abs(reversed$observed_gap[[1]]))
  )
} else {
  paste(
    "Here every observed segment gap points the same way as the expected one,",
    "but the intervals are wide enough that a different sample could reverse some of them."
  )
}

## ----segment-chart-data, include = FALSE, eval = have_ggplot2-----------------
segment_levels <- c("account", "technical", "billing", "All tickets")
chart_data <- bind_rows(
  mutate(pass_rates, segment = "All tickets"),
  by_segment
) |>
  mutate(
    segment = factor(segment, levels = segment_levels),
    y = as.integer(segment) + ifelse(target == "model-a", 0.16, -0.16)
  )
stopifnot(!anyNA(chart_data$segment))

series_colors <- c("model-a" = "#065D92", "model-b" = "#AD4D0E")

## ----segment-chart-title, include = FALSE, eval = have_ggplot2----------------
seg_wide <- inner_join(
  filter(by_segment, target == "model-a") |> select(segment, a_lower = lower, a_upper = upper),
  filter(by_segment, target == "model-b") |> select(segment, b_lower = lower, b_upper = upper),
  by = "segment"
)
separated <- seg_wide$segment[seg_wide$b_lower > seg_wide$a_upper | seg_wide$a_lower > seg_wide$b_upper]
gap_cases <- sum(paired$b) - sum(paired$a)
chart_title <- sprintf(
  "On %d tickets, model-b passed %d %s than model-a (paired 95%% interval %s to %s points)",
  n_cases,
  abs(gap_cases),
  if (gap_cases >= 0) "more" else "fewer",
  points(paired_ci[[1]]),
  points(paired_ci[[2]])
)
describe_row <- function(label) {
  rows <- chart_data[chart_data$segment == label, ]
  paste(
    sprintf("%s %s%% (%s to %s)", rows$target, points(rows$pass_rate), points(rows$lower), points(rows$upper)),
    collapse = ", "
  )
}
chart_alt <- paste0(
  chart_title, ". Dot and interval chart of pass rates with 95% Wilson intervals. ",
  paste(
    vapply(rev(segment_levels), function(label) paste0(label, ": ", describe_row(label)), character(1)),
    collapse = "; "
  ),
  "."
)
x_floor <- floor(min(chart_data$lower) * 10) / 10

## ----segment-chart, echo = FALSE, eval = have_ggplot2, fig.width = 7, fig.height = 3.9, out.width = "100%", fig.alt = get0("chart_alt", ifnotfound = "Pass rates with 95% intervals by segment for two targets.")----
top_row <- chart_data[chart_data$segment == "All tickets", ]
wrap_title <- function(text, width) paste(strwrap(text, width), collapse = "\n")

ggplot2::ggplot(chart_data, ggplot2::aes(y = y, colour = target)) +
  ggplot2::geom_segment(
    ggplot2::aes(x = lower, xend = upper, yend = y),
    linewidth = 0.9
  ) +
  ggplot2::geom_point(ggplot2::aes(x = pass_rate, shape = target), size = 2.8) +
  ggplot2::geom_text(
    data = top_row,
    ggplot2::aes(x = upper, label = target),
    hjust = -0.15,
    size = 3.6
  ) +
  ggplot2::scale_colour_manual(values = series_colors, guide = "none") +
  ggplot2::scale_shape_manual(values = c("model-a" = 16, "model-b" = 17), guide = "none") +
  ggplot2::scale_x_continuous(
    breaks = seq(x_floor, 1, by = 0.1),
    labels = function(x) paste0(round(100 * x), "%"),
    expand = ggplot2::expansion(mult = c(0.02, 0.02))
  ) +
  ggplot2::scale_y_continuous(
    breaks = seq_along(segment_levels),
    labels = segment_levels,
    expand = ggplot2::expansion(add = 0.45)
  ) +
  ggplot2::coord_cartesian(xlim = c(x_floor, 1.08), clip = "off") +
  ggplot2::labs(
    title = wrap_title(chart_title, 70),
    subtitle = "Pass rate with 95% Wilson interval, by ticket segment",
    x = "Pass rate",
    y = NULL,
    caption = wrap_title(
      paste(
        "Simulated results for 120 tickets. Segment intervals use only that segment's tickets;",
        "comparing many segments invites false positives. Overlap between two intervals is not a test of their difference."
      ),
      105
    )
  ) +
  ggplot2::theme_minimal(base_size = 12) +
  ggplot2::theme(
    plot.title = ggplot2::element_text(face = "bold"),
    plot.title.position = "plot",
    plot.subtitle = ggplot2::element_text(color = "grey30"),
    plot.caption = ggplot2::element_text(color = "grey40", hjust = 0),
    plot.caption.position = "plot",
    panel.grid.major.y = ggplot2::element_blank(),
    panel.grid.minor = ggplot2::element_blank(),
    axis.text.y = ggplot2::element_text(size = 11, color = "grey15"),
    plot.margin = ggplot2::margin(8, 40, 8, 8)
  )

## ----margins------------------------------------------------------------------
margin_of_error <- function(p, n) qnorm(0.975) * sqrt(p * (1 - p) / n)

tibble(cases = c(25, 50, 100, 200, 400, 800)) |>
  mutate(plus_or_minus_points = round(100 * margin_of_error(0.8, cases), 1))

## ----paired-sample-size-------------------------------------------------------
paired_cases_needed <- function(d, psi, alpha = 0.05, power = 0.8) {
  z_alpha <- qnorm(1 - alpha / 2)
  z_beta <- qnorm(power)
  ceiling((z_alpha * sqrt(psi) + z_beta * sqrt(psi - d^2))^2 / d^2)
}

observed_psi <- (only_a + only_b) / nrow(paired)
psi_ci <- wilson_interval(only_a + only_b, nrow(paired))
tibble(d = c(0.03, 0.05, 0.10)) |>
  filter(d <= observed_psi) |>
  mutate(cases_needed = paired_cases_needed(d, psi = observed_psi))

## ----judge-validation---------------------------------------------------------
validation <- tibble(human_pass = runif(80) < 0.7) |>
  mutate(
    judge_pass = ifelse(human_pass, runif(n()) < 0.92, runif(n()) < 0.30),
    human = ifelse(human_pass, "pass", "fail"),
    judge = ifelse(judge_pass, "pass", "fail")
  )

foundry_agreement(validation, estimate = "judge", truth = "human")
count(validation, human, judge)

## ----judge-counts, include = FALSE--------------------------------------------
false_pass <- sum(validation$judge == "pass" & validation$human == "fail")
false_fail <- sum(validation$judge == "fail" & validation$human == "pass")
human_passes <- sum(validation$human == "pass")
human_fails <- sum(validation$human == "fail")
judge_hits <- sum(validation$judge == "pass" & validation$human == "pass")
judge_correct_fails <- sum(validation$judge == "fail" & validation$human == "fail")
sensitivity <- judge_hits / human_passes
specificity <- judge_correct_fails / human_fails
sensitivity_ci <- wilson_interval(judge_hits, human_passes)
specificity_ci <- wilson_interval(judge_correct_fails, human_fails)
attenuation <- sensitivity + specificity - 1
judge_rate <- points(mean(validation$judge == "pass"))
human_rate <- points(mean(validation$human == "pass"))
judge_bias <- if (false_pass == false_fail) {
  sprintf(
    paste(
      "The errors cancel in the aggregate, and the judge and people both report a",
      "%s%% pass rate. Matching pass rates are not evidence of a good judge.",
      "The two kinds of error can stop cancelling when the mix of responses",
      "changes, for example between two targets, and then the judge distorts the",
      "comparison."
    ),
    judge_rate
  )
} else {
  sprintf(
    "The judge reports a %s%% pass rate against %s%% from people, so it makes targets look %s than they are.",
    judge_rate,
    human_rate,
    if (false_pass > false_fail) "better" else "worse"
  )
}

## ----cost-latency-------------------------------------------------------------
runs <- tibble(
  target = c("model-a", "model-b"),
  target_latency_p50_ms = c(640, 1420),
  target_latency_p95_ms = c(1900, 4100),
  target_cost = c(0.018, 0.071),
  target_cost_currency = "USD"
)

runs |>
  left_join(select(pass_rates, target, passed, cases), by = "target") |>
  mutate(cost_per_1000_cases = 1000 * target_cost / cases) |>
  select(target, target_latency_p95_ms, target_cost, passed, cost_per_1000_cases)

## ----cost-increment, include = FALSE------------------------------------------
extra_passes <- diff(pass_rates$passed[match(c("model-a", "model-b"), pass_rates$target)])
extra_cost <- diff(runs$target_cost)