Case-cohort hazard ratios and absolute risk

A case-cohort study draws a random subcohort from the full cohort at baseline and augments it with every case (regardless of subcohort membership). The key advantage over nested case-control sampling is that the subcohort serves as a control group simultaneously for all endpoints of interest — even endpoints ascertained years apart — because selection happens once, at the start. Controls are reused across failure times, which is what makes the pseudo-likelihood analysis necessary.

matchatr wraps survival::cch(), which implements the Prentice, Self-Prentice, Lin-Ying, and Borgan I/II pseudo-likelihoods with their correct asymptotic variances. The additional absolute-risk layer (IPW Breslow cumulative baseline hazard) is matchatr’s own.

Code

library(matchatr)
library(survival)

The Wilms tumour data

We use survival::nwtco, the standard case-cohort example from the Handbook of Statistical Methods for Case-Control Studies (Borgan et al. 2018). A random subcohort of 668 subjects was drawn from 4 028 children enrolled in the National Wilms Tumour Study. The endpoint is relapse (rel); the exposure of interest is histological type (histol: 1 = unfavourable). The subcohort indicator is in.subcohort.

Code

nwtco2 <- nwtco
nwtco2$subcohort <- as.logical(nwtco2$in.subcohort)
nwtco2$stage     <- factor(nwtco2$stage)

cat("Cohort size: ", nrow(nwtco2), "\n")
#> Cohort size:  4028
cat("Subcohort:   ", sum(nwtco2$subcohort), " (~",
    round(100 * mean(nwtco2$subcohort)), "%)\n", sep = "")
#> Subcohort:   668 (~17%)
cat("Relapses:    ", sum(nwtco2$rel), "\n")
#> Relapses:     571

matcha() takes the full cohort — it subsets internally to cases + subcohort members before calling cch(), so the denominators for the per-stratum cohort sizes are computed correctly.

Hazard ratio: Prentice pseudo-likelihood

The Prentice (1986) pseudo-likelihood augments each risk set with the case if the case falls outside the subcohort. contrast() reports the hazard ratio with the Self-Prentice corrected asymptotic variance (the naive information-matrix SE is not valid here because controls are reused across failure times).

Code

fit <- matcha(
  nwtco2,
  outcome     = "rel",
  exposure    = "histol",
  design      = case_cohort(subcohort = "subcohort", time = "edrel",
                            method = "Prentice"),
  confounders = ~ stage + age,
  estimator   = "cch"
)
fit
#> <matchatr_fit>
#>  Design:     Case-cohort
#>  Estimator:  cch  (engine: cch)
#>  Outcome:    rel
#>  Exposure:   histol
#>  Confounders: ~stage + age
#>  N:          4028  (cases: 571, controls: 3457)
contrast(fit)       # type = "hr" is the default for a case-cohort design
#> <matchatr_result>
#>  Estimator:  cch  (engine: cch)
#>  Estimand:   hazard ratio
#>  Contrast:   Hazard ratio
#>  CI method:  model
#>  N:          1239
#> 
#> Contrasts:
#>    comparison estimate        se ci_lower ci_upper
#>        <char>    <num>     <num>    <num>    <num>
#> 1:     histol 4.473017 0.7143638  3.27084 6.117047

Unfavourable histology roughly doubles the relapse hazard conditional on stage and age at diagnosis.

Alternative pseudo-likelihoods

Three further methods are supported. "SelfPrentice" uses the subcohort alone at each failure time (slightly more efficient); "LinYing" uses the full subcohort plus all failures (Lin and Ying 1993).

Code

fit_sp <- matcha(nwtco2, "rel", "histol",
                 design = case_cohort("subcohort", "edrel",
                                      method = "SelfPrentice"),
                 confounders = ~ stage + age, estimator = "cch")
fit_ly <- matcha(nwtco2, "rel", "histol",
                 design = case_cohort("subcohort", "edrel",
                                      method = "LinYing"),
                 confounders = ~ stage + age, estimator = "cch")

rbind(
  Prentice    = unlist(contrast(fit)$contrasts[1, c("estimate", "ci_lower", "ci_upper")]),
  SelfPrentice = unlist(contrast(fit_sp)$contrasts[1, c("estimate", "ci_lower", "ci_upper")]),
  LinYing     = unlist(contrast(fit_ly)$contrasts[1, c("estimate", "ci_lower", "ci_upper")])
)
#>              estimate ci_lower ci_upper
#> Prentice     4.473017 3.270840 6.117047
#> SelfPrentice 4.506659 3.295440 6.163054
#> LinYing      4.298614 3.239693 5.703652

All three methods agree closely; differences reflect the distinct variance estimators rather than the point estimates.

Stratified subcohorts: Borgan I/II

When the subcohort is drawn by stratified sampling (e.g. oversampling a rare exposure group), the Borgan et al. (2000) IPW estimators weight each risk-set contribution by the inverse of its stratum-specific subcohort sampling fraction N_s / n_sub_s. Provide the stratification column via stratum =.

Code

