Estimate lifetime risk of developing or dying from cancer using ltRISK

Overview of available approaches

The ltRISK package implements four approaches for estimating cancer risk. These methods differ in their assumptions, data requirements, and treatment of competing mortality.

Which method should I use?

• Use risk_func = "devcan" when incidence counts represent first primary cancers and you want estimates aligned with the DevCan current probability framework.
• Use risk_func = "wun" for historical comparisons with earlier DevCan approaches or when explicitly applying the Wun et al. method.
• Use risk_func = "amp" when incidence data may include multiple primary cancers and only aggregated registry data are available.
• Use risk_func = "cumulative" in calc_ltr(), or use cumulative() with get_risk(), when you need the traditional cumulative risk over a fixed age range and do not want to account for competing mortality.

The package also supports two approaches for confidence intervals in calc_ltr(): a gamma-based method and a delta-method approximation.

For routine analyses, start with calc_ltr() when the age-specific vectors are already available and calc_ltr_df() when the inputs are columns in a data frame or must be analyzed by population, sex, site, period, or another grouping variable. The lower-level devcan(), amp(), wun(), and cumulative() functions expose each point-estimation model directly, while get_risk() extracts specified age ranges from those fitted model objects. This separation keeps the unified confidence-interval workflow in calc_ltr() without hiding the underlying model calculations.

Literature and implementation status

Several related approaches have been proposed for estimating lifetime or age-conditional cancer risk from registry data. The methods differ in their target estimand, treatment of competing mortality, and handling of multiple primary cancers. The table below summarizes how the main methodological papers map to the current ltRISK implementation.

Overall, ltRISK currently implements the main probability estimators needed for routine registry-based lifetime risk analyses: Wun, DevCan/Fay, AMP, and traditional cumulative risk. Methods focused on trend modeling or alternative population-risk estimands can be combined with ltRISK outputs but are not implemented as separate core estimators.

Load the package

library(ltRISK)

Data preparation

The main functions in ltRISK use age-specific aggregated data. The required variables are:

• ages: lower bounds of the age intervals, such as 0, 1, 5, 10, ..., 85.
• cancer: number of incident cancer cases in each age interval. For DevCan and Wun, this usually represents first primary cancers. For AMP, routinely published incidence including multiple primary cancers may be used.
• cancer_death: number of deaths from the cancer or cancer group under study.
• death: number of deaths from all causes in each age interval.
• pys: person-years at risk or the corresponding population denominator for each age interval.

All vectors should have the same length and should refer to the same population, calendar period, sex, cancer site, and age grouping. The ages vector should be sorted in ascending order.

Example datasets

The package includes two example datasets.

• seer_example_data: data used to reproduce examples from DevCan-related publications, including breast cancer data from Fay et al. and all-sites SEER data embedded in DevCan.
• GCO_Today_2022: world-total, both-sexes-combined data assembled from Global Cancer Observatory estimates and World Population Prospects 2022 population and all-cause mortality data, stratified by age and cancer site. Seven site-age cells retain missing source estimates; the all-cancers-excluding-NMSC series is complete.

data("seer_example_data")
breast <- seer_example_data[seer_example_data$site == "Breast", ]
head(seer_example_data)
#>     site    period sex ages cancer death cancer_death     pys
#> 1 Breast 1996-1998   2    0      0  5893            0 4052953
#> 2 Breast 1996-1998   2    5      0   561            0 4032790
#> 3 Breast 1996-1998   2   10      1   628            1 3784789
#> 4 Breast 1996-1998   2   15      9  1367            0 3810986
#> 5 Breast 1996-1998   2   20     43  1547            6 3675646
#> 6 Breast 1996-1998   2   25    335  2064           35 4138795

data("GCO_Today_2022")
head(GCO_Today_2022)
#>   site_code  icd10        site_abbr ages cancer cancer_death   death       pys
#> 1         1 C00-06 Lip, oral cavity    0    329          114 4715187 654355571
#> 2         1 C00-06 Lip, oral cavity    5    319          107  709923 682607098
#> 3         1 C00-06 Lip, oral cavity   10    455          145  457714 664762157
#> 4         1 C00-06 Lip, oral cavity   15    648          261  683239 627243354
#> 5         1 C00-06 Lip, oral cavity   20   1023          683  895411 601601216
#> 6         1 C00-06 Lip, oral cavity   25   3762         1527  981907 589721628

The example datasets included in this package contain first primary cancer counts.

Data-frame and tidyverse-style interface

The core function calc_ltr() takes five age-specific vectors: ages, cancer, cancer_death, death, and pys. This interface is useful when the vectors have already been extracted from a registry table. In many applied workflows, however, cancer registry data are stored as a data frame with one row per age group. For this situation, ltRISK provides calc_ltr_df(), a data-frame interface to the same calculations.

