# library(fmsb) # includes the definition of clifetable() below # missing values have to be omitted before applying clifetable() clifetable <- function(qx) { nc <- length(qx) lx <- numeric(nc) dx <- numeric(nc) Lx <- numeric(nc) lx[1] <- 1e+05 for (i in 1:(nc - 1)) { dx[i] <- lx[i] * qx[i] lx[i + 1] <- lx[i] - dx[i] Lx[i] <- (lx[i] + lx[i + 1])/2 } dx[nc] <- lx[nc] Lx[nc] <- lx[nc]/2 Tx <- rev(cumsum(rev(Lx))) ex <- Tx/lx return(data.frame(qx, lx, dx, Lx, Tx, ex)) } clifetable(Jlife$qx2010M[!is.na(Jlife$qx2010M)]) clifetable(Jlife$qx2010F[!is.na(Jlife$qx2010F)]) # lifetable() function defined below # already included in fmsb package lifetable <- function (mx, ns = NULL, class = 5, mode = 1) { nc <- length(mx) qx <- numeric(nc) if (mode > 10) { mode <- mode%%10 grev <- TRUE } else { grev <- FALSE } if (is.null(ns)) { n <- rep(class, nc) ages <- c(0, cumsum(n)[1:(nc - 1)]) if (mode %in% 4:5) { mode <- mode - 2 } } else { n <- ns ages <- c(0, cumsum(n)[1:(nc - 1)]) if (mode %in% 2:3) { mode <- mode + 2 } } if (mode == 1) { ax <- c(rep(0.5, nc - 1), 1/mx[nc]) } else if (mode == 2) { ax <- c(0.1, rep(0.4, 4), rep(0.5, nc - 6), 1/mx[nc]) } else if (mode == 3) { ax <- c(0.3, rep(0.4, 4), rep(0.5, nc - 6), 1/mx[nc]) } else if (mode == 4) { ax <- c(0.1, 0.4, rep(0.5, nc - 3), 1/mx[nc]) } else if (mode == 5) { ax <- c(0.3, 0.4, rep(0.5, nc - 3), 1/mx[nc]) } else if (mode == 6) { if (mx[1] < 0.107) { ax <- c(0.045 + 2.684 * mx[1], (1.651 - 2.816 * mx[1])/4, rep(0.5, nc - 3), 1/mx[nc]) } else { ax <- c(0.33, 1.352/4, rep(0.5, nc - 3), 1/mx[nc]) } } else if (mode == 7) { if (mx[1] < 0.107) { ax <- c(0.053 + 2.8 * mx[1], (1.522 - 1.518 * mx[1])/4, rep(0.5, nc - 3), 1/mx[nc]) } else { ax <- c(0.35, 1.361/4, rep(0.5, nc - 3), 1/mx[nc]) } } else { ax <- c(rep(0.5, nc - 1), 1/mx[nc]) } if (!grev) { qx <- n * mx/(1 + n * (1 - ax) * mx) qx[nc] <- 1 } else { for (i in 1:(nc - 1)) { qx[i] <- mx[i]/(1/n[1] + mx[i] * (1/2 + n[i]/12 * (mx[i] - (log(mx[i + 1]) - log(mx[i]))/n[i])))/n[i] } qx[nc] <- 1} px <- dx <- lx <- Lx <- numeric(nc) lx[1] <- 1e+05 px <- 1 - qx for (i in 1:(nc - 1)) { dx[i] <- lx[i] * qx[i] lx[i + 1] <- lx[i] - dx[i] Lx[i] <- n[i] * (lx[i + 1] + ax[i] * dx[i]) } dx[nc] <- lx[nc] Lx[nc] <- lx[nc]/mx[nc] Tx <- rev(cumsum(rev(Lx))) ex <- Tx/lx return(data.frame(ages, n, ax, mx, qx, px, lx, dx, Lx, Tx, ex)) } # how to use lifetable() source("http://minato.sip21c.org/ldaR/mortality.R", encoding="CP932") options(digits=3) # avoiding too detailed digits lifetable(S60ASMR, class=5, mode=2)