# Reproduce from the repository root: Rscript site/evidence/reference.R
# Or download this script and difference.csv into one directory and run there.
# Required packages: survey, jsonlite. No Velocity code is used here.
if (dir.exists("validation/r/library")) .libPaths(c("validation/r/library", .libPaths()))
suppressPackageStartupMessages({ library(survey); library(jsonlite) })
folder <- if (file.exists("difference.csv")) "." else "site/evidence"
source <- read.csv(file.path(folder, "difference.csv"), na.strings = "")
stopifnot(nrow(source) == 16, !anyDuplicated(source$id))
results <- list()
for (mode in c("unweighted", "weighted", "extreme")) {
  for (base in c("answered", "everyone")) {
    for (missing in c("exclude", "category")) {
      d <- source
      d$w <- if (mode == "unweighted") 1 else d$weight
      if (mode == "extreme") d$w[d$id %in% c(5, 9)] <- 32
      valid <- !is.na(d$answer) & (missing == "category" | d$answer != 98)
      # The estimand is conditional on the chosen contributing population.
      # Everyone: non-affirmative and missing answers contribute 0, NOT imputation.
      # Answered: form the analysis design on valid respondents only.
      d$positive <- as.numeric(!is.na(d$answer) & d$answer == 1)
      if (base == "answered") d <- d[valid, ]
      d$group <- factor(d$group, levels = c("A", "B"))
      design <- svydesign(ids = ~1, strata = ~group, weights = ~w, data = d)
      estimates <- svyby(~positive, ~group, design, svymean, covmat = TRUE, keep.var = TRUE)
      contrast <- svycontrast(estimates, c(A = 1, B = -1))
      delta <- unname(coef(contrast)[1])
      se <- unname(SE(contrast)[1])
      df <- degf(design)
      critical <- qt(.975, df)
      groups <- lapply(c("A", "B"), function(g) {
        x <- d[d$group == g, ]
        estimate <- unname(coef(estimates)[g])
        list(group = g, rawN = nrow(x), weightSum = sum(x$w), weightSquares = sum(x$w^2),
             positiveWeight = unname(coef(svytotal(~positive, subset(design, group == g)))[1]),
             proportion = estimate, effectiveN = sum(x$w)^2 / sum(x$w^2))
      })
      kish_se <- sqrt(sum(vapply(groups, function(g) g$proportion * (1-g$proportion) / g$effectiveN, 0.0)))
      results[[paste(mode, base, missing, sep = "/")]] <- list(
        groups = groups,
        survey = list(difference = delta, standardError = se, degreesFreedom = df,
          pValue = 2 * pt(-abs(delta/se), df), normalPValue = 2 * pnorm(-abs(delta/se)),
          interval = c(delta-critical*se, delta+critical*se)),
        kishFormula = list(standardError = kish_se, z = delta/kish_se,
          pValue = 2 * pnorm(-abs(delta/kish_se)))
      )
    }
  }
}
reference <- list(
  schemaVersion = 1,
  runtime = list(R = R.version.string, survey = as.character(packageVersion("survey")),
                 jsonlite = as.character(packageVersion("jsonlite"))),
  design = "Independent respondents within two fixed group strata; sampling weights treated as fixed; no clusters, FPC, replicate weights or calibration variance. Design formed on the chosen contributing population. Conditional valid-case estimates do not adjust nonresponse bias.",
  inference = "svymean + svyby(covmat=TRUE) + svycontrast; Taylor linearization variance; two-sided t reference using degf(design), pointwise 95% t interval. This declared teaching reference is not a claim of Q/Displayr numerical parity.",
  scenarios = results
)
write_json(reference, file.path(folder, "reference.json"), auto_unbox = TRUE, pretty = TRUE, digits = 16)
cat("Wrote", length(results), "survey reference scenarios\n")