tests/testthat/test-published-nchs-us2019.R

In lifecontingencies: Financial and Actuarial Mathematics for Life Contingencies

## Published worked example: complete US period life tables 2019 (NCHS).
##
## Source: Arias E., Xu J. (2022). "United States Life Tables, 2019".
##   National Vital Statistics Reports, vol. 70, no. 19. Hyattsville, MD:
##   National Center for Health Statistics. Table 1 (total population) and
##   Table 2 (males): single-year ages 0..100+, columns qx, lx, dx, Lx, Tx, ex.
##   https://www.cdc.gov/nchs/data/nvsr/nvsr70/nvsr70-19.pdf
##   (These tables are the standard reference of US demography courses.)
##
## The published columns are rounded (lx, dx, Lx, Tx to integers, qx to 6 digits,
## ex to 1 decimal), so every comparison below uses a tolerance that reflects
## that rounding. The tables were transcribed and checked for internal
## consistency (dx = lx - lx+1, Tx = sum of Lx up to rounding, ex = Tx/lx).
##
## Convention differences between NCHS and lifecontingencies, documented here:
##  * Lxt() uses Lx = lx - fxt * dx with fxt = 0.5 by default. NCHS uses a_x = 0.5
##    for x >= 1 (so Lx matches to rounding) but a_0 ~ 0.126-0.128 (separation
##    factor of infants), so L0 matches only with fxt = 1 - a0.
##  * The last interval 100+ is closed by NCHS with L100 = l100 / m100 (open
##    interval), whereas Tx() closes it with l_omega / 2 (omega = 100). Hence
##    Tx_package(x) + (L100_published - l100 / 2) = Tx_published(x) up to rounding.

context("Published example: NCHS United States life tables 2019 (total, males)")

