Disease Survival Burden Assessment

Contents

1. US SSA Period Life Table
2. Background Population Survival
3. Disease Survival Model
4. Kaplan–Meier Survival Curves
5. Survival Table
6. Disease Burden Assessment

1. US SSA Period Life Table

The US Social Security Administration publishes annual Period Life Tables giving the conditional probability of dying within one year \(q(x)\) and remaining life expectancy \(e(x)\) at each exact age \(x\), separately for males and females. We use the survexp.us ratetable bundled with the survival package, which contains the SSA Period Life Table reformatted for survival analysis. Ages 50–99 are shown.

# survexp.us: SSA period life table as daily hazard rates (per survival package)
# Dimensions: age (one entry per year, stored in days), sex, calendar year
yr_use <- as.character(max(as.integer(dimnames(survexp.us)[[3]])))
n_ages <- dim(survexp.us)[1]

lt_all <- data.frame(
  age = seq_len(n_ages) - 1L,
  q_m = 1 - exp(-as.numeric(survexp.us[, "male",   yr_use]) * 365.25),
  q_f = 1 - exp(-as.numeric(survexp.us[, "female", yr_use]) * 365.25)
) |>
  mutate(
    l_m = round(100000 * c(1, cumprod(1 - head(q_m, -1)))),
    l_f = round(100000 * c(1, cumprod(1 - head(q_f, -1))))
  )

# Complete life expectancy ê(x) = Σ k_p_x + 0.5 (summed over remaining ages)
lt_ex_fn <- function(q) {
  n <- length(q)
  vapply(seq_len(n), function(i) round(sum(cumprod(1 - q[i:n])) + 0.5, 1), numeric(1L))
}

lt <- lt_all |>
  mutate(e_m = lt_ex_fn(q_m), e_f = lt_ex_fn(q_f)) |>
  filter(age >= 50, age <= 99) |>
  select(age, q_m, l_m, e_m, q_f, l_f, e_f)

kable(
  lt |>
    select(age, q_m, e_m, q_f, e_f) |>
    rename(
      Age           = age,
      `Male q(x)`   = q_m,
      `Male e(x)`   = e_m,
      `Female q(x)` = q_f,
      `Female e(x)` = e_f
    ),
  digits  = 4,
  caption = paste0("US SSA Period Life Table (", yr_use, ") — Ages 50–99")
)

2. Background Population Survival

From the life table we build the survival function \(S_{bg}(t)\) for a person of the specified age and sex. Starting from age \(a\), the probability of surviving at least \(t\) more years is:

\[S_{bg}(t) = \prod_{x=a}^{a+t-1}\bigl(1 - q(x)\bigr), \quad t = 0, 1, 2, \ldots\]

q_col <- if (pt_sex == "M") "q_m" else "q_f"
lt_pt <- lt |> filter(age >= pt_age)
q_x   <- lt_pt[[q_col]]

# Survival at each integer year from diagnosis
S_bg <- c(1, cumprod(1 - q_x))
t_bg <- 0L:(length(S_bg) - 1L)

# Median remaining years (first time S drops below 0.5)
bg_median <- t_bg[which(S_bg < 0.5)[1L]] - 1L

# Restricted mean survival time (RMST) over 30-year horizon
horizon <- min(30L, max(t_bg))
bg_rmst <- sum(S_bg[t_bg < horizon])   # ∫₀ᴴ S(t) dt for unit-step function

For a male aged 50, the SSA life table implies:

• Median remaining survival: 30 years (expected death near age 80)
• 30-year restricted mean survival (RMST): 24.7 years

3. Disease Survival Model

Given only a reported median and 95% CI, we fit a log-normal survival distribution, \(\log T \sim \mathcal{N}(\mu,\,\sigma^2)\).

The log-normal reproduces the stated median exactly. Its two parameters are derived as:

\[\mu = \ln(\text{median}), \qquad \sigma = \frac{\ln(\text{CI}_\text{high}) - \ln(\text{CI}_\text{low})}{2 \times 1.96}\]

The spread \(\sigma\) encodes how wide the survival distribution is — a large \(\sigma\) means high patient-to-patient variability in outcome.

mu_d    <- log(dis_median)
sigma_d <- (log(dis_ci_high) - log(dis_ci_low)) / (2 * 1.96)

