Multiple decrement tables in lifecontingencies"

In lifecontingencies: Financial and Actuarial Mathematics for Life Contingencies

knitr::opts_chunk$set(
  warning = FALSE,
  message = FALSE,
  collapse = TRUE,
  error = TRUE
)
library(lifecontingencies)

Introduction

Until now no R package provided a good tool to manage multiple decrement tables, even if @deshmukh2012multiple gives an R-based treatment of multiple decrement tables with applications. The topic is deeply related to multistate analysis of life histories, on which @willekens2014multistate provide a very good introduction. The mdt class of lifecontingencies has been specifically engineered to manage multiple decrement models in R.

This vignette walks through the whole workflow the package supports around multiple decrements: creating and validating a combined mdt table, computing decrement probabilities from it, deriving the Associated Single Decrement Table (ASDT) implied by each cause under the uniform distribution of decrements (UDD) assumption, using a mdt object in actuarial insurance and annuity calculations, and finally building a real, cause-specific mdt object end-to-end from an official government mortality report. For a shorter, self-contained example see the "Multiple Decrement Models" section of the main package vignette (vignette("an_introduction_to_lifecontingencies_package", package = "lifecontingencies")); this vignette develops the same material in full.

Following the notation in @finanMLC, let $l_x^{(\tau)} = \sum_{j=1\ldots m} l_x^{(j)}$ be the survivors to age $x$ that will, at future ages, be fully depleted by $m$ causes of decrement. $d_x^{(j)} = l_x^{(j)} - l_{x+1}^{(j)}$ is the expected number of lives exiting the population between ages $x$ and $x+1$ due to decrement $j$, so that ${}nd_x^{(j)} = \sum{t=0\ldots n-1} d_{x+t}^{(j)}$. The probability that a life aged $x$ will leave the group within one year as a result of decrement $j$ is ${}nq_x^{(j)} = \dfrac{{}_nd_x^{(j)}}{l_x^{(\tau)}}$, so that $q_x^{(\tau)} = \sum{j=1}^m q_x^{(j)}$ and ${}_tq_x^{(\tau)} = 1 - {}_tp_x^{(\tau)} = \sum_j {}_tq_x^{(j)}$.

The mdt class

Examples in this section are worked on slides provided by @valdezSlides (p. 4). We create a mdt class object:

valdezDf <- data.frame(
  x = c(50:54),
  lx = c(4832555, 4821937, 4810206, 4797185, 4782737),
  heart = c(5168, 5363, 5618, 5929, 6277),
  accidents = c(1157, 1206, 1443, 1679, 2152),
  other = c(4293, 5162, 5960, 6840, 7631)
)
valdezMdt <- new("mdt", name = "ValdezExample", table = valdezDf)

The mdt class is an S4 class object [@chambers2008software] comprised of a character slot name and a data.frame slot table, made of:

1. x: the age, from 0 to $\omega$.
2. lx: the subjects living (at risk) at the beginning of age $x$.
3. one or more columns for the different causes of decrement.

Values within table represent the absolute number of subjects at risk at the beginning of age $x$ dying of cause $j$ during the period $x$ to $x+1$.

setValidity performs consistency checks when a mdt object is created. In particular it verifies that:

1. x and lx exist and are consistent: x starts from 0 and increases by one, and $l_x = l_{x-1} - (d_{x-1,1} + d_{x-1,2} + \ldots + d_{x-1,k})$ for every $x$.
2. If the decrements (or x and lx) are supplied only for a partial age range, the table is completed below assuming a decrement rate of 0.01 for the first cause of death. This rate can be changed via the bottomCompletionSurvival argument of new("mdt", ...) (default 0.99, i.e. a 0.01 decrement rate, unchanged from previous releases); it only affects the synthetically reconstructed ages below the lowest age actually supplied and leaves all originally supplied rows untouched.
3. If the decrements at the last supplied age $\omega$ do not sum to $l_\omega$, the table is extended by one row such that $l_{\omega+1} = l_\omega - (d_{\omega,1} + d_{\omega,2} + \ldots + d_{\omega,j})$.

