args <- commandArgs(trailingOnly = FALSE)
file_arg <- grep("^--file=", args, value = TRUE)
analysis_dir <- if (length(file_arg)) dirname(normalizePath(sub("^--file=", "", file_arg[[1]]), mustWork = FALSE)) else file.path(getwd(), "analysis")
source(file.path(analysis_dir, "00_config.R"))
source(file.path(analysis_dir, "lib_timss.R"))

bind_rows <- function(rows) {
  rows <- Filter(Negate(is.null), rows)
  if (!length(rows)) return(data.frame())
  output <- do.call(rbind, rows)
  rownames(output) <- NULL
  output
}

core_components <- function(math_t0, math_t1, science_t0, science_t1, threshold) {
  state_t0 <- achievement_state(math_t0, science_t0, threshold)
  state_t1 <- achievement_state(math_t1, science_t1, threshold)
  list(
    baseline_bb = state_t0 == "BB",
    followup_aa = state_t1 == "AA",
    followup_bb = state_t1 == "BB",
    weakest_change = pmin(math_t1, science_t1) - pmin(math_t0, science_t0),
    math_change = math_t1 - math_t0,
    science_change = science_t1 - science_t0
  )
}

evaluate_core <- function(components, weight) {
  risk <- components$baseline_bb
  c(
    BB_to_AA = weighted_probability(components$followup_aa[risk], weight[risk]),
    BB_to_BB = weighted_probability(components$followup_bb[risk], weight[risk]),
    bb_exit = weighted_probability(!components$followup_bb[risk], weight[risk]),
    baseline_BB = weighted_probability(risk, weight),
    weakest_margin_change_BB = weighted_mean(components$weakest_change[risk], weight[risk]),
    math_change_BB = weighted_mean(components$math_change[risk], weight[risk]),
    science_change_BB = weighted_mean(components$science_change[risk], weight[risk])
  )
}

run_pv_analysis <- function(data, band, systems, estimator_builder, label) {
  observed_systems <- unique(data$system)
  if (!setequal(observed_systems, systems)) {
    stop(label, ": configured systems do not match data")
  }
  results <- setNames(vector("list", length(systems)), systems)
  for (system in systems) {
    message("  ", label, ": ", band, " / ", system)
    system_data <- data[data$system == system, , drop = FALSE]
    draw_results <- lapply(seq_len(config$pv_count), function(draw) {
      estimator <- estimator_builder(system_data, pv_variable_map(band, draw), draw)
      result <- jrr_estimate(system_data, estimator)
      result$replicates <- NULL
      result
    })
    results[[system]] <- combine_pv_results(draw_results)
  }
  results
}

result_row <- function(result, parameter, ...) {
  data.frame(
    ...,
    parameter = parameter,
    estimate = unname(result$estimate[parameter]),
    standard_error = unname(result$se[parameter]),
    ci_lower = unname(result$estimate[parameter] - 1.96 * result$se[parameter]),
    ci_upper = unname(result$estimate[parameter] + 1.96 * result$se[parameter]),
    stringsAsFactors = FALSE
  )
}

# Alternative thresholds and buffered (hysteresis) transitions -----------------
threshold_metrics <- c(
  "BB_to_AA", "BB_to_BB", "bb_exit", "baseline_BB",
  "weakest_margin_change_BB", "math_change_BB", "science_change_BB"
)
hysteresis_specs <- list(
  primary_450_to_500 = unname(config$primary_hysteresis),
  low_380_to_420 = unname(config$low_hysteresis)
)

threshold_hysteresis_builder <- function(system_data, map, draw) {
  threshold_components <- lapply(config$threshold_grid, function(threshold) {
    core_components(
      system_data[[map$math_t0]], system_data[[map$math_t1]],
      system_data[[map$science_t0]], system_data[[map$science_t1]], threshold
    )
  })
  names(threshold_components) <- as.character(config$threshold_grid)

  hysteresis_components <- lapply(hysteresis_specs, function(cutoffs) {
    baseline_risk <- system_data[[map$math_t0]] < cutoffs[1] & system_data[[map$science_t0]] < cutoffs[1]
    followup_recovery <- system_data[[map$math_t1]] >= cutoffs[2] & system_data[[map$science_t1]] >= cutoffs[2]
    followup_dual_deficit <- system_data[[map$math_t1]] < cutoffs[2] & system_data[[map$science_t1]] < cutoffs[2]
    list(
      risk = baseline_risk,
      recovery = followup_recovery,
      dual_deficit = followup_dual_deficit,
      weakest_change = pmin(system_data[[map$math_t1]], system_data[[map$science_t1]]) -
        pmin(system_data[[map$math_t0]], system_data[[map$science_t0]]),
      math_change = system_data[[map$math_t1]] - system_data[[map$math_t0]],
      science_change = system_data[[map$science_t1]] - system_data[[map$science_t0]]
    )
  })

  function(dataset, weight) {
    values <- numeric(0)
    for (threshold in names(threshold_components)) {
      metrics <- evaluate_core(threshold_components[[threshold]], weight)
      names(metrics) <- paste0("t", threshold, "__", names(metrics))
      values <- c(values, metrics)
    }
    for (specification in names(hysteresis_components)) {
      component <- hysteresis_components[[specification]]
      risk <- component$risk
      metrics <- c(
        dual_recovery = weighted_probability(component$recovery[risk], weight[risk]),
        dual_persistence = weighted_probability(component$dual_deficit[risk], weight[risk]),
        baseline_risk = weighted_probability(risk, weight),
        weakest_margin_change = weighted_mean(component$weakest_change[risk], weight[risk]),
        math_change = weighted_mean(component$math_change[risk], weight[risk]),
        science_change = weighted_mean(component$science_change[risk], weight[risk])
      )
      names(metrics) <- paste0("h_", specification, "__", names(metrics))
      values <- c(values, metrics)
    }
    values
  }
}

risk_count_summaries <- function(data, band, systems) {
  threshold_rows <- list()
  hysteresis_rows <- list()
  threshold_row <- hysteresis_row <- 1L
  for (system in systems) {
    system_data <- data[data$system == system, , drop = FALSE]
    for (threshold in config$threshold_grid) {
      counts <- vapply(seq_len(config$pv_count), function(draw) {
        map <- pv_variable_map(band, draw)
        sum(system_data[[map$math_t0]] < threshold & system_data[[map$science_t0]] < threshold)
      }, integer(1))
      threshold_rows[[threshold_row]] <- data.frame(
        band = band, system = system, threshold = threshold,
        risk_n_mean = mean(counts), risk_n_min = min(counts), risk_n_max = max(counts),
        stringsAsFactors = FALSE
      )
      threshold_row <- threshold_row + 1L
    }
    for (specification in names(hysteresis_specs)) {
      cutoffs <- hysteresis_specs[[specification]]
      counts <- vapply(seq_len(config$pv_count), function(draw) {
        map <- pv_variable_map(band, draw)
        sum(system_data[[map$math_t0]] < cutoffs[1] & system_data[[map$science_t0]] < cutoffs[1])
      }, integer(1))
      hysteresis_rows[[hysteresis_row]] <- data.frame(
        band = band, system = system, specification = specification,
        baseline_cutoff = cutoffs[1], followup_cutoff = cutoffs[2],
        risk_n_mean = mean(counts), risk_n_min = min(counts), risk_n_max = max(counts),
        stringsAsFactors = FALSE
      )
      hysteresis_row <- hysteresis_row + 1L
    }
  }
  list(threshold = bind_rows(threshold_rows), hysteresis = bind_rows(hysteresis_rows))
}

# Gate the complete standardization producer bundle before any expensive sensitivity work.
read_standardization_csv <- function(directory, filename) {
  path <- file.path(directory, filename)
  if (!file.exists(path)) stop("Required standardization artifact is missing: ", path)
  utils::read.csv(path, check.names = FALSE, stringsAsFactors = FALSE)
}

exact_schema <- function(data, expected, label) {
  identical(names(data), expected) || stop(
    label, " schema mismatch; expected ", paste(expected, collapse = ", "),
    "; observed ", paste(names(data), collapse = ", ")
  )
  TRUE
}

row_keys <- function(data, columns) {
  do.call(paste, c(lapply(data[columns], as.character), sep = "\r"))
}

exact_ordered_keys <- function(data, expected, columns, label) {
  observed_keys <- row_keys(data, columns)
  expected_keys <- row_keys(expected, columns)
  if (anyDuplicated(observed_keys) || !identical(observed_keys, expected_keys)) {
    stop(label, " keys are duplicated, missing, extra, or out of order")
  }
  TRUE
}

