Create simulated survival data with age and sex as predictors, where males have younger ages. Demonstrate that sex alone cannot differentiate survival, but it can once adjusted for age.

cox-age-confounding-sex.rmd

---
title: "Confounding/Suppression in Cox Regression: The Role of Age in the Sex–Survival Association"
author: "Your Name"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_float: true
    number_sections: true
    theme: flatly
    code_folding: show
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo    = TRUE,
  warning = FALSE,
  message = FALSE,
  fig.width  = 8,
  fig.height = 6
)
```

# Overview

This document demonstrates a classic **suppression** (a form of confounding)
scenario in survival analysis:

- When **Sex** is the only predictor, its effect is **not significant**.
- When **Sex + Age** are included together, the **Sex effect becomes
  significant**.

The mechanism: if males in the cohort happen to be **younger** on average than
females, the crude Kaplan–Meier comparison mixes two opposing effects:

1. The *protective* effect of younger age among males.
2. The *true harmful* biological effect of male sex on survival.

Because these background differences offset each other in the raw data, the
crude KM curves often lie nearly on top of each other and yield a
non-significant log-rank test. Adjusting for age unmasks the true sex effect.

```{r libraries}
library(survival)
library(survminer)
```

# Simulate the Data

We build age so that **males are younger** (the suppressor), while the true
hazard model makes **both** male sex and older age **harmful**.

```{r simulate}
set.seed(2024)

n <- 600

# Sex: 0 = Female, 1 = Male
sex <- rbinom(n, 1, 0.5)

# Age depends on sex: MALES are YOUNGER on average (the suppressor).
# This is what masks the true sex effect in the crude analysis.
age <- rnorm(n, mean = 65 - 8 * sex, sd = 8)   # males ~8 yrs younger

# True hazard model: both effects HARMFUL
beta_sex <- 0.55   # true log-HR for male sex (harmful)
beta_age <- 0.06   # log-HR per year of age (harmful)

# Linear predictor + exponential baseline hazard
lp            <- beta_sex * sex + beta_age * (age - 65)
baseline_rate <- 0.02
true_time     <- rexp(n, rate = baseline_rate * exp(lp))

# Administrative + random censoring
cens_time <- runif(n, 0, 60)
time      <- pmin(true_time, cens_time)
status    <- as.numeric(true_time <= cens_time)

dat <- data.frame(
  time   = time,
  status = status,
  sex    = factor(sex, levels = c(0, 1), labels = c("Female", "Male")),
  age    = age
)
```

Quick sanity check that males really are younger:

```{r check-age}
aggregate(age ~ sex, data = dat, mean)
```

# Cox Regression: Sex Only

Here the Sex effect is expected to be **non-significant** (masked by age).

```{r cox-sex}
cox_sex <- coxph(Surv(time, status) ~ sex, data = dat)
summary(cox_sex)
```

# Cox Regression: Sex + Age

Adding age holds it constant, so the Sex effect should now **become
significant** (HR for Male ≈ exp(0.55) ≈ 1.7).

```{r cox-both}
cox_both <- coxph(Surv(time, status) ~ sex + age, data = dat)
summary(cox_both)
```

# Crude Kaplan–Meier Curves by Sex

The crude curves should nearly overlap, with a **non-significant log-rank
test**.

```{r km-plot}
km_fit <- survfit(Surv(time, status) ~ sex, data = dat)

ggsurvplot(
  km_fit, data = dat,
  pval         = TRUE,
  risk.table   = TRUE,
  conf.int     = TRUE,
  title        = "Crude Kaplan-Meier by Sex (confounded by age)",
  legend.title = "Sex"
)
```

```{r logrank}
survdiff(Surv(time, status) ~ sex, data = dat)   # non-significant log-rank
```

# Age-Adjusted Survival Curves by Sex

Predicted survival from the two-predictor Cox model, holding **age at its mean**.
These curves should separate clearly, with males below females.

```{r adjusted-plot}
newdat <- data.frame(
  sex = factor(c("Female", "Male"), levels = c("Female", "Male")),
  age = mean(dat$age)          # adjust to overall mean age
)

adj_fit <- survfit(cox_both, newdata = newdat)

