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
}

publication_system <- function(system) {
  ifelse(system == "Korea, Rep. of", "Republic of Korea", system)
}

publication_band <- function(band) {
  ifelse(band == "g45", "Grades 4 to 5", "Grades 8 to 9")
}

publication_scope <- function(band) {
  ifelse(
    band == "g45",
    "Primary Grades 4-to-5 analysis",
    "Descriptive Grade 8-to-9 extension"
  )
}

format_percent_ci <- function(estimate, lower, upper, digits = 1) {
  sprintf(
    paste0("%.", digits, "f (95%% CI %.", digits, "f to %.", digits, "f)"),
    100 * estimate, 100 * lower, 100 * upper
  )
}

read_required_csv <- function(directory, filename) {
  path <- file.path(directory, filename)
  if (!file.exists(path)) stop("Required artifact is missing: ", path)
  utils::read.csv(path, check.names = FALSE, stringsAsFactors = FALSE)
}

require_gate <- function(condition, message) {
  if (length(condition) != 1L || is.na(condition) || !condition) stop(message)
  invisible(TRUE)
}

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

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)
  require_gate(!anyDuplicated(observed_keys), paste0(label, " has duplicate keys"))
  require_gate(identical(observed_keys, expected_keys), paste0(label, " key set or order is invalid"))
}

frames_agree <- function(left, right, label, tolerance = 1e-10) {
  require_gate(
    identical(names(left), names(right)) && nrow(left) == nrow(right),
    paste0(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_match <- is.finite(lhs) & is.finite(rhs) & abs(lhs - rhs) <= tolerance
      require_gate(all(matching_na | finite_match), paste0(label, " differs in ", column))
    } else {
      require_gate(identical(lhs, rhs), paste0(label, " differs in ", column))
    }
  }
  invisible(TRUE)
}

format_number <- function(value, digits = 1L, signed = FALSE, big_mark = TRUE) {
  require_gate(length(value) == 1L && is.finite(value), "Cannot format a nonfinite publication value")
  magnitude <- formatC(
    abs(value), format = "f", digits = digits,
    big.mark = if (big_mark) "," else ""
  )
  prefix <- if (value < 0) "−" else if (signed && value > 0) "+" else ""
  paste0(prefix, magnitude)
}

markdown_table <- function(headers, alignments, rows) {
  require_gate(length(headers) == length(alignments), "Markdown table header/alignment mismatch")
  rows <- as.matrix(rows)
  require_gate(ncol(rows) == length(headers), "Markdown table row width mismatch")
  escape_cell <- function(value) {
    value <- gsub("|", "&#124;", as.character(value), fixed = TRUE)
    gsub("[\r\n]+", " ", value)
  }
  row_line <- function(values) paste0("| ", paste(vapply(values, escape_cell, character(1)), collapse = " | "), " |")
  c(
    row_line(headers),
    row_line(alignments),
    apply(rows, 1L, row_line)
  )
}

write_utf8_lines <- function(lines, path) {
  dir.create(dirname(path), recursive = TRUE, showWarnings = FALSE)
  connection <- file(path, open = "wb")
  on.exit(close(connection))
  writeChar(paste0(paste(lines, collapse = "\n"), "\n"), connection, eos = NULL, useBytes = TRUE)
  invisible(path)
}

publish_staged_set <- function(staged, destinations, stale = character(), validate_installed) {
  require_gate(length(staged) == length(destinations) && length(staged) > 0L, "Publication install plan is empty or unbalanced")
  require_gate(all(file.exists(staged)) && !any(dir.exists(staged)), "A staged publication file is missing")
  require_gate(!anyDuplicated(destinations), "Publication install plan contains duplicate destinations")
  invisible(lapply(unique(dirname(destinations)), dir.create, recursive = TRUE, showWarnings = FALSE))

  stale <- unique(stale[file.exists(stale) & !stale %in% destinations])
  existing <- unique(c(destinations[file.exists(destinations)], stale))
  backup_root <- tempfile(pattern = "backups-", tmpdir = dirname(staged[[1L]]))
  if (!dir.create(backup_root, recursive = TRUE, showWarnings = FALSE) && !dir.exists(backup_root)) {
    stop("Could not create publication backup directory: ", backup_root)
  }
  backup_paths <- setNames(
    file.path(backup_root, sprintf("%04d-%s", seq_along(existing), basename(existing))),
    existing
  )
  moved <- character()
  installed <- character()
  failure <- NULL

  tryCatch({
    for (source in existing) {
      backup <- unname(backup_paths[[source]])
      if (!file.rename(source, backup)) stop("Could not back up publication file: ", source)
      moved <- c(moved, source)
    }
    for (index in seq_along(staged)) {
      destination <- destinations[[index]]
      if (!file.rename(staged[[index]], destination)) {
        stop("Could not install staged publication file: ", destination)
      }
      installed <- c(installed, destination)
    }
    validate_installed()
  }, error = function(error) {
    failure <<- error
  })

  if (!is.null(failure)) {
    rollback_errors <- character()
    for (destination in rev(installed)) {
      if (file.exists(destination) && unlink(destination) != 0L) {
        rollback_errors <- c(rollback_errors, paste0("remove ", destination))
      }
    }
    for (source in rev(moved)) {
      backup <- unname(backup_paths[[source]])
      if (file.exists(backup) && !file.rename(backup, source)) {
        rollback_errors <- c(rollback_errors, paste0("restore ", source))
      }
    }
    unlink(backup_root, recursive = TRUE, force = TRUE)
    if (length(rollback_errors)) {
      stop(
        "Publication install failed (", conditionMessage(failure), "); rollback was incomplete: ",
        paste(rollback_errors, collapse = "; ")
      )
    }
    stop(failure)
  }

  unlink(backup_root, recursive = TRUE, force = TRUE)
  invisible(destinations)
}

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

specification_labels <- c(
  primary = "Primary",
  population_target = "Pooled represented-BB target",
  leave_korea_out = "Leave Republic of Korea out",
  strict_common_cells = "Strict common cells",
  additive_linear = "Additive linear surface",
  linear_interaction = "Linear interaction surface"
)
standardization_status <- c(
  primary = "Primary",
  population_target = "Sensitivity",
  leave_korea_out = "Sensitivity",
  strict_common_cells = "Sensitivity",
  additive_linear = "Sensitivity",
  linear_interaction = "Sensitivity"
)
surface_labels <- c(
  quadratic_interaction = "Quadratic with interaction",
  additive_linear = "Additive linear",
  linear_interaction = "Linear with interaction"
)
target_labels <- c(
  equal_system = "Equal-system target",
  pooled_population = "Pooled represented-BB target"
)
support_labels <- c(
  rectangle_01_99 = "Joint marginal 1st-99th percentile rectangle",
  strict_common_cells = "Strict common-cell support"
)
allowed_analysis_status <- c("Primary", "Sensitivity", "Exploratory")
allowed_analysis_scope <- c(
  "Primary Grades 4-to-5 analysis",
  "Descriptive Grade 8-to-9 extension"
)

message("Loading and gating all validated analysis artifacts")
primary <- readRDS(file.path(paths$derived, "primary_transition_results.rds"))
standardization <- readRDS(file.path(paths$derived, "standardization_results.rds"))
sensitivity <- readRDS(file.path(paths$derived, "sensitivity_results.rds"))
g45 <- readRDS(file.path(paths$derived, "timss_l_g45_analysis.rds"))
g89 <- readRDS(file.path(paths$derived, "timss_l_g89_analysis.rds"))

inventory <- read_required_csv(paths$diagnostics, "data_inventory.csv")
official_validation <- read_required_csv(paths$diagnostics, "official_result_validation.csv")
sensitivity_validation <- read_required_csv(paths$diagnostics, "sensitivity_validation_checks.csv")
response_diagnostics <- read_required_csv(paths$diagnostics, "response_ipw_diagnostics.csv")

standardization_specification_results <- read_required_csv(paths$tables, "standardization_specification_results.csv")
standardization_alternative_contrasts <- read_required_csv(paths$tables, "standardization_alternative_contrasts.csv")
standardization_draws <- read_required_csv(paths$diagnostics, "standardization_support_by_draw.csv")
standardization_specification_diagnostics <- read_required_csv(
  paths$diagnostics, "standardization_specification_diagnostics_by_draw.csv"
)
standardization_model_diagnostics <- read_required_csv(
  paths$diagnostics, "standardization_model_diagnostics_by_draw.csv"
)
standardization_coefficient_diagnostics <- read_required_csv(
  paths$diagnostics, "standardization_coefficients_by_draw.csv"
)
standardization_replicate_diagnostics <- read_required_csv(
  paths$diagnostics, "standardization_replicate_diagnostics.csv"
)
standardization_covariance_validation <- read_required_csv(
  paths$diagnostics, "standardization_covariance_validation.csv"
)
pairwise_standardized <- read_required_csv(
  paths$tables, "standardized_pairwise_contrasts_exploratory.csv"
)

threshold_sensitivity <- read_required_csv(paths$tables, "threshold_sensitivity.csv")
hysteresis_sensitivity <- read_required_csv(paths$tables, "hysteresis_sensitivity.csv")
observed_ipw <- read_required_csv(paths$tables, "observed_response_ipw_sensitivity.csv")
mnar_sensitivity <- read_required_csv(paths$tables, "mnar_pattern_mixture_sensitivity.csv")
mnar_tipping <- read_required_csv(paths$tables, "mnar_tipping_points.csv")
measurement_crossing <- read_required_csv(paths$tables, "measurement_crossing_diagnostic.csv")
g89_domains <- read_required_csv(paths$tables, "g89_domain_extension.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"
  ),
  legacy_support = c(
    "draw", "system", "risk_n", "support_retained_weight", "support_effective_n",
    "math_lower", "math_upper", "science_lower", "science_upper", "model_iterations",
    "replicate_expected", "replicate_attempted", "replicate_completed"
  ),
  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"
  ),
  coefficient_diagnostics = c("draw", "surface", "system", "term", "estimate"),
  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"),
  pairwise = c("system_1", "system_2", "difference", "standard_error", "z", "p_value", "p_holm")
)