frames_agree <- function(left, right, label, tolerance = 1e-10) {
  if (!identical(names(left), names(right)) || nrow(left) != nrow(right)) {
    stop(label, " RDS/CSV dimensions or columns disagree")
  }
  for (column in names(left)) {
    lhs <- left[[column]]
    rhs <- right[[column]]
    if (is.numeric(lhs) && is.numeric(rhs)) {
      matching_na <- is.na(lhs) & is.na(rhs)
      finite <- is.finite(lhs) & is.finite(rhs)
      if (!all(matching_na | (finite & abs(lhs - rhs) <= tolerance))) {
        stop(label, " RDS/CSV values disagree in ", column)
      }
    } else if (!identical(lhs, rhs)) {
      stop(label, " RDS/CSV values disagree in ", column)
    }
  }
  TRUE
}

standardization_specification_systems <- function(specification) {
  system_set <- standardization_specifications[[specification]]$systems
  if (identical(system_set, "all")) return(config$g45_systems)
  if (identical(system_set, "without_korea")) return(setdiff(config$g45_systems, "Korea, Rep. of"))
  stop("Unknown configured standardization system set: ", system_set)
}

expected_standardization_keys <- bind_rows(lapply(names(standardization_specifications), function(specification) {
  data.frame(
    specification = specification,
    system = standardization_specification_systems(specification),
    stringsAsFactors = FALSE
  )
}))
expected_alternative_keys <- expected_standardization_keys[
  expected_standardization_keys$specification != "primary", , drop = FALSE
]
expected_support_keys <- bind_rows(lapply(seq_len(config$pv_count), function(draw) {
  transform(expected_standardization_keys, draw = draw)[, c("draw", "specification", "system")]
}))
expected_surfaces <- unique(vapply(standardization_specifications, `[[`, character(1), "surface"))
expected_model_keys <- bind_rows(lapply(seq_len(config$pv_count), function(draw) {
  bind_rows(lapply(expected_surfaces, function(surface) {
    data.frame(draw = draw, surface = surface, system = config$g45_systems, stringsAsFactors = FALSE)
  }))
}))
expected_replicate_keys <- bind_rows(lapply(seq_len(config$pv_count), function(draw) {
  data.frame(draw = draw, specification = names(standardization_specifications), stringsAsFactors = FALSE)
}))

standardization_bundle <- readRDS(file.path(paths$derived, "standardization_results.rds"))
standardization_specification_results <- read_standardization_csv(
  paths$tables, "standardization_specification_results.csv"
)
standardization_alternative_contrasts <- read_standardization_csv(
  paths$tables, "standardization_alternative_contrasts.csv"
)
standardization_specification_diagnostics <- read_standardization_csv(
  paths$diagnostics, "standardization_specification_diagnostics_by_draw.csv"
)
standardization_model_diagnostics <- read_standardization_csv(
  paths$diagnostics, "standardization_model_diagnostics_by_draw.csv"
)
standardization_replicate_diagnostics <- read_standardization_csv(
  paths$diagnostics, "standardization_replicate_diagnostics.csv"
)
standardization_producer_validation <- read_standardization_csv(
  paths$diagnostics, "standardization_covariance_validation.csv"
)

standardization_schemas <- list(
  specification_results = c(
    "specification", "system", "surface", "target", "support",
    "raw_full", "raw_full_se", "raw_support", "raw_support_se",
    "model_own_support", "model_own_support_se", "standardized_common",
    "standardized_common_se", "support_restriction_component", "support_restriction_se",
    "model_calibration_component", "model_calibration_se", "target_distribution_component",
    "target_distribution_se", "standardized_minus_raw_full",
    "standardized_minus_raw_full_se", "total_ci_lower", "total_ci_upper",
    "decomposition_residual"
  ),
  alternative_contrasts = c(
    "specification", "system", "alternative", "primary", "alternative_minus_primary",
    "standard_error", "ci_lower", "ci_upper"
  ),
  specification_diagnostics = c(
    "draw", "specification", "system", "surface", "target", "support", "risk_n",
    "supported_n", "support_retained_weight", "support_effective_n", "target_system_mass",
    "target_weight_sum", "target_mass_retained_before_renormalization", "math_lower",
    "math_upper", "science_lower", "science_upper", "grid_cells_retained",
    "grid_cells_total", "minimum_common_cell_n", "minimum_common_cell_ess"
  ),
  model_diagnostics = c(
    "draw", "surface", "system", "risk_n", "event_n", "design_columns", "model_rank",
    "condition_number", "model_iterations", "converged", "warning_count",
    "maximum_absolute_coefficient", "weighted_observed_rate", "weighted_fitted_rate",
    "weighted_calibration_difference", "weighted_brier_score", "linear_predictor_minimum",
    "linear_predictor_maximum", "predicted_probability_minimum",
    "predicted_probability_maximum", "pooled_clipped_prediction_n",
    "pooled_clipped_prediction_weighted_share", "own_system_clipped_prediction_n",
    "own_system_clipped_prediction_weighted_share", "common_target_clipped_prediction_n",
    "common_target_clipped_prediction_weighted_share"
  ),
  replicate_diagnostics = c(
    "draw", "specification", "surface", "stacked_replicates_expected",
    "stacked_replicates_attempted", "stacked_replicates_completed",
    "affected_replicates_expected", "affected_replicates_attempted",
    "affected_replicates_completed", "support_rebuilds_expected",
    "support_rebuilds_attempted", "support_rebuilds_completed",
    "shared_surface_refits_expected", "shared_surface_refits_attempted",
    "shared_surface_refits_completed", "refits_with_warnings",
    "warning_count_across_refits", "refits_with_any_pooled_clipping",
    "total_pooled_clipped_prediction_rows_across_refits",
    "maximum_pooled_clipped_prediction_weighted_share",
    "refits_with_any_own_system_clipping",
    "total_own_system_clipped_prediction_rows_across_refits",
    "maximum_own_system_clipped_prediction_weighted_share",
    "refits_with_any_common_target_clipping",
    "total_common_target_clipped_prediction_rows_across_refits",
    "maximum_common_target_clipped_prediction_weighted_share"
  ),
  validation = c("check", "passed", "value")
)

expected_bundle_names <- c(
  "version", "threshold", "systems", "specifications", "definitions", "combined",
  "combined_all", "diagnostics", "specification_diagnostics", "model_diagnostics",
  "replicate_diagnostics", "validation"
)
expected_combined_names <- c("estimate", "within", "between", "variance", "se", "draw_estimates")
expected_parameter_names <- c(
  paste0("raw_full::", config$g45_systems),
  unlist(lapply(names(standardization_specifications), function(specification) {
    unlist(lapply(standardization_specification_systems(specification), function(system) {
      paste(specification, c("raw_support", "model_own_support", "standardized_common"), system, sep = "::")
    }), use.names = FALSE)
  }), use.names = FALSE)
)

standardization_contract_checks <- c(
  identical(names(standardization_bundle), expected_bundle_names),
  identical(standardization_bundle$version, "2.1-auditability"),
  identical(standardization_bundle$threshold, config$primary_benchmark),
  identical(standardization_bundle$systems, config$g45_systems),
  identical(names(standardization_bundle$specifications), names(standardization_specifications)),
  identical(standardization_bundle$specifications, standardization_specifications),
  identical(names(standardization_bundle$combined_all), expected_combined_names),
  identical(names(standardization_bundle$combined_all$estimate), expected_parameter_names),
  length(standardization_bundle$combined_all$estimate) == 168L,
  identical(dim(standardization_bundle$combined_all$variance), c(168L, 168L)),
  identical(dim(standardization_bundle$combined_all$draw_estimates), c(config$pv_count, 168L)),
  exact_schema(standardization_specification_results, standardization_schemas$specification_results, "Standardization specification results"),
  exact_schema(standardization_alternative_contrasts, standardization_schemas$alternative_contrasts, "Standardization alternative contrasts"),
  exact_schema(standardization_specification_diagnostics, standardization_schemas$specification_diagnostics, "Standardization specification diagnostics"),
  exact_schema(standardization_model_diagnostics, standardization_schemas$model_diagnostics, "Standardization model diagnostics"),
  exact_schema(standardization_replicate_diagnostics, standardization_schemas$replicate_diagnostics, "Standardization replicate diagnostics"),
  exact_schema(standardization_producer_validation, standardization_schemas$validation, "Standardization producer validation"),
  nrow(standardization_specification_results) == 53L,
  nrow(standardization_alternative_contrasts) == 44L,
  nrow(standardization_specification_diagnostics) == 265L,
  nrow(standardization_model_diagnostics) == 135L,
  nrow(standardization_replicate_diagnostics) == 30L,
  nrow(standardization_producer_validation) == 11L,
  exact_ordered_keys(standardization_specification_results, expected_standardization_keys, c("specification", "system"), "Standardization specification results"),
  exact_ordered_keys(standardization_alternative_contrasts, expected_alternative_keys, c("specification", "system"), "Standardization alternative contrasts"),
  exact_ordered_keys(standardization_specification_diagnostics, expected_support_keys, c("draw", "specification", "system"), "Standardization specification diagnostics"),
  exact_ordered_keys(standardization_model_diagnostics, expected_model_keys, c("draw", "surface", "system"), "Standardization model diagnostics"),
  exact_ordered_keys(standardization_replicate_diagnostics, expected_replicate_keys, c("draw", "specification"), "Standardization replicate diagnostics"),
  frames_agree(standardization_bundle$specification_diagnostics, standardization_specification_diagnostics, "Standardization specification diagnostics"),
  frames_agree(standardization_bundle$model_diagnostics, standardization_model_diagnostics, "Standardization model diagnostics"),
  frames_agree(standardization_bundle$replicate_diagnostics, standardization_replicate_diagnostics, "Standardization replicate diagnostics"),
  frames_agree(standardization_bundle$validation, standardization_producer_validation, "Standardization validation")
)