# Theoretical mean of this log-normal
dis_mean_raw <- exp(mu_d + sigma_d^2 / 2)

# Disease RMST over the same horizon (fine-grid numerical integration)
t_fine     <- seq(0, horizon, by = 0.01)
S_dis_fine <- 1 - plnorm(t_fine, meanlog = mu_d, sdlog = sigma_d)
dis_rmst   <- sum(diff(t_fine) * S_dis_fine[-length(S_dis_fine)])

Note. \(\sigma\) is derived from the CI on the median estimate, not from patient-level heterogeneity data. It serves as a scale parameter for the shape of the modelled survival curve.

4. Kaplan–Meier Survival Curves

We generate \(n = 10{,}000\) synthetic survival times from each model and apply the standard Kaplan–Meier estimator.

• Background cohort: sampled via the inverse CDF of the SSA life table step function.
• Disease cohort: sampled from the fitted log-normal distribution; times are capped at the life-table horizon.

set.seed(2024)
n_sim <- 10000L

# Background: inverse-CDF of the step survival function
bg_times <- sapply(runif(n_sim), function(u) {
  idx <- which(S_bg <= u)
  if (length(idx) == 0L) max(t_bg) else t_bg[idx[1L]]
})

# Disease: log-normal sample, capped at life-table horizon
dis_times <- pmin(
  rlnorm(n_sim, meanlog = mu_d, sdlog = sigma_d),
  max(t_bg)
)

km_bg  <- survfit(Surv(bg_times,  rep(1L, n_sim)) ~ 1)
km_dis <- survfit(Surv(dis_times, rep(1L, n_sim)) ~ 1)

tidy_km <- function(km, label) {
  data.frame(
    time  = c(0, km$time),
    surv  = c(1, km$surv),
    lower = c(1, km$lower),
    upper = c(1, km$upper),
    group = label
  )
}

bg_label <- paste0(
  "Background (", ifelse(pt_sex == "M", "male", "female"), ", age ", pt_age, ")"
)

bg_df <- tidy_km(km_bg, bg_label)

# Disease CI: parametric bounds from reported median 95% CI.
# With n = 10,000 from a continuous distribution the KM CI is negligible;
# the meaningful uncertainty is in the reported median itself.
t_dis <- c(0, km_dis$time)
dis_df <- data.frame(
  time  = t_dis,
  surv  = c(1, km_dis$surv),
  lower = 1 - plnorm(t_dis, log(dis_ci_high), sigma_d),
  upper = 1 - plnorm(t_dis, log(dis_ci_low),  sigma_d),
  group = dis_name
)

plot_df <- rbind(bg_df, dis_df)

curve_colors <- setNames(c("#2470b3", "#902000"), c(bg_label, dis_name))
curve_types  <- setNames(c("solid",   "dashed"),  c(bg_label, dis_name))

ggplot(plot_df, aes(x = time, y = surv * 100, color = group, linetype = group)) +
  geom_ribbon(aes(ymin = lower * 100, ymax = upper * 100, fill = group),
              alpha = 0.15, color = NA) +
  geom_step(linewidth = 0.85) +
  scale_color_manual(values = curve_colors) +
  scale_fill_manual(values = curve_colors) +
  scale_linetype_manual(values = curve_types) +
  geom_hline(yintercept = 50, color = "gray55", linetype = "dotted", linewidth = 0.4) +
  annotate("text", x = horizon * 0.97, y = 53, label = "50%",
           hjust = 1, size = 3.2, color = "gray50") +
  scale_x_continuous(limits = c(0, horizon), breaks = seq(0, horizon, by = 5)) +
  scale_y_continuous(limits = c(0, 101),     breaks = seq(0, 100,    by = 20)) +
  labs(
    title    = paste0(dis_name, "  —  Survival vs Background Population"),
    subtitle = paste0(
      ifelse(pt_sex == "M", "Male", "Female"), ", age at diagnosis ", pt_age,
      "   |   Disease median ", dis_median,
      " yr  (95% CI ", dis_ci_low, "–", dis_ci_high, " yr)"
    ),
    x       = "Years from diagnosis",
    y       = "Survival (%)",
    color   = NULL, linetype = NULL,
    caption = paste0(
      "Background 95% CI: Greenwood formula (KM, n = ", format(n_sim, big.mark = ","), ").  ",
      "Disease 95% CI: log-normal survival at reported median CI bounds (",
      dis_ci_low, "–", dis_ci_high, " yr)."
    )
  ) +
  theme_gray() +
  theme(legend.position = "bottom") +
  guides(fill = "none")

