Add calibration measures

[bblodfon]

• Calibration slope from distribution-matching with KM (idea from Avati 2020, see Fig 5. => CDF vs Fraction of events)
  • Read tutorial/review paper - A tutorial on calibration measurements and calibration models for clinical prediction models
  • Expected Calibration Error (ECE) and MCE (maximum CE), this and this paper have the formulas

[bblodfon]

The KM-calibration statistic (point calibration) from the IPCW paper, section 4.1:

library(survival)
set.seed(123)

# --- 1. Simulate a simple dataset ---
n <- 20
dat <- data.frame(
  time  = rexp(n, 0.1),           # event or censor times
  status = rbinom(n, 1, 0.7),     # event indicator (1=event, 0=censor)
  pred_risk = runif(n)            # predicted risk (0-1)
)

# --- 2. Define binning ---
m <- 2 # number of bins
dat$bin <- cut(dat$pred_risk, breaks = quantile(dat$pred_risk, probs = 0:m/m),
               include.lowest = TRUE, labels = FALSE)

# --- 3. Choose a fixed time horizon tau ---
tau <- 5

# --- 4. Compute p_k^{KM} and its variance ---
results <- lapply(split(dat, dat$bin), function(df_bin) {
  fit <- survfit(Surv(time, status) ~ 1, data = df_bin)
  
  # Extract survival at tau
  # survfit gives survival probabilities; KM estimate of event prob = 1 - S(tau)
  surv_at_tau <- summary(fit, times = tau, extend = TRUE)
  
  S_tau <- surv_at_tau$surv
  var_S_tau <- surv_at_tau$std.err^2  # Greenwood variance of S(tau)
  
  p_k <- 1 - S_tau
  var_p_k <- var_S_tau                # var(1 - S) = var(S)
  
  data.frame(
    n_bin = nrow(df_bin),
    p_k = p_k,
    var_p_k = var_p_k
  )
})

results_df <- do.call(rbind, results)
results_df