# Statistics that compare a replication with the published estimate.

# Published estimate on the model's scale: log for ratio measures, as reported otherwise.
published_b <- function(published, ratio) {
  if (ratio) log(c(published$estimate, published$low, published$high)) else c(published$estimate, published$low, published$high)
}

# The published estimate's standard error on the analysis scale, from its 95% interval.
published_se <- function(published_b) (published_b[3] - published_b[2]) / (2 * stats::qnorm(0.975))

# Power of the replication's two-sided test at level alpha to detect the published effect,
# given the replication's standard error and degrees of freedom (noncentral t).
power_to_detect <- function(b_published, se, df, alpha = 0.05) {
  critical <- stats::qt(1 - alpha / 2, df)
  ncp <- b_published / se
  stats::pt(-critical, df, ncp) + 1 - stats::pt(critical, df, ncp)
}

# A distribution-free 95% confidence interval for a median: the order statistics whose ranks
# the binomial(n, 1/2) distribution gives (the exact interval for a median).
median_ci <- function(x, level = 0.95) {
  x <- sort(x)
  n <- length(x)
  lower <- stats::qbinom((1 - level) / 2, n, 0.5)
  upper <- n - lower + 1
  if (lower < 1) return(c(NA, NA))
  c(x[lower], x[upper])
}

# Before 2021-2023 is seen: the power its test would have if its standard error were the
# published one scaled by the square root of the paper's pooled years over the cycle's two, on
# the cycle's design degrees of freedom (30 PSUs in 15 strata). It ignores differences in
# subsample sizes, response, and design effects, so it is a projection only; no criterion uses it.
projected_power <- function(published_b, years, df = 15) {
  power_to_detect(published_b[1], published_se(published_b) * sqrt(years / 2), df)
}