5. Survival Table

KM estimates at clinically relevant time points. The 95% CI uses Greenwood’s formula (log-log transformation).

t_pts <- c(1, 2, 3, 5, 10, 15, 20)
t_pts <- t_pts[t_pts <= horizon]

extract_km <- function(km, times) {
  s <- summary(km, times = times, extend = TRUE)
  data.frame(
    n_risk = s$n.risk,
    surv   = round(s$surv  * 100, 1),
    lower  = round(s$lower * 100, 1),
    upper  = round(s$upper * 100, 1)
  )
}

bg_tbl  <- extract_km(km_bg,  t_pts)
dis_tbl <- extract_km(km_dis, t_pts)

km_out <- data.frame(
  `Year`                            = t_pts,
  `Bg n at risk`                    = bg_tbl$n_risk,
  `Background survival % (95% CI)`  = paste0(bg_tbl$surv,
                                              " (", bg_tbl$lower, "–", bg_tbl$upper, ")"),
  `Dis n at risk`                   = dis_tbl$n_risk,
  `Disease survival % (95% CI)`     = paste0(dis_tbl$surv,
                                             " (", dis_tbl$lower, "–", dis_tbl$upper, ")"),
  `Difference (pp)`                 = bg_tbl$surv - dis_tbl$surv,
  check.names = FALSE
)

kable(km_out,
      caption = paste0("Kaplan–Meier survival estimates — Background vs ", dis_name),
      align   = c("r", "r", "l", "r", "l", "r"))

6. Disease Burden Assessment

bg_km_tbl     <- summary(km_bg)$table
bg_med_km     <- round(bg_km_tbl[["median"]],  1)
bg_med_lower  <- round(bg_km_tbl[["0.95LCL"]], 1)
bg_med_upper  <- round(bg_km_tbl[["0.95UCL"]], 1)

dis_km_tbl    <- summary(km_dis)$table
dis_med_km    <- round(dis_km_tbl[["median"]],  1)
dis_med_lower <- round(dis_km_tbl[["0.95LCL"]], 1)
dis_med_upper <- round(dis_km_tbl[["0.95UCL"]], 1)

yll_med  <- round(bg_med_km - dis_median, 1)
yll_rmst <- round(bg_rmst   - dis_rmst,   1)

burden_tbl <- data.frame(
  Metric = c(
    paste0("Median remaining survival — background (",
           ifelse(pt_sex == "M", "male", "female"), ", age ", pt_age, ")"),
    paste0("Median overall survival — ", dis_name),
    "Years of life lost  (median comparison)",
    paste0(horizon, "-yr RMST — background"),
    paste0(horizon, "-yr RMST — ", dis_name),
    paste0("Years of life lost  (", horizon, "-yr RMST comparison)")
  ),
  Estimate = c(
    paste0(bg_med_km,  " years  (95% CI ", bg_med_lower, "–", bg_med_upper,
           " yr;  to age ≈", pt_age + as.integer(bg_med_km), ")"),
    paste0(dis_median, " years  (95% CI ", dis_ci_low, "–", dis_ci_high, " yr)"),
    paste0(yll_med,  " years"),
    paste0(round(bg_rmst,  1), " years"),
    paste0(round(dis_rmst, 1), " years"),
    paste0(yll_rmst, " years")
  ),
  stringsAsFactors = FALSE
)

kable(burden_tbl,
      col.names = c("Metric", "Estimate"),
      caption   = "Disease burden summary")

A male aged 50 without the disease has an expected median remaining survival of 31 years (95% CI 30–31 yr; to near age 81), per the US SSA Period Life Table.

With Disease X (median OS 3 yr, 95% CI 1–10 yr), the estimated years of life lost are:

• 28 years comparing medians
• 21.1 years comparing 30-year restricted mean survival times

Interpretation. The disease survival curve represents overall survival as typically reported from clinical trials, incorporating background mortality. The background curve is the expected survival of a matched general-population individual. The gap between the two curves is the attributable disease burden. To update this analysis, change the six parameters in the first code chunk (params) and re-knit.