The internal function .tableSanitizer implements this completion logic and reports on the console what it did. The table can be viewed with print and show (output omitted here for brevity), and coerced to a data.frame. Coercion to a markovchainList object needs the optional markovchain package: the chunk below is evaluated only when that package is installed, and is skipped otherwise.

print(valdezMdt)

valdezDf2 <- as(valdezMdt, "data.frame")
valdezMarkovChainList <- as(valdezMdt, "markovchainList")

Two specific methods are defined for mdt objects: getOmega, returning the maximum attainable age (as for the lifetable class), and getDecrements, returning the decrement names (the columns of table other than x and lx). A summary method is available as well.

getOmega(valdezMdt)
getDecrements(valdezMdt)
summary(valdezMdt)

Decrement probabilities calculation

lifecontingencies makes it easy to compute $d_x^{(j)}$, ${}_nd_x^{(j)}$ and ${}_nd_x^{(\tau)}$ via the dxt function, and the corresponding probabilities via qxt/pxt:

dxt(valdezMdt, x = 51, decrement = "other")
dxt(valdezMdt, x = 51, t = 2, decrement = "other")
dxt(valdezMdt, x = 51)
pxt(valdezMdt, x = 50, t = 3)
qxt(valdezMdt, x = 53, t = 2, decrement = 1)

It is also possible to generate random trajectories of a life subject to multiple causes of decrement:

rmdt(n = 2, object = valdezMdt, x = 50, t = 2)

Associated Single Decrement Table calculation

For each force of decrement $\mu_j(x+t)$, the Associated Single Decrement Table (ASDT) is a decrement model that assumes survivorship depends only on $j$. Within the ASDT the following identities hold:

$$ {}tp_x^{\prime(j)} = \exp!\left(\int_a^b \mu_j(x+s)\,ds\right), \qquad {}_tq_x^{\prime(j)} = 1 - {}_tp_x^{\prime(j)}, \qquad {}_tp_x^{(\tau)} = \prod{j=1}^m {}_tp_x^{\prime(j)}. $$

Assuming a uniform distribution of decrements (UDD), i.e. ${}_tq_x^{(j)} = t \cdot q_x^{(j)}$ for $t \in [0,1]$, then:

$$ p_x^{\prime(j)} = \left(1 - t \cdot q_x^{(\tau)}\right)^{\frac{q_x^{(j)}}{q_x^{(\tau)}}} = 1 - q_x^{\prime(j)}. $$

This is implemented in qxt.prime.fromMdt:

qxt.prime.fromMdt(object = valdezMdt, x = 53, decrement = "accidents")

If UDD also holds for each associated single decrement, then:

$$ {}tq_x^{(j)} = {}_tq_x^{\prime(j)} \cdot \int_0^t \prod{i \ne j}\left(1 - s \cdot q_x^{\prime(i)}\right) ds, $$

of which the case $m=2$, $t=1$ is a particular case:

$$ q_x^{(2)} = q_x^{\prime(2)}\left(1 - 0.5\, q_x^{\prime(1)}\right). $$

qxt.fromQxprime computes this; the example below replicates @finanMLC [Example 67.2]:

qxt.fromQxprime(qx.prime = 0.01, other.qx.prime = c(0.03, 0.06))

Extracting the full ASDT matrix: independentRatesFromMdt()

While qxt.prime.fromMdt returns a single scalar for one (age, decrement) pair, independentRatesFromMdt computes the full matrix of independent rates $q'^{(j)}_x$ for every age and every decrement in the table. This is convenient for inspecting and comparing rates across ages, and for the round-trip reconstruction described below.

qprime <- independentRatesFromMdt(valdezMdt, x = 50:54)
qprime

Each cell equals the value that qxt.prime.fromMdt would return for the same (age, decrement) pair:

qxt.prime.fromMdt(valdezMdt, x = 53, decrement = "accidents")
qprime["53", "accidents"]

A subset of ages may be specified:

independentRatesFromMdt(valdezMdt, x = 51:52)

Building an mdt from independent rates: buildMdtFromIndependentRates()