The following call gives the same type of result as calc_ltr(), but the input columns are selected directly from the data frame. Columns can be supplied as bare names, character strings, or character variables.

calc_ltr_df(
  breast,
  ages = ages,
  cancer = cancer,
  cancer_death = cancer_death,
  death = death,
  pys = pys,
  risk_func = "devcan",
  maj_method = "constant",
  ci_method = "gamma",
  age_start = c(0, 30, 50, 70),
  age_end = c(30, 50, 70, Inf),
  digits = 4
)
#>    start end    risk   lower   upper
#> 1      0  30  0.0470  0.0424  0.0519
#> 2      0  50  1.8995  1.8708  1.9286
#> 3      0  70  7.7861  7.7130  7.8598
#> 4      0 Inf 13.3198 13.2170 13.4235
#> 5     30  50  1.8817  1.8529  1.9108
#> 6     30  70  7.8609  7.7868  7.9355
#> 7     30 Inf 13.4816 13.3773 13.5868
#> 8     50  70  6.2505  6.1793  6.3224
#> 9     50 Inf 12.1264 12.0217 12.2320
#> 10    70 Inf  7.3149  7.2202  7.4109

This data-frame interface also works naturally with the base R pipe. For example, the breast cancer data can be filtered and passed directly to calc_ltr_df():

seer_example_data |>
  subset(site == "Breast") |>
  calc_ltr_df(
    ages = ages,
    cancer = cancer,
    cancer_death = cancer_death,
    death = death,
    pys = pys,
    risk_func = "devcan",
    maj_method = "constant",
    ci_method = "delta",
    age_start = c(0, 50, 70),
    age_end = Inf,
    digits = 4
  )
#>   start end    risk   lower   upper
#> 1     0 Inf 13.3198 13.2168 13.4228
#> 2    50 Inf 12.1264 12.0214 12.2313
#> 3    70 Inf  7.3149  7.2199  7.4100

For registry tables containing multiple populations, cancer sites, sexes, or regions, calc_ltr_df() can calculate risks separately within groups. The grouping variables are returned together with the risk estimates, producing a tidy result that can be used in downstream summaries, joins, or plots.

calc_ltr_df(
  seer_example_data,
  by = c("site", "sex"),
  ages = ages,
  cancer = cancer,
  cancer_death = cancer_death,
  death = death,
  pys = pys,
  risk_func = "devcan",
  maj_method = "constant",
  ci_method = "delta",
  age_start = 0,
  age_end = Inf,
  digits = 4
)
#>     site sex start end    risk   lower   upper
#> 1    All   0     0 Inf 39.1670 39.1072 39.2268
#> 2    All   1     0 Inf 39.4881 39.4018 39.5744
#> 3    All   2     0 Inf 38.9734 38.8899 39.0569
#> 4 Breast   2     0 Inf 13.3198 13.2168 13.4228

If dplyr is installed, grouped data frames are also supported. The grouping variables from dplyr::group_by() are detected automatically.

seer_example_data |>
  dplyr::group_by(site, sex) |>
  calc_ltr_df(
    ages = ages,
    cancer = cancer,
    cancer_death = cancer_death,
    death = death,
    pys = pys,
    risk_func = "devcan",
    maj_method = "constant",
    ci_method = "delta",
    age_start = 0,
    age_end = Inf,
    digits = 4
  )
#>     site sex start end    risk   lower   upper
#> 1    All   0     0 Inf 39.1670 39.1072 39.2268
#> 2    All   1     0 Inf 39.4881 39.4018 39.5744
#> 3    All   2     0 Inf 38.9734 38.8899 39.0569
#> 4 Breast   2     0 Inf 13.3198 13.2168 13.4228

PMAJ rate modeling

DevCan version 6.0 and later uses a piecewise mid-age group joinpoint (PMAJ) model (Fay 2004). In ltRISK, PMAJ is an age-specific rate-modeling option rather than a separate risk estimator. It is available for the DevCan, Wun, AMP, and cumulative-risk estimators by setting maj_method = "pmaj". The PMAJ step changes how rates are represented between reported age groups and the grid used for numerical integration; it does not change the definition of the selected risk estimator.

Wun retains maj_method = "constant" as its default so existing analyses continue to use the original Wun intervals. PMAJ must be requested explicitly for Wun. When surgery correction is enabled, its net surgery rate is represented on the same PMAJ grid; confidence intervals continue to use finite differences.

The following DevCan example uses PMAJ for both developing and dying risks.

# Extract all-sites cancer data for both sexes combined
all_sites <- seer_example_data[
  seer_example_data$site == "All" & seer_example_data$sex == 0,
]