expected_standardization_keys <- bind_rows(lapply(names(standardization_specifications), function(specification) {
  data.frame(
    specification = specification,
    system = 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_legacy_keys <- bind_rows(lapply(seq_len(config$pv_count), function(draw) {
  data.frame(draw = draw, system = config$g45_systems, stringsAsFactors = FALSE)
}))
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_terms <- list(
  quadratic_interaction = c(
    "intercept", "math_deficit", "science_deficit", "math_deficit_sq",
    "science_deficit_sq", "deficit_interaction"
  ),
  additive_linear = c("intercept", "math_deficit", "science_deficit"),
  linear_interaction = c("intercept", "math_deficit", "science_deficit", "deficit_interaction")
)
expected_coefficient_keys <- bind_rows(lapply(seq_len(config$pv_count), function(draw) {
  bind_rows(lapply(expected_surfaces, function(surface) {
    bind_rows(lapply(config$g45_systems, function(system) {
      data.frame(
        draw = draw, surface = surface, system = system, term = expected_terms[[surface]],
        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)
}))
expected_pairwise_keys <- bind_rows(lapply(seq_len(length(config$g45_systems) - 1L), function(first) {
  data.frame(
    system_1 = config$g45_systems[first],
    system_2 = config$g45_systems[(first + 1L):length(config$g45_systems)],
    stringsAsFactors = FALSE
  )
}))
expected_validation_checks <- c(
  "primary raw estimates reproduce transition pipeline",
  "primary raw SEs reproduce transition pipeline",
  "all component estimates are finite probabilities",
  "cross-specification covariance is finite and symmetric",
  "cross-specification covariance is positive semidefinite",
  "component decomposition identity holds",
  "all response surfaces converged at full rank",
  "strict common cells satisfy n and ESS criteria",
  "all targets sum to one",
  "all fail-fast replicate stages and diagnostics completed",
  "artifact dimensions match frozen specifications"
)

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_draws, standardization_schemas$legacy_support, "Legacy standardization support")
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_coefficient_diagnostics, standardization_schemas$coefficient_diagnostics, "Standardization coefficient diagnostics")
exact_schema(standardization_replicate_diagnostics, standardization_schemas$replicate_diagnostics, "Standardization replicate diagnostics")
exact_schema(standardization_covariance_validation, standardization_schemas$validation, "Standardization covariance validation")
exact_schema(pairwise_standardized, standardization_schemas$pairwise, "Exploratory standardized pairwise contrasts")

require_gate(nrow(standardization_specification_results) == 53L, "Standardization specification results must have 53 rows")
require_gate(nrow(standardization_alternative_contrasts) == 44L, "Standardization alternative contrasts must have 44 rows")
require_gate(nrow(standardization_draws) == 45L, "Legacy standardization support must have 45 rows")
require_gate(nrow(standardization_specification_diagnostics) == 265L, "Standardization specification diagnostics must have 265 rows")
require_gate(nrow(standardization_model_diagnostics) == 135L, "Standardization model diagnostics must have 135 rows")
require_gate(nrow(standardization_coefficient_diagnostics) == 585L, "Standardization coefficient diagnostics must have 585 rows")
require_gate(nrow(standardization_replicate_diagnostics) == 30L, "Standardization replicate diagnostics must have 30 rows")
require_gate(nrow(standardization_covariance_validation) == 11L, "Standardization covariance validation must have 11 rows")
require_gate(nrow(pairwise_standardized) == 36L, "Exploratory standardized pairwise table must have 36 rows")

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_draws, expected_legacy_keys, c("draw", "system"), "Legacy standardization support")
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_coefficient_diagnostics, expected_coefficient_keys, c("draw", "surface", "system", "term"), "Standardization coefficient diagnostics")
exact_ordered_keys(standardization_replicate_diagnostics, expected_replicate_keys, c("draw", "specification"), "Standardization replicate diagnostics")
exact_ordered_keys(pairwise_standardized, expected_pairwise_keys, c("system_1", "system_2"), "Exploratory standardized pairwise contrasts")
require_gate(
  !anyDuplicated(standardization_covariance_validation$check) &&
    identical(standardization_covariance_validation$check, expected_validation_checks),
  "Standardization producer validation keys are not exact or ordered"
)

expected_rds_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_full_parameter_names <- c(
  paste0("raw_full::", config$g45_systems),
  unlist(lapply(names(standardization_specifications), function(specification) {
    unlist(lapply(specification_systems(specification), function(system) {
      paste(specification, c("raw_support", "model_own_support", "standardized_common"), system, sep = "::")
    }), use.names = FALSE)
  }), use.names = FALSE)
)
expected_primary_parameter_names <- c(
  paste0("raw::", config$g45_systems), paste0("standardized::", config$g45_systems)
)

validate_combined_object <- function(object, expected_names, expected_draws, label) {
  require_gate(identical(names(object), expected_combined_names), paste0(label, " member names are invalid"))
  parameter_count <- length(expected_names)
  require_gate(identical(names(object$estimate), expected_names), paste0(label, " estimate names are invalid"))
  require_gate(identical(names(object$se), expected_names), paste0(label, " SE names are invalid"))
  for (matrix_name in c("within", "between", "variance")) {
    matrix_value <- object[[matrix_name]]
    require_gate(
      identical(dim(matrix_value), c(parameter_count, parameter_count)) &&
        identical(rownames(matrix_value), expected_names) &&
        identical(colnames(matrix_value), expected_names),
      paste0(label, " ", matrix_name, " dimensions or names are invalid")
    )
    require_gate(all(is.finite(matrix_value)), paste0(label, " ", matrix_name, " is nonfinite"))
  }
  require_gate(
    identical(dim(object$draw_estimates), c(expected_draws, parameter_count)) &&
      identical(colnames(object$draw_estimates), expected_names) &&
      all(is.finite(object$draw_estimates)),
    paste0(label, " draw estimates are invalid")
  )
  require_gate(all(is.finite(object$estimate)) && all(is.finite(object$se)), paste0(label, " estimates or SEs are nonfinite"))
}

require_gate(identical(names(standardization), expected_rds_names), "Standardization RDS top-level names are invalid")
require_gate(identical(standardization$version, "2.1-auditability"), "Standardization RDS version is not 2.1-auditability")
require_gate(identical(standardization$threshold, config$primary_benchmark), "Standardization threshold does not match configuration")
require_gate(identical(primary$threshold, config$primary_benchmark), "Primary threshold does not match configuration")
require_gate(identical(standardization$systems, config$g45_systems), "Standardization systems or order do not match configuration")
require_gate(
  identical(names(standardization$specifications), names(standardization_specifications)) &&
    identical(standardization$specifications, standardization_specifications),
  "Standardization specification definitions or order do not match configuration"
)
validate_combined_object(standardization$combined, expected_primary_parameter_names, config$pv_count, "Primary standardization combined object")
validate_combined_object(standardization$combined_all, expected_full_parameter_names, config$pv_count, "Cross-specification combined object")
require_gate(length(standardization$combined_all$estimate) == 168L, "Cross-specification combined object must have 168 parameters")

combined_all_variance <- standardization$combined_all$variance
combined_all_asymmetry <- max(abs(combined_all_variance - t(combined_all_variance)))
combined_all_minimum_eigenvalue <- min(eigen(
  (combined_all_variance + t(combined_all_variance)) / 2,
  symmetric = TRUE, only.values = TRUE
)$values)
combined_all_se_diagonal_difference <- max(abs(
  standardization$combined_all$se - sqrt(pmax(0, diag(combined_all_variance)))
))
require_gate(combined_all_asymmetry <= 1e-10, "Cross-specification covariance is asymmetric")
require_gate(combined_all_minimum_eigenvalue >= -1e-8, "Cross-specification covariance is not positive semidefinite")
require_gate(combined_all_se_diagonal_difference <= 1e-10, "Cross-specification SEs do not match the covariance diagonal")

full_primary_names <- c(
  paste0("raw_full::", config$g45_systems),
  paste("primary", "standardized_common", config$g45_systems, sep = "::")
)
require_gate(
  max(abs(standardization$combined$estimate - standardization$combined_all$estimate[full_primary_names])) <= 1e-10,
  "Primary combined estimates disagree with combined_all"
)
require_gate(
  max(abs(standardization$combined$variance - standardization$combined_all$variance[full_primary_names, full_primary_names])) <= 1e-10,
  "Primary combined covariance disagrees with combined_all"
)

frames_agree(standardization$diagnostics, standardization_draws, "Legacy standardization support")
frames_agree(standardization$specification_diagnostics, standardization_specification_diagnostics, "Standardization specification diagnostics")
frames_agree(standardization$model_diagnostics, standardization_model_diagnostics, "Standardization model diagnostics")
frames_agree(standardization$replicate_diagnostics, standardization_replicate_diagnostics, "Standardization replicate diagnostics")
frames_agree(standardization$validation, standardization_covariance_validation, "Standardization covariance validation")

for (row in seq_len(nrow(standardization_specification_results))) {
  item <- standardization_specification_results[row, ]
  definition <- standardization_specifications[[item$specification]]
  require_gate(
    identical(item$surface, definition$surface) && identical(item$target, definition$target) &&
      identical(item$support, definition$support),
    paste0("Standardization specification metadata mismatch at row ", row)
  )
  raw_name <- paste0("raw_full::", item$system)
  component_names <- paste(
    item$specification, c("raw_support", "model_own_support", "standardized_common"),
    item$system, sep = "::"
  )
  source_names <- c(raw_name, component_names)
  source_estimates <- c(item$raw_full, item$raw_support, item$model_own_support, item$standardized_common)
  source_ses <- c(item$raw_full_se, item$raw_support_se, item$model_own_support_se, item$standardized_common_se)
  require_gate(
    max(abs(source_estimates - standardization$combined_all$estimate[source_names])) <= 1e-10 &&
      max(abs(source_ses - standardization$combined_all$se[source_names])) <= 1e-10,
    paste0("Standardization specification results disagree with combined_all at row ", row)
  )
}
require_gate(
  all(is.finite(standardization_specification_results$decomposition_residual)) &&
    max(abs(standardization_specification_results$decomposition_residual)) <= 1e-12,
  "Standardization decomposition identity failed"
)

for (row in seq_len(nrow(standardization_alternative_contrasts))) {
  item <- standardization_alternative_contrasts[row, ]
  alternative_name <- paste(item$specification, "standardized_common", item$system, sep = "::")
  primary_name <- paste("primary", "standardized_common", item$system, sep = "::")
  require_gate(
    abs(item$alternative - standardization$combined_all$estimate[alternative_name]) <= 1e-10 &&
      abs(item$primary - standardization$combined_all$estimate[primary_name]) <= 1e-10 &&
      abs(item$alternative_minus_primary - (item$alternative - item$primary)) <= 1e-12,
    paste0("Alternative contrast values disagree with combined_all at row ", row)
  )
}
require_gate(
  all(is.finite(as.matrix(standardization_alternative_contrasts[, c(
    "alternative", "primary", "alternative_minus_primary", "standard_error", "ci_lower", "ci_upper"
  )]))) && all(standardization_alternative_contrasts$standard_error >= 0) &&
    all(standardization_alternative_contrasts$ci_lower <= standardization_alternative_contrasts$alternative_minus_primary) &&
    all(standardization_alternative_contrasts$ci_upper >= standardization_alternative_contrasts$alternative_minus_primary),
  "Covariance-aware alternative contrast intervals are invalid"
)

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)
require_gate(
  standardization_replicates_complete && standardization_replicate_diagnostics_valid,
  "A fail-fast standardization replicate stage or diagnostic is incomplete"
)
require_gate(all(standardization_covariance_validation$passed), "A standardization producer validation failed")
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"
)
require_gate(
  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)) &&
    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) &&
    all(standardization_model_diagnostics$predicted_probability_minimum >= 0) &&
    all(standardization_model_diagnostics$predicted_probability_maximum <= 1),
  "A full-sample standardization model or clipping diagnostic is invalid"
)
require_gate(
  all(is.finite(standardization_coefficient_diagnostics$estimate)),
  "A standardization coefficient diagnostic is nonfinite"
)

expected_surface_for_spec <- vapply(
  standardization_specification_diagnostics$specification,
  function(specification) standardization_specifications[[specification]]$surface,
  character(1)
)
expected_target_for_spec <- vapply(
  standardization_specification_diagnostics$specification,
  function(specification) standardization_specifications[[specification]]$target,
  character(1)
)
expected_support_for_spec <- vapply(
  standardization_specification_diagnostics$specification,
  function(specification) standardization_specifications[[specification]]$support,
  character(1)
)
require_gate(
  identical(standardization_specification_diagnostics$surface, unname(expected_surface_for_spec)) &&
    identical(standardization_specification_diagnostics$target, unname(expected_target_for_spec)) &&
    identical(standardization_specification_diagnostics$support, unname(expected_support_for_spec)),
  "Specification diagnostic definitions disagree with configuration"
)
require_gate(
  all(standardization_specification_diagnostics$risk_n > 0) &&
    all(standardization_specification_diagnostics$supported_n > 0) &&
    all(standardization_specification_diagnostics$supported_n <= standardization_specification_diagnostics$risk_n) &&
    all(standardization_specification_diagnostics$support_retained_weight > 0) &&
    all(standardization_specification_diagnostics$support_retained_weight <= 1) &&
    all(standardization_specification_diagnostics$support_effective_n > 0) &&
    all(abs(standardization_specification_diagnostics$target_weight_sum - 1) <= 1e-10) &&
    all(standardization_specification_diagnostics$target_system_mass > 0) &&
    all(standardization_specification_diagnostics$target_system_mass <= 1) &&
    all(standardization_specification_diagnostics$target_mass_retained_before_renormalization > 0) &&
    all(standardization_specification_diagnostics$target_mass_retained_before_renormalization <= 1) &&
    all(standardization_specification_diagnostics$math_lower < standardization_specification_diagnostics$math_upper) &&
    all(standardization_specification_diagnostics$science_lower < standardization_specification_diagnostics$science_upper),
  "A standardization support or target diagnostic violates its constraints"
)
target_mass_sums <- aggregate(
  target_system_mass ~ draw + specification,
  standardization_specification_diagnostics,
  sum
)
require_gate(all(abs(target_mass_sums$target_system_mass - 1) <= 1e-10), "Target system masses do not sum to one")
equal_target_rows <- standardization_specification_diagnostics$target == "equal_system"
equal_target_counts <- vapply(
  standardization_specification_diagnostics$specification[equal_target_rows],
  function(specification) length(specification_systems(specification)),
  integer(1)
)
require_gate(
  all(abs(
    standardization_specification_diagnostics$target_system_mass[equal_target_rows] -
      1 / equal_target_counts
  ) <= 1e-10),
  "Equal-system target masses are invalid"
)
strict_diagnostics <- standardization_specification_diagnostics[
  standardization_specification_diagnostics$specification == "strict_common_cells", , drop = FALSE
]
require_gate(
  nrow(strict_diagnostics) == config$pv_count * length(config$g45_systems) &&
    all(strict_diagnostics$grid_cells_total == config$strict_cell_grid^2) &&
    all(strict_diagnostics$grid_cells_retained > 0) &&
    all(strict_diagnostics$minimum_common_cell_n >= config$strict_cell_min_n) &&
    all(strict_diagnostics$minimum_common_cell_ess >= config$strict_cell_min_ess),
  "Strict common-cell diagnostics fail the frozen n or ESS criteria"
)
require_gate(
  all(standardization_draws$replicate_attempted == standardization_draws$replicate_expected) &&
    all(standardization_draws$replicate_completed == standardization_draws$replicate_expected),
  "Legacy standardization replicate completion gate failed"
)

require_gate(!is.null(sensitivity$standardization), "Sensitivity RDS lacks the compact standardization member")
require_gate(
  identical(sensitivity$standardization$version, standardization$version) &&
    identical(sensitivity$standardization$threshold, standardization$threshold),
  "Sensitivity RDS standardization metadata disagrees with the producer bundle"
)
require_gate(
  !any(names(sensitivity$standardization) %in% c(
    "combined", "combined_all", "covariance", "covariances", "draw_estimates",
    "coefficient_diagnostics", "full_diagnostics"
  )),
  "Sensitivity RDS duplicates prohibited full standardization objects"
)
frames_agree(
  sensitivity$standardization$specification_results,
  standardization_specification_results,
  "Sensitivity compact standardization results"
)
frames_agree(
  sensitivity$standardization$alternative_contrasts,
  standardization_alternative_contrasts,
  "Sensitivity compact standardization contrasts"
)
frames_agree(
  sensitivity$standardization$validation,
  standardization_covariance_validation,
  "Sensitivity compact standardization validation"
)

require_gate(nrow(sensitivity_validation) == 13L && all(sensitivity_validation$passed), "Sensitivity validation gate failed")
require_gate(!anyDuplicated(sensitivity_validation$check), "Sensitivity validation checks are duplicated")
require_gate(all(official_validation$passes_rounding_check), "Official-result reproduction gate failed")

required_official_columns <- c(
  "official_t0", "official_t1", "official_change", "official_at_or_above_400_t0",
  "official_at_or_above_400_t1", "official_t0_se", "official_t1_se", "official_change_se",
  "official_at_or_above_400_se_t0", "official_at_or_above_400_se_t1",
  "benchmark_se_suppressed_t0", "benchmark_se_suppressed_t1", "diff_t0", "diff_t1",
  "diff_change", "diff_benchmark_t0_pp", "diff_benchmark_t1_pp", "diff_t0_se",
  "diff_t1_se", "diff_change_se", "diff_benchmark_t0_se_pp", "diff_benchmark_t1_se_pp"
)
require_gate(all(required_official_columns %in% names(official_validation)), "Official validation schema lacks required fields")
component_pass <- function(official, difference, tolerance) {
  included <- is.finite(official)
  abs(difference[included]) <= tolerance
}
official_point_pass <- c(
  component_pass(official_validation$official_t0, official_validation$diff_t0, 0.5),
  component_pass(official_validation$official_t1, official_validation$diff_t1, 0.5),
  component_pass(official_validation$official_change, official_validation$diff_change, 0.5),
  component_pass(official_validation$official_at_or_above_400_t0, official_validation$diff_benchmark_t0_pp, 0.5),
  component_pass(official_validation$official_at_or_above_400_t1, official_validation$diff_benchmark_t1_pp, 0.5)
)
official_se_pass <- c(
  component_pass(official_validation$official_t0_se, official_validation$diff_t0_se, 0.05),
  component_pass(official_validation$official_t1_se, official_validation$diff_t1_se, 0.05),
  component_pass(official_validation$official_change_se, official_validation$diff_change_se, 0.05),
  component_pass(official_validation$official_at_or_above_400_se_t0, official_validation$diff_benchmark_t0_se_pp, 0.05),
  component_pass(official_validation$official_at_or_above_400_se_t1, official_validation$diff_benchmark_t1_se_pp, 0.05)
)
suppressed_benchmark_ses <- sum(
  official_validation$benchmark_se_suppressed_t0 %in% TRUE,
  official_validation$benchmark_se_suppressed_t1 %in% TRUE
)
benchmark_t0_available <- is.finite(official_validation$official_at_or_above_400_t0)
benchmark_t1_available <- is.finite(official_validation$official_at_or_above_400_t1)
benchmark_t0_se_handled <- !benchmark_t0_available |
  is.finite(official_validation$official_at_or_above_400_se_t0) |
  official_validation$benchmark_se_suppressed_t0 %in% TRUE
benchmark_t1_se_handled <- !benchmark_t1_available |
  is.finite(official_validation$official_at_or_above_400_se_t1) |
  official_validation$benchmark_se_suppressed_t1 %in% TRUE
require_gate(length(official_point_pass) == 108L && all(official_point_pass), "Official point-estimate validation must be 108 of 108")
require_gate(length(official_se_pass) == 102L && all(official_se_pass), "Official numeric-SE validation must be 102 of 102")
require_gate(
  suppressed_benchmark_ses == 6L && all(benchmark_t0_se_handled) && all(benchmark_t1_se_handled),
  "Suppressed official benchmark SEs are not handled exactly"
)

standardization_contract_passed <- TRUE
standardization_producer_passed <- all(standardization_covariance_validation$passed)
standardization_models_passed <- all(
  standardization_model_diagnostics$converged &
    standardization_model_diagnostics$model_rank == standardization_model_diagnostics$design_columns
)
strict_support_passed <- all(
  strict_diagnostics$minimum_common_cell_n >= config$strict_cell_min_n &
    strict_diagnostics$minimum_common_cell_ess >= config$strict_cell_min_ess
)

# Legacy primary overlap calculations are retained unchanged in substance.
points_in_convex_polygon <- function(point_x, point_y, polygon_x, polygon_y) {
  if (length(polygon_x) < 3L) return(rep(FALSE, length(point_x)))
  next_index <- c(seq.int(2L, length(polygon_x)), 1L)
  signed_area_twice <- sum(polygon_x * polygon_y[next_index] - polygon_y * polygon_x[next_index])
  orientation <- sign(signed_area_twice)
  if (orientation == 0) return(rep(FALSE, length(point_x)))
  coordinate_scale <- max(diff(range(polygon_x)), diff(range(polygon_y)), 1)
  tolerance <- 1e-9 * coordinate_scale^2
  inside <- rep(TRUE, length(point_x))
  for (edge in seq_along(polygon_x)) {
    following <- next_index[edge]
    cross_product <- (polygon_x[following] - polygon_x[edge]) * (point_y - polygon_y[edge]) -
      (polygon_y[following] - polygon_y[edge]) * (point_x - polygon_x[edge])
    if (orientation > 0) {
      inside <- inside & cross_product >= -tolerance
    } else {
      inside <- inside & cross_product <= tolerance
    }
  }
  inside
}

joint_overlap_for_draw <- function(data, draw, systems, bounds_row) {
  map <- pv_variable_map("g45", draw)
  state_t0 <- achievement_state(data[[map$math_t0]], data[[map$science_t0]], config$primary_benchmark)
  risk <- state_t0 == "BB"
  pooled <- data.frame(
    system = data$system[risk],
    math = data[[map$math_t0]][risk],
    science = data[[map$science_t0]][risk],
    weight = data$weight[risk],
    stringsAsFactors = FALSE
  )
  eligible <- pooled$math >= bounds_row$math_lower & pooled$math <= bounds_row$math_upper &
    pooled$science >= bounds_row$science_lower & pooled$science <= bounds_row$science_upper
  reference_weight <- numeric(nrow(pooled))
  for (system in systems) {
    included <- eligible & pooled$system == system
    included_weight <- sum(pooled$weight[included])
    if (included_weight <= 0) stop("Empty overlap support for ", system, " in draw ", draw)
    reference_weight[included] <- pooled$weight[included] / included_weight / length(systems)
  }
  if (abs(sum(reference_weight) - 1) > 1e-10) stop("Reference weights do not sum to one")
  reference_index <- which(reference_weight > 0)
  math_breaks <- weighted_quantile(
    pooled$math[reference_index], reference_weight[reference_index], seq(0, 1, length.out = 6)
  )
  science_breaks <- weighted_quantile(
    pooled$science[reference_index], reference_weight[reference_index], seq(0, 1, length.out = 6)
  )
  math_breaks[c(1L, length(math_breaks))] <- math_breaks[c(1L, length(math_breaks))] + c(-1e-8, 1e-8)
  science_breaks[c(1L, length(science_breaks))] <- science_breaks[c(1L, length(science_breaks))] + c(-1e-8, 1e-8)
  if (anyDuplicated(math_breaks) || anyDuplicated(science_breaks)) stop("Degenerate overlap grid in draw ", draw)
  math_bin <- cut(pooled$math, breaks = math_breaks, labels = FALSE, include.lowest = TRUE)
  science_bin <- cut(pooled$science, breaks = science_breaks, labels = FALSE, include.lowest = TRUE)
  cell <- (math_bin - 1L) * 5L + science_bin
  reference_mass_by_cell <- tapply(reference_weight, cell, sum)
  reference_cells <- as.integer(names(reference_mass_by_cell)[reference_mass_by_cell > 0])

  bind_rows(lapply(systems, function(system) {
    target <- eligible & pooled$system == system
    hull_index <- grDevices::chull(pooled$math[target], pooled$science[target])
    target_math <- pooled$math[target][hull_index]
    target_science <- pooled$science[target][hull_index]
    inside_hull <- points_in_convex_polygon(
      pooled$math[reference_index], pooled$science[reference_index], target_math, target_science
    )
    target_cell_counts <- table(factor(cell[target], levels = 1:25))
    supported_three <- as.integer(names(target_cell_counts)[target_cell_counts >= 3])
    supported_five <- as.integer(names(target_cell_counts)[target_cell_counts >= 5])
    mass_in_cells <- function(supported_cells) {
      if (!length(supported_cells)) return(0)
      pmin(1, pmax(0, sum(reference_weight[cell %in% supported_cells], na.rm = TRUE)))
    }
    data.frame(
      draw = draw,
      system = system,
      retained_student_n = sum(target),
      convex_hull_vertices = length(hull_index),
      reference_mass_within_convex_hull = pmin(1, pmax(0, sum(reference_weight[reference_index][inside_hull]))),
      reference_mass_in_cells_with_n_at_least_3 = mass_in_cells(supported_three),
      reference_mass_in_cells_with_n_at_least_5 = mass_in_cells(supported_five),
      reference_occupied_cells = length(reference_cells),
      target_cells_with_n_at_least_3 = length(intersect(reference_cells, supported_three)),
      target_cells_with_n_at_least_5 = length(intersect(reference_cells, supported_five)),
      stringsAsFactors = FALSE
    )
  }))
}

overlap_rows <- lapply(seq_len(config$pv_count), function(draw) {
  bounds_row <- standardization_draws[standardization_draws$draw == draw, ][1, ]
  joint_overlap_for_draw(g45, draw, config$g45_systems, bounds_row)
})
joint_overlap <- bind_rows(overlap_rows)
standardization_overlap_internal <- merge(
  standardization_draws, joint_overlap, by = c("draw", "system"), all = TRUE, sort = FALSE
)
standardization_overlap_internal <- standardization_overlap_internal[order(
  standardization_overlap_internal$draw,
  match(standardization_overlap_internal$system, config$g45_systems)
), ]
rownames(standardization_overlap_internal) <- NULL
require_gate(
  nrow(standardization_overlap_internal) == 45L &&
    all(is.finite(standardization_overlap_internal$reference_mass_within_convex_hull)) &&
    all(standardization_overlap_internal$reference_mass_within_convex_hull >= 0) &&
    all(standardization_overlap_internal$reference_mass_within_convex_hull <= 1),
  "Legacy joint-overlap diagnostics failed"
)

one_row <- function(data, condition, label) {
  result <- data[condition, , drop = FALSE]
  if (nrow(result) != 1L) stop("Expected exactly one row for ", label, "; found ", nrow(result))
  result
}

add_mixed_estimate_fields <- function(data, probability_metrics) {
  probability <- data$metric %in% probability_metrics
  data$metric_unit <- ifelse(probability, "percent / percentage points", "TIMSS points")
  data$estimate_percent <- ifelse(probability, 100 * data$estimate, NA_real_)
  data$standard_error_pp <- ifelse(probability, 100 * data$standard_error, NA_real_)
  data$ci_lower_percent <- ifelse(probability, 100 * data$ci_lower, NA_real_)
  data$ci_upper_percent <- ifelse(probability, 100 * data$ci_upper, NA_real_)
  data$estimate_points <- ifelse(probability, NA_real_, data$estimate)
  data$standard_error_points <- ifelse(probability, NA_real_, data$standard_error)
  data$ci_lower_points <- ifelse(probability, NA_real_, data$ci_lower)
  data$ci_upper_points <- ifelse(probability, NA_real_, data$ci_upper)
  data
}

message("Constructing the keyed 19-table publication catalog")
table_1 <- data.frame(
  grade_band_id = inventory$band,
  grade_band_label = publication_band(inventory$band),
  system_id = inventory$system,
  system_label = publication_system(inventory$system),
  students = inventory$students,
  schools = inventory$schools,
  jackknife_zones = inventory$zones,
  weighted_population = round(inventory$weighted_population),
  both_wave_students = inventory$both_wave_n,
  baseline_only_students = inventory$baseline_only_n,
  followup_only_students = inventory$followup_only_n,
  both_wave_weighted_percent = round(inventory$both_wave_weighted_pct, 1),
  missing_pv_cells = inventory$pv_missing_cells,
  stringsAsFactors = FALSE
)

primary_specification <- standardization_specification_results[
  standardization_specification_results$specification == "primary", , drop = FALSE
]
primary_specification <- primary_specification[match(config$g45_systems, primary_specification$system), ]
table_2 <- bind_rows(lapply(config$g45_systems, function(system) {
  source <- one_row(primary_specification, primary_specification$system == system, paste(system, "primary standardization"))
  support <- standardization_draws[standardization_draws$system == system, , drop = FALSE]
  raw <- source$raw_full
  raw_se <- source$raw_full_se
  standardized_value <- source$standardized_common
  standardized_se <- source$standardized_common_se
  data.frame(
    system_id = system,
    system_label = publication_system(system),
    analysis_status = "Primary",
    threshold = config$primary_benchmark,
    mean_unweighted_baseline_BB_n = mean(support$risk_n),
    baseline_BB_weighted_percent = 100 * primary$g45[[system]]$estimate["baseline_BB"],
    raw_recovery_percent = 100 * raw,
    raw_standard_error_pp = 100 * raw_se,
    raw_ci_lower_percent = 100 * (raw - 1.96 * raw_se),
    raw_ci_upper_percent = 100 * (raw + 1.96 * raw_se),
    standardized_recovery_percent = 100 * standardized_value,
    standardized_standard_error_pp = 100 * standardized_se,
    standardized_ci_lower_percent = 100 * (standardized_value - 1.96 * standardized_se),
    standardized_ci_upper_percent = 100 * (standardized_value + 1.96 * standardized_se),
    standardized_minus_raw_pp = 100 * source$standardized_minus_raw_full,
    difference_standard_error_pp = 100 * source$standardized_minus_raw_full_se,
    difference_ci_lower_pp = 100 * source$total_ci_lower,
    difference_ci_upper_pp = 100 * source$total_ci_upper,
    minimum_support_retained_percent = 100 * min(support$support_retained_weight),
    minimum_support_effective_n = min(support$support_effective_n),
    raw_formatted = format_percent_ci(raw, raw - 1.96 * raw_se, raw + 1.96 * raw_se),
    standardized_formatted = format_percent_ci(
      standardized_value,
      standardized_value - 1.96 * standardized_se,
      standardized_value + 1.96 * standardized_se
    ),
    stringsAsFactors = FALSE
  )
}))

robustness_rows <- lapply(config$g45_systems, function(system) {
  primary_result <- primary$g45[[system]]
  primary_standardized <- one_row(primary_specification, primary_specification$system == system, system)
  threshold_400 <- one_row(
    threshold_sensitivity,
    threshold_sensitivity$band == "g45" & threshold_sensitivity$system == system &
      threshold_sensitivity$threshold == 400 & threshold_sensitivity$metric == "BB_to_AA",
    paste(system, "threshold 400")
  )
  buffered <- one_row(
    hysteresis_sensitivity,
    hysteresis_sensitivity$band == "g45" & hysteresis_sensitivity$system == system &
      hysteresis_sensitivity$specification == "primary_450_to_500" &
      hysteresis_sensitivity$metric == "dual_recovery",
    paste(system, "hysteresis")
  )
  both <- one_row(
    observed_ipw,
    observed_ipw$band == "g45" & observed_ipw$system == system &
      observed_ipw$method == "both_wave_observed" & observed_ipw$metric == "BB_to_AA",
    paste(system, "both-wave")
  )
  ipw <- one_row(
    observed_ipw,
    observed_ipw$band == "g45" & observed_ipw$system == system &
      observed_ipw$method == "baseline_participant_response_ipw" & observed_ipw$metric == "BB_to_AA",
    paste(system, "IPW")
  )
  mnar <- one_row(
    mnar_sensitivity,
    mnar_sensitivity$band == "g45" & mnar_sensitivity$system == system &
      mnar_sensitivity$scenario == "combined_pattern_shift" &
      mnar_sensitivity$delta_points_per_affected_domain == -40 &
      mnar_sensitivity$metric == "BB_to_AA",
    paste(system, "MNAR")
  )
  crossing <- one_row(
    measurement_crossing,
    measurement_crossing$band == "g45" & measurement_crossing$system == system &
      measurement_crossing$comparison == "followup_same_wave_cross_pv",
    paste(system, "cross-PV")
  )
  data.frame(
    system_id = system,
    system_label = publication_system(system),
    analysis_status = "Sensitivity",
    primary_475_percent = 100 * primary_result$estimate["BB_to_AA"],
    primary_475_se_pp = 100 * primary_result$se["BB_to_AA"],
    standardized_475_percent = 100 * primary_standardized$standardized_common,
    standardized_475_se_pp = 100 * primary_standardized$standardized_common_se,
    threshold_400_percent = 100 * threshold_400$estimate,
    threshold_400_se_pp = 100 * threshold_400$standard_error,
    buffered_450_to_500_percent = 100 * buffered$estimate,
    buffered_450_to_500_se_pp = 100 * buffered$standard_error,
    both_wave_observed_percent = 100 * both$estimate,
    both_wave_observed_se_pp = 100 * both$standard_error,
    response_ipw_percent = 100 * ipw$estimate,
    response_ipw_se_pp = 100 * ipw$standard_error,
    combined_mnar_minus40_percent = 100 * mnar$estimate,
    combined_mnar_minus40_se_pp = 100 * mnar$standard_error,
    followup_same_wave_cross_pv_percent = 100 * crossing$estimate,
    followup_same_wave_cross_pv_sampling_se_pp = 100 * crossing$standard_error,
    stringsAsFactors = FALSE
  )
})
table_3 <- bind_rows(robustness_rows)

primary_long_rows <- list()
row <- 1L
for (band in c("g45", "g89")) {
  systems <- if (band == "g45") config$g45_systems else config$g89_systems
  for (system in systems) {
    result <- primary[[band]][[system]]
    for (parameter in names(result$estimate)) {
      primary_long_rows[[row]] <- data.frame(
        grade_band_id = band,
        grade_band_label = publication_band(band),
        system_id = system,
        system_label = publication_system(system),
        analysis_scope = publication_scope(band),
        analysis_status = "Primary",
        parameter = parameter,
        estimate = result$estimate[parameter],
        standard_error = result$se[parameter],
        ci_lower = result$estimate[parameter] - 1.96 * result$se[parameter],
        ci_upper = result$estimate[parameter] + 1.96 * result$se[parameter],
        stringsAsFactors = FALSE
      )
      row <- row + 1L
    }
  }
}
primary_long <- bind_rows(primary_long_rows)
transition_pattern <- "^(BB|MB|SB|AA)_to_(BB|MB|SB|AA)$"
table_s1 <- primary_long[grepl(transition_pattern, primary_long$parameter), ]
transition_parts <- strsplit(table_s1$parameter, "_to_", fixed = TRUE)
table_s1$origin_state <- vapply(transition_parts, `[[`, character(1), 1L)
table_s1$destination_state <- vapply(transition_parts, `[[`, character(1), 2L)
table_s1$estimate_percent <- 100 * table_s1$estimate
table_s1$standard_error_pp <- 100 * table_s1$standard_error
table_s1$ci_lower_percent <- 100 * table_s1$ci_lower
table_s1$ci_upper_percent <- 100 * table_s1$ci_upper
table_s1 <- table_s1[, c(
  "grade_band_id", "grade_band_label", "analysis_scope", "system_id", "system_label", "analysis_status",
  "origin_state", "destination_state", "estimate_percent", "standard_error_pp",
  "ci_lower_percent", "ci_upper_percent"
)]

state_pattern <- "^(baseline|followup)_(BB|MB|SB|AA)$"
table_s2 <- primary_long[grepl(state_pattern, primary_long$parameter), ]
state_parts <- strsplit(table_s2$parameter, "_", fixed = TRUE)
table_s2$wave <- vapply(state_parts, `[[`, character(1), 1L)
table_s2$state <- vapply(state_parts, `[[`, character(1), 2L)
table_s2$estimate_percent <- 100 * table_s2$estimate
table_s2$standard_error_pp <- 100 * table_s2$standard_error
table_s2 <- table_s2[, c(
  "grade_band_id", "grade_band_label", "analysis_scope", "system_id", "system_label", "analysis_status",
  "wave", "state", "estimate_percent", "standard_error_pp"
)]

continuous_parameters <- c("weakest_margin_change_BB", "math_change_BB", "science_change_BB")
table_s3 <- primary_long[primary_long$parameter %in% continuous_parameters, ]
names(table_s3)[names(table_s3) == "parameter"] <- "change_estimand"
table_s3 <- table_s3[, c(
  "grade_band_id", "grade_band_label", "analysis_scope", "system_id", "system_label", "analysis_status",
  "change_estimand", "estimate", "standard_error", "ci_lower", "ci_upper"
)]

specification_id <- standardization_specification_results$specification
table_s4 <- data.frame(
  specification_id = specification_id,
  specification_label = unname(specification_labels[specification_id]),
  analysis_status = unname(standardization_status[specification_id]),
  system_id = standardization_specification_results$system,
  system_label = publication_system(standardization_specification_results$system),
  surface_id = standardization_specification_results$surface,
  surface_label = unname(surface_labels[standardization_specification_results$surface]),
  target_id = standardization_specification_results$target,
  target_label = unname(target_labels[standardization_specification_results$target]),
  support_id = standardization_specification_results$support,
  support_label = unname(support_labels[standardization_specification_results$support]),
  raw_full_percent = 100 * standardization_specification_results$raw_full,
  raw_full_se_pp = 100 * standardization_specification_results$raw_full_se,
  raw_support_percent = 100 * standardization_specification_results$raw_support,
  raw_support_se_pp = 100 * standardization_specification_results$raw_support_se,
  model_own_support_percent = 100 * standardization_specification_results$model_own_support,
  model_own_support_se_pp = 100 * standardization_specification_results$model_own_support_se,
  standardized_common_percent = 100 * standardization_specification_results$standardized_common,
  standardized_common_se_pp = 100 * standardization_specification_results$standardized_common_se,
  support_restriction_component_pp = 100 * standardization_specification_results$support_restriction_component,
  support_restriction_se_pp = 100 * standardization_specification_results$support_restriction_se,
  model_calibration_component_pp = 100 * standardization_specification_results$model_calibration_component,
  model_calibration_se_pp = 100 * standardization_specification_results$model_calibration_se,
  target_distribution_component_pp = 100 * standardization_specification_results$target_distribution_component,
  target_distribution_se_pp = 100 * standardization_specification_results$target_distribution_se,
  standardized_minus_raw_full_pp = 100 * standardization_specification_results$standardized_minus_raw_full,
  standardized_minus_raw_full_se_pp = 100 * standardization_specification_results$standardized_minus_raw_full_se,
  total_ci_lower_pp = 100 * standardization_specification_results$total_ci_lower,
  total_ci_upper_pp = 100 * standardization_specification_results$total_ci_upper,
  decomposition_residual_pp = 100 * standardization_specification_results$decomposition_residual,
  stringsAsFactors = FALSE
)

alternative_id <- standardization_alternative_contrasts$specification
table_s5 <- data.frame(
  specification_id = alternative_id,
  specification_label = unname(specification_labels[alternative_id]),
  analysis_status = "Sensitivity",
  system_id = standardization_alternative_contrasts$system,
  system_label = publication_system(standardization_alternative_contrasts$system),
  alternative_standardized_percent = 100 * standardization_alternative_contrasts$alternative,
  primary_percent = 100 * standardization_alternative_contrasts$primary,
  alternative_minus_primary_pp = 100 * standardization_alternative_contrasts$alternative_minus_primary,
  covariance_aware_standard_error_pp = 100 * standardization_alternative_contrasts$standard_error,
  covariance_aware_ci_lower_pp = 100 * standardization_alternative_contrasts$ci_lower,
  covariance_aware_ci_upper_pp = 100 * standardization_alternative_contrasts$ci_upper,
  stringsAsFactors = FALSE
)

standardization_overlap_summary <- bind_rows(lapply(config$g45_systems, function(system) {
  subset <- standardization_overlap_internal[standardization_overlap_internal$system == system, ]
  data.frame(
    system_id = system,
    system_label = publication_system(system),
    analysis_status = "Primary",
    support_retained_percent_mean = 100 * mean(subset$support_retained_weight),
    support_retained_percent_minimum = 100 * min(subset$support_retained_weight),
    support_effective_n_mean = mean(subset$support_effective_n),
    support_effective_n_minimum = min(subset$support_effective_n),
    reference_mass_within_convex_hull_mean_percent = 100 * mean(subset$reference_mass_within_convex_hull),
    reference_mass_within_convex_hull_minimum_percent = 100 * min(subset$reference_mass_within_convex_hull),
    reference_mass_in_cells_n_at_least_3_mean_percent = 100 * mean(subset$reference_mass_in_cells_with_n_at_least_3),
    reference_mass_in_cells_n_at_least_5_mean_percent = 100 * mean(subset$reference_mass_in_cells_with_n_at_least_5),
    stringsAsFactors = FALSE
  )
}))
table_s6 <- data.frame(
  draw = standardization_overlap_internal$draw,
  system_id = standardization_overlap_internal$system,
  system_label = publication_system(standardization_overlap_internal$system),
  analysis_status = "Primary",
  risk_n = standardization_overlap_internal$risk_n,
  support_retained_percent = 100 * standardization_overlap_internal$support_retained_weight,
  support_effective_n = standardization_overlap_internal$support_effective_n,
  math_lower = standardization_overlap_internal$math_lower,
  math_upper = standardization_overlap_internal$math_upper,
  science_lower = standardization_overlap_internal$science_lower,
  science_upper = standardization_overlap_internal$science_upper,
  model_iterations = standardization_overlap_internal$model_iterations,
  replicate_expected = standardization_overlap_internal$replicate_expected,
  replicate_attempted = standardization_overlap_internal$replicate_attempted,
  replicate_completed = standardization_overlap_internal$replicate_completed,
  retained_student_n = standardization_overlap_internal$retained_student_n,
  convex_hull_vertices = standardization_overlap_internal$convex_hull_vertices,
  reference_mass_within_convex_hull_percent = 100 * standardization_overlap_internal$reference_mass_within_convex_hull,
  reference_mass_in_cells_with_n_at_least_3_percent = 100 * standardization_overlap_internal$reference_mass_in_cells_with_n_at_least_3,
  reference_mass_in_cells_with_n_at_least_5_percent = 100 * standardization_overlap_internal$reference_mass_in_cells_with_n_at_least_5,
  reference_occupied_cells = standardization_overlap_internal$reference_occupied_cells,
  target_cells_with_n_at_least_3 = standardization_overlap_internal$target_cells_with_n_at_least_3,
  target_cells_with_n_at_least_5 = standardization_overlap_internal$target_cells_with_n_at_least_5,
  stringsAsFactors = FALSE
)
table_s7 <- standardization_overlap_summary

probability_metrics <- c("BB_to_AA", "BB_to_BB", "bb_exit", "baseline_BB")
table_s8_source <- threshold_sensitivity
table_s8_source <- add_mixed_estimate_fields(table_s8_source, probability_metrics)
table_s8 <- data.frame(
  grade_band_id = table_s8_source$band,
  grade_band_label = publication_band(table_s8_source$band),
  system_id = table_s8_source$system,
  system_label = publication_system(table_s8_source$system),
  analysis_status = "Sensitivity",
  threshold = table_s8_source$threshold,
  metric = table_s8_source$metric,
  metric_unit = table_s8_source$metric_unit,
  risk_n_mean = table_s8_source$risk_n_mean,
  risk_n_min = table_s8_source$risk_n_min,
  risk_n_max = table_s8_source$risk_n_max,
  estimate_percent = table_s8_source$estimate_percent,
  standard_error_pp = table_s8_source$standard_error_pp,
  ci_lower_percent = table_s8_source$ci_lower_percent,
  ci_upper_percent = table_s8_source$ci_upper_percent,
  estimate_points = table_s8_source$estimate_points,
  standard_error_points = table_s8_source$standard_error_points,
  ci_lower_points = table_s8_source$ci_lower_points,
  ci_upper_points = table_s8_source$ci_upper_points,
  stringsAsFactors = FALSE
)

hysteresis_probability_metrics <- c("dual_recovery", "dual_persistence", "baseline_risk")
table_s9_source <- add_mixed_estimate_fields(hysteresis_sensitivity, hysteresis_probability_metrics)
table_s9 <- data.frame(
  grade_band_id = table_s9_source$band,
  grade_band_label = publication_band(table_s9_source$band),
  system_id = table_s9_source$system,
  system_label = publication_system(table_s9_source$system),
  analysis_status = "Sensitivity",
  specification_id = table_s9_source$specification,
  baseline_cutoff = table_s9_source$baseline_cutoff,
  followup_cutoff = table_s9_source$followup_cutoff,
  metric = table_s9_source$metric,
  metric_unit = table_s9_source$metric_unit,
  risk_n_mean = table_s9_source$risk_n_mean,
  risk_n_min = table_s9_source$risk_n_min,
  risk_n_max = table_s9_source$risk_n_max,
  estimate_percent = table_s9_source$estimate_percent,
  standard_error_pp = table_s9_source$standard_error_pp,
  ci_lower_percent = table_s9_source$ci_lower_percent,
  ci_upper_percent = table_s9_source$ci_upper_percent,
  estimate_points = table_s9_source$estimate_points,
  standard_error_points = table_s9_source$standard_error_points,
  ci_lower_points = table_s9_source$ci_lower_points,
  ci_upper_points = table_s9_source$ci_upper_points,
  stringsAsFactors = FALSE
)

observed_source <- add_mixed_estimate_fields(observed_ipw, probability_metrics)
observed_probability <- observed_source$metric %in% probability_metrics
table_s10 <- data.frame(
  grade_band_id = observed_source$band,
  grade_band_label = publication_band(observed_source$band),
  system_id = observed_source$system,
  system_label = publication_system(observed_source$system),
  analysis_status = "Sensitivity",
  method_id = observed_source$method,
  target_population = observed_source$target_population,
  metric = observed_source$metric,
  metric_unit = observed_source$metric_unit,
  estimate_percent = observed_source$estimate_percent,
  standard_error_pp = observed_source$standard_error_pp,
  ci_lower_percent = observed_source$ci_lower_percent,
  ci_upper_percent = observed_source$ci_upper_percent,
  estimate_points = observed_source$estimate_points,
  standard_error_points = observed_source$standard_error_points,
  ci_lower_points = observed_source$ci_lower_points,
  ci_upper_points = observed_source$ci_upper_points,
  method_minus_full_estimate_pp = ifelse(
    observed_probability, 100 * observed_source$descriptive_difference_from_full, NA_real_
  ),
  method_minus_full_standard_error_pp = ifelse(
    observed_probability, 100 * observed_source$difference_standard_error, NA_real_
  ),
  method_minus_full_ci_lower_pp = ifelse(
    observed_probability, 100 * observed_source$difference_ci_lower, NA_real_
  ),
  method_minus_full_ci_upper_pp = ifelse(
    observed_probability, 100 * observed_source$difference_ci_upper, NA_real_
  ),
  method_minus_full_estimate_points = ifelse(
    observed_probability, NA_real_, observed_source$descriptive_difference_from_full
  ),
  method_minus_full_standard_error_points = ifelse(
    observed_probability, NA_real_, observed_source$difference_standard_error
  ),
  method_minus_full_ci_lower_points = ifelse(
    observed_probability, NA_real_, observed_source$difference_ci_lower
  ),
  method_minus_full_ci_upper_points = ifelse(
    observed_probability, NA_real_, observed_source$difference_ci_upper
  ),
  stringsAsFactors = FALSE
)

mnar_source <- add_mixed_estimate_fields(mnar_sensitivity, probability_metrics)
table_s11 <- data.frame(
  grade_band_id = mnar_source$band,
  grade_band_label = publication_band(mnar_source$band),
  system_id = mnar_source$system,
  system_label = publication_system(mnar_source$system),
  analysis_status = "Sensitivity",
  scenario_id = mnar_source$scenario,
  affected_pattern = mnar_source$affected_pattern,
  delta_points_per_affected_domain = mnar_source$delta_points_per_affected_domain,
  metric = mnar_source$metric,
  metric_unit = mnar_source$metric_unit,
  estimate_percent = mnar_source$estimate_percent,
  standard_error_pp = mnar_source$standard_error_pp,
  ci_lower_percent = mnar_source$ci_lower_percent,
  ci_upper_percent = mnar_source$ci_upper_percent,
  estimate_points = mnar_source$estimate_points,
  standard_error_points = mnar_source$standard_error_points,
  ci_lower_points = mnar_source$ci_lower_points,
  ci_upper_points = mnar_source$ci_upper_points,
  stringsAsFactors = FALSE
)

table_s12 <- data.frame(
  grade_band_id = mnar_tipping$band,
  grade_band_label = publication_band(mnar_tipping$band),
  system_id = mnar_tipping$system,
  system_label = publication_system(mnar_tipping$system),
  analysis_status = "Sensitivity",
  scenario_id = mnar_tipping$scenario,
  mar_estimate_percent = 100 * mnar_tipping$mar_estimate,
  mar_standard_error_pp = 100 * mnar_tipping$mar_standard_error,
  estimate_at_minus_40_percent = 100 * mnar_tipping$estimate_at_minus_40,
  ratio_at_minus_40_to_mar = mnar_tipping$ratio_at_minus_40_to_mar,
  delta_crossing_half_mar = mnar_tipping$delta_crossing_half_mar,
  half_mar_status = mnar_tipping$half_mar_status,
  delta_crossing_mar_lower_95 = mnar_tipping$delta_crossing_mar_lower_95,
  mar_lower_95_status = mnar_tipping$mar_lower_95_status,
  interpretation = mnar_tipping$interpretation,
  stringsAsFactors = FALSE
)

table_s13 <- table_s1[table_s1$grade_band_id == "g89", ]
g89_probability_metrics <- c("below_to_at_or_above", "below_persistence", "baseline_below")
g89_domain_source <- add_mixed_estimate_fields(g89_domains, g89_probability_metrics)
table_s14 <- data.frame(
  grade_band_id = g89_domain_source$band,
  grade_band_label = publication_band(g89_domain_source$band),
  analysis_scope = publication_scope(g89_domain_source$band),
  system_id = g89_domain_source$system,
  system_label = publication_system(g89_domain_source$system),
  analysis_status = "Sensitivity",
  domain = g89_domain_source$domain,
  threshold = g89_domain_source$threshold,
  metric = g89_domain_source$metric,
  metric_unit = g89_domain_source$metric_unit,
  estimate_percent = g89_domain_source$estimate_percent,
  standard_error_pp = g89_domain_source$standard_error_pp,
  ci_lower_percent = g89_domain_source$ci_lower_percent,
  ci_upper_percent = g89_domain_source$ci_upper_percent,
  estimate_points = g89_domain_source$estimate_points,
  standard_error_points = g89_domain_source$standard_error_points,
  ci_lower_points = g89_domain_source$ci_lower_points,
  ci_upper_points = g89_domain_source$ci_upper_points,
  interpretation_caution = g89_domain_source$interpretation_caution,
  stringsAsFactors = FALSE
)

table_s15 <- data.frame(
  grade_band_id = measurement_crossing$band,
  grade_band_label = publication_band(measurement_crossing$band),
  analysis_scope = publication_scope(measurement_crossing$band),
  system_id = measurement_crossing$system,
  system_label = publication_system(measurement_crossing$system),
  analysis_status = ifelse(
    measurement_crossing$comparison == "actual_longitudinal_matched_pv",
    "Primary", "Exploratory"
  ),
  comparison_id = measurement_crossing$comparison,
  estimate_percent = 100 * measurement_crossing$estimate,
  standard_error_pp = 100 * measurement_crossing$standard_error,
  pair_standard_deviation_pp = 100 * measurement_crossing$pair_standard_deviation,
  pair_minimum_percent = 100 * measurement_crossing$pair_minimum,
  pair_maximum_percent = 100 * measurement_crossing$pair_maximum,
  variance_scope = measurement_crossing$variance_scope,
  interpretation = ifelse(
    measurement_crossing$comparison == "actual_longitudinal_matched_pv",
    "Observed longitudinal BB-to-AA transition to joint attainment",
    measurement_crossing$interpretation
  ),
  ratio_to_actual_longitudinal = measurement_crossing$ratio_to_actual_longitudinal,
  stringsAsFactors = FALSE
)

primary_raw_difference <- max(vapply(config$g45_systems, function(system) {
  abs(standardization$combined_all$estimate[paste0("raw_full::", system)] - primary$g45[[system]]$estimate["BB_to_AA"])
}, numeric(1)))
primary_raw_se_difference <- max(vapply(config$g45_systems, function(system) {
  abs(standardization$combined_all$se[paste0("raw_full::", system)] - primary$g45[[system]]$se["BB_to_AA"])
}, numeric(1)))
max_mean_growth_point_difference <- max(abs(c(
  official_validation$diff_t0, official_validation$diff_t1, official_validation$diff_change
)), na.rm = TRUE)
max_benchmark_point_difference <- max(abs(c(
  official_validation$diff_benchmark_t0_pp, official_validation$diff_benchmark_t1_pp
)), na.rm = TRUE)
max_mean_growth_se_difference <- max(abs(c(
  official_validation$diff_t0_se, official_validation$diff_t1_se, official_validation$diff_change_se
)), na.rm = TRUE)
max_benchmark_se_difference <- max(abs(c(
  official_validation$diff_benchmark_t0_se_pp, official_validation$diff_benchmark_t1_se_pp
)), na.rm = TRUE)
replicate_draw_completion <- unique(standardization_replicate_diagnostics[, c(
  "draw", "stacked_replicates_expected", "stacked_replicates_attempted",
  "stacked_replicates_completed", "support_rebuilds_expected",
  "support_rebuilds_attempted", "support_rebuilds_completed"
)])
require_gate(
  nrow(replicate_draw_completion) == config$pv_count &&
    !anyDuplicated(replicate_draw_completion$draw),
  "Stacked replicate draw-level completion counts are inconsistent"
)

table_s16 <- data.frame(
  validation_id = c(
    "official_point_estimates", "official_numeric_standard_errors", "official_mean_growth_point_tolerance",
    "official_benchmark_point_tolerance", "official_mean_growth_se_tolerance",
    "official_benchmark_se_tolerance", "missing_pv_cells", "sensitivity_checks",
    "standardization_producer_checks", "standardization_replicate_completion",
    "combined_all_covariance", "combined_all_se_diagonal", "primary_identity"
  ),
  validation = c(
    "Official point estimates reproduced within rounding tolerance",
    "Official numeric SEs reproduced and suppressed benchmark SEs handled",
    "Maximum absolute official mean/growth point discrepancy",
    "Maximum absolute official benchmark point discrepancy",
    "Maximum absolute official mean/growth SE discrepancy",
    "Maximum absolute official benchmark SE discrepancy",
    "Missing longitudinal plausible-value cells",
    "Sensitivity validation checks passed",
    "Standardization producer validation checks passed",
    "Completed fail-fast standardization replicate stages",
    "Combined covariance minimum eigenvalue / tolerance",
    "Maximum combined covariance-diagonal SE discrepancy",
    "Primary estimate and SE identity with transition pipeline"
  ),
  value = c(
    paste0(sum(official_point_pass), " of ", length(official_point_pass)),
    paste0(sum(official_se_pass), " of ", length(official_se_pass), "; ", suppressed_benchmark_ses, " suppressed handled"),
    sprintf("%.4f TIMSS points", max_mean_growth_point_difference),
    sprintf("%.4f percentage points", max_benchmark_point_difference),
    sprintf("%.4f TIMSS points", max_mean_growth_se_difference),
    sprintf("%.4f percentage points", max_benchmark_se_difference),
    as.character(sum(inventory$pv_missing_cells)),
    paste0(sum(sensitivity_validation$passed), " of ", nrow(sensitivity_validation)),
    paste0(sum(standardization_covariance_validation$passed), " of ", nrow(standardization_covariance_validation)),
    paste0(
      sum(replicate_draw_completion$stacked_replicates_completed), " of ",
      sum(replicate_draw_completion$stacked_replicates_expected),
      "; attempted=", sum(replicate_draw_completion$stacked_replicates_attempted)
    ),
    paste0(format(combined_all_minimum_eigenvalue, scientific = TRUE, digits = 7), " vs -1e-8"),
    format(combined_all_se_diagonal_difference, scientific = TRUE, digits = 5),
    paste0(
      "estimate=", format(primary_raw_difference, scientific = TRUE, digits = 4),
      "; SE=", format(primary_raw_se_difference, scientific = TRUE, digits = 4)
    )
  ),
  status = rep("PASS", 13L),
  stringsAsFactors = FALSE
)

# Data-driven interval scales use every displayed normal-approximation CI endpoint.
interval_scale <- function(endpoints, include_zero = FALSE, n = 8L, padding_fraction = 0.04) {
  values <- endpoints[is.finite(endpoints)]
  if (include_zero) values <- c(values, 0)
  if (!length(values)) stop("Cannot construct an interval scale from nonfinite values")
  limits <- range(values)
  span <- diff(limits)
  if (!is.finite(span) || span <= 0) span <- max(abs(limits), 1)
  limits <- limits + c(-1, 1) * span * padding_fraction
  if (all(values > 0) && limits[1] <= 0) limits[1] <- min(values) * (1 - padding_fraction)
  if (all(values < 0) && limits[2] >= 0) limits[2] <- max(values) * (1 - padding_fraction)
  ticks <- pretty(limits, n = n)
  ticks <- ticks[ticks >= limits[1] & ticks <= limits[2]]
  if (length(ticks) < 2L) ticks <- seq(limits[1], limits[2], length.out = max(2L, as.integer(n)))
  if (!all(is.finite(c(limits, ticks))) || any(diff(ticks) <= 0)) {
    stop("Constructed a nonfinite or unordered plotting scale")
  }
  list(limits = limits, ticks = ticks)
}

missingness_plot_data <- bind_rows(lapply(config$g45_systems, function(system) {
  full <- primary$g45[[system]]
  both <- one_row(observed_ipw, observed_ipw$band == "g45" & observed_ipw$system == system &
    observed_ipw$method == "both_wave_observed" & observed_ipw$metric == "BB_to_AA", system)
  ipw <- one_row(observed_ipw, observed_ipw$band == "g45" & observed_ipw$system == system &
    observed_ipw$method == "baseline_participant_response_ipw" & observed_ipw$metric == "BB_to_AA", system)
  mnar <- one_row(mnar_sensitivity, mnar_sensitivity$band == "g45" & mnar_sensitivity$system == system &
    mnar_sensitivity$scenario == "combined_pattern_shift" &
    mnar_sensitivity$delta_points_per_affected_domain == -40 & mnar_sensitivity$metric == "BB_to_AA", system)
  estimate <- 100 * c(full$estimate["BB_to_AA"], both$estimate, ipw$estimate, mnar$estimate)
  standard_error <- 100 * c(full$se["BB_to_AA"], both$standard_error, ipw$standard_error, mnar$standard_error)
  data.frame(
    system = system,
    method = c("Full PV", "Both", "IPW", "MNAR -40"),
    estimate = estimate,
    ci_lower = estimate - 1.96 * standard_error,
    ci_upper = estimate + 1.96 * standard_error,
    stringsAsFactors = FALSE
  )
}))
threshold_plot_data <- threshold_sensitivity[
  threshold_sensitivity$band == "g45" & threshold_sensitivity$metric == "BB_to_AA", , drop = FALSE
]
mnar_plot_data <- mnar_sensitivity[
  mnar_sensitivity$band == "g45" & mnar_sensitivity$scenario == "attriter_followup_shift" &
    mnar_sensitivity$metric == "BB_to_AA", , drop = FALSE
]
g89_plot_data <- data.frame(
  system = config$g89_systems,
  estimate = 100 * vapply(config$g89_systems, function(system) primary$g89[[system]]$estimate["BB_to_AA"], numeric(1)),
  standard_error = 100 * vapply(config$g89_systems, function(system) primary$g89[[system]]$se["BB_to_AA"], numeric(1)),
  stringsAsFactors = FALSE
)
g89_plot_data$ci_lower <- g89_plot_data$estimate - 1.96 * g89_plot_data$standard_error
g89_plot_data$ci_upper <- g89_plot_data$estimate + 1.96 * g89_plot_data$standard_error
crossing_plot_data <- measurement_crossing[measurement_crossing$band == "g45", , drop = FALSE]
crossing_plot_data$estimate_percent <- 100 * crossing_plot_data$estimate
crossing_plot_data$ci_lower_percent <- 100 * (crossing_plot_data$estimate - 1.96 * crossing_plot_data$standard_error)
crossing_plot_data$ci_upper_percent <- 100 * (crossing_plot_data$estimate + 1.96 * crossing_plot_data$standard_error)

plot_scales <- list(
  main_recovery = interval_scale(c(table_2$raw_ci_lower_percent, table_2$raw_ci_upper_percent)),
  raw_standardized = interval_scale(c(
    table_2$raw_ci_lower_percent, table_2$raw_ci_upper_percent,
    table_2$standardized_ci_lower_percent, table_2$standardized_ci_upper_percent
  )),
  threshold = interval_scale(100 * c(threshold_plot_data$ci_lower, threshold_plot_data$ci_upper)),
  missingness = interval_scale(c(missingness_plot_data$ci_lower, missingness_plot_data$ci_upper)),
  standardization_sensitivity = interval_scale(100 * c(
    standardization_alternative_contrasts$ci_lower,
    standardization_alternative_contrasts$ci_upper
  ), TRUE),
  mnar = interval_scale(100 * c(mnar_plot_data$ci_lower, mnar_plot_data$ci_upper)),
  grade89 = interval_scale(c(g89_plot_data$ci_lower, g89_plot_data$ci_upper)),
  crossing = interval_scale(c(crossing_plot_data$ci_lower_percent, crossing_plot_data$ci_upper_percent))
)
plot_endpoint_sets <- list(
  main_recovery = c(table_2$raw_ci_lower_percent, table_2$raw_ci_upper_percent),
  raw_standardized = c(
    table_2$raw_ci_lower_percent, table_2$raw_ci_upper_percent,
    table_2$standardized_ci_lower_percent, table_2$standardized_ci_upper_percent
  ),
  threshold = 100 * c(threshold_plot_data$ci_lower, threshold_plot_data$ci_upper),
  missingness = c(missingness_plot_data$ci_lower, missingness_plot_data$ci_upper),
  standardization_sensitivity = 100 * c(
    standardization_alternative_contrasts$ci_lower,
    standardization_alternative_contrasts$ci_upper
  ),
  mnar = 100 * c(mnar_plot_data$ci_lower, mnar_plot_data$ci_upper),
  grade89 = c(g89_plot_data$ci_lower, g89_plot_data$ci_upper),
  crossing = c(crossing_plot_data$ci_lower_percent, crossing_plot_data$ci_upper_percent)
)
plot_interval_validation <- bind_rows(lapply(names(plot_endpoint_sets), function(plot_id) {
  endpoints <- plot_endpoint_sets[[plot_id]]
  scale <- plot_scales[[plot_id]]
  data.frame(
    plot_id = plot_id,
    endpoint_minimum = min(endpoints),
    endpoint_maximum = max(endpoints),
    limit_lower = scale$limits[1],
    limit_upper = scale$limits[2],
    finite = all(is.finite(c(endpoints, scale$limits, scale$ticks))),
    all_ci_endpoints_visible = min(endpoints) >= scale$limits[1] && max(endpoints) <= scale$limits[2],
    stringsAsFactors = FALSE
  )
}))
require_gate(
  all(plot_interval_validation$finite & plot_interval_validation$all_ci_endpoints_visible),
  "At least one displayed CI endpoint would be clipped"
)

colors <- list(
  blue = "#0072B2", orange = "#D55E00", green = "#009E73",
  purple = "#CC79A7", sky = "#56B4E9", yellow = "#E69F00",
  gray = "#5B6573", light_gray = "#D9DEE5", black = "#222222"
)
publication_order <- publication_system(config$g45_systems)

plot_primary_recovery <- function() {
  old <- par(no.readonly = TRUE)
  on.exit(par(old))
  par(mar = c(4.5, 9.2, 3.2, 1.2), las = 1, family = "sans")
  data <- table_2[match(config$g45_systems, table_2$system_id), ]
  scale <- plot_scales$main_recovery
  y <- rev(seq_len(nrow(data)))
  plot(
    data$raw_recovery_percent, y, type = "n", xlim = scale$limits, ylim = c(0.5, nrow(data) + 0.5),
    yaxt = "n", xaxt = "n", xlab = "BB-to-AA transition to joint attainment (%)", ylab = "", bty = "n",
    main = "BB-to-AA transition to joint attainment, Grades 4 to 5"
  )
  abline(v = scale$ticks, col = "#EEF1F4", lwd = 0.8)
  segments(data$raw_ci_lower_percent, y, data$raw_ci_upper_percent, y, col = colors$blue, lwd = 2)
  points(data$raw_recovery_percent, y, pch = 21, bg = colors$blue, col = "white", cex = 1.25)
  axis(2, at = y, labels = data$system_label, tick = FALSE, las = 1)
  axis(1, at = scale$ticks)
}

plot_raw_standardized <- function() {
  old <- par(no.readonly = TRUE)
  on.exit(par(old))
  par(mar = c(4.7, 9.2, 3.2, 1.2), las = 1, family = "sans")
  data <- table_2[match(config$g45_systems, table_2$system_id), ]
  scale <- plot_scales$raw_standardized
  y <- rev(seq_len(nrow(data)))
  plot(
    data$raw_recovery_percent, y, type = "n", xlim = scale$limits, ylim = c(0.5, nrow(data) + 1.8),
    yaxt = "n", xaxt = "n", xlab = "BB-to-AA transition estimate (%)", ylab = "", bty = "n",
    main = "Observed and common-severity standardized BB-to-AA transition"
  )
  abline(v = scale$ticks, col = "#EEF1F4", lwd = 0.8)
  segments(data$raw_recovery_percent, y, data$standardized_recovery_percent, y, col = colors$light_gray, lwd = 2.5)
  segments(data$raw_ci_lower_percent, y + 0.11, data$raw_ci_upper_percent, y + 0.11, col = colors$black, lwd = 1.5)
  segments(data$standardized_ci_lower_percent, y - 0.11, data$standardized_ci_upper_percent, y - 0.11, col = colors$blue, lwd = 1.5)
  points(data$raw_recovery_percent, y + 0.11, pch = 21, bg = colors$black, col = "white", cex = 1.05)
  points(data$standardized_recovery_percent, y - 0.11, pch = 24, bg = colors$blue, col = "white", cex = 1.12)
  axis(2, at = y, labels = data$system_label, tick = FALSE, las = 1)
  axis(1, at = scale$ticks)
  legend(
    "top", inset = 0.01, horiz = TRUE,
    legend = c("Observed risk-set distribution", "Equal-system common-severity distribution"),
    pch = c(21, 24), pt.bg = c(colors$black, colors$blue), col = c(colors$black, colors$blue),
    bty = "n", cex = 0.68
  )
}

plot_threshold_sensitivity <- function() {
  old <- par(no.readonly = TRUE)
  on.exit(par(old))
  scale <- plot_scales$threshold
  par(mfrow = c(3, 3), mar = c(3.3, 3.8, 2.3, 0.8), oma = c(3.2, 4.5, 2.5, 0.5), family = "sans", las = 1)
  for (system in config$g45_systems) {
    data <- threshold_plot_data[threshold_plot_data$system == system, ]
    data <- data[order(data$threshold), ]
    plot(
      data$threshold, 100 * data$estimate, type = "n", xlim = range(config$threshold_grid), ylim = scale$limits,
      xlab = "", ylab = "", main = publication_system(system), bty = "n", cex.main = 0.92, axes = FALSE
    )
    abline(h = scale$ticks, col = "#EEF1F4", lwd = 0.8)
    axis(1, at = config$threshold_grid, labels = config$threshold_grid, cex.axis = 0.55)
    axis(2, at = scale$ticks, cex.axis = 0.72, las = 1)
    abline(v = config$primary_benchmark, lty = 2, col = colors$orange)
    segments(data$threshold, 100 * data$ci_lower, data$threshold, 100 * data$ci_upper,
      col = grDevices::adjustcolor(colors$blue, alpha.f = 0.45), lwd = 1)
    lines(data$threshold, 100 * data$estimate, col = colors$blue, lwd = 1.8)
    points(data$threshold, 100 * data$estimate, pch = 21, bg = colors$blue, col = "white", cex = 0.75)
  }
  mtext("Benchmark threshold (TIMSS points)", side = 1, outer = TRUE, line = 1.9, cex = 0.88)
  mtext("BB-to-AA transition estimate (%)", side = 2, outer = TRUE, line = 2.5, cex = 0.88, las = 0)
  mtext("Threshold sensitivity of the BB-to-AA transition", side = 3, outer = TRUE, line = 1.0, cex = 1.15, font = 2)
}

plot_missingness_robustness <- function() {
  old <- par(no.readonly = TRUE)
  on.exit(par(old))
  scale <- plot_scales$missingness
  method_ids <- c("Full PV", "Both", "IPW", "MNAR -40")
  method_labels <- c("PV-completed", "Both-wave", "Response-IPW", "MNAR −40")
  par(mfrow = c(3, 3), mar = c(5.2, 3.8, 2.3, 0.8), oma = c(2.8, 4.5, 2.5, 0.5), family = "sans", las = 1)
  for (system in config$g45_systems) {
    data <- missingness_plot_data[missingness_plot_data$system == system, ]
    data <- data[match(method_ids, data$method), ]
    plot(
      seq_len(nrow(data)), data$estimate, type = "n", xlim = c(0.7, 4.3), ylim = scale$limits,
      axes = FALSE, xlab = "", ylab = "", main = publication_system(system), bty = "n", cex.main = 0.92
    )
    abline(h = scale$ticks, col = "#EEF1F4")
    segments(seq_len(nrow(data)), data$ci_lower, seq_len(nrow(data)), data$ci_upper, col = colors$gray)
    lines(seq_len(nrow(data)), data$estimate, col = colors$blue, lwd = 1.5)
    points(seq_len(nrow(data)), data$estimate, pch = 21,
      bg = c(colors$black, colors$sky, colors$green, colors$orange), col = "white", cex = 0.9)
    axis(1, at = 1:4, labels = method_labels, las = 2, cex.axis = 0.58)
    axis(2, at = scale$ticks, cex.axis = 0.75)
  }
  mtext("BB-to-AA transition estimate (%)", side = 2, outer = TRUE, line = 2.5, cex = 0.88, las = 0)
  mtext("Missing-data and response-weighting sensitivity", side = 3, outer = TRUE, line = 1.0, cex = 1.15, font = 2)
}

plot_transition_heatmaps <- function() {
  old <- par(no.readonly = TRUE)
  on.exit(par(old))
  par(mfrow = c(3, 3), mar = c(2.7, 3.2, 2.3, 0.8), oma = c(2.8, 3.0, 2.5, 0.5), family = "sans")
  palette <- grDevices::colorRampPalette(c("#F7FBFF", "#6BAED6", "#08306B"))(101)
  states <- c("BB", "MB", "SB", "AA")
  for (system in config$g45_systems) {
    result <- primary$g45[[system]]
    matrix_values <- matrix(NA_real_, 4, 4, dimnames = list(states, states))
    for (origin in states) for (destination in states) {
      matrix_values[origin, destination] <- result$estimate[paste(origin, "to", destination, sep = "_")]
    }
    plot.new()
    plot.window(xlim = c(0.5, 4.5), ylim = c(0.5, 4.5), asp = 1)
    for (origin_index in seq_along(states)) {
      for (destination_index in seq_along(states)) {
        y <- 5 - origin_index
        value <- matrix_values[origin_index, destination_index]
        color_index <- max(1, min(101, round(value * 100) + 1))
        rect(destination_index - 0.5, y - 0.5, destination_index + 0.5, y + 0.5,
          col = palette[color_index], border = "white")
        text(destination_index, y, sprintf("%.1f", 100 * value),
          col = if (value > 0.5) "white" else colors$black, cex = 0.72, font = 2)
      }
    }
    axis(1, at = 1:4, labels = states, tick = FALSE, cex.axis = 0.8)
    axis(2, at = 4:1, labels = states, tick = FALSE, las = 1, cex.axis = 0.8)
    title(main = publication_system(system), cex.main = 0.9)
    box(col = "#C7CDD4")
  }
  mtext("Follow-up state", side = 1, outer = TRUE, line = 1.5)
  mtext("Baseline state", side = 2, outer = TRUE, line = 1.5)
  mtext("Four-state transition matrices at 475 (%)", side = 3, outer = TRUE, line = 1.0, cex = 1.15, font = 2)
}

plot_overlap_diagnostics <- function() {
  old <- par(no.readonly = TRUE)
  on.exit(par(old))
  par(
    mfrow = c(2, 2), mar = c(5.4, 4.5, 2.8, 0.8), oma = c(0.8, 0.8, 2.5, 0.5),
    family = "sans", las = 1
  )
  short_system_labels <- c(
    "Georgia", "Italy", "Jordan", "Korea", "Kosovo", "Montenegro",
    "N. Macedonia", "Slovenia", "Sweden"
  )
  system_factor <- factor(
    publication_system(standardization_overlap_internal$system),
    levels = publication_order,
    labels = short_system_labels
  )
  draw_boxplot <- function(values, color, ylab, title, ylim = NULL, log = "") {
    boxplot(
      values ~ system_factor,
      col = grDevices::adjustcolor(color, 0.55), border = color,
      las = 2, ylab = ylab, xlab = "", main = title, ylim = ylim,
      log = log, outline = FALSE, cex.axis = 0.68, cex.main = 0.92
    )
  }

  retained_values <- 100 * standardization_overlap_internal$support_retained_weight
  retained_scale <- interval_scale(retained_values)
  draw_boxplot(
    retained_values, colors$blue, "Retained risk-set weight (%)",
    "Marginal common-support rectangle", retained_scale$limits
  )

  draw_boxplot(
    standardization_overlap_internal$support_effective_n,
    colors$green, "Kish effective n (log scale)", "Effective sample size", log = "y"
  )

  hull_values <- 100 * standardization_overlap_internal$reference_mass_within_convex_hull
  hull_scale <- interval_scale(hull_values)
  draw_boxplot(
    hull_values, colors$purple, "Reference mass covered (%)",
    "Target-system convex hull", hull_scale$limits
  )
  abline(h = 100, col = colors$gray, lty = 3)

  cell_values <- 100 * standardization_overlap_internal$reference_mass_in_cells_with_n_at_least_5
  cell_scale <- interval_scale(cell_values)
  draw_boxplot(
    cell_values, colors$orange, "Reference mass covered (%)",
    "5 × 5 cells with target n ≥ 5", cell_scale$limits
  )
  abline(h = 100, col = colors$gray, lty = 3)

  mtext(
    "Common-support diagnostics across five plausible-value draws",
    side = 3, outer = TRUE, line = 0.8, cex = 1.08, font = 2
  )
}

plot_standardization_sensitivity <- function() {
  old <- par(no.readonly = TRUE)
  on.exit(par(old))
  scale <- plot_scales$standardization_sensitivity
  alternatives <- setdiff(names(standardization_specifications), "primary")
  short_system_labels <- c("Georgia", "Italy", "Jordan", "Korea", "Kosovo", "Montenegro", "N. Macedonia", "Slovenia", "Sweden")
  par(mfrow = c(3, 2), mar = c(3.4, 7.2, 2.5, 0.8), oma = c(3.0, 1.0, 2.6, 0.5), family = "sans", las = 1)
  for (specification in alternatives) {
    data <- standardization_alternative_contrasts[
      standardization_alternative_contrasts$specification == specification, , drop = FALSE
    ]
    systems <- specification_systems(specification)
    data <- data[match(systems, data$system), ]
    y <- rev(seq_len(nrow(data)))
    plot(
      100 * data$alternative_minus_primary, y, type = "n", xlim = scale$limits,
      ylim = c(0.5, nrow(data) + 0.5), axes = FALSE, xlab = "", ylab = "", bty = "n",
      main = specification_labels[specification], cex.main = 0.9
    )
    abline(v = scale$ticks, col = "#EEF1F4")
    abline(v = 0, col = colors$black, lty = 2, lwd = 1.2)
    segments(100 * data$ci_lower, y, 100 * data$ci_upper, y, col = colors$purple, lwd = 1.5)
    points(100 * data$alternative_minus_primary, y, pch = 21, bg = colors$purple, col = "white", cex = 0.82)
    labels <- short_system_labels[match(systems, config$g45_systems)]
    axis(2, at = y, labels = labels, tick = FALSE, las = 1, cex.axis = 0.66)
    axis(1, at = scale$ticks, cex.axis = 0.7)
  }
  plot.new()
  text(0.5, 0.62, "Five frozen one-factor alternatives", font = 2, cex = 1.0)
  text(0.5, 0.45, "Covariance-aware 95% CIs; dashed line = no difference", cex = 0.82)
  mtext("Alternative minus primary (percentage points)", side = 1, outer = TRUE, line = 1.7, cex = 0.88)
  mtext("Common-severity standardization sensitivity", side = 3, outer = TRUE, line = 1.0, cex = 1.12, font = 2)
}

plot_mnar_curves <- function() {
  old <- par(no.readonly = TRUE)
  on.exit(par(old))
  scale <- plot_scales$mnar
  par(mfrow = c(3, 3), mar = c(3.3, 3.8, 2.3, 0.8), oma = c(3.2, 4.5, 2.5, 0.5), family = "sans", las = 1)
  for (system in config$g45_systems) {
    data <- mnar_plot_data[mnar_plot_data$system == system, ]
    data <- data[order(data$delta_points_per_affected_domain), ]
    plot(
      data$delta_points_per_affected_domain, 100 * data$estimate, type = "n",
      xlim = range(data$delta_points_per_affected_domain), ylim = scale$limits,
      xlab = "", ylab = "", main = publication_system(system), bty = "n", cex.main = 0.92,
      axes = FALSE
    )
    abline(h = scale$ticks, col = "#EEF1F4")
    axis(1, at = sort(unique(data$delta_points_per_affected_domain)))
    axis(2, at = scale$ticks)
    segments(data$delta_points_per_affected_domain, 100 * data$ci_lower,
      data$delta_points_per_affected_domain, 100 * data$ci_upper,
      col = grDevices::adjustcolor(colors$orange, 0.5))
    lines(data$delta_points_per_affected_domain, 100 * data$estimate, col = colors$orange, lwd = 1.8)
    points(data$delta_points_per_affected_domain, 100 * data$estimate,
      pch = 21, bg = colors$orange, col = "white", cex = 0.78)
  }
  mtext("Imposed follow-up PV shift for baseline-only students (points per domain)", side = 1, outer = TRUE, line = 1.9, cex = 0.86)
  mtext("BB-to-AA transition estimate (%)", side = 2, outer = TRUE, line = 2.5, cex = 0.88, las = 0)
  mtext("Pattern-mixture sensitivity for follow-up attrition", side = 3, outer = TRUE, line = 1.0, cex = 1.15, font = 2)
}

plot_grade89_recovery <- function() {
  old <- par(no.readonly = TRUE)
  on.exit(par(old))
  scale <- plot_scales$grade89
  par(mar = c(4.5, 8.0, 3.2, 1.2), las = 1, family = "sans")
  y <- rev(seq_len(nrow(g89_plot_data)))
  plot(g89_plot_data$estimate, y, type = "n", xlim = scale$limits, ylim = c(0.5, nrow(g89_plot_data) + 0.5),
    yaxt = "n", xaxt = "n", xlab = "BB-to-AA transition to joint attainment (%)", ylab = "", bty = "n",
    main = "Grade 8-to-9 extension")
  abline(v = scale$ticks, col = "#EEF1F4")
  segments(g89_plot_data$ci_lower, y, g89_plot_data$ci_upper, y, col = colors$green, lwd = 2)
  points(g89_plot_data$estimate, y, pch = 21, bg = colors$green, col = "white", cex = 1.25)
  axis(2, at = y, labels = publication_system(g89_plot_data$system), tick = FALSE)
  axis(1, at = scale$ticks)
}

plot_measurement_crossing <- function() {
  old <- par(no.readonly = TRUE)
  on.exit(par(old))
  scale <- plot_scales$crossing
  labels <- c("Actual longitudinal", "Baseline discordance", "Follow-up discordance")
  comparisons <- c("actual_longitudinal_matched_pv", "baseline_same_wave_cross_pv", "followup_same_wave_cross_pv")
  base_y <- rev(seq_along(config$g45_systems))
  offsets <- c(0.18, 0, -0.18)
  par(mar = c(4.7, 9.2, 3.2, 1.2), las = 1, family = "sans")
  plot(NA, xlim = scale$limits, ylim = c(0.5, length(config$g45_systems) + 0.7), yaxt = "n", xaxt = "n",
    xlab = "BB-to-AA classification (%)", ylab = "", bty = "n",
    main = "Longitudinal transition versus same-wave classification discordance")
  abline(v = scale$ticks, col = "#EEF1F4")
  for (comparison_index in seq_along(comparisons)) {
    subset <- crossing_plot_data[crossing_plot_data$comparison == comparisons[comparison_index], ]
    subset <- subset[match(config$g45_systems, subset$system), ]
    y <- base_y + offsets[comparison_index]
    segments(subset$ci_lower_percent, y, subset$ci_upper_percent, y,
      col = c(colors$black, colors$purple, colors$orange)[comparison_index], lwd = 1.3)
    points(subset$estimate_percent, y, pch = c(21, 22, 24)[comparison_index],
      bg = c(colors$black, colors$purple, colors$orange)[comparison_index], col = "white", cex = 0.9)
  }
  axis(2, at = base_y, labels = publication_order, tick = FALSE)
  axis(1, at = scale$ticks)
  legend("bottomright", legend = labels, pch = c(21, 22, 24),
    pt.bg = c(colors$black, colors$purple, colors$orange), bty = "n", cex = 0.78)
}

catalog_table_row <- function(table_id, role, number, title, note, caption, declared_rows) {
  data.frame(
    table_id = table_id,
    file = paste0(table_id, ".csv"),
    role = role,
    number = number,
    title = title,
    note = note,
    caption = caption,
    declared_rows = as.integer(declared_rows),
    stringsAsFactors = FALSE
  )
}

catalog_figure_row <- function(figure_id, role, number, title, note, caption, width, height) {
  data.frame(
    figure_id = figure_id,
    pdf_path = file.path("outputs", "figures", paste0(figure_id, ".pdf")),
    png_path = file.path("outputs", "figures", paste0(figure_id, ".png")),
    role = role,
    number = number,
    title = title,
    note = note,
    caption = caption,
    width_inches = width,
    height_inches = height,
    stringsAsFactors = FALSE
  )
}

table_files <- list(
  table_1_sample_participation = table_1,
  table_2_primary_standardized_recovery = table_2,
  table_3_robustness_summary = table_3,
  table_s1_full_transition_matrices = table_s1,
  table_s2_state_distributions = table_s2,
  table_s3_continuous_changes = table_s3,
  table_s4_standardization_specification_decomposition = table_s4,
  table_s5_standardization_alternative_contrasts = table_s5,
  table_s6_standardization_overlap_by_draw = table_s6,
  table_s7_standardization_overlap_summary = table_s7,
  table_s8_threshold_sensitivity = table_s8,
  table_s9_hysteresis_sensitivity = table_s9,
  table_s10_observed_ipw_sensitivity = table_s10,
  table_s11_mnar_pattern_mixture = table_s11,
  table_s12_mnar_tipping_diagnostics = table_s12,
  table_s13_grade89_joint_transitions = table_s13,
  table_s14_grade89_domain_extension = table_s14,
  table_s15_measurement_crossing_diagnostic = table_s15,
  table_s16_validation_summary = table_s16
)

table_catalog <- bind_rows(list(
  catalog_table_row("table_1_sample_participation", "Main", "Table 1", "Sample composition and participation", "Machine IDs and publication labels are retained.", "System and grade-band sample composition, design records, participation patterns, and PV completeness.", 12),
  catalog_table_row("table_2_primary_standardized_recovery", "Main", "Table 2", "Primary raw and common-severity standardized BB-to-AA transition", "Primary 475-point estimand with covariance-aware standardized-minus-raw contrasts.", "Observed and common-severity standardized Grade 4-to-5 BB-to-AA transitions to joint attainment.", 9),
  catalog_table_row("table_3_robustness_summary", "Main", "Table 3", "Robustness and sensitivity summary", "Same-wave unequal-index quantities remain exploratory diagnostics only.", "Selected threshold, buffer, missingness, MNAR, and classification-discordance sensitivity results.", 9),
  catalog_table_row("table_s1_full_transition_matrices", "Supplement", "Table S1", "Full transition matrices", "Both grade bands; percentages and percentage-point SEs.", "Complete 4 by 4 baseline-to-follow-up transition matrices.", 192),
  catalog_table_row("table_s2_state_distributions", "Supplement", "Table S2", "State distributions", "Baseline and follow-up four-state distributions.", "Joint mathematics-science state shares by wave.", 96),
  catalog_table_row("table_s3_continuous_changes", "Supplement", "Table S3", "Continuous score changes", "TIMSS-point changes in the baseline-BB risk set.", "Weakest-domain and domain-specific changes among baseline-BB students.", 36),
  catalog_table_row("table_s4_standardization_specification_decomposition", "Supplement", "Table S4", "Standardization specifications and decomposition", "Primary specification; five frozen sensitivity alternatives.", "Four estimand components and covariance-aware decomposition for all 53 specification-system combinations.", 53),
  catalog_table_row("table_s5_standardization_alternative_contrasts", "Supplement", "Table S5", "Alternative versus primary standardization contrasts", "SEs and intervals are consumed from the joint cross-specification covariance.", "Five frozen alternative-minus-primary sensitivity contrasts by system.", 44),
  catalog_table_row("table_s6_standardization_overlap_by_draw", "Supplement", "Table S6", "Standardization overlap by draw", "Legacy primary overlap calculations are preserved.", "Draw-level marginal support, ESS, convex-hull, and occupied-cell diagnostics.", 45),
  catalog_table_row("table_s7_standardization_overlap_summary", "Supplement", "Table S7", "Standardization overlap summary", "Across-draw summary of the primary standardization support.", "System-level summary of support retention and joint overlap.", 9),
  catalog_table_row("table_s8_threshold_sensitivity", "Supplement", "Table S8", "Threshold sensitivity", "Probability metrics are percentages; score changes are TIMSS points.", "All six planned thresholds and seven metrics in both grade bands.", 504),
  catalog_table_row("table_s9_hysteresis_sensitivity", "Supplement", "Table S9", "Buffered crossing sensitivity", "Both planned baseline-to-follow-up buffers are included.", "Buffered BB-to-AA transitions, joint-deficit persistence, risk shares, and score changes.", 144),
  catalog_table_row("table_s10_observed_ipw_sensitivity", "Supplement", "Table S10", "Observed and response-IPW sensitivity", "Targets remain explicitly distinguished; method-minus-primary intervals use the joint covariance.", "PV-completed, both-wave observed, and response-IPW estimates with covariance-aware method contrasts.", 252),
  catalog_table_row("table_s11_mnar_pattern_mixture", "Supplement", "Table S11", "MNAR pattern-mixture sensitivity", "Imposed score shifts are bounded sensitivity analyses, not corrections.", "Three participation-pattern scenarios across four score shifts.", 1008),
  catalog_table_row("table_s12_mnar_tipping_diagnostics", "Supplement", "Table S12", "MNAR interpolation diagnostics", "Diagnostic interpolation is not a causal tipping-point analysis.", "Bounded scenario-level interpolation diagnostics.", 36),
  catalog_table_row("table_s13_grade89_joint_transitions", "Supplement", "Table S13", "Grade 8-to-9 joint transitions", "Separate extension; never pooled with Grades 4 to 5.", "Complete Grade 8-to-9 joint transition matrices.", 48),
  catalog_table_row("table_s14_grade89_domain_extension", "Supplement", "Table S14", "Grade 8-to-9 domain extension", "Jordan science retains its documented interpretation caution.", "Domain-specific Grade 8-to-9 benchmark crossing and score change.", 24),
  catalog_table_row("table_s15_measurement_crossing_diagnostic", "Supplement", "Table S15", "Exploratory same-wave classification-discordance diagnostic", "Unequal-index intervals are JRR sampling-only conditional on five released PVs and are not corrections.", "Actual longitudinal BB-to-AA transitions and same-wave unequal-index classification-discordance means.", 36),
  catalog_table_row("table_s16_validation_summary", "Supplement", "Table S16", "Validation summary", "Uses the 168-parameter combined_all covariance contract.", "Official reproduction, sensitivity, standardization, covariance, and identity validation.", 13)
))

expected_table_ids <- c(
  "table_1_sample_participation", "table_2_primary_standardized_recovery", "table_3_robustness_summary",
  paste0("table_s", 1:16, c(
    "_full_transition_matrices", "_state_distributions", "_continuous_changes",
    "_standardization_specification_decomposition", "_standardization_alternative_contrasts",
    "_standardization_overlap_by_draw", "_standardization_overlap_summary", "_threshold_sensitivity",
    "_hysteresis_sensitivity", "_observed_ipw_sensitivity", "_mnar_pattern_mixture",
    "_mnar_tipping_diagnostics", "_grade89_joint_transitions", "_grade89_domain_extension",
    "_measurement_crossing_diagnostic", "_validation_summary"
  ))
)
require_gate(identical(names(table_files), expected_table_ids), "Publication table IDs or order are invalid")
require_gate(identical(table_catalog$table_id, expected_table_ids), "Table catalog IDs or order are invalid")
require_gate(!anyDuplicated(table_catalog$table_id) && !anyDuplicated(table_catalog$file), "Table catalog contains duplicate IDs or files")
require_gate(
  identical(vapply(table_files, nrow, integer(1)), setNames(table_catalog$declared_rows, table_catalog$table_id)),
  "Publication table row counts do not match the declared catalog"
)
require_gate(sum(table_catalog$role == "Main") == 3L && sum(table_catalog$role == "Supplement") == 16L, "Table role counts are not 3 main and 16 supplement")
analysis_status_values <- unlist(lapply(table_files, function(data) {
  if ("analysis_status" %in% names(data)) unique(data$analysis_status) else character(0)
}), use.names = FALSE)
require_gate(all(analysis_status_values %in% allowed_analysis_status), "A publication table contains an invalid analysis_status")
analysis_scope_values <- unlist(lapply(table_files, function(data) {
  if ("analysis_scope" %in% names(data)) unique(data$analysis_scope) else character(0)
}), use.names = FALSE)
require_gate(
  length(analysis_scope_values) > 0L && all(analysis_scope_values %in% allowed_analysis_scope),
  "A publication table contains an invalid analysis_scope"
)
scope_consistency <- unlist(lapply(table_files, function(data) {
  if (all(c("grade_band_id", "analysis_scope") %in% names(data))) {
    data$analysis_scope == publication_scope(data$grade_band_id)
  } else {
    logical(0)
  }
}), use.names = FALSE)
require_gate(length(scope_consistency) > 0L && all(scope_consistency), "Grade-band analysis_scope labels are inconsistent")

figure_catalog <- bind_rows(list(
  catalog_figure_row(
    "figure_1_primary_recovery", "Main", "Figure 1",
    "BB-to-AA transition to joint attainment, Grades 4 to 5",
    "Primary observed transition at the 475-point threshold; systems remain in the planned reporting order.",
    "Among students below both mathematics and science thresholds (BB) at baseline, the percentage at or above both (AA) at follow-up. Points are matched plausible-value estimates and bars are jackknife repeated replication 95% confidence intervals.",
    8.2, 5.8
  ),
  catalog_figure_row(
    "figure_2_raw_vs_standardized", "Main", "Figure 2",
    "Observed and common-severity standardized BB-to-AA transition",
    "The synthetic common-severity target gives each of the nine systems equal mass over shared support.",
    "Observed baseline-BB risk-set estimates and modeled estimates integrated over the equal-system common-severity target. Bars are covariance-aware 95% confidence intervals; connecting lines identify each system rather than a temporal change.",
    8.2, 5.8
  ),
  catalog_figure_row(
    "figure_3_threshold_sensitivity", "Main", "Figure 3",
    "Threshold sensitivity of the BB-to-AA transition",
    "Each threshold changes both the baseline-BB denominator and the follow-up-AA event; the dashed vertical line marks the primary threshold of 475.",
    "System-specific Grade 4-to-5 BB-to-AA transition percentages across six planned TIMSS thresholds. Vertical bars are 95% confidence intervals, including endpoints below zero from the untruncated normal approximation where present.",
    9.0, 7.2
  ),
  catalog_figure_row(
    "figure_4_missingness_robustness", "Main", "Figure 4",
    "Missing-data and response-weighting sensitivity",
    "PV-completed, both-wave observed, response-IPW, and combined minus-40-point pattern-mixture analyses have explicitly different targets or assumptions.",
    "Grade 4-to-5 BB-to-AA transition percentages under plausible-value (PV) completion, restriction to both-wave participants, inverse-probability weighting (IPW) for response, and a combined missing-not-at-random (MNAR) minus-40-point scenario. Bars are 95% confidence intervals.",
    9.0, 7.2
  ),
  catalog_figure_row(
    "figure_s1_transition_matrices", "Supplement", "Figure S1",
    "Four-state transition matrices at 475",
    "Grade 4-to-5 rows condition on the baseline state and sum to 100% apart from rounding; darker cells indicate larger percentages.",
    "Matched plausible-value/jackknife repeated replication percentages for all transitions among below-both (BB), mathematics-below-only (MB), science-below-only (SB), and at/above-both (AA) states. Rows are baseline states and columns are follow-up states.",
    9.0, 7.2
  ),
  catalog_figure_row(
    "figure_s2_overlap_diagnostics", "Supplement", "Figure S2",
    "Common-support diagnostics",
    "Boxplots summarize five plausible-value draws; Korea has the weakest effective sample size and occupied-cell coverage.",
    "Four draw-level diagnostics for the primary common-severity standardization: retained baseline-BB weight, Kish effective sample size, equal-system reference mass inside each target system's convex hull, and reference mass in 5 by 5 cells containing at least five target-system cases.",
    10.0, 8.4
  ),
  catalog_figure_row(
    "figure_s3_standardization_sensitivity", "Supplement", "Figure S3",
    "Standardization sensitivity",
    "The Republic of Korea is absent only from the leave-Korea-out panel; the dashed line denotes no alternative-minus-primary difference.",
    "Five internally frozen one-factor standardization alternatives minus the primary equal-system, quadratic-surface estimate. Points are system-specific percentage-point contrasts and bars are covariance-aware 95% confidence intervals.",
    10.0, 8.4
  ),
  catalog_figure_row(
    "figure_s4_mnar_sensitivity", "Supplement", "Figure S4",
    "Pattern-mixture sensitivity for follow-up attrition",
    "Negative shifts are imposed sensitivity scenarios for baseline-only students, not estimated corrections for nonresponse.",
    "Grade 4-to-5 BB-to-AA transition percentages after shifting both follow-up-domain plausible values by 0, minus 10, minus 20, or minus 40 TIMSS points for baseline-only students. Vertical bars are 95% confidence intervals.",
    9.0, 7.2
  ),
  catalog_figure_row(
    "figure_s5_grade89_extension", "Supplement", "Figure S5",
    "Grade 8-to-9 extension",
    "This three-system analysis is a separate descriptive extension, not an independent replication; Jordan science has a documented cross-wave organization change.",
    "Separate Grade 8-to-9 percentages transitioning from below both 475-point thresholds at baseline to at or above both at follow-up. Points are matched plausible-value estimates and bars are jackknife repeated replication 95% confidence intervals.",
    7.2, 4.3
  ),
  catalog_figure_row(
    "figure_s6_measurement_crossing", "Supplement", "Figure S6",
    "Longitudinal BB-to-AA transition versus same-wave classification discordance",
    "Unequal-index means are exploratory; their intervals are JRR sampling-only conditional on the five released plausible values.",
    "Actual matched-plausible-value longitudinal BB-to-AA transitions compared with baseline and follow-up same-wave unequal-index classification-discordance means. Diagnostic bars are JRR sampling-only 95% confidence intervals conditional on five released plausible values.",
    8.5, 5.8
  )
))
expected_figure_ids <- c(
  "figure_1_primary_recovery", "figure_2_raw_vs_standardized", "figure_3_threshold_sensitivity",
  "figure_4_missingness_robustness", "figure_s1_transition_matrices", "figure_s2_overlap_diagnostics",
  "figure_s3_standardization_sensitivity", "figure_s4_mnar_sensitivity",
  "figure_s5_grade89_extension", "figure_s6_measurement_crossing"
)
require_gate(identical(figure_catalog$figure_id, expected_figure_ids), "Figure catalog IDs or order are invalid")
require_gate(!anyDuplicated(figure_catalog$figure_id) && !anyDuplicated(figure_catalog$pdf_path) && !anyDuplicated(figure_catalog$png_path), "Figure catalog contains duplicate IDs or paths")
require_gate(sum(figure_catalog$role == "Main") == 4L && sum(figure_catalog$role == "Supplement") == 6L, "Figure role counts are not 4 main and 6 supplement")
require_gate(
  all(basename(figure_catalog$pdf_path) == paste0(figure_catalog$figure_id, ".pdf")) &&
    all(basename(figure_catalog$png_path) == paste0(figure_catalog$figure_id, ".png")),
  "Figure catalog path pairs are invalid"
)

figure_functions <- list(
  figure_1_primary_recovery = plot_primary_recovery,
  figure_2_raw_vs_standardized = plot_raw_standardized,
  figure_3_threshold_sensitivity = plot_threshold_sensitivity,
  figure_4_missingness_robustness = plot_missingness_robustness,
  figure_s1_transition_matrices = plot_transition_heatmaps,
  figure_s2_overlap_diagnostics = plot_overlap_diagnostics,
  figure_s3_standardization_sensitivity = plot_standardization_sensitivity,
  figure_s4_mnar_sensitivity = plot_mnar_curves,
  figure_s5_grade89_extension = plot_grade89_recovery,
  figure_s6_measurement_crossing = plot_measurement_crossing
)
require_gate(identical(names(figure_functions), expected_figure_ids), "Figure function map is not keyed exactly")
require_gate(
  all(nzchar(figure_catalog$title)) && all(nchar(figure_catalog$caption) >= 140L) &&
    !any(grepl("\\b(confirmatory|prespecified|recovery)\\b", paste(figure_catalog$title, figure_catalog$note, figure_catalog$caption), ignore.case = TRUE)) &&
    !any(grepl("Grade 8-to-9 replication", paste(figure_catalog$title, figure_catalog$note, figure_catalog$caption), ignore.case = TRUE)),
  "Figure catalog captions are incomplete or use obsolete analysis-status wording"
)
require_gate(
  all(nzchar(table_catalog$title)) && all(nzchar(table_catalog$note)) && all(nzchar(table_catalog$caption)) &&
    !any(grepl("\\b(confirmatory|prespecified)\\b", paste(table_catalog$title, table_catalog$note, table_catalog$caption), ignore.case = TRUE)),
  "Table catalog metadata are incomplete or use obsolete analysis-status wording"
)

format_ci_display <- function(estimate, lower, upper, signed_estimate = FALSE) {
  paste0(
    format_number(estimate, 1L, signed = signed_estimate), " (",
    format_number(lower, 1L), "–", format_number(upper, 1L), ")"
  )
}

main_table_1_rows <- cbind(
  gsub(" to ", "–", table_1$grade_band_label, fixed = TRUE),
  table_1$system_label,
  vapply(table_1$students, format_number, character(1), digits = 0L),
  vapply(table_1$schools, format_number, character(1), digits = 0L),
  vapply(table_1$jackknife_zones, format_number, character(1), digits = 0L),
  vapply(table_1$both_wave_students, format_number, character(1), digits = 0L),
  vapply(table_1$baseline_only_students, format_number, character(1), digits = 0L),
  vapply(table_1$followup_only_students, format_number, character(1), digits = 0L),
  vapply(table_1$both_wave_weighted_percent, format_number, character(1), digits = 1L)
)
main_table_1_fragment <- c(
  "**Table 1. Analytic samples and participation patterns**",
  "",
  markdown_table(
    c("Grade band", "System", "Students", "Schools", "JRR zones", "Both-wave *n*", "Baseline-only *n*", "Follow-up-only *n*", "Both-wave weighted %"),
    c("---", "---", "---:", "---:", "---:", "---:", "---:", "---:", "---:"),
    main_table_1_rows
  ),
  "",
  "*Note.* JRR = jackknife repeated replication. A wave is marked participated when either relevant process flag equals 1; source codes were gated to {1, 2}. No students were classified as participating in neither wave in the linked analytic files. Weighted percentages use the final student weight (TOTWGT)."
)

main_table_2_rows <- cbind(
  table_2$system_label,
  vapply(table_2$mean_unweighted_baseline_BB_n, format_number, character(1), digits = 1L),
  vapply(table_2$baseline_BB_weighted_percent, format_number, character(1), digits = 1L),
  mapply(format_ci_display, table_2$raw_recovery_percent, table_2$raw_ci_lower_percent, table_2$raw_ci_upper_percent, USE.NAMES = FALSE),
  mapply(format_ci_display, table_2$standardized_recovery_percent, table_2$standardized_ci_lower_percent, table_2$standardized_ci_upper_percent, USE.NAMES = FALSE),
  mapply(
    format_ci_display,
    table_2$standardized_minus_raw_pp,
    table_2$difference_ci_lower_pp,
    table_2$difference_ci_upper_pp,
    MoreArgs = list(signed_estimate = TRUE), USE.NAMES = FALSE
  ),
  vapply(table_2$minimum_support_effective_n, format_number, character(1), digits = 1L)
)
main_table_2_fragment <- c(
  "**Table 2. Raw and common-severity standardized Grade 4-to-5 BB-to-AA transition at 475**",
  "",
  markdown_table(
    c("System", "Mean BB *n*", "Baseline BB %", "Raw BB-to-AA, % (95% CI)", "Standardized BB-to-AA, % (95% CI)", "Standardized − raw, pp (95% CI)", "Minimum support ESS"),
    c("---", "---:", "---:", "---:", "---:", "---:", "---:"),
    main_table_2_rows
  ),
  "",
  "*Note.* BB = below both mathematics and science thresholds; ESS = Kish effective sample size; pp = percentage points. Mean BB *n* averages the unweighted risk-set count across five plausible-value draws. Confidence intervals use total matched-PV/JRR covariance."
)

main_table_3_columns <- c(
  "primary_475_percent", "standardized_475_percent", "threshold_400_percent",
  "buffered_450_to_500_percent", "both_wave_observed_percent", "response_ipw_percent",
  "combined_mnar_minus40_percent", "followup_same_wave_cross_pv_percent"
)
main_table_3_rows <- cbind(
  table_3$system_label,
  do.call(cbind, lapply(table_3[main_table_3_columns], function(values) {
    vapply(values, format_number, character(1), digits = 1L)
  }))
)
main_table_3_fragment <- c(
  "**Table 3. Selected robustness results for the Grade 4-to-5 BB-to-AA transition (percent)**",
  "",
  markdown_table(
    c("System", "Primary 475", "Std. 475", "Cut 400", "450→500 buffer", "Both wave", "Response IPW", "MNAR −40", "Follow-up discordance"),
    c("---", rep("---:", 8L)),
    main_table_3_rows
  ),
  "",
  "*Note.* Std. 475 is common-severity standardized at 475. Cut 400 and the 450→500 buffer define different denominators and events. Both wave denotes observed both-wave participants; response IPW denotes the inverse-probability-weighted response analysis. MNAR −40 is an imposed missing-not-at-random score-shift scenario. Follow-up discordance is the exploratory follow-up same-wave unequal-index classification-discordance diagnostic; its interval is JRR sampling-only conditional on five released plausible values, not an adjusted longitudinal estimate. Complete standard errors and intervals are in Additional files 1 and 2."
)

grade_band_order <- c("g45", "g89")
supplement_sample_rows <- do.call(rbind, lapply(grade_band_order, function(band) {
  subset <- table_1[table_1$grade_band_id == band, , drop = FALSE]
  c(
    gsub(" to ", "–", subset$grade_band_label[[1L]], fixed = TRUE),
    format_number(nrow(subset), 0L),
    format_number(sum(subset$students), 0L),
    format_number(sum(subset$schools), 0L),
    paste0(
      format_number(min(subset$both_wave_weighted_percent), 1L), "%–",
      format_number(max(subset$both_wave_weighted_percent), 1L), "%"
    )
  )
}))
supplement_sample_fragment <- c(
  markdown_table(
    c("Grade band", "Systems", "Students", "System-specific school records", "Weighted both-wave participation range"),
    c("---", "---:", "---:", "---:", "---:"),
    supplement_sample_rows
  ),
  "",
  "*Note.* Schools are counted within systems. Complete system-level counts, jackknife zones, participation patterns, and plausible-value completeness are in Table 1 and `outputs/diagnostics/data_inventory.csv`."
)

expected_publication_gate_count <- 17L
supplement_validation_rows <- rbind(
  cbind(table_s16$validation, table_s16$value, table_s16$status),
  c("Publication-output gates", paste0(expected_publication_gate_count, " of ", expected_publication_gate_count), "PASS")
)
supplement_validation_fragment <- markdown_table(
  c("Validation", "Value", "Status"),
  c("------------------------------", "---------------:", "-------"),
  supplement_validation_rows
)

supplement_overlap_rows <- cbind(
  table_s7$system_label,
  vapply(table_s7$support_retained_percent_mean, format_number, character(1), digits = 1L),
  vapply(table_s7$support_retained_percent_minimum, format_number, character(1), digits = 1L),
  vapply(table_s7$support_effective_n_minimum, format_number, character(1), digits = 1L),
  vapply(table_s7$reference_mass_within_convex_hull_mean_percent, format_number, character(1), digits = 1L),
  vapply(table_s7$reference_mass_in_cells_n_at_least_5_mean_percent, format_number, character(1), digits = 1L)
)
supplement_overlap_fragment <- markdown_table(
  c("System", "Mean retained weight, %", "Minimum retained weight, %", "Minimum support ESS", "Mean reference mass inside hull, %", "Mean reference mass in cells with at least 5 cases, %"),
  c("---", "---:", "---:", "---:", "---:", "---:"),
  supplement_overlap_rows
)

supplement_grade89_rows <- cbind(
  publication_system(g89_plot_data$system),
  vapply(g89_plot_data$estimate, format_number, character(1), digits = 1L),
  vapply(g89_plot_data$standard_error, format_number, character(1), digits = 1L),
  mapply(
    function(lower, upper) paste0(format_number(lower, 1L), "–", format_number(upper, 1L)),
    g89_plot_data$ci_lower, g89_plot_data$ci_upper, USE.NAMES = FALSE
  )
)
supplement_grade89_fragment <- markdown_table(
  c("System", "BB-to-AA, %", "SE, pp", "95% CI, %"),
  c("---", "---:", "---:", "---:"),
  supplement_grade89_rows
)

supplement_catalog <- table_catalog[table_catalog$role == "Supplement", , drop = FALSE]
supplement_inventory_rows <- cbind(
  sub("Table ", "", supplement_catalog$number, fixed = TRUE),
  vapply(supplement_catalog$declared_rows, format_number, character(1), digits = 0L),
  supplement_catalog$title
)
supplement_inventory_fragment <- markdown_table(
  c("Table", "Rows", "Scope"),
  c("---", "---:", "---"),
  supplement_inventory_rows
)

manuscript_fragments <- setNames(
  list(
    main_table_1_fragment,
    main_table_2_fragment,
    main_table_3_fragment,
    supplement_sample_fragment,
    supplement_validation_fragment,
    supplement_overlap_fragment,
    supplement_grade89_fragment,
    supplement_inventory_fragment
  ),
  c(
    "main_table_1.md", "main_table_2.md", "main_table_3.md",
    "supplement_sample_summary.md", "supplement_validation_summary.md",
    "supplement_overlap_summary.md", "supplement_grade89_summary.md",
    "supplement_table_inventory.md"
  )
)
expected_fragment_basenames <- names(manuscript_fragments)
require_gate(
  length(manuscript_fragments) == 8L && !anyDuplicated(expected_fragment_basenames) &&
    all(vapply(manuscript_fragments, function(lines) length(lines) >= 4L && all(nzchar(lines[nzchar(lines)])), logical(1))),
  "Generated manuscript table-fragment catalog is invalid"
)

# All source gates have passed. Build publication artifacts in a temporary stage.
(function() {
  stage_root <- tempfile(pattern = ".publication-stage-", tmpdir = project_root)
  created <- dir.create(stage_root, recursive = TRUE, showWarnings = FALSE)
  if (!created && !dir.exists(stage_root)) {
    stop("Could not create publication staging directory: ", stage_root)
  }
  on.exit(if (dir.exists(stage_root)) unlink(stage_root, recursive = TRUE, force = TRUE), add = TRUE)

stage_tables <- file.path(stage_root, "tables")
stage_figures <- file.path(stage_root, "figures")
stage_fragments <- file.path(stage_root, "manuscript-generated")
stage_diagnostics <- file.path(stage_root, "diagnostics")
for (directory in c(stage_tables, stage_figures, stage_fragments, stage_diagnostics)) {
  if (!dir.create(directory, recursive = TRUE, showWarnings = FALSE) && !dir.exists(directory)) {
    stop("Could not create staged publication directory: ", directory)
  }
}

for (table_id in expected_table_ids) {
  utils::write.csv(table_files[[table_id]], file.path(stage_tables, paste0(table_id, ".csv")), row.names = FALSE)
}
utils::write.csv(table_catalog, file.path(stage_tables, "table_manifest.csv"), row.names = FALSE)
for (fragment_name in expected_fragment_basenames) {
  write_utf8_lines(manuscript_fragments[[fragment_name]], file.path(stage_fragments, fragment_name))
}

save_figure <- function(directory, stem, width, height, plot_function) {
  require_gate(capabilities("cairo"), "R was built without Cairo support required for publication figures")
  pdf_path <- file.path(directory, paste0(stem, ".pdf"))
  png_path <- file.path(directory, paste0(stem, ".png"))
  render_once <- function(open_device) {
    open_device()
    opened_device <- grDevices::dev.cur()
    tryCatch(
      plot_function(),
      finally = if (grDevices::dev.cur() == opened_device) grDevices::dev.off()
    )
  }
  render_once(function() grDevices::cairo_pdf(pdf_path, width = width, height = height, family = "sans"))
  render_once(function() {
    grDevices::png(
      png_path, width = width, height = height, units = "in", res = 300,
      type = "cairo", bg = "white"
    )
  })
  c(pdf = pdf_path, png = png_path)
}

message("Rendering 10 staged vector and 300-dpi raster figures")
figure_paths <- setNames(vector("list", length(expected_figure_ids)), expected_figure_ids)
for (figure_id in expected_figure_ids) {
  metadata <- figure_catalog[figure_catalog$figure_id == figure_id, ]
  figure_paths[[figure_id]] <- save_figure(
    stage_figures, figure_id, metadata$width_inches, metadata$height_inches,
    figure_functions[[figure_id]]
  )
}
utils::write.csv(figure_catalog, file.path(stage_figures, "figure_manifest.csv"), row.names = FALSE)

png_dimensions <- function(path) {
  connection <- file(path, "rb")
  on.exit(close(connection))
  header <- readBin(connection, what = "raw", n = 24)
  png_signature <- as.raw(c(0x89, 0x50, 0x4e, 0x47, 0x0d, 0x0a, 0x1a, 0x0a))
  ihdr_signature <- as.raw(c(0x49, 0x48, 0x44, 0x52))
  if (
    length(header) < 24 || !identical(header[1:8], png_signature) ||
      !identical(header[13:16], ihdr_signature)
  ) return(c(width = NA_integer_, height = NA_integer_))
  raw_to_integer <- function(bytes) sum(as.integer(bytes) * 256^(rev(seq_along(bytes)) - 1L))
  c(width = raw_to_integer(header[17:20]), height = raw_to_integer(header[21:24]))
}

validate_table_artifacts <- function(directory) {
  bind_rows(lapply(expected_table_ids, function(table_id) {
    path <- file.path(directory, paste0(table_id, ".csv"))
    require_gate(file.exists(path), paste0("Publication table is missing: ", table_id))
    reread <- utils::read.csv(path, check.names = FALSE, stringsAsFactors = FALSE)
    expected <- table_files[[table_id]]
    frames_agree(expected, reread, paste0("Publication table round trip: ", table_id))
    data.frame(
      artifact = table_id,
      expected_rows = nrow(expected),
      observed_rows = nrow(reread),
      expected_columns = ncol(expected),
      observed_columns = ncol(reread),
      file_bytes = file.info(path)$size,
      passed = file.info(path)$size > 20,
      stringsAsFactors = FALSE
    )
  }))
}

validate_figure_artifacts <- function(directory) {
  bind_rows(lapply(expected_figure_ids, function(figure_id) {
    pdf_path <- file.path(directory, paste0(figure_id, ".pdf"))
    png_path <- file.path(directory, paste0(figure_id, ".png"))
    pdf_exists <- file.exists(pdf_path)
    png_exists <- file.exists(png_path)
    png_size <- if (png_exists) file.info(png_path)$size else NA_real_
    pdf_size <- if (pdf_exists) file.info(pdf_path)$size else NA_real_
    dimensions <- if (png_exists) png_dimensions(png_path) else c(width = NA_integer_, height = NA_integer_)
    pdf_signature_valid <- FALSE
    if (pdf_exists) {
      connection <- file(pdf_path, "rb")
      signature <- readBin(connection, what = "raw", n = 5L)
      close(connection)
      pdf_signature_valid <- identical(signature, charToRaw("%PDF-"))
    }
    metadata <- figure_catalog[figure_catalog$figure_id == figure_id, , drop = FALSE]
    expected_width <- as.integer(round(300 * metadata$width_inches))
    expected_height <- as.integer(round(300 * metadata$height_inches))
    png_signature_valid <- all(is.finite(dimensions))
    data.frame(
      artifact = figure_id,
      pdf_exists = pdf_exists,
      png_exists = png_exists,
      pdf_signature_valid = pdf_signature_valid,
      png_signature_valid = png_signature_valid,
      pdf_bytes = pdf_size,
      png_bytes = png_size,
      expected_png_width_pixels = expected_width,
      expected_png_height_pixels = expected_height,
      png_width_pixels = dimensions["width"],
      png_height_pixels = dimensions["height"],
      passed = pdf_exists && png_exists && pdf_signature_valid && png_signature_valid &&
        is.finite(pdf_size) && pdf_size > 5000 && is.finite(png_size) && png_size > 10000 &&
        dimensions["width"] == expected_width && dimensions["height"] == expected_height &&
        dimensions["width"] >= 2000 && dimensions["height"] >= 1200,
      stringsAsFactors = FALSE
    )
  }))
}

validate_fragment_artifacts <- function(directory) {
  bind_rows(lapply(expected_fragment_basenames, function(fragment_name) {
    path <- file.path(directory, fragment_name)
    observed <- if (file.exists(path)) readLines(path, encoding = "UTF-8", warn = FALSE) else character()
    expected <- manuscript_fragments[[fragment_name]]
    size <- if (file.exists(path)) file.info(path)$size else NA_real_
    data.frame(
      artifact = fragment_name,
      file_bytes = size,
      passed = identical(observed, expected) && is.finite(size) && size > 20,
      stringsAsFactors = FALSE
    )
  }))
}

stage_table_checks <- validate_table_artifacts(stage_tables)
stage_figure_checks <- validate_figure_artifacts(stage_figures)
stage_fragment_checks <- validate_fragment_artifacts(stage_fragments)
require_gate(all(stage_table_checks$passed), "A staged publication table failed validation")
require_gate(all(stage_figure_checks$passed), "A staged publication figure failed validation")
require_gate(all(stage_fragment_checks$passed), "A staged manuscript table fragment failed validation")
frames_agree(
  table_catalog,
  utils::read.csv(file.path(stage_tables, "table_manifest.csv"), check.names = FALSE, stringsAsFactors = FALSE),
  "Staged table catalog"
)
frames_agree(
  figure_catalog,
  utils::read.csv(file.path(stage_figures, "figure_manifest.csv"), check.names = FALSE, stringsAsFactors = FALSE),
  "Staged figure catalog"
)

expected_table_basenames <- paste0(expected_table_ids, ".csv")
expected_figure_basenames <- c(paste0(expected_figure_ids, ".pdf"), paste0(expected_figure_ids, ".png"))
staged_table_basenames <- sort(list.files(stage_tables, pattern = "^table_([123]|s[0-9]+[a-z]?)_.*[.]csv$"))
staged_figure_basenames <- sort(list.files(stage_figures, pattern = "^figure_([1-9][0-9]*|s[0-9]+)_.*[.](pdf|png)$"))
staged_fragment_basenames <- sort(list.files(stage_fragments, pattern = "[.]md$"))
require_gate(identical(staged_table_basenames, sort(expected_table_basenames)), "Staged numbered table file set is not exact")
require_gate(identical(staged_figure_basenames, sort(expected_figure_basenames)), "Staged numbered figure file set is not exact")
require_gate(identical(staged_fragment_basenames, sort(expected_fragment_basenames)), "Staged manuscript fragment file set is not exact")

output_validation <- data.frame(
  check = c(
    "official point estimates: 108 of 108",
    "official numeric SEs: 102 of 102 and suppressed values handled",
    "all 13 sensitivity checks passed",
    "exact version-2.1 standardization bundle contract passed",
    "all 11 standardization producer checks passed",
    "all fail-fast standardization replicate stages and diagnostics completed",
    "combined_all covariance and SE diagonal passed",
    "all standardization models converged at full rank",
    "strict common support satisfies frozen n and ESS criteria",
    "table catalog has exact ordered IDs, declared rows, and 3/16 roles",
    "numbered table file set is exact with no stale artifacts",
    "all table artifacts passed schema and row round-trip checks",
    "figure catalog has exact ordered IDs, paths, and 4/6 roles",
    "numbered figure pair set is exact with no stale artifacts",
    "all figure artifacts passed existence, size, and resolution checks",
    "all displayed CI endpoints lie inside finite data-driven limits",
    "all eight manuscript table fragments match canonical analysis objects"
  ),
  passed = c(
    length(official_point_pass) == 108L && all(official_point_pass),
    length(official_se_pass) == 102L && all(official_se_pass) && suppressed_benchmark_ses == 6L,
    nrow(sensitivity_validation) == 13L && all(sensitivity_validation$passed),
    standardization_contract_passed,
    standardization_producer_passed,
    standardization_replicates_complete && standardization_replicate_diagnostics_valid,
    combined_all_asymmetry <= 1e-10 && combined_all_minimum_eigenvalue >= -1e-8 &&
      combined_all_se_diagonal_difference <= 1e-10,
    standardization_models_passed,
    strict_support_passed,
    identical(table_catalog$table_id, expected_table_ids) &&
      identical(table_catalog$declared_rows, unname(vapply(table_files, nrow, integer(1)))) &&
      sum(table_catalog$role == "Main") == 3L && sum(table_catalog$role == "Supplement") == 16L,
    identical(staged_table_basenames, sort(expected_table_basenames)),
    all(stage_table_checks$passed),
    identical(figure_catalog$figure_id, expected_figure_ids) &&
      sum(figure_catalog$role == "Main") == 4L && sum(figure_catalog$role == "Supplement") == 6L,
    identical(staged_figure_basenames, sort(expected_figure_basenames)),
    all(stage_figure_checks$passed),
    all(plot_interval_validation$finite & plot_interval_validation$all_ci_endpoints_visible),
    identical(staged_fragment_basenames, sort(expected_fragment_basenames)) && all(stage_fragment_checks$passed)
  ),
  stringsAsFactors = FALSE
)
require_gate(
  nrow(output_validation) == expected_publication_gate_count &&
    !anyDuplicated(output_validation$check),
  "Publication output validation has an unexpected gate count or duplicate gates"
)
if (!all(output_validation$passed)) {
  stop("Publication output validation failed: ", paste(output_validation$check[!output_validation$passed], collapse = "; "))
}

utils::write.csv(stage_figure_checks, file.path(stage_diagnostics, "figure_artifact_validation.csv"), row.names = FALSE)
utils::write.csv(stage_table_checks, file.path(stage_diagnostics, "table_artifact_validation.csv"), row.names = FALSE)
utils::write.csv(plot_interval_validation, file.path(stage_diagnostics, "figure_interval_validation.csv"), row.names = FALSE)
# This set-level commit marker is deliberately staged and installed last.
utils::write.csv(output_validation, file.path(stage_diagnostics, "publication_output_validation.csv"), row.names = FALSE)

# Superseded numbered outputs and generated fragments participate in the same
# rollback transaction rather than being deleted before the new set commits.
existing_numbered_tables <- list.files(paths$tables, pattern = "^table_([123]|s[0-9]+[a-z]?)_.*[.]csv$", full.names = TRUE)
existing_numbered_figures <- list.files(paths$figures, pattern = "^figure_([1-9][0-9]*|s[0-9]+)_.*[.](pdf|png)$", full.names = TRUE)
generated_manuscript_directory <- file.path(paths$manuscript, "generated")
existing_fragments <- if (dir.exists(generated_manuscript_directory)) {
  list.files(generated_manuscript_directory, pattern = "[.]md$", full.names = TRUE)
} else {
  character()
}
stale_paths <- c(
  existing_numbered_tables[!basename(existing_numbered_tables) %in% expected_table_basenames],
  existing_numbered_figures[!basename(existing_numbered_figures) %in% expected_figure_basenames],
  existing_fragments[!basename(existing_fragments) %in% expected_fragment_basenames]
)

table_staged <- c(
  file.path(stage_tables, expected_table_basenames),
  file.path(stage_tables, "table_manifest.csv")
)
table_destinations <- c(
  file.path(paths$tables, expected_table_basenames),
  file.path(paths$tables, "table_manifest.csv")
)
figure_staged <- c(
  file.path(stage_figures, expected_figure_basenames),
  file.path(stage_figures, "figure_manifest.csv")
)
figure_destinations <- c(
  file.path(paths$figures, expected_figure_basenames),
  file.path(paths$figures, "figure_manifest.csv")
)
fragment_staged <- file.path(stage_fragments, expected_fragment_basenames)
fragment_destinations <- file.path(generated_manuscript_directory, expected_fragment_basenames)
diagnostic_basenames <- c(
  "figure_artifact_validation.csv",
  "table_artifact_validation.csv",
  "figure_interval_validation.csv",
  "publication_output_validation.csv"
)
diagnostic_staged <- file.path(stage_diagnostics, diagnostic_basenames)
diagnostic_destinations <- file.path(paths$diagnostics, diagnostic_basenames)
staged_install_plan <- c(table_staged, figure_staged, fragment_staged, diagnostic_staged)
destination_install_plan <- c(
  table_destinations, figure_destinations, fragment_destinations, diagnostic_destinations
)
require_gate(
  basename(tail(staged_install_plan, 1L)) == "publication_output_validation.csv" &&
    basename(tail(destination_install_plan, 1L)) == "publication_output_validation.csv",
  "Publication validation commit marker is not installed last"
)

final_table_checks <- NULL
final_figure_checks <- NULL
final_fragment_checks <- NULL
publish_staged_set(
  staged_install_plan,
  destination_install_plan,
  stale_paths,
  validate_installed = function() {
    final_numbered_tables <- sort(list.files(paths$tables, pattern = "^table_([123]|s[0-9]+[a-z]?)_.*[.]csv$"))
    final_numbered_figures <- sort(list.files(paths$figures, pattern = "^figure_([1-9][0-9]*|s[0-9]+)_.*[.](pdf|png)$"))
    final_fragments <- sort(list.files(generated_manuscript_directory, pattern = "[.]md$"))
    final_table_checks <<- validate_table_artifacts(paths$tables)
    final_figure_checks <<- validate_figure_artifacts(paths$figures)
    final_fragment_checks <<- validate_fragment_artifacts(generated_manuscript_directory)
    final_table_manifest <- read_required_csv(paths$tables, "table_manifest.csv")
    final_figure_manifest <- read_required_csv(paths$figures, "figure_manifest.csv")
    final_output_validation <- read_required_csv(paths$diagnostics, "publication_output_validation.csv")
    final_figure_validation <- read_required_csv(paths$diagnostics, "figure_artifact_validation.csv")
    final_table_validation <- read_required_csv(paths$diagnostics, "table_artifact_validation.csv")
    final_interval_validation <- read_required_csv(paths$diagnostics, "figure_interval_validation.csv")

    frames_agree(table_catalog, final_table_manifest, "Installed table catalog")
    frames_agree(figure_catalog, final_figure_manifest, "Installed figure catalog")
    frames_agree(stage_table_checks, final_table_checks, "Installed table artifact checks")
    frames_agree(stage_figure_checks, final_figure_checks, "Installed figure artifact checks")
    frames_agree(stage_table_checks, final_table_validation, "Installed table validation record")
    frames_agree(stage_figure_checks, final_figure_validation, "Installed figure validation record")
    frames_agree(plot_interval_validation, final_interval_validation, "Installed interval validation record")
    frames_agree(output_validation, final_output_validation, "Installed publication commit marker")
    require_gate(identical(final_numbered_tables, sort(expected_table_basenames)), "Installed numbered table file set is not exact")
    require_gate(identical(final_numbered_figures, sort(expected_figure_basenames)), "Installed numbered figure file set is not exact")
    require_gate(identical(final_fragments, sort(expected_fragment_basenames)), "Installed manuscript fragment file set is not exact")
    require_gate(all(final_table_checks$passed), "An installed publication table failed round-trip validation")
    require_gate(all(final_figure_checks$passed), "An installed publication figure failed artifact validation")
    require_gate(all(final_fragment_checks$passed), "An installed manuscript fragment differs from its canonical data")
    require_gate(
      all(file.exists(file.path(project_root, final_figure_manifest$pdf_path))) &&
        all(file.exists(file.path(project_root, final_figure_manifest$png_path))),
      "A figure catalog path does not resolve to its declared pair"
    )
    invisible(TRUE)
  }
)

message(
  "Publication outputs complete: ", length(table_files), " tables (3 main, 16 supplement), ",
  length(figure_functions), " figures (4 main, 6 supplement), and all ", nrow(output_validation),
  " publication gates passed"
)
})()