nchs_total <- read.csv(text = "
age,qx,lx,dx,Lx,Tx,ex
0,0.005575,100000,557,99513,7884823,78.8
1,0.000379,99443,38,99424,7785310,78.3
2,0.000234,99405,23,99393,7685886,77.3
3,0.000180,99382,18,99373,7586493,76.3
4,0.000145,99364,14,99357,7487120,75.4
5,0.000137,99349,14,99343,7387764,74.4
6,0.000124,99336,12,99330,7288421,73.4
7,0.000114,99323,11,99318,7189091,72.4
8,0.000105,99312,10,99307,7089774,71.4
9,0.000097,99302,10,99297,6990467,70.4
10,0.000095,99292,9,99287,6891170,69.4
11,0.000104,99282,10,99277,6791883,68.4
12,0.000131,99272,13,99266,6692606,67.4
13,0.000181,99259,18,99250,6593340,66.4
14,0.000248,99241,25,99229,6494090,65.4
15,0.000324,99217,32,99200,6394861,64.5
16,0.000401,99184,40,99165,6295661,63.5
17,0.000483,99145,48,99121,6196496,62.5
18,0.000569,99097,56,99069,6097375,61.5
19,0.000657,99040,65,99008,5998307,60.6
20,0.000748,98975,74,98938,5899299,59.6
21,0.000839,98901,83,98860,5800361,58.6
22,0.000917,98818,91,98773,5701501,57.7
23,0.000976,98728,96,98680,5602728,56.7
24,0.001020,98631,101,98581,5504048,55.8
25,0.001057,98531,104,98479,5405467,54.9
26,0.001097,98427,108,98373,5306988,53.9
27,0.001139,98319,112,98263,5208616,53.0
28,0.001190,98207,117,98148,5110353,52.0
29,0.001248,98090,122,98029,5012205,51.1
30,0.001309,97967,128,97903,4914176,50.2
31,0.001371,97839,134,97772,4816273,49.2
32,0.001436,97705,140,97635,4718501,48.3
33,0.001500,97565,146,97492,4620866,47.4
34,0.001564,97418,152,97342,4523374,46.4
35,0.001636,97266,159,97186,4426032,45.5
36,0.001711,97107,166,97024,4328846,44.6
37,0.001782,96941,173,96854,4231822,43.7
38,0.001847,96768,179,96679,4134967,42.7
39,0.001912,96589,185,96497,4038289,41.8
40,0.001990,96404,192,96309,3941792,40.9
41,0.002087,96213,201,96112,3845483,40.0
42,0.002205,96012,212,95906,3749371,39.1
43,0.002344,95800,225,95688,3653465,38.1
44,0.002501,95576,239,95456,3557777,37.2
45,0.002679,95336,255,95209,3462321,36.3
46,0.002877,95081,274,94944,3367112,35.4
47,0.003094,94808,293,94661,3272168,34.5
48,0.003334,94514,315,94357,3177507,33.6
49,0.003608,94199,340,94029,3083151,32.7
50,0.003904,93859,366,93676,2989121,31.8
51,0.004238,93493,396,93295,2895445,31.0
52,0.004639,93097,432,92881,2802151,30.1
53,0.005107,92665,473,92428,2709270,29.2
54,0.005615,92191,518,91933,2616842,28.4
55,0.006123,91674,561,91393,2524910,27.5
56,0.006630,91112,604,90810,2433517,26.7
57,0.007167,90508,649,90184,2342706,25.9
58,0.007754,89860,697,89511,2252522,25.1
59,0.008394,89163,748,88789,2163011,24.3
60,0.009090,88415,804,88013,2074222,23.5
61,0.009808,87611,859,87181,1986209,22.7
62,0.010518,86751,912,86295,1899028,21.9
63,0.011202,85839,962,85358,1812733,21.1
64,0.011890,84877,1009,84373,1727375,20.4
65,0.012621,83868,1059,83339,1643002,19.6
66,0.013499,82810,1118,82251,1559663,18.8
67,0.014468,81692,1182,81101,1477412,18.1
68,0.015571,80510,1254,79883,1396311,17.3
69,0.016804,79256,1332,78590,1316428,16.6
70,0.018143,77924,1414,77218,1237837,15.9
71,0.019683,76511,1506,75758,1160620,15.2
72,0.021193,75005,1590,74210,1084862,14.5
73,0.023529,73415,1727,72551,1010652,13.8
74,0.025583,71688,1834,70771,938101,13.1
75,0.028288,69854,1976,68866,867330,12.4
76,0.030996,67878,2104,66826,798464,11.8
77,0.034485,65774,2268,64640,731638,11.1
78,0.038071,63506,2418,62297,666999,10.5
79,0.042153,61088,2575,59800,604702,9.9
80,0.046589,58513,2726,57150,544902,9.3
81,0.051588,55787,2878,54348,487752,8.7
82,0.057742,52909,3055,51381,433404,8.2
83,0.064540,49854,3218,48245,382023,7.7
84,0.071641,46636,3341,44966,333778,7.2
85,0.080529,43295,3487,41552,288812,6.7
86,0.089744,39809,3573,38022,247260,6.2
87,0.100480,36236,3641,34416,209238,5.8
88,0.112250,32595,3659,30766,174822,5.4
89,0.125094,28936,3620,27126,144057,5.0
90,0.139038,25317,3520,23557,116930,4.6
91,0.154096,21797,3359,20117,93374,4.3
92,0.170257,18438,3139,16868,73256,4.0
93,0.187494,15299,2868,13864,56388,3.7
94,0.205752,12430,2558,11151,42524,3.4
95,0.224954,9873,2221,8762,31372,3.2
96,0.244995,7652,1875,6714,22610,3.0
97,0.265747,5777,1535,5010,15896,2.8
98,0.287060,4242,1218,3633,10886,2.6
99,0.308767,3024,934,2557,7253,2.4
100,1.000000,2090,2090,4696,4696,2.2
")

nchs_male <- read.csv(text = "
age,qx,lx,dx,Lx,Tx,ex
0,0.006080,100000,608,99470,7631130,76.3
1,0.000415,99392,41,99371,7531660,75.8
2,0.000257,99351,26,99338,7432289,74.8
3,0.000193,99325,19,99316,7332951,73.8
4,0.000153,99306,15,99298,7233635,72.8
5,0.000149,99291,15,99283,7134337,71.9
6,0.000137,99276,14,99269,7035053,70.9
7,0.000126,99263,12,99256,6935784,69.9
8,0.000113,99250,11,99244,6836528,68.9
9,0.000101,99239,10,99234,6737283,67.9
10,0.000095,99229,9,99224,6638049,66.9
11,0.000106,99219,11,99214,6538825,65.9
12,0.000145,99209,14,99202,6439611,64.9
13,0.000221,99194,22,99183,6340410,63.9
14,0.000324,99172,32,99156,6241226,62.9
15,0.000438,99140,43,99119,6142070,62.0
16,0.000554,99097,55,99069,6042951,61.0
17,0.000679,99042,67,99008,5943882,60.0
18,0.000809,98975,80,98935,5844874,59.1
19,0.000940,98895,93,98848,5745939,58.1
20,0.001077,98802,106,98748,5647091,57.2
21,0.001211,98695,120,98635,5548342,56.2
22,0.001324,98576,131,98510,5449707,55.3
23,0.001407,98445,138,98376,5351196,54.4
24,0.001465,98307,144,98235,5252820,53.4
25,0.001514,98163,149,98088,5154586,52.5
26,0.001564,98014,153,97937,5056497,51.6
27,0.001614,97861,158,97782,4958560,50.7
28,0.001669,97703,163,97621,4860778,49.8
29,0.001730,97540,169,97455,4763156,48.8
30,0.001794,97371,175,97284,4665701,47.9
31,0.001858,97196,181,97106,4568417,47.0
32,0.001926,97016,187,96922,4471311,46.1
33,0.001998,96829,193,96732,4374389,45.2
34,0.002072,96636,200,96535,4277656,44.3
35,0.002154,96435,208,96332,4181121,43.4
36,0.002241,96228,216,96120,4084789,42.4
37,0.002323,96012,223,95900,3988669,41.5
38,0.002399,95789,230,95674,3892769,40.6
39,0.002476,95559,237,95441,3797095,39.7
40,0.002569,95323,245,95200,3701654,38.8
41,0.002685,95078,255,94950,3606454,37.9
42,0.002824,94822,268,94688,3511504,37.0
43,0.002983,94555,282,94414,3416815,36.1
44,0.003162,94273,298,94124,3322402,35.2
45,0.003365,93974,316,93816,3228278,34.4
46,0.003596,93658,337,93490,3134462,33.5
47,0.003856,93321,360,93142,3040972,32.6
48,0.004154,92962,386,92769,2947831,31.7
49,0.004497,92575,416,92367,2855062,30.8
50,0.004867,92159,449,91935,2762695,30.0
51,0.005284,91711,485,91468,2670760,29.1
52,0.005789,91226,528,90962,2579292,28.3
53,0.006385,90698,579,90408,2488330,27.4
54,0.007034,90119,634,89802,2397921,26.6
55,0.007687,89485,688,89141,2308120,25.8
56,0.008338,88797,740,88427,2218979,25.0
57,0.009024,88057,795,87659,2130552,24.2
58,0.009768,87262,852,86836,2042893,23.4
59,0.010575,86410,914,85953,1956057,22.6
60,0.011445,85496,979,85007,1870104,21.9
61,0.012342,84517,1043,83996,1785097,21.1
62,0.013242,83474,1105,82922,1701101,20.4
63,0.014130,82369,1164,81787,1618180,19.6
64,0.015033,81205,1221,80595,1536393,18.9
65,0.016003,79984,1280,79344,1455798,18.2
66,0.017154,78704,1350,78029,1376454,17.5
67,0.018348,77354,1419,76644,1298425,16.8
68,0.019601,75935,1488,75191,1221780,16.1
69,0.020915,74446,1557,73668,1146590,15.4
70,0.022303,72889,1626,72077,1072922,14.7
71,0.023930,71264,1705,70411,1000845,14.0
72,0.025553,69558,1777,68670,930434,13.4
73,0.028257,67781,1915,66823,861764,12.7
74,0.030698,65866,2022,64855,794941,12.1
75,0.033891,63844,2164,62762,730086,11.4
76,0.037078,61680,2287,60537,667324,10.8
77,0.041154,59393,2444,58171,606788,10.2
78,0.045300,56949,2580,55659,548617,9.6
79,0.049951,54369,2716,53011,492958,9.1
80,0.054931,51653,2837,50235,439947,8.5
81,0.060709,48816,2964,47334,389712,8.0
82,0.067466,45852,3093,44306,342378,7.5
83,0.074945,42759,3205,41157,298073,7.0
84,0.083495,39554,3303,37903,256916,6.5
85,0.093585,36252,3393,34555,219013,6.0
86,0.103957,32859,3416,31151,184458,5.6
87,0.116315,29443,3425,27731,153307,5.2
88,0.129805,26018,3377,24330,125576,4.8
89,0.144448,22641,3270,21006,101246,4.5
90,0.160250,19371,3104,17819,80240,4.1
91,0.177194,16267,2882,14825,62422,3.8
92,0.195239,13384,2613,12078,47596,3.6
93,0.214316,10771,2308,9617,35519,3.3
94,0.234332,8463,1983,7471,25902,3.1
95,0.255162,6480,1653,5653,18431,2.8
96,0.276659,4826,1335,4159,12778,2.6
97,0.298653,3491,1043,2970,8619,2.5
98,0.320956,2448,786,2056,5649,2.3
99,0.343371,1663,571,1377,3594,2.2
100,1.000000,1092,1092,2217,2217,2.0
")

nchs_tables <- list(total = nchs_total, male = nchs_male)

nchs_lifetable <- function(d, name) {
  new("lifetable", x = d$age, lx = d$lx, name = name)
}

for (nm in names(nchs_tables)) {
  local({
    d <- nchs_tables[[nm]]
    lt <- nchs_lifetable(d, paste("US 2019", nm))
    omega <- 100
    ages <- 0:99
    tail_gap <- d$Lx[101] - d$lx[101] / 2   # NCHS open-interval L100 vs package l_omega/2

    test_that(paste("US 2019", nm, ": table structure and radix"), {
      expect_equal(getOmega(lt), omega)
      expect_equal(lt@lx[1], 100000)
      expect_equal(length(lt@x), 101)
    })

    test_that(paste("US 2019", nm, ": dxt reproduces published dx"), {
      dx <- sapply(ages, function(a) dxt(lt, a))
      expect_true(max(abs(dx - d$dx[1:100])) <= 1)  # integer rounding of both columns
    })

    test_that(paste("US 2019", nm, ": qxt/pxt reproduce published qx and px"), {
      qx <- sapply(0:89, function(a) qxt(lt, a, 1))
      px <- sapply(0:89, function(a) pxt(lt, a, 1))
      expect_true(max(abs(qx - d$qx[1:90])) < 5e-5)
      expect_equal(px, 1 - qx, tolerance = 1e-12)
      # multi-year survival is the ratio of lx
      expect_equal(pxt(lt, 20, 40), d$lx[61] / d$lx[21], tolerance = 1e-12)
      expect_equal(pxt(lt, 65, 25), d$lx[91] / d$lx[66], tolerance = 1e-12)
      # 20|10 q40 = (l60 - l70)/l40 from the published lx column
      expect_equal(pxt(lt, 40, 20) * qxt(lt, 60, 10),
                   (d$lx[61] - d$lx[71]) / d$lx[41], tolerance = 1e-12)
    })

    test_that(paste("US 2019", nm, ": Lxt reproduces published Lx (ages >= 1)"), {
      L <- sapply(1:99, function(a) Lxt(lt, a))
      expect_true(max(abs(L - d$Lx[2:100])) <= 1)
    })

    test_that(paste("US 2019", nm, ": Lxt age 0 depends on fxt (infant separation factor)"), {
      # package default (fxt = 0.5): L0 = l0 - 0.5 * d0
      expect_equal(Lxt(lt, 0), lt@lx[1] - 0.5 * (lt@lx[1] - lt@lx[2]), tolerance = 1e-12)
      # a0 implied by the published L0: L0 = l1 + a0 * d0  (approx. 0.13 in 2019)
      a0 <- (d$Lx[1] - d$lx[2]) / (d$lx[1] - d$lx[2])
      expect_true(a0 > 0.10 && a0 < 0.15)
      # Lxt uses fxt as the fraction of the year lived by those who die: fxt = 1 - a0
      expect_true(abs(Lxt(lt, 0, fxt = 1 - a0) - d$Lx[1]) <= 1)
    })

    test_that(paste("US 2019", nm, ": mxt = dx / Lx and matches published ratio"), {
      expect_equal(mxt(lt, 40, 1), dxt(lt, 40, 1) / Lxt(lt, 40, 1), tolerance = 1e-12)
      m <- sapply(1:90, function(a) mxt(lt, a, 1))
      expect_equal(m, d$dx[2:91] / d$Lx[2:91], tolerance = 2e-3)
    })

    # tolerance 8: published Tx accumulates unrounded Lx, so it drifts a few units from sum(rounded Lx)
    test_that(paste("US 2019", nm, ": Tx matches published Tx up to the open-interval convention"), {
      for (a in c(1, 5, 20, 40, 65, 80, 90, 99)) {
        expect_true(abs(Tx(lt, a) + tail_gap - d$Tx[a + 1]) <= 8,
                    info = paste("age", a))
      }
      # Tx(0) additionally differs by the age-0 separation factor
      gap0 <- d$Lx[1] - (d$lx[1] - 0.5 * d$dx[1])
      expect_true(abs(Tx(lt, 0) + tail_gap + gap0 - d$Tx[1]) <= 8)
      # differences of Tx give one-year Lx
      expect_equal(Tx(lt, 30) - Tx(lt, 31), Lxt(lt, 30), tolerance = 1e-9)
    })

    test_that(paste("US 2019", nm, ": complete life expectancy vs published ex"), {
      for (a in c(1, 20, 40, 65, 80, 90)) {
        exPkg <- exn(lt, a, type = "complete")
        expect_equal(exPkg, Tx(lt, a) / lt@lx[a + 1], tolerance = 1e-10)
        # correct for the open-interval convention, then compare to the 1-decimal ex
        expect_true(abs(exPkg + tail_gap / d$lx[a + 1] - d$ex[a + 1]) <= 0.051,
                    info = paste("age", a))
      }
    })

    test_that(paste("US 2019", nm, ": curtate expectation equals direct sum of published lx"), {
      for (a in c(0, 25, 50, 75)) {
        direct <- sum(d$lx[(a + 2):101]) / d$lx[a + 1]
        expect_equal(exn(lt, a, type = "curtate"), direct, tolerance = 1e-10)
      }
    })

    test_that(paste("US 2019", nm, ": probs2lifetable rebuilds lx from published qx"), {
      rebuilt <- probs2lifetable(probs = d$qx[1:100], radix = 100000, type = "qx",
                                 name = "rebuilt")
      expect_equal(rebuilt@lx[1], 100000)
      expect_equal(rebuilt@lx[1:96], d$lx[1:96], tolerance = 2e-4)  # qx rounded to 6 digits
    })
  })
}