fit_b2 <- matcha(
  nwtco2,
  outcome     = "rel",
  exposure    = "histol",
  design      = case_cohort(subcohort = "subcohort", time = "edrel",
                            method = "II.Borgan", stratum = "stage"),
  confounders = ~ age,
  estimator   = "cch"
)
contrast(fit_b2)
#> <matchatr_result>
#>  Estimator:  cch  (engine: cch)
#>  Estimand:   hazard ratio
#>  Contrast:   Hazard ratio
#>  CI method:  model
#>  N:          1154
#> 
#> Contrasts:
#>    comparison estimate        se ci_lower ci_upper
#>        <char>    <num>     <num>    <num>    <num>
#> 1:     histol 4.437016 0.6305204 3.358391 5.862065

Requesting a Borgan method without supplying stratum is declined:

Code

matcha(nwtco2, "rel", "histol",
       design = case_cohort("subcohort", "edrel", method = "I.Borgan"),
       confounders = ~ stage + age, estimator = "cch")
#> Error in `fit_cch()` at matchatr/R/dispatch.R:236:3:
#> ! `method = "I.Borgan"` requires a subcohort stratification column.
#> ℹ Supply `stratum = "<column>"` in `case_cohort()` to specify which column defines the subcohort sampling strata.

Absolute risk: IPW Breslow cumulative baseline hazard

absolute_risk() evaluates the cumulative incidence F_x(t) = 1 − exp(−exp(β̂ᵀ x) Λ̂₀(t)) at specified times, where Λ̂₀(t) is the IPW Breslow estimator. The subcohort-only risk-set denominator at each event time is scaled by the inverse sampling fraction N / n_sub (or per-stratum N_s / n_sub_s for Borgan methods) so the estimate targets the full-cohort baseline hazard. Pointwise confidence intervals use the delta method on the complementary log-log scale (Borgan and Liestøl 1990).

Code

# Two covariate profiles: favourable vs unfavourable histology, median age
nd <- data.frame(
  histol = c(0L, 1L),
  stage  = factor(c(2, 2), levels = 1:4),
  age    = c(3, 3)
)

ar <- absolute_risk(fit, newdata = nd, times = c(500, 1000, 2000))
tidy(ar)
#>      row  time   estimate   ci_lower   ci_upper
#>    <int> <num>      <num>      <num>      <num>
#> 1:     1   500 0.01838153 0.01290246 0.02615629
#> 2:     1  1000 0.02437886 0.01726042 0.03438128
#> 3:     1  2000 0.02600946 0.01844491 0.03661796
#> 4:     2   500 0.07963593 0.05018924 0.12518352
#> 5:     2  1000 0.10452260 0.06657944 0.16213006
#> 6:     2  2000 0.11119772 0.07099865 0.17195364

We can compare the case-cohort absolute risk against the full-cohort Kaplan-Meier / Cox curves. Because the subcohort is only ~17% of the cohort, individual-run sampling variability is visible, but the estimates track the full-cohort curve closely.

Code

cox_full <- coxph(Surv(edrel, rel) ~ histol + stage + age, data = nwtco2)
sf_full  <- survfit(cox_full, newdata = nd)
f_full   <- summary(sf_full, times = c(500, 1000, 2000))$surv
f_full   <- 1 - f_full

cat("Full-cohort F_x(t) [histol=0, histol=1]:\n")
#> Full-cohort F_x(t) [histol=0, histol=1]:
print(round(f_full, 3))
#>          1     2
#> [1,] 0.014 0.067
#> [2,] 0.019 0.089
#> [3,] 0.020 0.095

cat("\nIPW Breslow F_x(t) [histol=0, histol=1]:\n")
#> 
#> IPW Breslow F_x(t) [histol=0, histol=1]:
est_wide <- matrix(ar$estimates$estimate, nrow = 2, byrow = FALSE)
print(round(est_wide, 3))
#>       [,1]  [,2]  [,3]
#> [1,] 0.018 0.026 0.105
#> [2,] 0.024 0.080 0.111

One scale per design

A case-cohort fit identifies the hazard ratio. Requesting an odds ratio (appropriate for matched case-control data) is declined:

Code

contrast(fit, type = "or")
#> Error in `contrast()`:
#> ! A case-cohort design is reported on the hazard-ratio scale.
#> ℹ The Prentice / Self-Prentice / Lin-Ying pseudo-likelihood identifies the hazard ratio. Use `type = "hr"` (the default).

A marginal risk difference / risk ratio from the classical case-cohort analysis would require a source-population prevalence q0 that estimator = "cch" does not supply — but the design-weighted g-computation below gets there a different way.

Marginal causal survival

The hazard ratio above is a conditional effect on the hazard scale. For a marginal causal contrast — the absolute risk, risk difference, risk ratio, or restricted mean survival time (RMST) the whole cohort would show under one exposure versus the other — use estimator = "surv_gcomp". It g-computes on the same design-weighted Cox: it standardizes the model-based absolute risk \(F(t \mid x, W) = 1 - \exp(-\hat\Lambda_0(t)\,e^{\hat\beta^\top x})\) over the subcohort (a random sample of the cohort, weighted by the inverse subcohort sampling fraction) to the treat-all and treat-none marginal risks, then contrasts them. No q0 is needed: a case-cohort samples the cohort directly, so the subcohort already represents the source-population covariate distribution. Variance is a design-preserving bootstrap.