The inverse operation — constructing a combined multiple-decrement table from a matrix of ASDT independent rates — is provided by buildMdtFromIndependentRates. This is useful when independent rates are estimated separately (e.g. from cause-specific mortality studies or from different data sources) and need to be combined into a single mdt object.

For each age, the absolute (combined) rate of decrement $j$ is obtained by the UDD integration formula: $$ q_x^{(j)} = q'^{(j)}x \int_0^1 \prod{i \ne j} \bigl(1 - s\, q'^{(i)}_x\bigr)\,ds, $$ and the survivorship column is built recursively from $p_x^{(\tau)} = \prod_j (1 - q'^{(j)}_x)$.

The following example replicates @finanMLC [Example 67.4]: three decrements (death, disability, retirement) over two ages, with 1000 initial lives:

qp <- matrix(c(0.010, 0.030, 0.100,
                0.013, 0.050, 0.200),
             nrow = 2, byrow = TRUE,
             dimnames = list(NULL, c("death", "disability", "retirement")))
qp

mdt674 <- buildMdtFromIndependentRates(
  x = 60:61, qx.primes = qp,
  radix = 1000, name = "Finan 67.4"
)
print(mdt674)

We can verify the key numerical results from the example. The combined survival probability at age 60 is $p^{(\tau)}{60} = (1-0.010)(1-0.030)(1-0.100) = 0.86427$, giving $l^{(\tau)}{61} = 864.27$:

ptau60 <- prod(1 - qp[1, ])
round(ptau60, 5)
tbl674 <- mdt674@table
round(tbl674$lx[tbl674$x == 61], 2)

Round-trip: mdt $\to$ ASDT $\to$ mdt

The two functions compose naturally: extracting independent rates from an mdt and then rebuilding should recover the original table (up to floating-point rounding). We verify this on the Valdez example:

qprime <- independentRatesFromMdt(valdezMdt, x = 50:54)
rebuilt <- buildMdtFromIndependentRates(
  x = 50:54, qx.primes = qprime,
  radix = valdezMdt@table$lx[valdezMdt@table$x == 50],
  name = "Roundtrip"
)

## Compare lx
data.frame(
  age = 50:54,
  original = valdezMdt@table$lx[valdezMdt@table$x %in% 50:54],
  rebuilt  = round(rebuilt@table$lx[rebuilt@table$x %in% 50:54], 0)
)

Visualising mdt objects: plot()

The plot S4 method for mdt objects produces a ggplot2 visualisation of the decrement structure. Three views are available:

• "area" (default): stacked-area chart of decrement counts $d^{(j)}_x$.
• "bar": stacked bar chart of decrement counts.
• "probability": line chart of decrement-specific probabilities $q^{(j)}_x = d^{(j)}_x / l^{(\tau)}_x$.

plot(valdezMdt)

plot(valdezMdt, type = "bar")

plot(valdezMdt, type = "probability")

The plot method also works on larger tables; here is the CDC-based mdt object built later in this vignette, aggregated across its 18 causes:

plot(cdcMdt, type = "probability")

Actuarial applications

Two functions value contracts on a mdt object. Axn.mdt gives the actuarial present value (APV) of a term insurance paying, at the end of the year of decrement, a benefit $b_j$ that may depend on the cause $j$: $$ \sum_{j \in J} b_j \sum_{h=0}^{n-1} v^{h+1}\,{}hp^{(\tau)}_x\, q^{(j)}{x+h}, $$ where $J$ is the set of covered decrements (all of them if decrement is omitted). axn.mdt gives the APV of an annuity payable while the insured has not yet left the table, with optional deferment m, k payments per year and payments in advance or in arrears. Benefit premiums and reserves follow from the equivalence principle.

The example below, from @finanMLC [Problem 68.1], considers a 3-year term issued to (16) paying 20,000 at the end of the year of death if death results from an accident:

myTable <- data.frame(
  x = c(16, 17, 18),
  lx = c(20000, 17600, 14520),
  da = c(1300, 1870, 2380),
  doc = c(1100, 1210, 1331)
)
myMdt <- new("mdt", table = myTable, name = "Sample")

The value of $A^{1}_{16:\overline{3}|}$, restricted to decrement da (accidental death), is computed by Axn.mdt:

20000 * Axn.mdt(object = myMdt, x = 16, n = 3, i = .1, decrement = "da")

If n is omitted the cover runs to the end of the table, which here would include the synthetic closing row (age 19) added by new("mdt", ...): for term covers it is safer to always pass n.

@finanMLC [Example 69.1] uses the same table to price a level annual premium. Finan's printed solution uses the column of non-accidental deaths (1100, 1210, 1331), giving $A = 3000/20000$, $\ddot a^{(\tau)}_{16:\overline{3}|} = 2.4$ and a premium of 1250:

A691 <- Axn.mdt(myMdt, x = 16, n = 3, i = 0.10, decrement = "doc")
a691 <- axn.mdt(myMdt, x = 16, n = 3, i = 0.10)
c(APV = 20000 * A691, annuity = a691, premium = 20000 * A691 / a691)

Benefits that depend on the cause are passed through benefits. @finanMLC [Example 68.1] pays 1 for cause 1 and 2 for cause 2 over two years at $i = 50\%$ (printed answer 0.6852), and [Example 69.2] computes the benefit reserve ${}_2V = 11.091$ of a 4-year term paying 2000 for decrement 1 and 1000 for decrement 2, with annual premium 34 and $v = 0.95$:

m681 <- new("mdt", table = data.frame(x = 50:51, lx = c(1200, 800),
                                      d1 = c(100, 200), d2 = c(300, 300)))
Axn.mdt(m681, x = 50, n = 2, i = 0.5, decrement = c("d1", "d2"),
        benefits = c(1, 2))

m692 <- new("mdt", table = data.frame(x = 41:43, lx = c(800, 776, 752),
                                      d1 = 8, d2 = 16))
i692 <- 1 / 0.95 - 1
Axn.mdt(m692, x = 42, n = 2, i = i692, decrement = c("d1", "d2"),
        benefits = c(2000, 1000)) - 34 * axn.mdt(m692, x = 42, n = 2, i = i692)

Another example is inspired by data in @de2016assicurazioni. We use life tables from that source and a small helper function (still based on single decrement tables) to compute the APV of a constant or variable annuity with multiple decrement, where benefits are payable while the annuitant is in the current state and cease upon transition to another state:

axnmdt.firsttype <- function(object, x, n, i, payment = "advance", delta = 0) {
  # delta is the annuity indexing
  out <- numeric(1)
  if (!(class(object) %in% c("lifetable", "actuarialtable", "mdt")))
    stop("Error! Only lifetable, actuarialtable or mdt classes are accepted")
  if (missing(object))
    stop("Error! Need a Multiple decrement table")
  if (missing(x))
    stop("Error! Need age!")
  if (x > getOmega(object)) {
    stop("Age greater than Omega")
  }
  if (class(object) == "mdt") {
    if (x < min(object@table$x)) {
      stop("Age lower than minimum age")
    }
  }
  if (class(object) == "actuarialtable") {
    if (x < min(object@x)) {
      stop("Age lower than minimum age")
    }
  }
  if (!(missing(i))) {
    interest <- i
  } else {
    if (class(object) == "actuarialtable") {
      interest <- object@interest
    } else {
      stop("Needed Interest Rate ")
    }
  }
  if (missing(n))
    n <- (getOmega(object) - x)
  if (n == 0) {
    stop("Contract duration equal to zero")
  }
  probs <- numeric(n)
  times <- seq(from = 0, to = n - 1, by = 1)
  if (payment == "arrears") times <- times + 1
  for (j in 1:length(times)) probs[j] <- pxt(object, x, times[j])
  out <- sum(apply(cbind(probs, ((1 + interest) / (1 + delta))^-times), 1, prod))
  return(out)
}

Data and examples are taken from @de2016assicurazioni:

data("de_angelis_di_falco")
HealthyMaleTable2013 <- de_angelis_di_falco$HealthyMaleTable2013
DAT <- new("actuarialtable",
  x = de_angelis_di_falco$DisabledMaleLifeTable$age,
  lx = de_angelis_di_falco$DisabledMaleLifeTable$'2013',
  name = "DisabledTable", i = 0.03
)
axnmdt.firsttype(DAT, x = 65, n = 10, i = 0.03, payment = "arrears", delta = 0.02)
axnmdt.firsttype(DAT, 65, 10, payment = "arrears", delta = 0.02)
axnmdt.firsttype(DAT, 65, 10, payment = "arrears", i = 0.03, delta = 0.02)
# Last case equal to axn
axnmdt.firsttype(DAT, 65, 10, payment = "arrears", delta = 0)
axn(DAT, 65, 10, payment = "arrears")

A worked example on real cause-of-death data, with a note on graduation

The examples above use small, illustrative tables. This section works through the full pipeline — from an official, publicly available cause-specific mortality report to a validated mdt object — using real data: Table 8 of @ariasNVSR2013, "United States Life Tables Eliminating Certain Causes of Death, 1999-2001", published by the U.S. National Center for Health Statistics. Table 8 tabulates, for 22 five-year (and one open-ended) age bands from birth to 100+, the number of decrements attributable to each of 33 causes of death out of a radix of $l_0 = 10{,}000{,}000$.

From an official report to a mutually exclusive partition

A mdt object requires its non-x/lx columns to be a partition of total mortality: at every age the columns must sum to exactly the decrement implied by lx. The 33 causes tabulated by @ariasNVSR2013, however, are not mutually exclusive: several rows are parent categories that contain other rows as subtypes (e.g. "Malignant neoplasms" contains site-specific cancers such as colon, lung, breast and prostate), and others are cross-cutting aggregates that overlap several disease groups (e.g. "Alcohol-induced causes"). Summing all 33 columns as if they were disjoint would double count decrements.

We therefore curate a mutually exclusive partition made of the 17 top-level (non-nested, non-aggregate) causes of the original 33, plus a residual other category that closes the table to total mortality at every age. The columns kept are septicemia, hiv, cancer, diabetes, alzheimer, heart_disease, hypertension, cerebrovascular, influenza_pneumonia, copd, pneumonitis, liver_disease, nephritis, congenital, accidents, suicide and homicide; other is obtained, for each age band, as the total decrements minus the sum of these 17 columns, and is non-negative at every age band, confirming no double counting remains.

We embed the resulting curated, mutually exclusive 5-year table below (partition_sum is the sum of the 17 named causes, kept only as a diagnostic column to show that partition_sum + other == total_dx at every age; it is dropped before building the mdt):

bandStart <- c(0, 1, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65,
               70, 75, 80, 85, 90, 95, 100)
cdc5yr <- data.frame(
  age_group = c("0-1", "1-5", "5-10", "10-15", "15-20", "20-25", "25-30",
    "30-35", "35-40", "40-45", "45-50", "50-55", "55-60", "60-65",
    "65-70", "70-75", "75-80", "80-85", "85-90", "90-95", "95-100",
    "100 and over"),
  lx = c(10000000.0, 9930523.0, 9917600.0, 9909689.0, 9899801.0, 9866403.0,
    9820252.0, 9775095.0, 9720086.0, 9642151.0, 9527362.0, 9360142.0,
    9123237.0, 8764226.0, 8233035.0, 7489085.0, 6464402.0, 5088466.0,
    3451511.0, 1849572.0, 687944.7, 147917.1),
  total_dx = c(69477.0, 12923.0, 7911.0, 9888.0, 33398.0, 46151.0, 45157.0,
    55009.0, 77935.0, 114789.0, 167220.0, 236905.0, 359011.0, 531191.0,
    743950.0, 1024683.0, 1375936.0, 1636955.0, 1601939.0, 1161627.3,
    540027.6, 147917.1),
  partition_sum = c(20749.0, 9384.9, 6159.0, 7847.6, 29571.0, 40140.9,
    38078.3, 45424.2, 63928.2, 94335.2, 138926.7, 202948.7, 315395.2,
    470524.7, 658241.3, 900367.8, 1191618.0, 1395849.7, 1344478.8,
    958568.3, 436904.4, 116411.8),
  other = c(48728.0, 3538.1, 1752.0, 2040.4, 3827.0, 6010.1, 7078.7,
    9584.8, 14006.8, 20453.8, 28293.3, 33956.3, 43615.8, 60666.3,
    85708.7, 124315.2, 184318.0, 241105.3, 257460.2, 203059.0, 103123.2,
    31505.3),
  septicemia = c(722.6, 247.7, 91.8, 76.5, 105.3, 194.7, 257.8, 418.8,
    671.8, 1120.7, 1847.0, 2783.5, 4280.8, 6452.9, 9611.9, 13677.9,
    19662.0, 23350.0, 23250.0, 15851.9, 6872.4, 1626.4),
  hiv = c(28.3, 43.0, 58.0, 62.8, 88.2, 427.6, 1680.0, 4167.0, 6114.4,
    6706.2, 5899.3, 4113.8, 2677.1, 1727.4, 1046.5, 539.9, 257.5, 81.0,
    25.5, 19.5, 10.6, 2.7),
  cancer = c(187.9, 1060.0, 1201.0, 1242.0, 1814.4, 2526.2, 3518.9,
    6114.5, 11915.0, 23559.7, 44358.6, 77442.7, 131212.1, 198969.2,
    267824.2, 332678.3, 368548.5, 334200.1, 233736.2, 113426.7, 35785.0,
    6122.5),
  diabetes = c(6.7, 16.9, 19.3, 58.8, 120.0, 277.8, 513.9, 940.3, 1524.6,
    2623.1, 4656.5, 8029.6, 13462.8, 20843.6, 29206.0, 37986.0, 46899.8,
    48482.2, 38608.7, 21835.1, 8005.9, 1545.5),
  alzheimer = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.9, 0.8, 2.4, 2.8, 8.6, 42.5,
    126.0, 482.8, 1366.0, 3477.9, 9876.6, 25418.6, 49458.4, 63722.0,
    54057.0, 25433.7, 5666.6),
  heart_disease = c(1244.1, 496.3, 257.8, 402.5, 990.5, 1632.3, 2575.5,
    4841.8, 9546.1, 19013.9, 33655.2, 56478.6, 92602.3, 141365.9,
    202096.5, 288960.8, 415183.0, 537740.9, 577672.7, 457905.7, 227285.0,
    66431.5),
  hypertension = c(3.3, 4.2, 1.6, 3.2, 11.4, 33.8, 64.9, 138.3, 220.6,
    485.8, 838.2, 1359.1, 1972.7, 3237.0, 4635.8, 6594.0, 10279.8,
    14091.1, 15639.9, 12697.0, 6581.1, 1933.8),
  cerebrovascular = c(278.6, 123.9, 71.7, 106.2, 164.9, 324.5, 513.9,
    926.9, 1870.1, 3586.8, 5965.7, 8723.2, 13848.4, 22185.8, 34745.3,
    58994.0, 100084.4, 145206.2, 161264.4, 122436.7, 54220.7, 12630.6),
  influenza_pneumonia = c(755.1, 290.7, 112.0, 107.1, 166.6, 298.6, 337.8,
    506.5, 893.9, 1322.6, 1880.1, 2539.4, 3857.3, 6272.5, 9849.1,
    17726.2, 31608.2, 51735.7, 64665.0, 58841.1, 32133.0, 10956.1),
  copd = c(94.0, 124.7, 112.0, 195.6, 218.0, 262.2, 313.4, 429.9, 686.0,
    1299.8, 2657.7, 5557.3, 13792.5, 27507.5, 49137.5, 76948.9, 100455.3,
    105576.7, 81889.6, 41622.6, 14325.5, 2877.8),
  pneumonitis = c(31.6, 31.2, 13.7, 21.7, 25.3, 62.3, 72.4, 105.1, 131.2,
    243.3, 336.5, 528.8, 866.2, 1408.2, 2466.9, 4832.3, 8971.9, 14562.1,
    17864.6, 14771.2, 7068.0, 1928.4),
  liver_disease = c(9.1, 4.2, 1.6, 4.0, 18.0, 57.1, 198.0, 769.7, 2337.2,
    4790.8, 7989.0, 8922.6, 9780.0, 11092.9, 11344.1, 10895.8, 9658.5,
    6664.8, 3434.6, 1153.5, 225.4, 24.3),
  nephritis = c(383.3, 37.9, 26.6, 33.0, 62.9, 151.5, 218.2, 365.1, 599.5,
    977.3, 1607.2, 2564.0, 4302.1, 7118.9, 11608.2, 16664.5, 24123.1,
    29663.4, 29579.4, 21013.8, 9269.0, 2227.8),
  congenital = c(13912.0, 1349.1, 472.9, 495.9, 573.3, 586.8, 543.4,
    596.6, 578.2, 636.4, 732.9, 848.9, 1003.8, 1066.8, 970.3, 879.6,
    1001.5, 1063.5, 870.6, 504.0, 236.0, 56.6),
  accidents = c(2249.1, 4586.7, 3332.5, 3845.8, 16420.4, 19138.2,
    14938.4, 14067.2, 15926.3, 17254.8, 16556.9, 14316.4, 13565.0,
    13626.6, 14401.4, 17356.4, 23649.9, 29109.6, 29211.3, 21255.5,
    9191.2, 2357.0),
  suicide = c(0.0, 0.0, 12.9, 655.5, 3959.4, 6077.0, 6007.4, 6197.6,
    6794.9, 7200.1, 7130.0, 6521.5, 6005.7, 4995.4, 4726.3, 4904.1,
    5078.7, 4297.4, 2719.2, 1042.0, 214.9, 18.9),
  homicide = c(843.2, 968.4, 373.7, 537.2, 4832.4, 8089.4, 6323.6,
    4836.4, 4115.3, 3505.5, 2773.4, 2093.2, 1683.3, 1288.2, 1093.5,
    852.3, 737.2, 566.6, 325.3, 135.1, 46.9, 5.4)
)
## sanity check: the curated partition closes exactly to total mortality
stopifnot(all(abs(cdc5yr$partition_sum + cdc5yr$other - cdc5yr$total_dx) < 1e-6))