ggsurvplot(
  adj_fit, data = dat,
  conf.int     = TRUE,
  legend.labs  = c("Female", "Male"),
  legend.title = "Sex (adjusted, age = mean)",
  title        = "Age-adjusted survival by Sex (Cox model)"
)
```

# Age-Adjusted Curves via `ggadjustedcurves()`

`survminer::ggadjustedcurves()` provides a convenient alternative that computes
adjusted survival curves directly from the fitted Cox model. It supports several
adjustment methods; the two most common are shown below.

## Conditional method (`method = "conditional"`)

Despite the name, this method does **not** plot a single "mean-age" subject.
Instead it uses the whole sample: every subject is assigned to *each* level of
`sex` (keeping their own age), an individual survival curve is predicted for
each, and these are averaged within each group. Because both groups are averaged
over the **same** (full-sample) age distribution, between-group age imbalance is
removed. This standardization to a common covariate distribution is what makes
it the appropriate choice for a confounding-adjusted comparison — and it differs
from the manual `survfit(newdata = ...)` plot above, which instead evaluates one
representative subject at the mean age.

```{r ggadj-conditional}
ggadjustedcurves(
  cox_both,
  data     = dat,
  variable = "sex",
  method   = "conditional",
  palette  = c("#F8766D", "#00BFC4"),   # match ggsurvplot default (Female, Male)
  ggtheme  = theme_bw()
) +
  ggplot2::labs(
    title = "Age-adjusted survival by Sex (conditional method)",
    x     = "Time",
    y     = "Survival probability"
  )
```

## Marginal / average method (`method = "average"`)

The marginal (population-averaged) method also averages individual predicted
survival curves, but each group is averaged over its **own** observed covariate
distribution rather than the full sample's. Because males and females have
different age distributions here, this method does **not** fully remove the age
imbalance — it reflects each group as it actually occurs. Contrast this with the
`"conditional"` method above, which standardizes both groups to the same
(full-sample) age distribution and is therefore the cleaner adjusted comparison.

```{r ggadj-average}
ggadjustedcurves(
  cox_both,
  data     = dat,
  variable = "sex",
  method   = "average",
  palette  = c("#F8766D", "#00BFC4"),   # match ggsurvplot default (Female, Male)
  ggtheme  = theme_bw()
) +
  ggplot2::labs(
    title = "Age-adjusted survival by Sex (marginal/average method)",
    x     = "Time",
    y     = "Survival probability"
  )
```

Both methods should show males with worse (lower) survival than females once age
is accounted for, matching the significant Sex effect from the two-predictor
Cox model.

# Adjusted Curves Using Base `survival` (Standardization)

**Key fact (from `survminer`'s source):** `ggadjustedcurves()` does **not** use
`survfit.coxph()`. Internally it calls `survexp(~ variable, ratetable = fit)`,
which computes *expected* survival from the fitted Cox model. For the
`"conditional"` method it first replicates every subject once per level of the
grouping variable (so each subject appears as both Female and Male), then calls
`survexp` on that expanded data. That replication is what standardizes both
groups to the full-sample age distribution.

So to reproduce `ggadjustedcurves(method = "conditional")` in base `survival`,
use `survexp()` on the replicated data — not `survfit()`.

```{r survexp-standardization}
# Replicate every subject once per sex level (this is what survminer does)
lev   <- levels(dat$sex)
ndata <- dat[rep(seq_len(nrow(dat)), each = length(lev)),
             setdiff(names(dat), "sex")]
ndata$sex <- factor(rep(lev, times = nrow(dat)), levels = lev)

# survexp with ratetable = the Cox fit -> expected (standardized) curves
se_fit <- survexp(~ sex, data = ndata, ratetable = cox_both)

adj_survexp <- do.call(rbind, lapply(seq_along(lev), function(i) {
  data.frame(time = c(0, se_fit$time),
             surv = c(1, se_fit$surv[, i]),
             sex  = lev[i])
}))
adj_survexp$sex <- factor(adj_survexp$sex, levels = lev)
```

```{r survexp-plot}
library(ggplot2)

ggplot(adj_survexp, aes(x = time, y = surv, color = sex)) +
  geom_step(linewidth = 0.9) +
  scale_color_manual(values = c("#F8766D", "#00BFC4")) +   # match ggsurvplot
  coord_cartesian(ylim = c(0, 1)) +
  labs(
    title = "Age-standardized survival by Sex (survexp, matches ggadjustedcurves conditional)",
    x = "Time", y = "Survival probability", color = "Sex"
  ) +
  theme_bw()
```

This curve should overlay `ggadjustedcurves(method = "conditional")`, since it
runs the same computation. (The `survminer` source also applies a small
monotonicity correction that forces the curve to be non-increasing; if you see a
tiny discrepancy at a step, that is the likely cause.)

# Interpretation

- **Sex-only Cox:** HR near 1 with a non-significant p-value. Younger males
  (protected by age) offset the harmful biological male effect.
- **Sex + Age Cox:** recovers a significant male HR (~1.7), because age is now
  held constant.
- **Crude KM curves:** nearly overlap, non-significant log-rank test.
- **Age-adjusted curves:** separate clearly, with males showing worse survival.

# Tuning the Effect

- Increase the age gap (`- 8 * sex` → e.g. `- 12 * sex`) for stronger
  suppression (the crude effect can even flip sign).
- Adjust `beta_sex` / `beta_age` to control how dramatic the reversal is.
- Change `n` for more or less statistical power.

# Session Info

```{r session-info}
sessionInfo()
```