# Estimate age-conditional and lifetime risks of developing cancer
calc_ltr(
  ages = all_sites$ages,
  cancer = all_sites$cancer,
  cancer_death = all_sites$cancer_death,
  death = all_sites$death,
  pys = all_sites$pys,
  risk_func = "devcan",
  maj_method = "pmaj",
  age_start = c(0, seq(40, 90, 10)),
  age_end = Inf,
  digits = 4
)
#>   start end    risk   lower   upper
#> 1     0 Inf 38.9045 38.8450 38.9640
#> 2    40 Inf 39.1821 39.1203 39.2440
#> 3    50 Inf 38.3877 38.3241 38.4513
#> 4    60 Inf 36.1041 36.0367 36.1715
#> 5    70 Inf 30.3517 30.2772 30.4264
#> 6    80 Inf 20.9244 20.8408 21.0083
#> 7    90 Inf 11.5426 11.4408 11.6452

# Estimate age-conditional and lifetime risks of dying from cancer
calc_ltr(
  ages = all_sites$ages,
  cancer = all_sites$cancer,
  cancer_death = all_sites$cancer_death,
  death = all_sites$death,
  pys = all_sites$pys,
  risk_func = "devcan",
  maj_method = "pmaj",
  type = "dying",
  age_start = c(0, seq(40, 90, 10)),
  age_end = Inf,
  digits = 4
)
#>   start end    risk   lower   upper
#> 1     0 Inf 17.2576 17.2195 17.2958
#> 2    40 Inf 17.7570 17.7175 17.7965
#> 3    50 Inf 17.8328 17.7925 17.8732
#> 4    60 Inf 17.5314 17.4895 17.5733
#> 5    70 Inf 16.1641 16.1192 16.2092
#> 6    80 Inf 13.2571 13.2067 13.3076
#> 7    90 Inf  9.0968  9.0317  9.1623

Using these example data, the estimated lifetime risk of developing cancer is approximately 38.9045%, and the estimated lifetime risk of dying from cancer is approximately 17.2576%.

The same interface applies PMAJ to the Wun, AMP, and cumulative-risk estimators. Wun and AMP can use an open-ended lifetime target, whereas traditional cumulative risk is normally reported to a finite upper age.

# Wun with PMAJ-smoothed event and mortality rates
calc_ltr(
  ages = all_sites$ages,
  cancer = all_sites$cancer,
  cancer_death = all_sites$cancer_death,
  death = all_sites$death,
  pys = all_sites$pys,
  risk_func = "wun",
  maj_method = "pmaj",
  pmaj_sub_interval = 1,
  ci_method = "none",
  age_start = 0,
  age_end = Inf,
  digits = 4
)
#>   start end    risk
#> 1     0 Inf 38.7774

# AMP with PMAJ-smoothed incidence and removal rates
calc_ltr(
  ages = all_sites$ages,
  cancer = all_sites$cancer,
  cancer_death = all_sites$cancer_death,
  death = all_sites$death,
  pys = all_sites$pys,
  risk_func = "amp",
  maj_method = "pmaj",
  pmaj_sub_interval = 1,
  ci_method = "none",
  age_start = 0,
  age_end = Inf,
  digits = 4
)
#>   start end   risk
#> 1     0 Inf 33.229

# Traditional cumulative risk with a PMAJ-smoothed event rate
calc_ltr(
  ages = all_sites$ages,
  cancer = all_sites$cancer,
  cancer_death = all_sites$cancer_death,
  death = all_sites$death,
  pys = all_sites$pys,
  risk_func = "cumulative",
  maj_method = "pmaj",
  pmaj_sub_interval = 1,
  ci_method = "none",
  age_start = 0,
  age_end = 75,
  digits = 4
)
#>       start end    risk
#> [77,]     0  75 27.3233

Why interpolate age-specific rates?

Registry rates are commonly reported for age bands rather than for every single year of age. A piecewise-constant model assigns one rate to an entire band and changes it abruptly at the next boundary. Those jumps reflect the chosen grouping as well as the underlying age pattern. PMAJ instead linearly interpolates between the rates at adjacent age-group midpoints and represents the interpolated curve on short, piecewise-constant sub-intervals. This gives a more gradual age profile when the underlying incidence or mortality rate is expected to change continuously with age.

The following synthetic example makes the contrast visible. The observed 5-year rates are identical under both approaches; only the way rates are represented between their midpoints changes.

example_ages <- seq(40, 70, by = 5)
example_rates <- c(5, 8, 14, 24, 38, 55, 70) / 100000
age_width <- 5