Graduating from 5-year bands to single years of age

The mdt class requires one row per single year of age, while Table 8 is tabulated in 5-year bands. A modern, purpose-built way to "ungroup" banded counts into single-year counts is Eilers' penalized composite link model (PCLM, @eilers2007pclm), implemented in the ungroup and DemoTools packages; it lets the within-band shape — including how the cause-of-death mix evolves within a band — be estimated smoothly, and is generally the preferred approach for this task.

This vignette does not use PCLM. lifecontingencies does not depend on ungroup or DemoTools, and adding either as a Suggests purely to compile one example was judged not worth the extra dependency surface. Instead, we graduate lx with a monotone cubic Hermite interpolant (the Fritsch-Carlson method, @fritschCarlson1980), available in base R as stats::splinefun(..., method = "monoH.FC"), and then split each single-year total decrement across causes using each band's own, fixed cause mix (i.e. the cause shares are constant within a band and only the total number of decrements is graduated). This is a materially weaker assumption than PCLM: real cause-of-death mixes drift gradually within a 5-year band (e.g. the balance between accidents and degenerative disease shifts within the 15-20 age band), and a constant-share assumption cannot reproduce that drift. The single-year figures produced below should therefore be read as a reasonable, fully reproducible, dependency-free approximation, not as a substitute for a proper PCLM ungrouping. Because the split is proportional to a graduated total that itself sums exactly across each band, the construction below reaggregates exactly back to the original 5-year, cause-specific totals (checked below), which a naive graduation might not guarantee.

