Nothing
# -------------------------------------------------------------- #
# Author: Marius D. PASCARIU
# Last Update: Thu Jul 20 21:27:19 2023
# -------------------------------------------------------------- #
#' Compute Life Tables from Mortality Data
#'
#' Construct either a full or abridged life table with various input choices
#' like: death counts and mid-interval population estimates \code{(Dx, Ex)} or
#' age-specific death rates \code{(mx)} or death probabilities \code{(qx)}
#' or survivorship curve \code{(lx)} or a distribution of deaths \code{(dx)}.
#' If one of these options is specified, the other can be ignored. The input
#' data can be an object of class: numerical \code{vector}, \code{matrix} or
#' \code{data.frame}.
#'
#' @details
#' The "life table" is also called "mortality table" or "actuarial table".
#' This shows, for each age, what the probability is that a person of that
#' age will die before his or her next birthday, the expectation of life across
#' different age ranges or the survivorship of people from a certain population.
#' @usage
#' LifeTable(x, Dx = NULL, Ex = NULL,
#' mx = NULL,
#' qx = NULL,
#' lx = NULL,
#' dx = NULL,
#' sex = NULL,
#' lx0 = 1e5,
#' ax = NULL)
#' @param x Vector of ages at the beginning of the age interval.
#' @param Dx Object containing death counts. An element of the \code{Dx} object
#' represents the number of deaths during the year to persons aged x to x+n.
#' @param Ex Exposure in the period. \code{Ex} can be approximated by the
#' mid-year population aged x to x+n.
#' @param mx Life table death rate in age interval [x, x+n).
#' @param qx Probability of dying in age interval [x, x+n).
#' @param lx Probability of survival up until exact age x (if l(0) = 1), or
#' the number of survivors at exact age x, assuming l(0) > 1.
#' @param dx Deaths by life-table population in the age interval [x, x+n).
#' @param sex Sex of the population considered here. Default: \code{NULL}.
#' This argument affects the first two values in the life table ax column.
#' If sex is specified the values are computed based on the Coale-Demeny method
#' and are slightly different for males than for females.
#' Options: \code{NULL, male, female, total}.
#' @param lx0 Radix. Default: 100 000.
#' @param ax Numeric scalar. Subject-time alive in age-interval for those who
#' die in the same interval. If \code{NULL} this will be estimated. A common
#' assumption is \code{ax = 0.5}, i.e. the deaths occur in the middle of
#' the interval. Default: \code{NULL}.
#' @return The output is of the \code{"LifeTable"} class with the components:
#' \item{lt}{Computed life table;}
#' \item{call}{\code{Call} in which all of the specified arguments are
#' specified by their full names;}
#' \item{process_date}{Time stamp.}
#' @seealso
#' \code{\link{LawTable}}
#' \code{\link{convertFx}}
#' @author Marius D. Pascariu
#' @examples
#' # Example 1 --- Full life tables with different inputs ---
#'
#' y <- 1900
#' x <- as.numeric(rownames(ahmd$mx))
#' Dx <- ahmd$Dx[, paste(y)]
#' Ex <- ahmd$Ex[, paste(y)]
#'
#' LT1 <- LifeTable(x, Dx = Dx, Ex = Ex)
#' LT2 <- LifeTable(x, mx = LT1$lt$mx)
#' LT3 <- LifeTable(x, qx = LT1$lt$qx)
#' LT4 <- LifeTable(x, lx = LT1$lt$lx)
#' LT5 <- LifeTable(x, dx = LT1$lt$dx)
#'
#' LT1
#' LT5
#' ls(LT5)
#'
#' # Example 2 --- Compute multiple life tables at once ---
#'
#' LTs = LifeTable(x, mx = ahmd$mx)
#' LTs
#' # A warning is printed if the input contains missing values.
#' # Some of the missing values can be handled by the function.
#'
#' # Example 3 --- Abridged life table ------------
#'
#' x <- c(0, 1, seq(5, 110, by = 5))
#' mx <- c(.053, .005, .001, .0012, .0018, .002, .003, .004,
#' .004, .005, .006, .0093, .0129, .019, .031, .049,
#' .084, .129, .180, .2354, .3085, .390, .478, .551)
#' LT6 <- LifeTable(x, mx = mx, sex = "female")
#' LT6
#'
#' # Example 4 --- Abridged life table w using my own 'ax' ------------
#' # In this examples we are using the ages (x) and death rates (mx) from
#' # example 3. Note that 'ax' must have the same length as the 'x' vector
#' # otherwise an error message will be returned.
#'
#' my_ax <- c(0.1, 1.5, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
#' 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1)
#'
#' LT7 <- LifeTable(x = x, mx = mx, ax = my_ax)
#'
#' @export
LifeTable <- function(x,
Dx = NULL,
Ex = NULL,
mx = NULL,
qx = NULL,
lx = NULL,
dx = NULL,
sex = NULL,
lx0 = 1e5,
ax = NULL){
input <- c(as.list(environment()))
X <- LifeTable.check(input)
LT <- NULL
if (any(X$iclass == "numeric")) {
LT <- with(X, LifeTable.core(x, Dx, Ex, mx, qx, lx, dx, sex, lx0, ax))
} else {
for (i in 1:X$nLT) {
LTi <- with(
X,
LifeTable.core(
x,
Dx = Dx[, i],
Ex = Ex[, i],
mx = mx[, i],
qx = qx[, i],
lx = lx[, i],
dx = dx[, i],
sex = sex,
lx0 = lx0,
ax = ax
)
)
N <- X$LTnames
LTn <- if (is.na(N[i])) i else N[i]
LTi <- cbind(LT = LTn, LTi)
LT <- rbind(LT, LTi)
}
}
# Exit
out <- list(
lt = LT,
call = match.call(),
process_date = date()
)
out <- structure(class = "LifeTable", out)
return(out)
}
#' LifeTable.core
#' @inheritParams LifeTable
#' @return A data.frame containing life table results
#' @keywords internal
LifeTable.core <- function(x, Dx, Ex, mx, qx, lx, dx, sex, lx0, ax){
my.case <- find.my.case(Dx, Ex, mx, qx, lx, dx)$case
gr_names <- paste0("[", x,",", c(x[-1], "+"), ")")
N <- length(x)
df <- diff(x)
nx <- c(df, df[N - 1])
if (my.case == "C1_DxEx") {
Dx <- as.numeric(Dx)
Ex <- as.numeric(Ex)
mx <- Dx/Ex
mx <- uxAbove100(x, mx)
qx <- mx_qx(x, nx, mx, out = "qx")
lx <- lx0 * c(1, cumprod(1 - qx)[1:(N - 1)])
dx <- dx_lx(lx, out = "dx")
}
if (my.case == "C2_mx") {
mx <- as.numeric(mx)
qx <- mx_qx(x, nx, mx, out = "qx")
lx <- lx0 * c(1, cumprod(1 - qx)[1:(N - 1)])
dx <- dx_lx(lx, out = "dx")
}
if (my.case == "C3_qx") {
qx <- as.numeric(qx)
mx <- mx_qx(x, nx, qx, out = "mx")
lx <- lx0 * c(1, cumprod(1 - qx)[1:(N - 1)])
dx <- dx_lx(lx, out = "dx")
}
if (my.case == "C4_lx") {
lx <- as.numeric(lx)
lx <- lx * lx0/lx[1]
dx <- dx_lx(lx, out = "dx")
qx <- dx/lx
qx <- uxAbove100(x, qx)
mx <- mx_qx(x, nx, qx, out = "mx")
}
if (my.case == "C5_dx") {
dx <- as.numeric(dx)
dx <- dx * lx0/sum(dx)
lx <- dx_lx(dx, out = "lx")
qx <- dx/lx
qx <- uxAbove100(x, qx)
mx <- mx_qx(x, nx, qx, out = "mx")
}
if (is.null(ax)) {
ax <- compute.ax(x, mx, qx)
if (!is.null(sex)) ax <- coale.demeny.ax(x, mx, ax, sex)
} else if (length(ax) == 1){
ax <- rep(ax, N)
} else if (length(ax) == N){
ax <- as.numeric(ax)
}
Lx <- nx * lx - (nx - ax) * dx
Lx[N] <- ax[N] * dx[N]
Lx[is.na(Lx)] <- 0
Tx <- rev(cumsum(rev(Lx)))
ex <- Tx/lx
ex[is.na(ex)] <- 0
ex[N] <- if (ex[N - 1] == 0) 0 else ax[N]
last_check = all(is.na(mx)) | all(is.nan(mx)) | all(is.infinite(mx)) | all(mx == 0)
if (last_check) mx = qx = ax = lx = dx = Lx = Tx = ex <- NA
out <- data.frame(x.int = gr_names,
x = x,
mx = mx,
qx = qx,
ax = ax,
lx = lx,
dx = dx,
Lx = Lx,
Tx = Tx,
ex = ex)
return(out)
}
#' Function that identifies the case/problem we have to solve
#' @inheritParams LifeTable
#' @return A list containing problem solving details
#' @keywords internal
find.my.case <- function(Dx = NULL,
Ex = NULL,
mx = NULL,
qx = NULL,
lx = NULL,
dx = NULL) {
input <- c(as.list(environment()))
# Matrix of possible cases --------------------
rn <- c("C1_DxEx", "C2_mx", "C3_qx", "C4_lx", "C5_dx")
cn <- c("Dx", "Ex", "mx", "qx", "lx", "dx")
mat <- matrix(
ncol = 6,
byrow = TRUE,
dimnames = list(rn, cn),
data = c(T,T,F,F,F,F,
F,F,T,F,F,F,
F,F,F,T,F,F,
F,F,F,F,T,F,
F,F,F,F,F,T)
)
# ----------------------------------------------
L1 <- !unlist(lapply(input, is.null))
L2 <- apply(mat, 1, function(x) all(L1 == x))
my_case <- rn[L2]
if (sum(L1[c(1, 2)]) == 1) {
stop("If you input 'Dx' you must input 'Ex' as well, and viceversa",
call. = FALSE)
}
if (!any(L2)) {
stop("The input is not specified correctly. Check again the function ",
"arguments and make sure the input data is added properly.",
call. = FALSE)
}
X <- input[L1][[1]]
if (length(dim(X)) == 1){
# TR changed here:
# The following input isn't detected with is.vector()
# X = structure(c(0.0036542739116619, 0.000150960486092765, 2.77881983521598e-05,
# 0.000136941279579316, 0.00018946827083136, 0.00026606712873658,
# 0.000258838755220895, 0.000557913403869528, 0.000665917509468515,
# 0.00114979916336982, 0.0021040113433655, 0.00384681333263752,
# 0.00632225868307671, 0.00937214337262448, 0.0154497258185624,
# 0.023117866694913, 0.0363617070194365, 0.0609184108563059, 0.120398389987367,
# 0.221461187214612, 0.420152946468736), .Dim = 21L)
# nLT would come as NA
# and iclass would be array
X <- c(X)
}
nLT <- 1
LTnames <- NA
# TR: change from !is.vector
if (length(dim(X)) == 2 ) {
nLT <- ncol(X) # number of LTs to be created
LTnames <- colnames(X) # the names to be assigned to LTs
}
out <- list(case = my_case,
iclass = class(X), # TR: if inputs are matrix,
# then this is two elements
nLT = nLT,
LTnames = LTnames)
return(out)
}
#' mx to qx
#'
#' Function to convert mx into qx and back, using the constant force of
#' mortality assumption (CFM).
#' @inheritParams LifeTable
#' @param nx Length of the age-intervals.
#' @param ux A vector of mx or qx.
#' @param out Type of the output: mx or qx.
#' @return A vector of rates
#' @keywords internal
mx_qx <- function(x, nx, ux, out = c("qx", "mx")){
out <- match.arg(out)
if (out == "qx") {
eta <- 1 - exp(-nx * ux)
eta[length(nx)] <- 1 # The life table should always close with q[x] = 1
} else {
eta <- suppressWarnings(-log(1 - ux)/nx)
# If qx[last-age] = 1 then mx[last-age] = Inf. Not nice to have Inf's;
# they distort the results in the subsequent processes.
# We apply a simple extrapolation method of the last mx.
N <- length(x)
eta[N] <- eta[N - 1]^2 / eta[N - 2]
}
eta <- uxAbove100(x, eta)
return(eta)
}
#' Educate mx or qx on how to behave above age 100 if it gets in trouble
#' (with NA's, zero's and Inf)
#' @inheritParams LifeTable
#' @inheritParams mx_qx
#' @param omega Threshold age. Default: 100.
#' @param verbose A logical value. Set \code{verbose = FALSE} to silent
#' the process that take place inside the function and avoid progress messages.
#' @return A vector of rates
#' @keywords internal
uxAbove100 <- function(x,
ux,
omega = 100,
verbose = FALSE) {
if (is.vector(ux)) {
L <- x >= 100 & (is.na(ux) | is.infinite(ux) | ux == 0)
if (any(L)) {
mux <- max(ux[!L])
ux[L] <- mux
if (verbose)
warning("The input data contains NA's, Inf or zero's over the age of ",
"100. These have been replaced with maximum observed value: ",
round(mux, 4), call. = FALSE)
}
} else {
for (i in 1:ncol(ux)) {
ux[, i] = uxAbove100(x, ux[, i], omega, verbose)
}
}
return(ux)
}
#' dx to lx
#'
#' Function to convert dx into lx and back
#' @param ux A vector of dx or lx data.
#' @param out Type of the output: dx or lx.
#' @return A vector containing dx or lx values
#' @keywords internal
dx_lx <- function(ux, out = c("dx", "lx")) {
out <- match.arg(out)
if (out == "dx") {
ux_ <- rev(diff(rev(ux)))
d <- ux[1] - sum(ux_)
eta <- c(ux_, d)
} else {
eta <- rev(cumsum(rev(ux)))
}
return(eta)
}
#' Find ax indicator
#'
#' @inheritParams LifeTable
#' @return \code{ax} - the point in the age internal where 50% of the deaths
#' have already occurred
#' @keywords internal
compute.ax <- function(x, mx, qx) {
nx <- c(diff(x), Inf)
N <- length(x)
ax <- nx + 1/mx - nx/qx
for (i in 1:(N - 1)) {
if (is.infinite(ax[i + 1]) | is.na(ax[i + 1])) ax[i + 1] = ax[i]
}
return(ax)
}
#' Find ax[1:2] indicators using Coale-Demeny coefficients
#' Here we adjust the first two values of ax to account for infant
#' mortality more accurately
#' @inheritParams LifeTable
#' @return A vector of coefficients
#' @keywords internal
coale.demeny.ax <- function(x, mx, ax, sex) {
if (mx[1] < 0) stop("'m[1]' must be greater than 0", call. = FALSE)
nx <- c(diff(x), Inf)
m0 <- mx[1]
a0M <- ifelse(m0 >= 0.107, 0.330, 0.045 + 2.684 * m0)
a1M <- ifelse(m0 >= 0.107, 0.330, 1.651 - 2.816 * m0)
a0F <- ifelse(m0 >= 0.107, 0.350, 0.053 + 2.800 * m0)
a1F <- ifelse(m0 >= 0.107, 0.350, 1.522 - 1.518 * m0)
a0T <- (a0M + a0F)/2
a1T <- (a1M + a1F)/2
f <- nx[1:2] / c(1, 4)
if (sex == "male") ax[1:2] <- c(a0M, a1M) * f
if (sex == "female") ax[1:2] <- c(a0F, a1F) * f
if (sex == "total") ax[1:2] <- c(a0T, a1T) * f
return(ax)
}
#' Check LifeTable input
#' @param input A list containing the input arguments of the LifeTable functions.
#' @return A list of life table validated data
#' @keywords internal
LifeTable.check <- function(input) {
with(input, {
# ----------------------------------------------
K <- find.my.case(Dx, Ex, mx, qx, lx, dx)
C <- K$case
valid_classes <- c("numeric", "matrix", "data.frame", NULL)
if (!any(K$iclass %in% valid_classes)) {
stop(paste0("The class of the input should be: ",
paste(valid_classes, collapse = ", ")), call. = FALSE)
}
# ----------------------------------------------
SMS <- "contains missing values. These have been replaced with "
if (!is.null(sex)) {
if (!any(sex %in% c("male", "female", "total")))
stop("'sex' should be: 'male', 'female', 'total' or 'NULL'.",
call. = FALSE)
}
if (C == "C1_DxEx") {
if (any(is.na(Dx))) warning("'Dx'", SMS, 0, call. = FALSE)
if (any(is.na(Ex))) warning("'Ex'", SMS, 0.01, call. = FALSE)
Dx[is.na(Dx)] <- 0
Ex[is.na(Ex) | Ex == 0] <- 0.01
}
if (C == "C2_mx") {
mx <- uxAbove100(x, mx)
}
if (C == "C3_qx") {
qx <- uxAbove100(x, qx)
}
if (C == "C4_lx") {
if (any(is.na(lx))) warning("'lx'", SMS, 0, call. = FALSE)
lx[is.na(lx) & x >= 100] <- 0
}
if (C == "C5_dx") {
if (any(is.na(dx))) warning("'dx'", SMS, 0, call. = FALSE)
dx[is.na(dx)] <- 0
}
if (!is.null(ax)) {
if (!is.numeric(ax))
stop("'ax' must be a numeric scalar (or NULL)", call. = FALSE)
if (!any(length(ax) %in% c(1, length(x))))
stop("'ax' must be a scalar of lenght 1 or a ",
"vector of the same dimension as 'x'",
call. = FALSE)
}
# Exit
out <- list(x = x,
Dx = Dx,
Ex = Ex,
mx = mx,
qx = qx,
lx = lx,
dx = dx,
sex = sex,
lx0 = lx0,
ax = ax,
iclass = K$iclass,
nLT = K$nLT,
LTnames = K$LTnames)
return(out)
})
}
#' Print LifeTable
#' @param x An object of class \code{"LifeTable"}
#' @param ... Further arguments passed to or from other methods.
#' @return Print data on the console
#' @keywords internal
#' @export
print.LifeTable <- function(x, ...){
LT <- x$lt
lt <- with(
LT,
data.frame(
x.int = x.int,
x = x,
mx = round(mx, 6),
qx = round(qx, 6),
ax = round(ax, 2),
lx = round(lx),
dx = round(dx),
Lx = round(Lx),
Tx = round(Tx),
ex = round(ex, 2)
)
)
if (colnames(LT)[1] == "LT") lt <- data.frame(LT = LT$LT, lt)
dimnames(lt) <- dimnames(LT)
nx <- length(unique(LT$x))
nlt <- nrow(LT) / nx
out <- head_tail(lt, hlength = 6, tlength = 3, ...)
step <- diff(LT$x)
step <- step[step > 0]
type1 <- if (all(step == 1)) "Full" else "Abridged"
type2 <- if (nlt == 1) "Life Table" else "Life Tables"
cat("\n", type1, " ", type2, "\n\n", sep = "")
cat("Number of life tables:", nlt, "\n")
cat("Dimension:", nrow(LT), "x", ncol(LT), "\n")
cat("Age intervals:", head_tail(lt$x.int, hlength = 3, tlength = 3), "\n\n")
print(out, row.names = FALSE)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.