Age, sex and race matched US reference survival — the dashed comparison line on a clinical survival figure — computed in R.
Status: working. us_matched(), us_cohort_curve(),
us_lifetable_vintages() and us_lifetable_model() are implemented and tested
against the fitted blocks.
Production SAS acceptance is complete: all 32 stored %usmatchd answers from
two studies matched to 7.9e-15 worst case, and one study carries a repeatable
acceptance document outside this package. The remaining gap is narrower: the
AVR/LV-function study's 31 table84 answers have not yet been installed as a
repeatable test in that study. See
the design spec.
An R replacement for the Cleveland Clinic SAS macro %usmatchd.
The surprise that shapes the whole package: %usmatchd is not a life-table
lookup. It is PROC HAZPRED evaluating a stored three-phase parametric
hazard fit on the age axis — time origin is birth, not surgery. So this is
not a data package carrying life tables. It ships ~27 small fitted-parameter
blocks and a thin evaluator over
TemporalHazard.
survival::survexp.usr cannot substitute. Twelve vintage × interpolation
combinations were tested; the best individual-curve error in every one was
0.09–0.11, with an age tilt no vintage removes. The reason is structural —
there is no table on the R side of that comparison.
Signatures, not runnable calls — vintage and stratum appear bare because
they have no defaults:
us_matched(age, male, other, times,
id = seq_along(age),
vintage,
table = c("sexrace", "race", "sex", "overall"),
scale = c("years", "months", "days"),
individual = TRUE)
us_cohort_curve(x, by = NULL) # the cohort curve, or one per group
us_lifetable_vintages() # vintages, provenance, race semantics
us_lifetable_model(vintage, stratum) # the raw parameter set
Argument names mirror the macro so a SAS job translates by inspection.
A 62-year-old white male and a 71-year-old white female, ten years of follow-up, against the 1984 fits:
library(hvtiRlifetables)
us_matched(age = c(62, 71),
male = c(1, 0),
other = c(0, 0),
times = seq(0, 10, by = 5),
vintage = "table84")
#> id time agesurv smatched hmatched
#> 1 1 0 0.7925575 1.0000000 0.01826492
#> 2 1 5 0.7925575 0.8914976 0.02843656
#> 3 1 10 0.7925575 0.7453229 0.04439022
#> 4 2 0 0.7773741 1.0000000 0.02241349
#> 5 2 5 0.7773741 0.8640873 0.03736907
#> 6 2 10 0.7773741 0.6758045 0.06334817
agesurv is survival from birth to the patient's current age; smatched is
conditional on having reached it, so it starts at 1. Add
individual = FALSE for the cohort mean curve, which is what a figure's
dashed line usually needs.
us_cohort_curve() reduces the per-patient frame to that dashed line, and with
by gives one line per group of your own choosing:
x <- us_matched(age = c(62, 71, 84, 88), male = c(1, 0, 1, 0),
other = c(0, 0, 0, 0), times = seq(0, 10, by = 5),
vintage = "table84")
band <- factor(c("<80", "<80", "80+", "80+")[x$id], levels = c("<80", "80+"))
us_cohort_curve(x, by = band)
#> group time smatched hmatched
#> 1 <80 0 1.0000000 0.02033920
#> 2 <80 5 0.8777925 0.03283308
#> 3 <80 10 0.7105637 0.05340550
#> 4 80+ 0 1.0000000 0.13316718
#> 5 80+ 5 0.4201996 0.22110436
#> 6 80+ 10 0.1001218 0.36277801
x holds one row per patient per time, so a patient-level grouping is expanded
to it by id — that is what band[x$id] is doing above.
by is not table. table chooses which strata the life table itself is
built from; by chooses how the resulting curves are grouped for reporting.
The two axes are independent, and a report broken down by age band is normally
still matched on "sexrace".
The averaging is the macro's, and it is made in one place so that two studies
cannot average differently without noticing: survival and the density
hmatched * smatched are averaged across patients at each time, unweighted,
and the mean hazard is recovered by dividing. It is not the mean of the
individual hazards.
Omitting vintage is an error, deliberately:
us_matched(70, 1, 0, times = 5)
#> Error: `vintage` has no default and must be given.
#> Available: table84, table2008, table2023.
vintage has no defaultOmitting it is an error.
%usmatchd's own default silently moved twice: from table84 to table2008
between 2010 and 2024, then to table2023 on 2025-12-23. Every job re-run
across either move got different numbers with no signal. An
analysis that does not state its reference vintage is not reproducible, so the
package refuses to guess. This is deliberate friction and it is the package's
main value over calling the macro.
table84 names its non-white strata hzico*, honestly "other".
table2008 and table2023 name the same strata hzicb*, and the macro's own
comment claims these are Black-race estimates.
For 2023 that comment is false. The 2023 "other" category is a
risk-weighted average of Black, Asian, American Indian and Hispanic death
rates, weighted by number at risk, stored under B only to keep the macro's
naming consistent. This package does not propagate the macro's error.
hmatched is per year, alwaysWhatever scale= you pass, time is in those units but hmatched is a
hazard per year.
This matches the macro's source. usmatchd.sas:227 assigns HMATCHED = _HAZARD with no SCALEF multiplication, and PROC HAZPRED is evaluated at
AGE_YR in years. The macro's own header comment at line 47 claims
HMATCHED is "scaled by SCALE" and is wrong — but it has never mattered,
because every job that consumes it uses scale="years", where SCALEF = 1.
This package reproduces the code, not the comment, because bit-fidelity to
the macro is the acceptance criterion. If you need a per-month or per-day
hazard, multiply by the corresponding SCALEF yourself.
hvtiRlifetables is part of the HVTI R family. Install it with the rest of the family through hvtiR:
# install.packages("pak")
pak::pak("ehrlinger/hvtiR")
hvtiR::install()
hvtiR::install() installs every member in one call, so the versions the
members require of one another resolve together. Update later with
hvtiR::update(), from a fresh R session: it will not overwrite a package
that is already loaded.
To install hvtiRlifetables on its own:
# pak::pak("ehrlinger/hvtiRlifetables")
install.packages("hvtiRlifetables", repos = c('https://ehrlinger.r-universe.dev', 'https://cloud.r-project.org'))
The package requires TemporalHazard (>= 1.2.0), and CRAN currently carries
1.1.0. The bound is real, not defensive: the evaluator calls
hzr_decompos_g3(), which 1.1.0 does not provide. Both commands above take
TemporalHazard from GitHub instead: hvtiR::install() installs it with the
family, and the Remotes: field in DESCRIPTION points pak at it when you
install hvtiRlifetables on its own.
vignette("hvtiRlifetables"), takes a synthetic cohort from us_lifetable_vintages() through us_matched() to us_cohort_curve(), with and without a reporting group.GPL-3. The fitted US life-table hazard parameters are the work of The Cleveland Clinic Foundation.