amvillegas / StMoMo

Stochastic Mortality Modelling
21 stars 6 forks source link

How to create Confidence Interval for bootstrapped parameters #28

Closed andersmaxDK closed 4 years ago

andersmaxDK commented 4 years ago

Hi Andres,

I am a happy user of your package, so thank you for this contribution. Is there any easy way to create confidence intervals for the bootstrapped parameters? And furthermore, are there any ways to estimate life expectancy on all the bootstrap samples?

Kind regards, Anders

amvillegas commented 4 years ago

here is no simple way to extract coefficients or rates from the boostrap.fitStMoMo, definitely something we may want to add in the package in one of the next releases.

For the moment you can use a calculation like the one below (which works for the LC model). In your case, it depends on the model you used, but once you have the parameters you can just reconstruct the rates.

library("StMoMo")

LCfit <- fit(lc(), data = EWMaleData)

#number of bootsrap samples
nBoot <- 5

LCResBoot <- bootstrap(LCfit, nBoot = nBoot, type = "residual")

#just for the first bootstrap sample
ax1 <- LCResBoot$bootParameters[[1]]$ax
bx1 <- as.numeric(LCResBoot$bootParameters[[1]]$bx)
kt1 <- as.numeric(LCResBoot$bootParameters[[1]]$kt)
mxt1 <- exp(ax1 + outer(bx1, kt1))

#for all bootstrap samples
axBoot <- sapply(LCResBoot$bootParameters, function(x) x$ax)
bxBoot <- sapply(LCResBoot$bootParameters, function(x) x$bx)
ktBoot <- sapply(LCResBoot$bootParameters, function(x) x$kt)

mxtBoot <- array(data = NA, dim = c(nrow(axBoot), nrow(ktBoot), nBoot))

for (i in 1:nBoot)
{
mxtBoot[, , i] <- exp(axBoot[, i] + outer(bxBoot[, i], ktBoot[, i]))
}

#age 65
matplot(mxtBoot[66, , ], type = "l", col = "black")