--- title: "Multiple decrement tables in lifecontingencies" author: "Giorgio Alfredo Spedicato" date: "`r format(Sys.time(), '%d %B, %Y')`" output: rmarkdown::html_vignette bibliography: lifecontingenciesBiblio.bib vignette: > %\VignetteEngine{knitr::rmarkdown} %\VignetteIndexEntry{Multiple decrement tables in lifecontingencies} %!\VignetteEncoding{UTF-8} \usepackage[utf8]{inputenc} --- ```{r echo=FALSE} 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: ```{r mdt1} 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. ```{r md3, eval=FALSE} print(valdezMdt) ``` ```{r md3b, eval = requireNamespace("markovchain", quietly = TRUE)} 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. ```{r mdt4} 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`: ```{r mdt.dx1} 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: ```{r mdt.randomSamples} 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`: ```{r mdt.udd1} 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]: ```{r mdt.udd2} 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. ```{r asdt.matrix1} qprime <- independentRatesFromMdt(valdezMdt, x = 50:54) qprime ``` Each cell equals the value that `qxt.prime.fromMdt` would return for the same (age, decrement) pair: ```{r asdt.matrix2} qxt.prime.fromMdt(valdezMdt, x = 53, decrement = "accidents") qprime["53", "accidents"] ``` A subset of ages may be specified: ```{r asdt.matrix3} 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: ```{r asdt.build1} 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$: ```{r asdt.build2} 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: ```{r asdt.roundtrip} 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$. ```{r plot.area, fig.width=7, fig.height=4} plot(valdezMdt) ``` ```{r plot.bar, fig.width=7, fig.height=4} plot(valdezMdt, type = "bar") ``` ```{r plot.prob, fig.width=7, fig.height=4} 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: ```{r plot.cdc, fig.width=7, fig.height=5, eval=FALSE} 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: ```{r mdt.act1} 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`: ```{r mdt.act2} 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: ```{r mdt.act3} 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$: ```{r mdt.act4} 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: ```{r deadifalco.1a} 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: ```{r deadifalco.1b} 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`): ```{r cdc.data} 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. ```{r cdc.graduate} 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. ```{r cdc.build} 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: ```{r cdc.compute} 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. ```{r long1} 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: ```{r long2} 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: ```{r long3} cdcLong <- mdtToLong(cdcMdt, x = 70) tapply(cdcLong$count, cdcLong$status, sum)[c("cancer", "heart_disease")] / cdcMdt@table$lx[cdcMdt@table$x == 70] ``` ## References