# Analysis Clinic Case 003: complete separation in logistic regression
#
# Purpose: create a deterministic binary-outcome data set in which x perfectly
# separates y, show why ordinary maximum-likelihood logistic regression does
# not provide a finite estimand, and verify that Firth's bias-reduced fit gives
# finite penalized-likelihood estimates. This is a teaching example, not a rule
# that Firth regression is the right response for every sparse-data problem.

suppressPackageStartupMessages(library(logistf))

# Twenty non-events all have x = 0; twenty events all have x = 1. The predictor
# therefore separates the outcome classes completely.
dat <- data.frame(
  y = c(rep(0L, 20), rep(1L, 20)),
  x = c(rep(0L, 20), rep(1L, 20))
)

glm_early_warnings <- character()
fit_mle_early <- withCallingHandlers(
  glm(y ~ x, family = binomial(), data = dat, control = glm.control(maxit = 10)),
  warning = function(w) {
    glm_early_warnings <<- c(glm_early_warnings, conditionMessage(w))
    invokeRestart("muffleWarning")
  }
)
fit_mle <- glm(y ~ x, family = binomial(), data = dat)

# logistf uses Firth's penalized likelihood by default and computes penalized
# profile-likelihood confidence intervals. The finite estimate is conditional
# on this method and model; it does not make the original MLE finite.
fit_firth <- logistf(
  y ~ x,
  data = dat,
  plcontrol = logistpl.control(maxit = 1000)
)

outcome_by_x <- table(dat$y, dat$x)
mle_coef <- coef(fit_mle)
mle_se <- sqrt(diag(vcov(fit_mle)))
firth_coef <- coef(fit_firth)
firth_ci <- confint(fit_firth)
firth_probability_at_midpoint <- unname(predict(fit_firth, newdata = data.frame(x = 0.5), type = "response"))

cat("R version:", R.version.string, "\n")
cat("logistf version:", as.character(packageVersion("logistf")), "\n")
cat("total observations:", nrow(dat), "\n")
cat("events:", sum(dat$y == 1), "\n")
cat("non-events:", sum(dat$y == 0), "\n\n")

cat("OUTCOME BY x\n")
print(outcome_by_x)
cat("\n")

cat("ORDINARY MAXIMUM-LIKELIHOOD GLM\n")
cat("10-iteration warning:", paste(glm_early_warnings, collapse = " | "), "\n")
cat("10-iteration converged flag:", fit_mle_early$converged, "\n")
cat("10-iteration slope estimate:", format(unname(coef(fit_mle_early)["x"]), scientific = TRUE), "\n")
cat("default converged flag:", fit_mle$converged, "\n")
cat("default iterations:", fit_mle$iter, "\n")
cat("intercept estimate:", format(unname(mle_coef["(Intercept)"]), scientific = TRUE), "\n")
cat("slope estimate:", format(unname(mle_coef["x"]), scientific = TRUE), "\n")
cat("slope standard error:", format(unname(mle_se["x"]), scientific = TRUE), "\n")
cat("fitted probability range:", paste(format(range(fitted(fit_mle)), scientific = TRUE), collapse = " to "), "\n\n")

cat("FIRTH BIAS-REDUCED LOGISTIC REGRESSION\n")
cat("intercept estimate:", round(unname(firth_coef["(Intercept)"]), 6), "\n")
cat("slope estimate:", round(unname(firth_coef["x"]), 6), "\n")
cat("slope profile-likelihood 95% CI:", paste(round(firth_ci["x", ], 6), collapse = " to "), "\n")
cat("predicted probability at x = 0.5:", round(firth_probability_at_midpoint, 6), "\n\n")

cat("INTERPRETATION\n")
cat("The binary predictor x perfectly separates every event from every non-event.\n")
cat("The ordinary logistic maximum-likelihood slope is not finite; its printed coefficient changes as the algorithm advances and is not a reportable MLE.\n")
cat("The default glm run reports convergence here, showing that its convergence flag alone cannot rule out separation.\n")
cat("Firth's penalized-likelihood method returns finite estimates and a profile-likelihood interval for this controlled example.\n")
cat("An empirical analysis must still verify coding, sparse cells, the prespecified estimand, and whether bias reduction, exact methods, redesign, or qualified reporting fits the study.\n")

stopifnot(
  identical(unname(outcome_by_x["0", "0"]), 20L),
  identical(unname(outcome_by_x["0", "1"]), 0L),
  identical(unname(outcome_by_x["1", "0"]), 0L),
  identical(unname(outcome_by_x["1", "1"]), 20L),
  any(grepl("algorithm did not converge", glm_early_warnings)),
  identical(fit_mle_early$converged, FALSE),
  identical(fit_mle$converged, TRUE),
  unname(mle_coef["x"]) - unname(coef(fit_mle_early)["x"]) > 20,
  abs(unname(mle_coef["x"])) > 50,
  unname(mle_se["x"]) > 10000,
  all(is.finite(firth_coef)),
  all(is.finite(firth_ci)),
  firth_ci["x", 1] > 0,
  abs(firth_probability_at_midpoint - 0.5) < 1e-8
)