R/extra_mortality.R
lt_rule_m_extrapolate.RdExtrapolate old-age human mortality curve using mortality laws
lt_rule_m_extrapolate( mx, x, x_fit = x, x_extr, law = "kannisto", opt.method = "LF2", ... )
| mx | Vector or matrix of age specific death-rates. |
|---|---|
| x | Vector of ages at the beginning of the age interval. |
| x_fit | Ages to be considered in estimating the mortality model parameters.
|
| x_extr | Ages for which to extrapolate the death-rates. |
| law | The name of the mortality law/model to be used. The following options are available:
|
| opt.method | character. Default |
| ... | Other arguments to be passed on to the
|
An object of class lt_rule_m_extrapolate with the following components:
List with arguments provided in input. Saved for convenience.
An unevaluated function call, that is, an unevaluated expression that consists of the named function applied to the given arguments.
An object of class MortalityLaw. Here one can find fitted values, residuals, goodness of fit measures etc.
A vector or matrix containing the complete mortality data, that is the modified input data following the extrapolation procedure.
If fitting fails to converge, then we refit assuming Gompertz mortality with explicit starting parameters of parS = c(A = 0.005, B = 0.13) and a warning is issued.
Marius D. Pascariu rpascariu@outlook.com
# Example 1 - abridged data # Age-specific death rates mx <- c(.0859, .0034, .0009, .0007, .0016, .0029, .0036, .0054, .0053, .0146, .0127, .0269, .0170, .0433, .0371, .0784, .0930, .1399, .1875, .2250, .2500, .3000) # Vector of ages x <- c(0, 1, seq(5, 100, by = 5)) names(mx) <- x # Fit the models / Extrapolate the mortality curve x_fit = c(80, 85, 90, 95, 100) x_extr = 90:110 f1 <- lt_rule_m_extrapolate(mx, x, x_fit, x_extr, law = "kannisto") f2 <- lt_rule_m_extrapolate(mx, x, x_fit, x_extr, law = "kannisto_makeham") f3 <- lt_rule_m_extrapolate(mx, x, x_fit, x_extr, law = "gompertz") f4 <- lt_rule_m_extrapolate(mx, x, x_fit, x_extr, law = "ggompertz") f5 <- lt_rule_m_extrapolate(mx, x, x_fit, x_extr, law = "makeham")#> Warning: Extrapolation failed to converge #> Falling back to Gompertz with starting parameters: #> parS = c(A = 0.005, B = 0.13))f6 <- lt_rule_m_extrapolate(mx, x, x_fit, x_extr, law = "beard") f7 <- lt_rule_m_extrapolate(mx, x, x_fit, x_extr, law = "beard_makeham") f8 <- lt_rule_m_extrapolate(mx, x, x_fit, x_extr, law = "quadratic") # Plot the results if (FALSE) { par(mfrow = c(1, 2)) plot(x, mx, pch = 16, xlim = c(60, 110), ylim = c(0, 0.6), cex = 1.5) points(x_fit, mx[paste(x_fit)], pch = 16, col = 4, cex = 1.5) lines(x_extr, f1$values[paste(x_extr)], lty = 1, col = 2, lwd = 2) lines(x_extr, f2$values[paste(x_extr)], lty = 2, col = 3, lwd = 2) lines(x_extr, f3$values[paste(x_extr)], lty = 3, col = 4, lwd = 2) lines(x_extr, f4$values[paste(x_extr)], lty = 4, col = 5, lwd = 2) lines(x_extr, f5$values[paste(x_extr)], lty = 5, col = 6, lwd = 2) lines(x_extr, f6$values[paste(x_extr)], lty = 6, col = 7, lwd = 2) lines(x_extr, f7$values[paste(x_extr)], lty = 7, col = 8, lwd = 2) lines(x_extr, f8$values[paste(x_extr)], lty = 8, col = 9, lwd = 2) legend("topleft", bty = "n", legend = c("Obs. Values", "Obs. Values used in fitting", "Kannisto", "Kannisto-Makeham", "Gompertz", "Gamma-Gompertz", "Makeham", "Beard", "Beard-Makeham", "Quadratic"), lty = c(NA, NA, 1:8), pch = c(16, 16, rep(NA, 8)), col = c(1, 4, 2:9), lwd = 2, pt.cex = 2) } # ---------------------------------------------- # Example 2 - 1-year age data # Age-specific death rates mx1 <- c(.0070, .0082, .0091, .0096, .0108, .0122, .0141, .0150, .0165, .0186, .0205, .0229, .0259, .0294, .0334, .0379, .0426, .0482, .0550, .0628, .0716, .0806, .0897, .1003, .1149, .1264, .1558, .1563, .1812, .2084, .2298, .2536, .2813, .3143, .3352, .3651, .4128) # Vector of ages x1 <- 65:101 names(mx1) <- x1 # Fit the models / Extrapolate the mortality curve x_fit = 80:95 x_extr = 80:125 g1 <- lt_rule_m_extrapolate(mx1, x1, x_fit, x_extr, law = "kannisto") g2 <- lt_rule_m_extrapolate(mx1, x1, x_fit, x_extr, law = "kannisto_makeham") g3 <- lt_rule_m_extrapolate(mx1, x1, x_fit, x_extr, law = "gompertz") g4 <- lt_rule_m_extrapolate(mx1, x1, x_fit, x_extr, law = "ggompertz") g5 <- lt_rule_m_extrapolate(mx1, x1, x_fit, x_extr, law = "makeham") g6 <- lt_rule_m_extrapolate(mx1, x1, x_fit, x_extr, law = "beard") g7 <- lt_rule_m_extrapolate(mx1, x1, x_fit, x_extr, law = "beard_makeham")#> Warning: Extrapolation failed to converge #> Falling back to Gompertz with starting parameters: #> parS = c(A = 0.005, B = 0.13))g8 <- lt_rule_m_extrapolate(mx1, x1, x_fit, x_extr, law = "quadratic") # Plot if (FALSE) { plot(x1, mx1, log = "y", ylim = c(0.001, 5), pch = 16, xlim = c(65, 125), cex = 1.3) points(x_fit, mx1[paste(x_fit)], pch = 16, col = 4, cex = 1.5) lines(x_extr, g1$values[paste(x_extr)], lty = 1, col = 2, lwd = 2) lines(x_extr, g2$values[paste(x_extr)], lty = 2, col = 3, lwd = 2) lines(x_extr, g3$values[paste(x_extr)], lty = 3, col = 4, lwd = 2) lines(x_extr, g4$values[paste(x_extr)], lty = 4, col = 5, lwd = 2) lines(x_extr, g5$values[paste(x_extr)], lty = 5, col = 6, lwd = 2) lines(x_extr, g6$values[paste(x_extr)], lty = 6, col = 7, lwd = 2) lines(x_extr, g7$values[paste(x_extr)], lty = 7, col = 8, lwd = 2) lines(x_extr, g8$values[paste(x_extr)], lty = 8, col = 9, lwd = 2) legend("topleft", bty = "n", legend = c("Obs. Values", "Obs. Values used in fitting", "Kannisto", "Kannisto-Makeham", "Gompertz", "Gamma-Gompertz", "Makeham", "Beard", "Beard-Makeham", "Quadratic"), lty = c(NA, NA, 1:8), pch = c(16, 16, rep(NA, 8)), col = c(1, 4, 2:9), lwd = 2, pt.cex = 2) } # ---------------------------------------------- # Example 3 - Extrapolate mortality for multiple years at once # Create some data mx_matrix <- matrix(rep(mx1, 3), ncol = 3) %*% diag(c(1, 1.05, 1.1)) dimnames(mx_matrix) <- list(age = x1, year = c("year1", "year2", "year3")) # TR: temporary warning suppression until case handling is fixed # in MortalityLaws package F1 <- suppressWarnings(lt_rule_m_extrapolate(mx_matrix, x = x1, x_fit, x_extr, law = "kannisto")) F1#> $input #> $input$opt.choices #> [1] "poissonL" "LF2" "LF1" "LF3" "LF4" "LF5" #> [7] "LF6" "binomialL" #> #> $input$all_the_laws_we_care_about #> [1] "kannisto" "kannisto_makeham" "makeham" "gompertz" #> [5] "ggompertz" "beard" "beard_makeham" "quadratic" #> #> $input$dm #> [1] 37 3 #> #> $input$mx #> year #> age year1 year2 year3 #> 65 0.0070 0.007350 0.00770 #> 66 0.0082 0.008610 0.00902 #> 67 0.0091 0.009555 0.01001 #> 68 0.0096 0.010080 0.01056 #> 69 0.0108 0.011340 0.01188 #> 70 0.0122 0.012810 0.01342 #> 71 0.0141 0.014805 0.01551 #> 72 0.0150 0.015750 0.01650 #> 73 0.0165 0.017325 0.01815 #> 74 0.0186 0.019530 0.02046 #> 75 0.0205 0.021525 0.02255 #> 76 0.0229 0.024045 0.02519 #> 77 0.0259 0.027195 0.02849 #> 78 0.0294 0.030870 0.03234 #> 79 0.0334 0.035070 0.03674 #> 80 0.0379 0.039795 0.04169 #> 81 0.0426 0.044730 0.04686 #> 82 0.0482 0.050610 0.05302 #> 83 0.0550 0.057750 0.06050 #> 84 0.0628 0.065940 0.06908 #> 85 0.0716 0.075180 0.07876 #> 86 0.0806 0.084630 0.08866 #> 87 0.0897 0.094185 0.09867 #> 88 0.1003 0.105315 0.11033 #> 89 0.1149 0.120645 0.12639 #> 90 0.1264 0.132720 0.13904 #> 91 0.1558 0.163590 0.17138 #> 92 0.1563 0.164115 0.17193 #> 93 0.1812 0.190260 0.19932 #> 94 0.2084 0.218820 0.22924 #> 95 0.2298 0.241290 0.25278 #> 96 0.2536 0.266280 0.27896 #> 97 0.2813 0.295365 0.30943 #> 98 0.3143 0.330015 0.34573 #> 99 0.3352 0.351960 0.36872 #> 100 0.3651 0.383355 0.40161 #> 101 0.4128 0.433440 0.45408 #> #> $input$x #> [1] 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 #> [20] 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 #> #> $input$x_fit #> [1] 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 #> #> $input$x_extr #> [1] 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 #> [20] 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 #> [39] 118 119 120 121 122 123 124 125 #> #> $input$law #> [1] "kannisto" #> #> $input$opt.method #> [1] "LF2" #> #> #> $call #> lt_rule_m_extrapolate(mx = mx_matrix, x = x1, x_fit = x_fit, #> x_extr = x_extr, law = "kannisto") #> #> $fitted.model #> Kannisto model: mu[x] = A exp(Bx) / [1 + A exp(Bx)] #> Fitted values: mx #> #> $values #> year1 year2 year3 #> 65 0.00700000 0.00735000 0.00770000 #> 66 0.00820000 0.00861000 0.00902000 #> 67 0.00910000 0.00955500 0.01001000 #> 68 0.00960000 0.01008000 0.01056000 #> 69 0.01080000 0.01134000 0.01188000 #> 70 0.01220000 0.01281000 0.01342000 #> 71 0.01410000 0.01480500 0.01551000 #> 72 0.01500000 0.01575000 0.01650000 #> 73 0.01650000 0.01732500 0.01815000 #> 74 0.01860000 0.01953000 0.02046000 #> 75 0.02050000 0.02152500 0.02255000 #> 76 0.02290000 0.02404500 0.02519000 #> 77 0.02590000 0.02719500 0.02849000 #> 78 0.02940000 0.03087000 0.03234000 #> 79 0.03340000 0.03507000 0.03674000 #> 80 0.03745784 0.03927167 0.04109451 #> 81 0.04262432 0.04471125 0.04680677 #> 82 0.04846753 0.05086438 0.05326894 #> 83 0.05506568 0.05781305 0.06056660 #> 84 0.06250308 0.06564534 0.06879136 #> 85 0.07086966 0.07445486 0.07804023 #> 86 0.08026030 0.08433990 0.08841454 #> 87 0.09077369 0.09540206 0.10001835 #> 88 0.10251075 0.10774441 0.11295636 #> 89 0.11557251 0.12146911 0.12733121 #> 90 0.13005740 0.13667429 0.14324000 #> 91 0.14605789 0.15345036 0.16077026 #> 92 0.16365656 0.17187565 0.17999521 #> 93 0.18292148 0.19201153 0.20096857 #> 94 0.20390129 0.21389709 0.22371903 #> 95 0.22661997 0.23754380 0.24824479 #> 96 0.25107163 0.26293033 0.27450847 #> 97 0.27721589 0.28999802 0.30243293 #> 98 0.30497395 0.31864757 0.33189849 #> 99 0.33422622 0.34873727 0.36274199 #> 100 0.36481154 0.38008331 0.39475819 #> 101 0.39652877 0.41246246 0.42770363 #> 102 0.42914050 0.44561716 0.46130296 #> 103 0.46237932 0.47926300 0.49525760 #> 104 0.49595593 0.51309802 0.52925604 #> 105 0.52956907 0.54681342 0.56298517 #> 106 0.56291612 0.58010475 0.59614170 #> 107 0.59570398 0.61268274 0.62844287 #> 108 0.62765916 0.64428304 0.65963558 #> 109 0.65853639 0.67467413 0.68950343 #> 110 0.68812536 0.70366298 0.71787131 #> 111 0.71625509 0.73109825 0.74460740 #> 112 0.74279598 0.75687115 0.76962293 #> 113 0.76765960 0.78091408 0.79286980 #> 114 0.79079662 0.80319764 0.81433678 #> 115 0.81219328 0.82372632 0.83404460 #> 116 0.83186683 0.84253343 0.85204056 #> 117 0.84986047 0.85967580 0.86839299 #> 118 0.86623806 0.87522835 0.88318594 #> 119 0.88107914 0.88927917 0.89651427 #> 120 0.89447417 0.90192500 0.90847944 #> 121 0.90652048 0.91326738 0.91918589 #> 122 0.91731876 0.92340939 0.92873810 #> 123 0.92697017 0.93245318 0.93723833 #> 124 0.93557409 0.94049789 0.94478492 #> 125 0.94322643 0.94763829 0.95147108 #> #> attr(,"class") #> [1] "lt_rule_m_extrapolate"#> [1] "call" "fitted.model" "input" "values"#> NULL