standardization_replicates_complete <- all(
  standardization_replicate_diagnostics$stacked_replicates_attempted ==
    standardization_replicate_diagnostics$stacked_replicates_expected &
    standardization_replicate_diagnostics$stacked_replicates_completed ==
      standardization_replicate_diagnostics$stacked_replicates_expected &
    standardization_replicate_diagnostics$affected_replicates_attempted ==
      standardization_replicate_diagnostics$affected_replicates_expected &
    standardization_replicate_diagnostics$affected_replicates_completed ==
      standardization_replicate_diagnostics$affected_replicates_expected &
    standardization_replicate_diagnostics$support_rebuilds_attempted ==
      standardization_replicate_diagnostics$support_rebuilds_expected &
    standardization_replicate_diagnostics$support_rebuilds_completed ==
      standardization_replicate_diagnostics$support_rebuilds_expected &
    standardization_replicate_diagnostics$shared_surface_refits_attempted ==
      standardization_replicate_diagnostics$shared_surface_refits_expected &
    standardization_replicate_diagnostics$shared_surface_refits_completed ==
      standardization_replicate_diagnostics$shared_surface_refits_expected
)
replicate_clipping_share_columns <- c(
  "maximum_pooled_clipped_prediction_weighted_share",
  "maximum_own_system_clipped_prediction_weighted_share",
  "maximum_common_target_clipped_prediction_weighted_share"
)
replicate_count_columns <- setdiff(
  names(standardization_replicate_diagnostics),
  c("draw", "specification", "surface", replicate_clipping_share_columns)
)
standardization_replicate_diagnostics_valid <-
  all(is.finite(as.matrix(standardization_replicate_diagnostics[replicate_count_columns]))) &&
  all(as.matrix(standardization_replicate_diagnostics[replicate_count_columns]) >= 0) &&
  all(is.finite(as.matrix(standardization_replicate_diagnostics[replicate_clipping_share_columns]))) &&
  all(as.matrix(standardization_replicate_diagnostics[replicate_clipping_share_columns]) >= 0 &
    as.matrix(standardization_replicate_diagnostics[replicate_clipping_share_columns]) <= 1)
model_clipping_share_columns <- c(
  "pooled_clipped_prediction_weighted_share",
  "own_system_clipped_prediction_weighted_share",
  "common_target_clipped_prediction_weighted_share"
)
model_clipping_count_columns <- c(
  "pooled_clipped_prediction_n", "own_system_clipped_prediction_n",
  "common_target_clipped_prediction_n"
)
standardization_model_diagnostics_valid <-
  all(standardization_model_diagnostics$warning_count >= 0) &&
  all(as.matrix(standardization_model_diagnostics[model_clipping_count_columns]) >= 0) &&
  all(is.finite(as.matrix(standardization_model_diagnostics[model_clipping_share_columns]))) &&
  all(as.matrix(standardization_model_diagnostics[model_clipping_share_columns]) >= 0 &
    as.matrix(standardization_model_diagnostics[model_clipping_share_columns]) <= 1)
standardization_producer_checks <- c(
  all(standardization_producer_validation$passed),
  standardization_replicates_complete,
  standardization_replicate_diagnostics_valid,
  all(standardization_model_diagnostics$converged),
  all(standardization_model_diagnostics$model_rank == standardization_model_diagnostics$design_columns),
  all(is.finite(standardization_model_diagnostics$condition_number)),
  all(is.finite(standardization_model_diagnostics$maximum_absolute_coefficient)),
  standardization_model_diagnostics_valid,
  all(standardization_replicate_diagnostics$warning_count_across_refits >= 0)
)
standardization_contract_valid <- all(standardization_contract_checks)
standardization_producer_validations <- all(standardization_producer_checks)
if (!standardization_contract_valid) stop("Standardization bundle contract gate failed")
if (!standardization_producer_validations) stop("Standardization producer validation gate failed")

standardization_definition_map <- data.frame(
  specification = names(standardization_specifications),
  analysis_status = c("Primary", rep("Sensitivity", length(standardization_specifications) - 1L)),
  systems = vapply(standardization_specifications, `[[`, character(1), "systems"),
  surface = vapply(standardization_specifications, `[[`, character(1), "surface"),
  target = vapply(standardization_specifications, `[[`, character(1), "target"),
  support = vapply(standardization_specifications, `[[`, character(1), "support"),
  stringsAsFactors = FALSE
)
standardization_support_summary <- bind_rows(lapply(names(standardization_specifications), function(specification) {
  rows <- standardization_specification_diagnostics[
    standardization_specification_diagnostics$specification == specification, , drop = FALSE
  ]
  data.frame(
    specification = specification,
    draw_system_rows = nrow(rows),
    minimum_retained_weight = min(rows$support_retained_weight),
    minimum_effective_n = min(rows$support_effective_n),
    minimum_target_mass_before_renormalization = min(rows$target_mass_retained_before_renormalization),
    minimum_grid_cells_retained = if (all(is.na(rows$grid_cells_retained))) NA_integer_ else min(rows$grid_cells_retained, na.rm = TRUE),
    maximum_grid_cells_retained = if (all(is.na(rows$grid_cells_retained))) NA_integer_ else max(rows$grid_cells_retained, na.rm = TRUE),
    stringsAsFactors = FALSE
  )
}))
standardization_model_summary <- bind_rows(lapply(expected_surfaces, function(surface) {
  rows <- standardization_model_diagnostics[standardization_model_diagnostics$surface == surface, ]
  data.frame(
    surface = surface,
    fits = nrow(rows),
    converged = all(rows$converged),
    full_rank = all(rows$model_rank == rows$design_columns),
    maximum_condition_number = max(rows$condition_number),
    maximum_pooled_clipped_weighted_share = max(rows$pooled_clipped_prediction_weighted_share),
    maximum_own_system_clipped_weighted_share = max(rows$own_system_clipped_prediction_weighted_share),
    maximum_common_target_clipped_weighted_share = max(rows$common_target_clipped_prediction_weighted_share),
    stringsAsFactors = FALSE
  )
}))
standardization_replicate_summary <- bind_rows(lapply(names(standardization_specifications), function(specification) {
  rows <- standardization_replicate_diagnostics[
    standardization_replicate_diagnostics$specification == specification, , drop = FALSE
  ]
  data.frame(
    specification = specification,
    draws = nrow(rows),
    minimum_affected_replicates_expected = min(rows$affected_replicates_expected),
    maximum_affected_replicates_expected = max(rows$affected_replicates_expected),
    stacked_replicates_expected = sum(rows$stacked_replicates_expected),
    stacked_replicates_attempted = sum(rows$stacked_replicates_attempted),
    stacked_replicates_completed = sum(rows$stacked_replicates_completed),
    support_rebuilds_expected = sum(rows$support_rebuilds_expected),
    support_rebuilds_attempted = sum(rows$support_rebuilds_attempted),
    support_rebuilds_completed = sum(rows$support_rebuilds_completed),
    shared_surface_refits_expected = sum(rows$shared_surface_refits_expected),
    shared_surface_refits_attempted = sum(rows$shared_surface_refits_attempted),
    shared_surface_refits_completed = sum(rows$shared_surface_refits_completed),
    refits_with_warnings = sum(rows$refits_with_warnings),
    warning_count_across_refits = sum(rows$warning_count_across_refits),
    maximum_pooled_clipped_weighted_share = max(rows$maximum_pooled_clipped_prediction_weighted_share),
    maximum_own_system_clipped_weighted_share = max(rows$maximum_own_system_clipped_prediction_weighted_share),
    maximum_common_target_clipped_weighted_share = max(rows$maximum_common_target_clipped_prediction_weighted_share),
    stringsAsFactors = FALSE
  )
}))