causeCols <- setdiff(names(cdc5yr),
                      c("age_group", "lx", "total_dx", "partition_sum", "other"))
causeCols <- c(causeCols, "other")   # 17 named causes + residual

## monotone Hermite (Fritsch-Carlson) graduation of lx, ages 0..100
lxFun <- splinefun(x = bandStart, y = cdc5yr$lx, method = "monoH.FC")
ages <- 0:100
lxSingle <- c(lxFun(ages), 0)          # append the closing lx=0 at age 101
dxTotalSingle <- -diff(lxSingle)       # single-year total decrements, ages 0..100
band <- findInterval(ages, bandStart)  # which 5-year band each single age falls in

cdcSingle <- data.frame(x = ages, lx = lxSingle[1:101])
for (cs in causeCols) {
  shareBand <- cdc5yr[[cs]] / cdc5yr$total_dx
  cdcSingle[[cs]] <- dxTotalSingle * shareBand[band]
}

## verify exact reaggregation back to the original 5-year cause totals
recheck <- aggregate(cdcSingle[causeCols],
                      by = list(age_group = cdc5yr$age_group[band]), FUN = sum)
recheck <- recheck[match(cdc5yr$age_group, recheck$age_group), ]
max(abs(as.matrix(recheck[causeCols]) - as.matrix(cdc5yr[causeCols])))