The exposure must be binary, so we code unfavourable histology as a 0/1 indicator:

Code

nwtco2$histol_uh <- as.integer(nwtco2$histol == 2L)

fit_g <- matcha(
  nwtco2,
  outcome     = "rel",
  exposure    = "histol_uh",
  design      = case_cohort(subcohort = "subcohort", time = "edrel"),
  confounders = ~ stage + age,
  estimator   = "surv_gcomp"
)

# Marginal risk difference at two relapse-free follow-up times (days).
contrast(fit_g, type = "difference", times = c(1500, 3000), n_boot = 200)
#> <matchatr_result>
#>  Estimator:  surv_gcomp  (engine: surv_gcomp)
#>  Estimand:   marginal risk difference
#>  Contrast:   Risk difference
#>  CI method:  bootstrap
#>  N:          668
#> 
#> Contrasts:
#>                 comparison  time  estimate         se  ci_lower  ci_upper
#>                     <char> <num>     <num>      <num>     <num>     <num>
#> 1: treat-all vs treat-none  1500 0.2782414 0.04199591 0.2068035 0.3698705
#> 2: treat-all vs treat-none  3000 0.2830524 0.04267281 0.2087195 0.3775254

Unfavourable histology raises the marginal relapse risk by roughly 28 percentage points — a population-level statement the hazard ratio cannot make directly. The restricted mean survival time difference up to a horizon summarises the gap as expected relapse-free days lost:

Code

contrast(fit_g, type = "rmst", times = 3000, n_boot = 200)
#> <matchatr_result>
#>  Estimator:  surv_gcomp  (engine: surv_gcomp)
#>  Estimand:   marginal RMST difference
#>  Contrast:   rmst
#>  CI method:  bootstrap
#>  N:          668
#> 
#> Contrasts:
#>                 comparison  time  estimate       se  ci_lower  ci_upper
#>                     <char> <num>     <num>    <num>     <num>     <num>
#> 1: treat-all vs treat-none  3000 -746.4646 112.9313 -1014.731 -569.4708

The risk ratio (type = "ratio") is available on the same fit.

Covered combinations

Legend. ✅ truth-pinned in tests · ⛔ rejected with an informative error.

Method	Subcohort	Estimand	Variance	Status
Prentice	simple	HR	Self-Prentice asymptotic	✅
SelfPrentice / LinYing	simple	HR	asymptotic	✅
I.Borgan	stratified	HR	Borgan asymptotic	✅
II.Borgan	stratified	HR	plug-in asymptotic	✅
any method	—	absolute risk F_x(t)	IPW Breslow + delta log-log CI	✅
surv_gcomp	simple / stratified	marginal RD(t) / RR(t) / RMST	design-preserving bootstrap	✅
Borgan I/II, no `stratum`	—	—	⛔ `matchatr_bad_design`
cch	—	odds ratio	⛔ use `type = "hr"`
cch	—	RD / RR	⛔ need `q0` (use `surv_gcomp`)
cch	—	sandwich / bootstrap CI	⛔ `matchatr_unsupported_variance`
surv_gcomp	—	conditional OR / HR	⛔ `matchatr_unidentified_estimand`
surv_gcomp	—	non-binary exposure	⛔ `matchatr_bad_input`

See FEATURE_COVERAGE_MATRIX.md for the full combination table.

References

Borgan, Ørnulf, Norman Breslow, Nilanjan Chatterjee, Mitchell H. Gail, Alastair Scott, and Chris J. Wild, eds. 2018. Handbook of Statistical Methods for Case-Control Studies. Chapman & Hall/CRC. https://doi.org/10.1201/9781315154084.

Borgan, Ørnulf, Bryan Langholz, Sven Ove Samuelsen, Larry Goldstein, and Janice Pogoda. 2000. “Exposure Stratified Case-Cohort Designs.” Lifetime Data Analysis 6 (1): 39–58. https://doi.org/10.1023/A:1009661900674.

Borgan, Ørnulf, and Knut Liestøl. 1990. “A Note on Confidence Intervals and Bands for the Survival Function Based on Transformations.” Scandinavian Journal of Statistics 17 (1): 35–41.

Lin, Danyu Y., and Zhiliang Ying. 1993. “Cox Regression with Incomplete Covariate Measurements.” Journal of the American Statistical Association 88 (424): 1341–49. https://doi.org/10.1080/01621459.1993.10476416.

Prentice, Ross L. 1986. “A Case-Cohort Design for Epidemiologic Cohort Studies and Disease Prevention Trials.” Biometrika 73 (1): 1–11. https://doi.org/10.1093/biomet/73.1.1.

Self, Steven G., and Ross L. Prentice. 1988. “Asymptotic Distribution Theory and Efficiency Results for Case-Cohort Studies.” The Annals of Statistics 16 (1): 64–81. https://doi.org/10.1214/aos/1176350691.