message("Standardization bundle gate passed; running threshold and hysteresis sensitivity analyses")
g45 <- readRDS(file.path(paths$derived, "timss_l_g45_analysis.rds"))
g89 <- readRDS(file.path(paths$derived, "timss_l_g89_analysis.rds"))
primary <- readRDS(file.path(paths$derived, "primary_transition_results.rds"))
band_data <- list(g45 = g45, g89 = g89)
band_systems <- list(g45 = config$g45_systems, g89 = config$g89_systems)

threshold_hysteresis_results <- lapply(names(band_data), function(band) {
  run_pv_analysis(
    band_data[[band]], band, band_systems[[band]],
    threshold_hysteresis_builder, "threshold/hysteresis"
  )
})
names(threshold_hysteresis_results) <- names(band_data)
risk_counts <- lapply(names(band_data), function(band) {
  risk_count_summaries(band_data[[band]], band, band_systems[[band]])
})
names(risk_counts) <- names(band_data)
threshold_risk <- bind_rows(lapply(risk_counts, `[[`, "threshold"))
hysteresis_risk <- bind_rows(lapply(risk_counts, `[[`, "hysteresis"))

threshold_rows <- list()
row <- 1L
for (band in names(threshold_hysteresis_results)) {
  for (system in band_systems[[band]]) {
    result <- threshold_hysteresis_results[[band]][[system]]
    for (threshold in config$threshold_grid) {
      risk_row <- threshold_risk[
        threshold_risk$band == band & threshold_risk$system == system & threshold_risk$threshold == threshold,
      ]
      for (metric in threshold_metrics) {
        parameter <- paste0("t", threshold, "__", metric)
        threshold_rows[[row]] <- result_row(
          result, parameter,
          band = band, system = system, threshold = threshold, metric = metric,
          risk_n_mean = risk_row$risk_n_mean,
          risk_n_min = risk_row$risk_n_min,
          risk_n_max = risk_row$risk_n_max
        )
        row <- row + 1L
      }
    }
  }
}
threshold_table <- bind_rows(threshold_rows)
threshold_table$parameter <- NULL

hysteresis_rows <- list()
row <- 1L
hysteresis_metrics <- c(
  "dual_recovery", "dual_persistence", "baseline_risk",
  "weakest_margin_change", "math_change", "science_change"
)
for (band in names(threshold_hysteresis_results)) {
  for (system in band_systems[[band]]) {
    result <- threshold_hysteresis_results[[band]][[system]]
    for (specification in names(hysteresis_specs)) {
      risk_row <- hysteresis_risk[
        hysteresis_risk$band == band & hysteresis_risk$system == system &
          hysteresis_risk$specification == specification,
      ]
      for (metric in hysteresis_metrics) {
        parameter <- paste0("h_", specification, "__", metric)
        hysteresis_rows[[row]] <- result_row(
          result, parameter,
          band = band, system = system, specification = specification,
          baseline_cutoff = risk_row$baseline_cutoff,
          followup_cutoff = risk_row$followup_cutoff,
          metric = metric,
          risk_n_mean = risk_row$risk_n_mean,
          risk_n_min = risk_row$risk_n_min,
          risk_n_max = risk_row$risk_n_max
        )
        row <- row + 1L
      }
    }
  }
}
hysteresis_table <- bind_rows(hysteresis_rows)
hysteresis_table$parameter <- NULL

# Both-wave observed and response-IPW analyses ---------------------------------
core_builder <- function(system_data, map, draw) {
  components <- core_components(
    system_data[[map$math_t0]], system_data[[map$math_t1]],
    system_data[[map$science_t0]], system_data[[map$science_t1]],
    config$primary_benchmark
  )
  function(dataset, weight) evaluate_core(components, weight)
}

response_design <- function(math_t0, science_t0) {
  math_centered <- (math_t0 - config$primary_benchmark) / 100
  science_centered <- (science_t0 - config$primary_benchmark) / 100
  cbind(
    intercept = 1,
    math_centered = math_centered,
    science_centered = science_centered,
    math_centered_sq = math_centered^2,
    science_centered_sq = science_centered^2,
    centered_interaction = math_centered * science_centered
  )
}

fit_response_model <- function(
  design, response, weight,
  propensity_floor = config$response_propensity_floor
) {
  valid <- is.finite(response) & is.finite(weight) & weight > 0 & rowSums(!is.finite(design)) == 0
  if (!any(valid) || !all(response[valid] %in% c(0, 1)) || length(unique(response[valid])) != 2L) {
    stop("Response-IPW propensity model requires finite binary outcomes with both levels")
  }
  model_weight <- weight[valid]
  model_weight <- model_weight / mean(model_weight)
  fit <- suppressWarnings(stats::glm.fit(
    x = design[valid, , drop = FALSE], y = response[valid], weights = model_weight,
    family = stats::quasibinomial(), intercept = FALSE,
    control = stats::glm.control(maxit = 200, epsilon = 1e-9)
  ))
  if (!isTRUE(fit$converged) || fit$rank != ncol(design) || any(!is.finite(fit$coefficients))) {
    stop("Response-IPW propensity model failed convergence, rank, or coefficient gate")
  }
  propensity <- rep(NA_real_, length(response))
  linear_predictor <- as.numeric(design[valid, , drop = FALSE] %*% fit$coefficients)
  clipped_predictor <- pmax(
    -config$response_eta_limit,
    pmin(config$response_eta_limit, linear_predictor)
  )
  propensity[valid] <- stats::plogis(clipped_predictor)
  propensity_floored <- pmax(propensity_floor, propensity)
  clipped <- rep(FALSE, length(response))
  clipped[valid] <- linear_predictor != clipped_predictor
  list(
    propensity = propensity,
    propensity_used = propensity_floored,
    floored = is.finite(propensity) & propensity < propensity_floor,
    clipped = clipped,
    coefficients = fit$coefficients,
    rank = fit$rank,
    converged = fit$converged,
    iterations = fit$iter
  )
}

prefixed_parameters <- function(values, prefix) {
  setNames(as.numeric(values), paste(prefix, names(values), sep = "::"))
}

missingness_joint_builder <- function(system_data, map, draw) {
  baseline_participant <- system_data$participated_2023 == 1
  response <- system_data$participated_2024
  both <- baseline_participant & response == 1
  design <- response_design(system_data[[map$math_t0]], system_data[[map$science_t0]])
  components <- core_components(
    system_data[[map$math_t0]], system_data[[map$math_t1]],
    system_data[[map$science_t0]], system_data[[map$science_t1]],
    config$primary_benchmark
  )
  function(dataset, weight) {
    full_result <- evaluate_core(components, weight)
    both_weight <- ifelse(both, weight, 0)
    both_result <- evaluate_core(components, both_weight)
    fit <- fit_response_model(
      design[baseline_participant, , drop = FALSE],
      response[baseline_participant], weight[baseline_participant]
    )
    propensity_used <- rep(NA_real_, nrow(system_data))
    propensity_used[baseline_participant] <- fit$propensity_used
    stabilization <- weighted_probability(
      response[baseline_participant] == 1,
      weight[baseline_participant]
    )
    if (!is.finite(stabilization) || stabilization <= 0 || stabilization >= 1) {
      stop("Response-IPW stabilization factor is outside (0, 1)")
    }
    ipw_weight <- numeric(length(weight))
    ipw_weight[both] <- weight[both] * stabilization / propensity_used[both]
    ipw_result <- evaluate_core(components, ipw_weight)
    c(
      prefixed_parameters(full_result, "full_pv_completed"),
      prefixed_parameters(both_result, "both_wave_observed"),
      prefixed_parameters(ipw_result, "baseline_participant_response_ipw")
    )
  }
}

extract_joint_method <- function(result, method) {
  prefix <- paste0(method, "::")
  selected <- startsWith(names(result$estimate), prefix)
  if (sum(selected) != length(threshold_metrics)) {
    stop("Joint missingness result lacks the expected method block: ", method)
  }
  parameter_names <- substring(names(result$estimate)[selected], nchar(prefix) + 1L)
  subset_matrix <- function(value) {
    output <- value[selected, selected, drop = FALSE]
    dimnames(output) <- list(parameter_names, parameter_names)
    output
  }
  estimate <- result$estimate[selected]
  names(estimate) <- parameter_names
  standard_error <- result$se[selected]
  names(standard_error) <- parameter_names
  draw_estimates <- result$draw_estimates[, selected, drop = FALSE]
  colnames(draw_estimates) <- parameter_names
  list(
    estimate = estimate,
    within = subset_matrix(result$within),
    between = subset_matrix(result$between),
    variance = subset_matrix(result$variance),
    se = standard_error,
    draw_estimates = draw_estimates
  )
}