Building and validating the mdt object

The single-year table cdcSingle (columns x, lx and the 18 mutually exclusive causes) is now a valid input for new("mdt", ...); because it already spans exactly $x = 0, \ldots, 100$ and closes exactly to lx[0], no meaningful synthetic bottom- or top-completion is needed from .tableSanitizer (described earlier in this vignette). At this scale, floating-point arithmetic can still make the very last row's decrements sum to a value that differs from lx by a fraction of a unit, in which case .tableSanitizer adds one further, economically negligible completion row at $x = 101$ — a harmless artifact of exact floating-point comparison, not a data or modelling error.

cdcMdt <- new("mdt",
              name = "NCHS US Life Table 1999-2001, cause-specific (Table 8)",
              table = cdcSingle)
getOmega(cdcMdt)
getDecrements(cdcMdt)
isTRUE(validObject(cdcMdt))

We can now use the same dxt/qxt/pxt functions used throughout this vignette on a real, cause-specific national life table. For instance, at age 70:

qxt(cdcMdt, x = 70, t = 1)                          # total decrement rate
qxt(cdcMdt, x = 70, t = 1, decrement = "cancer")     # cancer only
qxt(cdcMdt, x = 70, t = 1, decrement = "heart_disease")
qxt(cdcMdt, x = 70, t = 1, decrement = "cancer") /
  qxt(cdcMdt, x = 70, t = 1)                          # cancer's share of q70