fine_breaks <- seq(
  min(example_ages),
  max(example_ages) + age_width / 2,
  by = 0.5
)
fine_midages <- head(fine_breaks, -1) + diff(fine_breaks) / 2
fine_rates <- head(
  pmaj(
    ages = example_ages,
    rate_hat = example_rates,
    pmaj_sub_interval = 0.5
  ),
  -1
)

constant_breaks <- c(
  example_ages,
  max(example_ages) + age_width / 2
)
constant_rates <- c(example_rates, tail(example_rates, 1))

plot(
  constant_breaks,
  constant_rates * 100000,
  type = "s",
  lty = 2,
  lwd = 2,
  col = "grey45",
  xlab = "Age (years)",
  ylab = "Rate per 100,000"
)
lines(fine_midages, fine_rates * 100000, lwd = 2, col = "#0072B2")
points(
  example_ages + age_width / 2,
  example_rates * 100000,
  pch = 19
)
legend(
  "topleft",
  legend = c("Piecewise constant", "PMAJ", "Age-group rate"),
  lty = c(2, 1, NA),
  lwd = c(2, 2, NA),
  pch = c(NA, NA, 19),
  col = c("grey45", "#0072B2", "black"),
  bty = "n"
)

Piecewise-constant and PMAJ representations of the same grouped age-specific rates.

This interpolation can reduce artifacts caused by arbitrary age-group boundaries and gives the risk calculation a finer grid for numerical integration. It is especially useful when rates rise or fall rapidly across adjacent age bands or when risks are needed for endpoints inside an original band. PMAJ does not create new observations, however, and it does not by itself remove sampling noise from sparse age-specific rates. Decreasing pmaj_sub_interval improves the numerical approximation to the continuous MAJ calculation; it should not be interpreted as increasing the resolution of the source data. Smaller sub-intervals create more integration steps and therefore increase computation time, particularly when confidence intervals require many risk evaluations.

Because the PMAJ model evaluates age-specific rates on a finer age grid, it can also be used to estimate risks over age ranges whose endpoints are not original age-group boundaries. For example, if the input data are grouped in 5-year age bands, a piecewise-constant calculation is naturally tied to those 5-year boundaries. By contrast, maj_method = "pmaj" with pmaj_sub_interval = 1 creates an approximately 1-year grid, allowing endpoints such as age 78 or age 66. Smaller values of pmaj_sub_interval provide finer age resolution and more closely approximate the continuous mid-age group joinpoint calculation.

# Estimate risks for age ranges such as 0--78 and 40--66 years.
# age_start and age_end are crossed by calc_ltr(), so the rows below include
# all combinations of the supplied starts and ends.
calc_ltr(
  ages = all_sites$ages,
  cancer = all_sites$cancer,
  cancer_death = all_sites$cancer_death,
  death = all_sites$death,
  pys = all_sites$pys,
  risk_func = "devcan",
  maj_method = "pmaj",
  pmaj_sub_interval = 1,
  ci_method = "delta",
  age_start = c(0, 40),
  age_end = c(78, 66),
  digits = 4
)
#>   start end    risk   lower   upper
#> 1     0  66 15.5343 15.5046 15.5641
#> 2     0  78 28.4377 28.3941 28.4812
#> 3    40  66 14.3424 14.3134 14.3714
#> 4    40  78 28.0513 28.0068 28.0957

The same endpoint behavior applies to DevCan, Wun, AMP, and cumulative risk. With maj_method = "constant", finite endpoints generally need to coincide with the original age-group boundaries. With PMAJ, endpoints must lie on the finer grid created by pmaj_sub_interval.

Confidence intervals

The calc_ltr() function reports point estimates together with confidence intervals for lifetime and age-conditional risks. The inference framework follows the approach described by Fay et al. (2003) for age-conditional probabilities of developing cancer. In that work, the risk estimator is treated as a smooth function of the observed age-specific counts, including incident cancer cases, cancer deaths, and deaths from other causes.

Let \(\hat{A}\) denote the estimated risk for a specified age range, and let \(\widehat{V}(\hat{A})\) denote its estimated variance. ltRISK propagates count variability with the delta method. By default, variance_method = "auto" uses analytic gradients for DevCan, AMP, Wun without surgery correction, and cumulative risk. These gradients avoid repeatedly refitting the risk model and are much faster for confidence-interval calculations. variance_method = "finite_difference" remains available as an independent implementation check. Wun calculations with surgery correction currently use finite differences because the API does not specify sampling variances or covariances for the surgery counts H and Ch.

calc_ltr(..., risk_func = "devcan", ci_method = "gamma")
calc_ltr(..., risk_func = "devcan", ci_method = "delta")
calc_ltr(..., risk_func = "amp", ci_method = "gamma")
calc_ltr(..., risk_func = "wun", ci_method = "gamma")
calc_ltr(..., risk_func = "cumulative", ci_method = "delta")