missingness_contrast <- function(joint_result, method, parameter) {
  if (method == "full_pv_completed") {
    return(c(estimate = 0, standard_error = 0, ci_lower = 0, ci_upper = 0))
  }
  method_parameter <- paste0(method, "::", parameter)
  full_parameter <- paste0("full_pv_completed::", parameter)
  if (!all(c(method_parameter, full_parameter) %in% names(joint_result$estimate))) {
    stop("Joint missingness result lacks a requested contrast parameter")
  }
  contrast <- setNames(numeric(length(joint_result$estimate)), names(joint_result$estimate))
  contrast[method_parameter] <- 1
  contrast[full_parameter] <- -1
  estimate <- unname(crossprod(contrast, joint_result$estimate))
  variance <- unname(crossprod(contrast, joint_result$variance %*% contrast))
  if (!is.finite(variance) || variance < -1e-10) stop("Invalid missingness contrast variance")
  standard_error <- sqrt(max(0, variance))
  c(
    estimate = estimate,
    standard_error = standard_error,
    ci_lower = estimate - 1.96 * standard_error,
    ci_upper = estimate + 1.96 * standard_error
  )
}

response_diagnostic_rows <- function(data, band, systems) {
  rows <- list()
  row <- 1L
  probabilities <- c(0.01, 0.05, 0.5, 0.95, 0.99)
  for (system in systems) {
    system_data <- data[data$system == system, , drop = FALSE]
    for (draw in seq_len(config$pv_count)) {
      map <- pv_variable_map(band, draw)
      baseline <- system_data$participated_2023 == 1
      response <- system_data$participated_2024
      design <- response_design(system_data[[map$math_t0]], system_data[[map$science_t0]])
      fit <- fit_response_model(
        design[baseline, , drop = FALSE],
        response[baseline], system_data$weight[baseline]
      )
      propensity <- rep(NA_real_, nrow(system_data))
      propensity[baseline] <- fit$propensity
      propensity_used <- rep(NA_real_, nrow(system_data))
      propensity_used[baseline] <- fit$propensity_used
      both <- baseline & response == 1
      response_rate <- weighted_probability(response[baseline] == 1, system_data$weight[baseline])
      multipliers <- response_rate / propensity_used[both]
      ipw_weight <- system_data$weight[both] * multipliers
      propensity_quantiles <- weighted_quantile(
        propensity[baseline], system_data$weight[baseline], probabilities
      )
      multiplier_quantiles <- weighted_quantile(
        multipliers, system_data$weight[both], probabilities
      )

      baseline_design <- design[baseline, -1L, drop = FALSE]
      baseline_weight <- system_data$weight[baseline]
      respondent <- response[baseline] == 1
      respondent_weight <- baseline_weight[respondent]
      respondent_ipw_weight <- respondent_weight * response_rate / fit$propensity_used[respondent]
      balance_before <- balance_after <- setNames(
        numeric(ncol(baseline_design)), colnames(baseline_design)
      )
      for (term in colnames(baseline_design)) {
        value <- baseline_design[, term]
        target_mean <- weighted_mean(value, baseline_weight)
        target_sd <- sqrt(weighted_mean((value - target_mean)^2, baseline_weight))
        if (!is.finite(target_sd) || target_sd <= 0) {
          stop("Response-IPW balance diagnostic has a nonpositive target SD")
        }
        balance_before[term] <- (
          weighted_mean(value[respondent], respondent_weight) - target_mean
        ) / target_sd
        balance_after[term] <- (
          weighted_mean(value[respondent], respondent_ipw_weight) - target_mean
        ) / target_sd
      }

      original_ess <- weighted_effective_n(system_data$weight[both])
      adjusted_ess <- weighted_effective_n(ipw_weight)
      rows[[row]] <- data.frame(
        band = band, system = system, draw = draw,
        baseline_participant_n = sum(baseline), both_wave_n = sum(both),
        weighted_followup_response = response_rate,
        propensity_min = min(propensity[baseline]), propensity_p01 = propensity_quantiles[1],
        propensity_p05 = propensity_quantiles[2], propensity_median = propensity_quantiles[3],
        propensity_p95 = propensity_quantiles[4], propensity_p99 = propensity_quantiles[5],
        propensity_max = max(propensity[baseline]),
        propensity_floor = config$response_propensity_floor,
        weighted_fraction_floored = weighted_probability(fit$floored, system_data$weight[baseline]),
        weighted_fraction_linear_predictor_clipped = weighted_probability(
          fit$clipped, system_data$weight[baseline]
        ),
        ipw_multiplier_min = min(multipliers),
        ipw_multiplier_p01 = multiplier_quantiles[1],
        ipw_multiplier_p05 = multiplier_quantiles[2],
        ipw_multiplier_median = multiplier_quantiles[3],
        ipw_multiplier_p95 = multiplier_quantiles[4],
        ipw_multiplier_p99 = multiplier_quantiles[5],
        ipw_multiplier_max = max(multipliers),
        maximum_ipw_multiplier = max(multipliers),
        original_complete_case_ess = original_ess,
        response_ipw_ess = adjusted_ess,
        response_ipw_to_original_ess_ratio = adjusted_ess / original_ess,
        unadjusted_max_abs_standardized_mean_difference = max(abs(balance_before)),
        response_ipw_max_abs_standardized_mean_difference = max(abs(balance_after)),
        model_rank = fit$rank,
        model_converged = fit$converged, model_iterations = fit$iterations,
        stringsAsFactors = FALSE
      )
      row <- row + 1L
    }
  }
  bind_rows(rows)
}

message("Running joint PV-completed, both-wave observed, and response-IPW analyses")
missingness_joint_results <- list()
full_pv_missingness_results <- list()
both_wave_results <- list()
response_ipw_results <- list()
response_diagnostics <- list()
for (band in names(band_data)) {
  missingness_joint_results[[band]] <- run_pv_analysis(
    band_data[[band]], band, band_systems[[band]],
    missingness_joint_builder, "joint missingness/response"
  )
  full_pv_missingness_results[[band]] <- lapply(
    missingness_joint_results[[band]], extract_joint_method, method = "full_pv_completed"
  )
  both_wave_results[[band]] <- lapply(
    missingness_joint_results[[band]], extract_joint_method, method = "both_wave_observed"
  )
  response_ipw_results[[band]] <- lapply(
    missingness_joint_results[[band]], extract_joint_method,
    method = "baseline_participant_response_ipw"
  )
  response_diagnostics[[band]] <- response_diagnostic_rows(
    band_data[[band]], band, band_systems[[band]]
  )
}
response_diagnostics <- bind_rows(response_diagnostics)

observed_response_rows <- list()
row <- 1L
for (band in names(band_data)) {
  for (system in band_systems[[band]]) {
    method_results <- list(
      full_pv_completed = full_pv_missingness_results[[band]][[system]],
      both_wave_observed = both_wave_results[[band]][[system]],
      baseline_participant_response_ipw = response_ipw_results[[band]][[system]]
    )
    joint_result <- missingness_joint_results[[band]][[system]]
    for (method in names(method_results)) {
      result <- method_results[[method]]
      for (metric in threshold_metrics) {
        target <- switch(
          method,
          full_pv_completed = "all sampled longitudinal students represented by PV completion",
          both_wave_observed = "students participating in both waves",
          baseline_participant_response_ipw = "baseline participants, adjusted for follow-up nonresponse"
        )
        contrast <- missingness_contrast(joint_result, method, metric)
        output_row <- result_row(
          result, metric,
          band = band, system = system, method = method, target_population = target,
          metric = metric
        )
        output_row$full_pv_estimate <- unname(
          full_pv_missingness_results[[band]][[system]]$estimate[metric]
        )
        output_row$descriptive_difference_from_full <- unname(contrast["estimate"])
        output_row$difference_standard_error <- unname(contrast["standard_error"])
        output_row$difference_ci_lower <- unname(contrast["ci_lower"])
        output_row$difference_ci_upper <- unname(contrast["ci_upper"])
        observed_response_rows[[row]] <- output_row
        row <- row + 1L
      }
    }
  }
}
observed_response_table <- bind_rows(observed_response_rows)
observed_response_table$parameter <- NULL
if (max(abs(
  observed_response_table$descriptive_difference_from_full -
    (observed_response_table$estimate - observed_response_table$full_pv_estimate)
)) > 1e-12) {
  stop("Joint missingness contrast point estimates do not equal method-minus-primary differences")
}