At age 70 in this table, cancer alone accounts for roughly a third of all decrements, illustrating how a mdt object built this way lets standard lifecontingencies functions (dxt, qxt, pxt, Axn.mdt, ...) be applied directly to official, cause-specific national mortality data, once the table has been reduced to a mutually exclusive partition and graduated to single years of age.

Competing-risks analysis: mdtToLong()

A mdt object is a life-table summary of a competing-risks process. For analyses with standard survival tools, mdtToLong turns it into an aggregated long data set: one row per (exit time, cause), a status factor whose first level is "censored" (the layout expected by survival::Surv for multi-state data) and a count column to be used as case weights. Lives still in the table at the end of the follow-up are censored there.

valdezLong <- mdtToLong(valdezMdt, x = 50, t = 5)
head(valdezLong)

With exits placed at the end of the year of age (the default), the Aalen-Johansen estimate of the cumulative incidence of each cause computed by survival coincides with the decrement probabilities of the table:

if (requireNamespace("survival", quietly = TRUE)) {
  fit <- survival::survfit(survival::Surv(time, status) ~ 1,
                           data = valdezLong, weights = count)
  aj <- summary(fit, times = 1:5)$pstate[, -1]
  colnames(aj) <- getDecrements(valdezMdt)
  byTable <- sapply(getDecrements(valdezMdt), function(d)
    qxt(valdezMdt, x = 50, t = 1:5, decrement = d))
  max(abs(aj - byTable))
}

The same conversion applies to the cause-specific national table built above, e.g. to obtain the probability that a life aged 70 eventually dies of cancer or of heart disease:

cdcLong <- mdtToLong(cdcMdt, x = 70)
tapply(cdcLong$count, cdcLong$status, sum)[c("cancer", "heart_disease")] /
  cdcMdt@table$lx[cdcMdt@table$x == 70]

References

lifecontingencies documentation built on Oct. 10, 2026, 9:06 a.m.