library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.4 v dplyr 1.0.7
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.0.2 v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(HMDHFDplus)
# select data and download from HMD ----
# HMD username and password must be in session memory as 'un' and 'pw'
sel_country = 'USA'
sel_year = 1980
HMDdata = readHMDweb(sel_country,item='fltper_1x1', username=un, password = pw) %>%
filter(Year %in% sel_year) %>%
select(Year, Age, mx, dx, lx, ex,qx)
regr_ages = 70:89 # ages from which to estimate extrapolation coeffs
extrap_ages = 90:120 # ages to which we want to extrapolate mx
# mx values for regression
mx_regr = HMDdata %>%
filter(Age %in% regr_ages) %>%
pull(mx)
logit = function(p) log( p/(1-p) )
Yx = logit( mx_regr )
ab = coef( lm(Yx ~ I(regr_ages))) # regress logit(mx) on x for regression ages
ahat = exp(ab[1]) # ahat = e^(intercept)
bhat = ab[2] # bhat = slope
# Kannisto extrapolation is m(x) = a*e^(bx) / (1+a*e^(bx))
num = ahat * exp(bhat * extrap_ages)
mx_extrap = num/(1+num)
mx_hat = c(mx_regr, mx_extrap)
# illustrate results when extrapolating from (70:89) to (90:120)
ggplot() +
aes(x = 0:89, y= log10(HMDdata$mx[1+0:89])) +
geom_point() +
theme_bw() +
labs(x='Age', y='log10[mx]', title='Kannisto Extrapolation') +
geom_polygon(aes(x=c(70,70,90,90), y=c(-4,0,0,-4)), fill='red', alpha=.10) +
geom_line(aes(x=c(regr_ages,extrap_ages),
y=log10(mx_hat)), color='red', lwd=1.5, alpha=.60)
ggsave(filename='Kannisto.png')
## Saving 7 x 5 in image