# Pattern-mixture MNAR score shifts --------------------------------------------
mnar_scenarios <- c("attriter_followup_shift", "late_entrant_baseline_shift", "combined_pattern_shift")
delta_tag <- function(delta) if (delta == 0) "0" else paste0("m", abs(delta))

mnar_builder <- function(system_data, map, draw) {
  original <- list(
    math_t0 = system_data[[map$math_t0]], math_t1 = system_data[[map$math_t1]],
    science_t0 = system_data[[map$science_t0]], science_t1 = system_data[[map$science_t1]]
  )
  attriter <- system_data$participation_pattern == "2023_only"
  late_entrant <- system_data$participation_pattern == "2024_only"
  components <- list()
  for (scenario in mnar_scenarios) {
    for (delta in config$mnar_deltas) {
      shifted <- original
      if (scenario %in% c("attriter_followup_shift", "combined_pattern_shift")) {
        shifted$math_t1[attriter] <- shifted$math_t1[attriter] + delta
        shifted$science_t1[attriter] <- shifted$science_t1[attriter] + delta
      }
      if (scenario %in% c("late_entrant_baseline_shift", "combined_pattern_shift")) {
        shifted$math_t0[late_entrant] <- shifted$math_t0[late_entrant] + delta
        shifted$science_t0[late_entrant] <- shifted$science_t0[late_entrant] + delta
      }
      name <- paste0(scenario, "__d", delta_tag(delta))
      components[[name]] <- core_components(
        shifted$math_t0, shifted$math_t1, shifted$science_t0, shifted$science_t1,
        config$primary_benchmark
      )
    }
  }
  function(dataset, weight) {
    values <- numeric(0)
    for (name in names(components)) {
      metrics <- evaluate_core(components[[name]], weight)
      names(metrics) <- paste0(name, "__", names(metrics))
      values <- c(values, metrics)
    }
    values
  }
}

message("Running MNAR pattern-mixture sensitivity analyses")
mnar_results <- lapply(names(band_data), function(band) {
  run_pv_analysis(band_data[[band]], band, band_systems[[band]], mnar_builder, "MNAR shifts")
})
names(mnar_results) <- names(band_data)

mnar_rows <- list()
row <- 1L
for (band in names(mnar_results)) {
  for (system in band_systems[[band]]) {
    result <- mnar_results[[band]][[system]]
    for (scenario in mnar_scenarios) {
      for (delta in config$mnar_deltas) {
        for (metric in threshold_metrics) {
          parameter <- paste0(scenario, "__d", delta_tag(delta), "__", metric)
          mnar_rows[[row]] <- result_row(
            result, parameter,
            band = band, system = system, scenario = scenario,
            delta_points_per_affected_domain = delta, metric = metric,
            affected_pattern = switch(
              scenario,
              attriter_followup_shift = "2023_only: both follow-up domain PVs",
              late_entrant_baseline_shift = "2024_only: both baseline domain PVs",
              combined_pattern_shift = "2023_only follow-up and 2024_only baseline PVs"
            )
          )
          row <- row + 1L
        }
      }
    }
  }
}
mnar_table <- bind_rows(mnar_rows)
mnar_table$parameter <- NULL

first_crossing <- function(delta, estimate, target) {
  order_index <- order(delta, decreasing = TRUE)
  delta <- delta[order_index]
  estimate <- estimate[order_index]
  if (!all(is.finite(c(delta, estimate, target)))) return(NA_real_)
  reached <- which(estimate <= target)
  if (!length(reached)) return(NA_real_)
  index <- reached[1]
  if (index == 1L) return(delta[1])
  x1 <- delta[index - 1L]
  x2 <- delta[index]
  y1 <- estimate[index - 1L]
  y2 <- estimate[index]
  if (y1 == y2) return(x2)
  x1 + (target - y1) * (x2 - x1) / (y2 - y1)
}

mnar_tipping_rows <- list()
row <- 1L
mnar_recovery <- mnar_table[mnar_table$metric == "BB_to_AA", ]
for (band in names(band_data)) {
  for (system in band_systems[[band]]) {
    for (scenario in mnar_scenarios) {
      subset <- mnar_recovery[
        mnar_recovery$band == band & mnar_recovery$system == system &
          mnar_recovery$scenario == scenario,
      ]
      mar <- subset[subset$delta_points_per_affected_domain == 0, ]
      half_target <- mar$estimate / 2
      mar_lower_target <- mar$estimate - 1.96 * mar$standard_error
      minus_40 <- subset[subset$delta_points_per_affected_domain == -40, ]
      half_delta <- first_crossing(subset$delta_points_per_affected_domain, subset$estimate, half_target)
      lower_delta <- first_crossing(subset$delta_points_per_affected_domain, subset$estimate, mar_lower_target)
      mnar_tipping_rows[[row]] <- data.frame(
        band = band, system = system, scenario = scenario,
        mar_estimate = mar$estimate, mar_standard_error = mar$standard_error,
        estimate_at_minus_40 = minus_40$estimate,
        ratio_at_minus_40_to_mar = minus_40$estimate / mar$estimate,
        delta_crossing_half_mar = half_delta,
        half_mar_status = ifelse(is.na(half_delta), "not_reached_within_0_to_minus40", "crossed_within_grid"),
        delta_crossing_mar_lower_95 = lower_delta,
        mar_lower_95_status = ifelse(is.na(lower_delta), "not_reached_within_0_to_minus40", "crossed_within_grid"),
        interpretation = "Diagnostic interpolation across planned shifts; not a causal missing-data correction",
        stringsAsFactors = FALSE
      )
      row <- row + 1L
    }
  }
}
mnar_tipping_table <- bind_rows(mnar_tipping_rows)

# Exploratory same-wave unequal-index classification-discordance diagnostic ----
precompute_cross_pv <- function(system_data, band, threshold) {
  maps <- lapply(seq_len(config$pv_count), function(draw) pv_variable_map(band, draw))
  waves <- list(
    baseline = list(math = "math_t0", science = "science_t0"),
    followup = list(math = "math_t1", science = "science_t1")
  )
  pairs <- list()
  for (wave in names(waves)) {
    wave_pairs <- list()
    pair <- 1L
    for (origin_draw in seq_len(config$pv_count)) {
      for (destination_draw in seq_len(config$pv_count)) {
        if (origin_draw == destination_draw) next
        origin_map <- maps[[origin_draw]]
        destination_map <- maps[[destination_draw]]
        origin_state <- achievement_state(
          system_data[[origin_map[[waves[[wave]]$math]]]],
          system_data[[origin_map[[waves[[wave]]$science]]]], threshold
        )
        destination_state <- achievement_state(
          system_data[[destination_map[[waves[[wave]]$math]]]],
          system_data[[destination_map[[waves[[wave]]$science]]]], threshold
        )
        wave_pairs[[pair]] <- list(
          risk = origin_state == "BB", destination = destination_state == "AA",
          origin_draw = origin_draw, destination_draw = destination_draw
        )
        pair <- pair + 1L
      }
    }
    pairs[[wave]] <- wave_pairs
  }
  pairs
}

evaluate_cross_pairs <- function(pairs, weight) {
  vapply(pairs, function(pair) {
    weighted_probability(pair$destination[pair$risk], weight[pair$risk])
  }, numeric(1))
}

run_cross_pv_diagnostic <- function(data, band, systems) {
  results <- setNames(vector("list", length(systems)), systems)
  pair_summaries <- list()
  row <- 1L
  for (system in systems) {
    message("  same-wave unequal-index classification-discordance diagnostic: ", band, " / ", system)
    system_data <- data[data$system == system, , drop = FALSE]
    pairs <- precompute_cross_pv(system_data, band, config$primary_benchmark)
    estimator <- function(dataset, weight) {
      c(
        baseline_cross_pv_BB_to_AA = mean(evaluate_cross_pairs(pairs$baseline, weight)),
        followup_cross_pv_BB_to_AA = mean(evaluate_cross_pairs(pairs$followup, weight))
      )
    }
    result <- jrr_estimate(system_data, estimator)
    result$replicates <- NULL
    result$se <- sqrt(diag(result$variance))
    results[[system]] <- result
    for (wave in names(pairs)) {
      pair_values <- evaluate_cross_pairs(pairs[[wave]], system_data$weight)
      pair_summaries[[row]] <- data.frame(
        band = band, system = system, wave = wave,
        ordered_pv_pairs = length(pair_values), pair_mean = mean(pair_values),
        pair_standard_deviation = stats::sd(pair_values), pair_minimum = min(pair_values),
        pair_maximum = max(pair_values),
        stringsAsFactors = FALSE
      )
      row <- row + 1L
    }
  }
  list(results = results, pair_summaries = bind_rows(pair_summaries))
}

message("Running exploratory same-wave unequal-index classification-discordance diagnostics")
cross_pv <- lapply(names(band_data), function(band) {
  run_cross_pv_diagnostic(band_data[[band]], band, band_systems[[band]])
})
names(cross_pv) <- names(band_data)

cross_pv_rows <- list()
row <- 1L
for (band in names(cross_pv)) {
  for (system in band_systems[[band]]) {
    actual <- primary[[band]][[system]]
    cross_result <- cross_pv[[band]]$results[[system]]
    pair_summary <- cross_pv[[band]]$pair_summaries
    cross_pv_rows[[row]] <- data.frame(
      band = band, system = system, comparison = "actual_longitudinal_matched_pv",
      estimate = unname(actual$estimate["BB_to_AA"]),
      standard_error = unname(actual$se["BB_to_AA"]),
      pair_standard_deviation = NA_real_, pair_minimum = NA_real_, pair_maximum = NA_real_,
      variance_scope = "PV plus JRR total variance",
      interpretation = "Observed matched-PV longitudinal BB-to-AA transition to joint attainment",
      stringsAsFactors = FALSE
    )
    row <- row + 1L
    for (wave in c("baseline", "followup")) {
      parameter <- paste0(wave, "_cross_pv_BB_to_AA")
      summary_row <- pair_summary[pair_summary$system == system & pair_summary$wave == wave, ]
      cross_pv_rows[[row]] <- data.frame(
        band = band, system = system, comparison = paste0(wave, "_same_wave_cross_pv"),
        estimate = unname(cross_result$estimate[parameter]),
        standard_error = unname(cross_result$se[parameter]),
        pair_standard_deviation = summary_row$pair_standard_deviation,
        pair_minimum = summary_row$pair_minimum,
        pair_maximum = summary_row$pair_maximum,
        variance_scope = "JRR sampling variance of the 20-pair mean conditional on five released PVs; pair dispersion separate",
        interpretation = "Exploratory same-wave unequal-index classification discordance; not a longitudinal or measurement-error correction",
        stringsAsFactors = FALSE
      )
      row <- row + 1L
    }
  }
}
cross_pv_table <- bind_rows(cross_pv_rows)
actual_lookup <- cross_pv_table[cross_pv_table$comparison == "actual_longitudinal_matched_pv", c("band", "system", "estimate")]
names(actual_lookup)[3] <- "actual_longitudinal_estimate"
cross_pv_table <- merge(cross_pv_table, actual_lookup, by = c("band", "system"), all.x = TRUE, sort = FALSE)
cross_pv_table$ratio_to_actual_longitudinal <- cross_pv_table$estimate / cross_pv_table$actual_longitudinal_estimate

# Separate Grade 8-to-9 domain extension ------------------------------------
domain_builder <- function(system_data, map, draw) {
  domain_components <- list(
    mathematics = list(
      t0 = system_data[[map$math_t0]], t1 = system_data[[map$math_t1]]
    ),
    science = list(
      t0 = system_data[[map$science_t0]], t1 = system_data[[map$science_t1]]
    )
  )
  function(dataset, weight) {
    values <- numeric(0)
    for (domain in names(domain_components)) {
      component <- domain_components[[domain]]
      below <- component$t0 < config$primary_benchmark
      metrics <- c(
        below_to_at_or_above = weighted_probability(component$t1[below] >= config$primary_benchmark, weight[below]),
        below_persistence = weighted_probability(component$t1[below] < config$primary_benchmark, weight[below]),
        baseline_below = weighted_probability(below, weight),
        score_change_if_below = weighted_mean((component$t1 - component$t0)[below], weight[below])
      )
      names(metrics) <- paste0(domain, "__", names(metrics))
      values <- c(values, metrics)
    }
    values
  }
}

message("Running separate Grade 8-to-9 domain extension")
g89_domain_results <- run_pv_analysis(
  g89, "g89", config$g89_systems, domain_builder, "Grade 8-to-9 domains"
)
g89_domain_rows <- list()
row <- 1L
for (system in config$g89_systems) {
  result <- g89_domain_results[[system]]
  for (domain in c("mathematics", "science")) {
    for (metric in c("below_to_at_or_above", "below_persistence", "baseline_below", "score_change_if_below")) {
      parameter <- paste0(domain, "__", metric)
      caution <- if (system == "Jordan" && domain == "science") {
        "Interpret separately: TIMSS documents a change in Jordan's science instructional/assessment organization across waves."
      } else {
        "None specific to this system-domain beyond the general noncausal design."
      }
      g89_domain_rows[[row]] <- result_row(
        result, parameter,
        band = "g89", system = system, domain = domain,
        threshold = config$primary_benchmark, metric = metric,
        interpretation_caution = caution
      )
      row <- row + 1L
    }
  }
}
g89_domain_table <- bind_rows(g89_domain_rows)
g89_domain_table$parameter <- NULL

# Validation and output gate ---------------------------------------------------
validate_combined_collection <- function(collection, label) {
  rows <- list()
  row <- 1L
  for (band in names(collection)) {
    for (system in names(collection[[band]])) {
      result <- collection[[band]][[system]]
      variance <- result$variance
      eigenvalues <- eigen((variance + t(variance)) / 2, symmetric = TRUE, only.values = TRUE)$values
      rows[[row]] <- data.frame(
        analysis = label, band = band, system = system,
        estimates_finite = all(is.finite(result$estimate)),
        covariance_finite = all(is.finite(variance)),
        covariance_symmetric = max(abs(variance - t(variance))) < 1e-10,
        covariance_psd = min(eigenvalues) > -1e-8,
        se_matches_diagonal = max(abs(result$se - sqrt(pmax(0, diag(variance))))) < 1e-10,
        minimum_eigenvalue = min(eigenvalues),
        stringsAsFactors = FALSE
      )
      row <- row + 1L
    }
  }
  bind_rows(rows)
}

validation_rows <- bind_rows(list(
  validate_combined_collection(threshold_hysteresis_results, "threshold_hysteresis"),
  validate_combined_collection(missingness_joint_results, "missingness_joint"),
  validate_combined_collection(both_wave_results, "both_wave_observed"),
  validate_combined_collection(response_ipw_results, "response_ipw"),
  validate_combined_collection(mnar_results, "mnar_pattern_mixture"),
  validate_combined_collection(list(g89 = g89_domain_results), "g89_domain_extension")
))

primary_differences <- list()
row <- 1L
for (band in names(band_data)) {
  for (system in band_systems[[band]]) {
    sensitivity <- threshold_hysteresis_results[[band]][[system]]
    for (metric in threshold_metrics) {
      parameter <- paste0("t", config$primary_benchmark, "__", metric)
      primary_differences[[row]] <- data.frame(
        band = band, system = system, metric = metric,
        estimate_difference = unname(sensitivity$estimate[parameter] - primary[[band]][[system]]$estimate[metric]),
        se_difference = unname(sensitivity$se[parameter] - primary[[band]][[system]]$se[metric]),
        stringsAsFactors = FALSE
      )
      row <- row + 1L
    }
  }
}
primary_differences <- bind_rows(primary_differences)

missingness_full_differences <- list()
row <- 1L
for (band in names(band_data)) {
  for (system in band_systems[[band]]) {
    result <- full_pv_missingness_results[[band]][[system]]
    for (metric in threshold_metrics) {
      missingness_full_differences[[row]] <- data.frame(
        band = band, system = system, metric = metric,
        estimate_difference = unname(
          result$estimate[metric] - primary[[band]][[system]]$estimate[metric]
        ),
        se_difference = unname(
          result$se[metric] - primary[[band]][[system]]$se[metric]
        ),
        stringsAsFactors = FALSE
      )
      row <- row + 1L
    }
  }
}
missingness_full_differences <- bind_rows(missingness_full_differences)

mnar_zero_differences <- list()
row <- 1L
for (band in names(band_data)) {
  for (system in band_systems[[band]]) {
    result <- mnar_results[[band]][[system]]
    for (scenario in mnar_scenarios) {
      for (metric in threshold_metrics) {
        parameter <- paste0(scenario, "__d0__", metric)
        mnar_zero_differences[[row]] <- data.frame(
          band = band, system = system, scenario = scenario, metric = metric,
          estimate_difference = unname(result$estimate[parameter] - primary[[band]][[system]]$estimate[metric]),
          se_difference = unname(result$se[parameter] - primary[[band]][[system]]$se[metric]),
          stringsAsFactors = FALSE
        )
        row <- row + 1L
      }
    }
  }
}
mnar_zero_differences <- bind_rows(mnar_zero_differences)

attriter_recovery <- mnar_table[
  mnar_table$scenario == "attriter_followup_shift" & mnar_table$metric == "BB_to_AA",
]
monotonicity_groups <- split(
  attriter_recovery,
  interaction(attriter_recovery$band, attriter_recovery$system, drop = TRUE)
)
monotonicity <- vapply(monotonicity_groups, function(group) {
  ordered <- group[order(group$delta_points_per_affected_domain, decreasing = TRUE), ]
  all(diff(ordered$estimate) <= 1e-12)
}, logical(1))

validation_checks <- data.frame(
  check = c(
    "all covariance collections valid",
    "threshold 475 and joint missingness full method exactly reproduce primary estimates and SEs",
    "MNAR delta zero exactly reproduces primary estimates and SEs",
    "response models converged at full rank",
    "response propensities satisfy finite positivity and floor gates",
    "response-IPW weights, balance diagnostics, contrasts, and effective sample sizes are valid",
    "all probability estimates within bounds",
    "attriter follow-up shifts are monotone nonincreasing",
    "same-wave unequal-index diagnostics use all 20 ordered PV pairs",
    "configured threshold grid executed",
    "Grade 8-to-9 system set is exact"
  ),
  passed = c(
    all(validation_rows$estimates_finite & validation_rows$covariance_finite &
      validation_rows$covariance_symmetric & validation_rows$covariance_psd &
      validation_rows$se_matches_diagonal),
    max(abs(c(
      primary_differences$estimate_difference, primary_differences$se_difference,
      missingness_full_differences$estimate_difference, missingness_full_differences$se_difference
    ))) < 1e-10,
    max(abs(c(mnar_zero_differences$estimate_difference, mnar_zero_differences$se_difference))) < 1e-10,
    all(response_diagnostics$model_converged & response_diagnostics$model_rank == 6L),
    all(is.finite(response_diagnostics$propensity_min) & response_diagnostics$propensity_min > 0 &
      is.finite(response_diagnostics$propensity_max) & response_diagnostics$propensity_max < 1 &
      response_diagnostics$weighted_fraction_floored == 0),
    all(is.finite(as.matrix(response_diagnostics[c(
      "ipw_multiplier_min", "ipw_multiplier_p01", "ipw_multiplier_p05",
      "ipw_multiplier_median", "ipw_multiplier_p95", "ipw_multiplier_p99",
      "ipw_multiplier_max", "original_complete_case_ess", "response_ipw_ess",
      "response_ipw_to_original_ess_ratio",
      "unadjusted_max_abs_standardized_mean_difference",
      "response_ipw_max_abs_standardized_mean_difference"
    )]))) &&
      all(response_diagnostics$ipw_multiplier_min > 0 &
        response_diagnostics$response_ipw_ess > 0 &
        response_diagnostics$original_complete_case_ess > 0 &
        response_diagnostics$response_ipw_to_original_ess_ratio > 0) &&
      all(is.finite(observed_response_table$difference_standard_error) &
        observed_response_table$difference_standard_error >= 0 &
        is.finite(observed_response_table$difference_ci_lower) &
        is.finite(observed_response_table$difference_ci_upper)),
    all(threshold_table$estimate[threshold_table$metric %in% c("BB_to_AA", "BB_to_BB", "bb_exit", "baseline_BB")] >= 0 &
      threshold_table$estimate[threshold_table$metric %in% c("BB_to_AA", "BB_to_BB", "bb_exit", "baseline_BB")] <= 1) &&
      all(hysteresis_table$estimate[hysteresis_table$metric %in% c("dual_recovery", "dual_persistence", "baseline_risk")] >= 0 &
        hysteresis_table$estimate[hysteresis_table$metric %in% c("dual_recovery", "dual_persistence", "baseline_risk")] <= 1) &&
      all(observed_response_table$estimate[observed_response_table$metric %in% c("BB_to_AA", "BB_to_BB", "bb_exit", "baseline_BB")] >= 0 &
        observed_response_table$estimate[observed_response_table$metric %in% c("BB_to_AA", "BB_to_BB", "bb_exit", "baseline_BB")] <= 1) &&
      all(mnar_table$estimate[mnar_table$metric %in% c("BB_to_AA", "BB_to_BB", "bb_exit", "baseline_BB")] >= 0 &
        mnar_table$estimate[mnar_table$metric %in% c("BB_to_AA", "BB_to_BB", "bb_exit", "baseline_BB")] <= 1),
    all(as.logical(monotonicity)),
    all(bind_rows(lapply(cross_pv, `[[`, "pair_summaries"))$ordered_pv_pairs == 20),
    setequal(unique(threshold_table$threshold), config$threshold_grid),
    setequal(names(g89_domain_results), config$g89_systems)
  ),
  stringsAsFactors = FALSE
)

validation_checks <- rbind(
  validation_checks,
  data.frame(
    check = c(
      "standardization bundle contract is exact",
      "standardization producer validations passed"
    ),
    passed = c(standardization_contract_valid, standardization_producer_validations),
    stringsAsFactors = FALSE
  )
)
validation_checks$passed[is.na(validation_checks$passed)] <- FALSE
if (!all(validation_checks$passed)) {
  failed <- validation_checks$check[!validation_checks$passed]
  stop("Sensitivity validation failed: ", paste(failed, collapse = "; "))
}

sensitivity_results <- list(
  definitions = list(
    threshold_grid = config$threshold_grid,
    hysteresis = hysteresis_specs,
    primary_threshold = config$primary_benchmark,
    response_ipw = paste(
      "System- and PV-specific weighted quadratic-logit follow-up response model among baseline participants;",
      "stabilized inverse-probability weights with propensity floor",
      config$response_propensity_floor,
      "and covariance-aware method-minus-primary contrasts from a joint JRR/PV stack."
    ),
    mnar = paste(
      "Pattern-mixture offsets add the listed delta to both affected domain PVs for 2023-only follow-up",
      "scores, 2024-only baseline scores, or both; delta zero is the PV-completed primary analysis."
    ),
    cross_pv = paste(
      "Exploratory mean of all 20 ordered unequal-index BB-to-AA classifications within the same wave;",
      "JRR sampling-only conditional on five released PVs, not a longitudinal or measurement-error correction."
    )
  ),
  threshold_hysteresis = threshold_hysteresis_results,
  missingness_joint = missingness_joint_results,
  both_wave_observed = both_wave_results,
  response_ipw = response_ipw_results,
  mnar = mnar_results,
  cross_pv = cross_pv,
  g89_domains = g89_domain_results,
  standardization = list(
    version = standardization_bundle$version,
    threshold = standardization_bundle$threshold,
    definitions = list(
      producer = standardization_bundle$definitions,
      specification_status = standardization_definition_map
    ),
    specification_results = standardization_specification_results,
    alternative_contrasts = standardization_alternative_contrasts,
    validation = standardization_producer_validation,
    diagnostic_summaries = list(
      support = standardization_support_summary,
      models = standardization_model_summary,
      replicates = standardization_replicate_summary
    )
  ),
  validation = validation_checks
)

saveRDS(sensitivity_results, file.path(paths$derived, "sensitivity_results.rds"), compress = "xz")
write.csv(threshold_table, file.path(paths$tables, "threshold_sensitivity.csv"), row.names = FALSE)
write.csv(hysteresis_table, file.path(paths$tables, "hysteresis_sensitivity.csv"), row.names = FALSE)
write.csv(observed_response_table, file.path(paths$tables, "observed_response_ipw_sensitivity.csv"), row.names = FALSE)
write.csv(mnar_table, file.path(paths$tables, "mnar_pattern_mixture_sensitivity.csv"), row.names = FALSE)
write.csv(mnar_tipping_table, file.path(paths$tables, "mnar_tipping_points.csv"), row.names = FALSE)
write.csv(cross_pv_table, file.path(paths$tables, "measurement_crossing_diagnostic.csv"), row.names = FALSE)
write.csv(g89_domain_table, file.path(paths$tables, "g89_domain_extension.csv"), row.names = FALSE)
stale_g89_domain_path <- file.path(paths$tables, "g89_domain_replication.csv")
if (file.exists(stale_g89_domain_path) && unlink(stale_g89_domain_path) != 0L) {
  stop("Could not remove superseded legacy Grade 8-to-9 domain artifact")
}
write.csv(response_diagnostics, file.path(paths$diagnostics, "response_ipw_diagnostics.csv"), row.names = FALSE)
write.csv(validation_rows, file.path(paths$diagnostics, "sensitivity_covariance_validation.csv"), row.names = FALSE)
write.csv(validation_checks, file.path(paths$diagnostics, "sensitivity_validation_checks.csv"), row.names = FALSE)

message(
  "Sensitivity analyses complete: ", nrow(threshold_table), " threshold rows, ",
  nrow(mnar_table), " MNAR rows, and all ", nrow(validation_checks), " validation checks passed"
)
