R/LC.R
In PPgp/MortCast: Estimation and Projection of Age-Specific Mortality Rates
Defines functions
.kranges
ultimate.bx
rotate.leecarter
lileecarter.estimate
.finish.bx
bx.estimate
leecarter.estimate
Documented in
leecarter.estimate
lileecarter.estimate
rotate.leecarter
ultimate.bx
#' @title Lee-Carter Estimation
#' @description Estimate Lee-Carter parameters (Lee and Carter 1992).
#' @details
#' The function estimates parameters of \eqn{\log(m_x(t)) = a_x + b_x k(t) + \epsilon_x(t)} (Lee and Carter 1992).
#' The argument \code{ax.index} determines which time periods to use to
#' estimate the \eqn{a_x} parameter, while \code{ax.smooth} controls if
#' the resulting \eqn{a_x} should be smoothened over ages (see Sevcikova et al. 2016 for details).
#'
#' @param mx A matrix of age-specific mortality rates where rows correspond to age groups
#' and columns correspond to time periods. Rownames define the starting ages of the age groups.
#' @param ax.index A vector of column indices of \code{mx} to be used to estimate the \eqn{a_x} parameter.
#' By default all time periods are used.
#' @param ax.smooth Logical allowing to smooth the \eqn{a_x} over ages.
#' @param ax.smooth.df Degree of freedom for smoothing if \code{ax.smooth} is \code{TRUE}.
#' Default is half the length of \eqn{a_x}.
#' @param bx.postprocess Logical determining if numerical anomalies in \eqn{b_x} should be dealt with.
#' @param nx Size of age groups. By default ages are determined by rownames of \code{mx}. This argument is only used if
#' \code{mx} has no rownames. If \code{nx} is 5, the age groups are interpreted as 0, 1, 5, 10, \dots. For \code{nx} equals 1,
#' the age groups are interpreted as 0, 1, 2, 3, \dots.
#' @return List with elements \code{ax}, \code{bx} and \code{kt} corresponding to the estimated parameters.
#' @export
#' @seealso \code{\link{mortcast}}, \code{\link{lileecarter.estimate}}
#' @references
#' Lee, R. D. and Carter, L. (1992). Modeling and forecasting the time series of
#' US mortality. Journal of the American Statistical Association, 87, 659-671.
#'
#' Sevcikova H., Li N., Kantorova V., Gerland P., Raftery A.E. (2016).
#' Age-Specific Mortality and Fertility Rates for Probabilistic Population Projections.
#' In: Schoen R. (eds) Dynamic Demographic Analysis. The Springer Series on Demographic Methods
#' and Population Analysis, vol 39. Springer, Cham
#'
#' @examples
#' data(mxM, package = "wpp2017")
#' mx <- subset(mxM, name == "Netherlands")[,4:16]
#' rownames(mx) <- c(0,1, seq(5, 100, by=5))
#' lc.ax.avg <- leecarter.estimate(mx)
#' lc.ax.last <- leecarter.estimate(mx, ax.index=ncol(mx))
#' plot(lc.ax.avg$ax, type="l")
#' lines(lc.ax.last$ax, col="blue")
#'
leecarter.estimate <- function(mx, ax.index = NULL, ax.smooth = FALSE,
ax.smooth.df = NULL, bx.postprocess = TRUE, nx = 5) {
if(length(dim(mx)) < 2 || ncol(mx) < 2) stop("mx must be a matrix (age x time) of at least two columns.")
lmx <- log(mx)
if(length(dim(lmx))==0)
lmx <- as.matrix(lmx)
if(is.null(rownames(lmx))) # default ages
rownames(lmx) <- if(nx > 1) c(0,1, seq(nx, by=nx, length=nrow(lmx)-2)) else 0:(nrow(lmx)-1)
nest <- ncol(lmx)
if(is.null(ax.index)) ax.index <- 1:nest
ax <- apply(lmx[,ax.index, drop=FALSE], 1, sum, na.rm=TRUE) / length(ax.index)
if(ax.smooth) {
if (is.null(ax.smooth.df)) ax.smooth.df <- ceiling(length(ax)/2)
ax.sm <- smooth.spline(ax, df = ax.smooth.df)$y
ax[-1] <- ax.sm[-1] # keep value of the first age group
}
kt <- apply(lmx, 2, sum) - sum(ax)
bx <- bx.estimate(lmx, ax, kt, postprocess = bx.postprocess)
return(list(ax = ax, bx = bx, kt = kt, nx = nx))
}
bx.estimate <- function(lmx, ax, kt, postprocess=TRUE) {
x2 <- sum(kt*kt)
x1 <- rep(NA, nrow(lmx))
for (i in 1:nrow(lmx))
x1[i] <- sum((lmx[i,]-ax[i])*kt)
bx <- x1/x2
if(postprocess)
bx <- .finish.bx(bx)
names(bx) <- rownames(lmx)
return(bx)
}
.finish.bx <- function(bx) {
negbx <- which(bx <= 0)
lnegbx <- length(negbx)
if(lnegbx > 0 && negbx[1] == 1) {
bx[1] <- 0
negbx <- if(lnegbx > 1) negbx[2:lnegbx] else c()
lnegbx <- length(negbx)
}
while(lnegbx > 0) {
bx[negbx] <- 0.5 * bx[negbx-1]
negbx <- which(bx < 0)
lnegbx <- length(negbx)
}
lbx <- length(bx)
for (i in 1:(lbx-1)) {
if (bx[lbx+1 - i] == 0) bx[lbx+1 - i] <- bx[lbx - i]
}
bx <- bx/sum(bx) # must sum to 1
return(bx)
}
#' @title Coherent Lee-Carter Estimation
#' @description Estimate coherent Lee-Carter parameters (Li and Lee 2005).
#' @details The coherent Lee-Carter parameters for male and female mortality rates
#' share the same \eqn{b_x} which is the average of the age-specific
#' \eqn{b_x} parameters.
#'
#' The function in addition computes the ultimate \eqn{b^u_x} as defined in
#' Li et al. (2013) based on the coherent \eqn{b_x}.
#' @param mxM A matrix of male age-specific mortality rates where rows correspond to age groups
#' and columns correspond to time periods. For 5-year age groups, the first and second rows corresponds to
#' 0-1 and 1-5 age groups, respectively. Rownames define the starting ages of the respective groups.
#' @param mxF A matrix of female mortality rates of the same shape as \code{mxM}.
#' @param nx Size of age groups. Should be either 5 or 1.
#' @param ... Additional arguments passed to \code{\link{leecarter.estimate}}.
#' @return List containing elements \code{bx} (coherent \eqn{b_x} parameter),
#' \code{ultimate.bx} (ultimate \eqn{b^u_x} parameter), \code{ages} (age groups), \code{nx} (age group interval), and
#' lists \code{female} and \code{male}, each with the Lee-Carter parameters.
#' @export
#'
#' @references
#' Li, N. and Lee, R. D. (2005). Coherent mortality forecasts for a group of populations:
#' An extension of the Lee-Carter method. Demography, 42, 575-594.
#'
#' Li, N., Lee, R. D. and Gerland, P. (2013). Extending the Lee-Carter method to model the rotation
#' of age patterns of mortality decline for long-term projections. Demography, 50, 2037-2051.
#'
#' @examples
#' data(mxM, mxF, package = "wpp2017")
#' country <- "Germany"
#' mxm <- subset(mxM, name == country)[,4:16]
#' mxf <- subset(mxF, name == country)[,4:16]
#' rownames(mxm) <- rownames(mxf) <- c(0,1, seq(5, 100, by=5))
#' lc <- lileecarter.estimate(mxm, mxf)
#' plot(lc$bx, type="l")
#' lines(lc$ultimate.bx, lty=2)
#'
lileecarter.estimate <- function(mxM, mxF, nx = 5, ...) {
lc.male <- leecarter.estimate(mxM, nx = nx, ...)
lc.female <- leecarter.estimate(mxF, nx = nx, ...)
if(length(lc.female$ax) != length(lc.male$ax) || length(lc.female$kt) != length(lc.male$kt))
stop("Mismatch in dimensions of male and female mortality. Check the mxM and mxF arguments.")
bx <- (lc.male$bx + lc.female$bx)/2
ages <- as.integer(names(lc.female$bx))
# correct nx if wrong input
nx <- ages[4] - ages[3]
return(list(bx = bx, ultimate.bx = ultimate.bx(bx), ages = ages, nx = nx,
female=list(ax=lc.female$ax, bx=bx, kt=lc.female$kt, sex.bx = lc.female$bx),
male=list(ax=lc.male$ax, bx=bx, kt=lc.male$kt, sex.bx = lc.male$bx)
))
}
#' @title Rotated Lee-Carter
#' @description Rotate the Lee-Carter parameter \eqn{b_x} over time to reach an ultimate \eqn{b^u_x},
#' as described in Li et al. (2013).
#'
#' @param bx A vector of the Lee-Carter \eqn{b_x} parameter
#' (from e.g. \code{\link{lileecarter.estimate}} or \code{\link{leecarter.estimate}}).
#' @param ultimate.bx A vector of the ultimate \eqn{b^u_x} parameter as defined in Li, Lee, Gerland (2013)
#' (obtained using \code{\link{lileecarter.estimate}} or \code{\link{ultimate.bx}}).
#' @param e0 A time series of life expectancies.
#' @param e0l Level of life expectancy at which the rotation starts.
#' @param e0u Level of life expectancy at which the rotation finishes.
#' @param p Exponent of the smooth function.
#' @return Function \code{rotate.leecarter} returns a matrix of rotated \eqn{B_x(t)} where rows correspond to age groups and columns
#' correspond to time periods (given by the vector \code{e0}).
#' @export
#'
#' @references
#' Li, N., Lee, R. D. and Gerland, P. (2013). Extending the Lee-Carter method to model the rotation
#' of age patterns of mortality decline for long-term projections. Demography, 50, 2037-2051.
#'
#' @examples
#' data(mxF, mxM, e0Fproj, e0Mproj, package = "wpp2017")
#' country <- "Japan"
#' mxm <- subset(mxM, name == country)[,4:16]
#' mxf <- subset(mxF, name == country)[,4:16]
#' e0f <- as.numeric(subset(e0Fproj, name == country)[-(1:2)])
#' e0m <- as.numeric(subset(e0Mproj, name == country)[-(1:2)])
#' rownames(mxm) <- rownames(mxf) <- c(0,1, seq(5, 100, by=5))
#' lc <- lileecarter.estimate(mxm, mxf)
#' rotlc <- rotate.leecarter(lc$bx, lc$ultimate.bx, (e0f + e0m)/2)
#' plot(lc$bx, type="l")
#' lines(lc$ultimate.bx, col="red")
#' for(i in 1:ncol(rotlc)) lines(rotlc[,i], col="grey")
#'
rotate.leecarter <- function(bx, ultimate.bx, e0, e0l = 80, e0u = 102, p = 0.5) {
npred <- length(e0)
wt <- (e0 - e0l)/(e0u-e0l)
wst <- (0.5*(1+(sin(pi/2*(2*wt-1)))))^p
Bxt <- matrix(NA, nrow=length(bx), ncol=npred)
for(t in 1:npred) {
Bxt[,t] <- switch(cut(e0[t], c(0, e0l, e0u, 9999), labels=FALSE, right=FALSE),
bx,
(1-wst[t])*bx + wst[t]*ultimate.bx,
ultimate.bx)
}
dimnames(Bxt) <- list(names(bx), names(e0))
return(Bxt)
}
#' @rdname rotate.leecarter
#' @return Function \code{ultimate.bx} returns a vector of the ultimate \eqn{b^u_x}.
#' @export
ultimate.bx <- function(bx) {
# ultimate bx (Li, Lee, Gerland 2013)
bux <- bx
lbx <- length(bx)
idx15 <- which(names(bx) == 15)
idx60 <- which(names(bx) == 60)
bux[1:idx60] <- mean(bux[idx15:idx60])
bux[(idx60+1):lbx] <- bux[(idx60+1):lbx] * (bux[idx60]/bux[idx60+1]) # adjust so that b(65-69)=b(60-64)
return(bux/sum(bux)) # must sum to 1
}
.kranges <- function(Bxt, axm, axf) {
lmin <- -12
lmax <- 0
machine.max <- log(.Machine$double.xmax)
machine.min <- log(.Machine$double.xmin)
npred <- ncol(Bxt)
kranges <- list(male=list(), female=list())
ax <- list(male=axm, female=axf)
bx.lim <- rbind(apply(Bxt, 2, function(x) min(x[x>0])),
apply(Bxt, 2, function(x) max(x[x>0])))
for(sex in c("male", "female")) {
this.lmin <- floor(min(lmin, min(ax[[sex]])-0.1))
kranges[[sex]]$kl <- pmin((lmax - apply(ax[[sex]], 2, max))/bx.lim[2,], machine.max)
kranges[[sex]]$ku <- pmax((this.lmin - apply(ax[[sex]], 2, min))/bx.lim[1,], machine.min)
}
return(kranges)
}
#' @title Coherent Rotated Lee-Carter Prediction
#' @description Predict age-specific mortality rates using the coherent rotated Lee-Carter method.
#' @details This function implements Steps 6-9 of Algorithm 2 in Sevcikova et al. (2016).
#' It uses the abridged or unabridged life table function to find the level of mortality that coresponds to the given
#' life expectancy. Thus, it can be used for both, mortality for 5- or 1-year age groups.
#' @param e0m A time series of future male life expectancy.
#' @param e0f A time series of future female life expectancy.
#' @param lc.pars A list of coherent Lee-Carter parameters with elements \code{bx}, \code{ultimate.bx},
#' \code{ages}, \code{nx},
#' \code{female} and \code{male} as returned by \code{\link{lileecarter.estimate}}.
#' The \code{female} and \code{male} objects are again lists that should contain a vector
#' \code{ax} and optionally a matrix \code{axt} if the \eqn{a_x} parameter
#' needs to be defined as time dependent. In such a case, rows are age groups and columns are
#' time periods corresponding to the length of the \code{e0f} and \code{e0m} vectors.
#' @param rotate If \code{TRUE} the rotation method of \eqn{b_x} is used as described in Li et al. (2013).
#' @param keep.lt Logical. If \code{TRUE} additional life table columns are kept in the
#' resulting object.
#' @param constrain.all.ages By default the method constrains the male mortality to be above female
#' mortality for old ages if the male life expectancy is below the female life expectancy. Setting
#' this argument to \code{TRUE} causes this constraint to be applied to all ages.
#' @param \dots Additional life table arguments.
#' @return List with elements \code{female} and \code{male}, each of which contains a matrix \code{mx}
#' with the predicted mortality rates. If \code{keep.lt} is \code{TRUE}, it also
#' contains matrices \code{sr} (survival rates), and life table quantities \code{Lx} and \code{lx}.
#' @export
#'
#' @seealso \code{\link{rotate.leecarter}}, \code{\link{leecarter.estimate}}, \code{\link{lileecarter.estimate}},
#' \code{\link{mortcast.blend}}
#' @references
#' Li, N., Lee, R. D. and Gerland, P. (2013). Extending the Lee-Carter method to model the rotation
#' of age patterns of mortality decline for long-term projections. Demography, 50, 2037-2051.
#'
#' Sevcikova H., Li N., Kantorova V., Gerland P., Raftery A.E. (2016).
#' Age-Specific Mortality and Fertility Rates for Probabilistic Population Projections.
#' In: Schoen R. (eds) Dynamic Demographic Analysis. The Springer Series on Demographic Methods
#' and Population Analysis, vol 39. Springer, Cham
#'
#' @examples
#' # estimate parameters from historical mortality data (5-year age groups)
#' data(mxM, mxF, e0Fproj, e0Mproj, package = "wpp2017")
#' country <- "Brazil"
#' mxm <- subset(mxM, name == country)[,4:16]
#' mxf <- subset(mxF, name == country)[,4:16]
#' rownames(mxm) <- rownames(mxf) <- c(0,1, seq(5, 100, by=5))
#' lc <- lileecarter.estimate(mxm, mxf)
#'
#' # project into future for given levels of life expectancy
#' e0f <- as.numeric(subset(e0Fproj, name == country)[-(1:2)])
#' e0m <- as.numeric(subset(e0Mproj, name == country)[-(1:2)])
#' pred <- mortcast(e0m, e0f, lc)
#'
#' # plot first projection in black and the remaining ones in grey
#' plot(lc$ages, pred$female$mx[,1], type="b", log="y", ylim=range(pred$female$mx),
#' ylab="female mx", xlab="Age", main=paste(country, "(5-year age groups)"), cex=0.5)
#' for(i in 2:ncol(pred$female$mx)) lines(lc$ages, pred$female$mx[,i], col="grey")
#'
#' # similarly for 1-year age groups
#' # derive toy 1-year mx using model life tables at given level of e0
#' mxm1y <- mlt(seq(65, 71, length = 4), sex = "male", nx = 1)
#' mxf1y <- mlt(seq(73, 78, length = 4), sex = "female", nx = 1)
#'
#' # estimate parameters
#' lc1y <- lileecarter.estimate(mxm1y, mxf1y, nx = 1)
#'
#' # project into the future
#' pred1y <- mortcast(e0m, e0f, lc1y)
#'
#' # plot first projection in black and the remaining ones in grey
#' plot(lc1y$ages, pred1y$female$mx[,1], type="b", log="y", ylim=range(pred1y$female$mx),
#' ylab="female mx", xlab="Age", main="1-year age groups", cex=0.5)
#' for(i in 2:ncol(pred1y$female$mx)) lines(lc1y$ages, pred1y$female$mx[,i], col="grey")
#'
mortcast <- function (e0m, e0f, lc.pars, rotate = TRUE, keep.lt = FALSE,
constrain.all.ages = FALSE, ...) {
get.a0rule <- function(a0rule = c("ak", "cd"), ...)
list(ak = 1, cd = 2)[[match.arg(a0rule)]]
# if e0 is a data.frame, convert to vector (it would not drop dimension without as.matrix)
if(length(dim(e0m)) > 0) e0m <- drop(as.matrix(e0m))
if(length(dim(e0f)) > 0) e0f <- drop(as.matrix(e0f))
# prepare for computation
e0 <- list(female=e0f, male=e0m)
npred <- length(e0f)
nage <- length(lc.pars$ages)
if(lc.pars$nx > 1) {
resnage <- nage-1 # number of age groups of the resulting matrices
age.groups <- lc.pars$ages[-2] # group 0-1 collapsed into 0-5
} else {
resnage <- nage # all ages
age.groups <- lc.pars$ages
}
a0cat <- get.a0rule(...)
zeromatsr <- matrix(0, nrow=resnage, ncol=npred)
zeromatmx <- matrix(0, nrow=nage, ncol=npred)
# in an abridged case, lx, Lx and sr will be returned for 5-year intervals
ressex <- list(mx=zeromatmx, lx=zeromatsr, sr=zeromatsr, Lx=zeromatsr)
result <- list(female = ressex, male = ressex)
# rotate bx if needed
if(rotate)
Bxt <- rotate.leecarter(lc.pars$bx, lc.pars$ultimate.bx, (e0f + e0m)/2)
else
Bxt <- matrix(lc.pars$bx, nrow=nage, ncol=npred)
# allow for time-dependent ax
for(sex in c("female", "male"))
if(is.null(lc.pars[[sex]]$axt))
lc.pars[[sex]]$axt <- matrix(lc.pars[[sex]]$ax, nrow=nrow(Bxt), ncol=npred)
# compute ranges for k(t)
kranges <- .kranges(Bxt, lc.pars$male$axt, lc.pars$female$axt)
#Get the projected kt from e0, and make projection of Mx
for (sex in c("female", "male")) { # iterate over female and male (order matters because of the constrain)
LCres <- .C("LC", as.integer(npred), as.integer(c(female=2, male=1)[sex]), as.integer(nage), as.integer(lc.pars$nx),
as.numeric(lc.pars[[sex]]$axt), as.numeric(Bxt), as.numeric(e0[[sex]]),
Kl=as.numeric(kranges[[sex]]$kl), Ku=as.numeric(kranges[[sex]]$ku),
# 1 for constraining old ages only; 2 for constraining all ages
constrain=as.integer((sex == "male") * ((sex == "male") + (constrain.all.ages == TRUE))),
FMx=as.numeric(result$female$mx), FEop=as.numeric(e0$female), a0rule = as.integer(a0cat),
LLm = as.numeric(result[[sex]]$Lx), Sr=as.numeric(result[[sex]]$sr),
lx=as.numeric(result[[sex]]$lx), Mx=as.numeric(result[[sex]]$mx))
result[[sex]]$mx <- matrix(LCres$Mx, nrow=nage,
dimnames=list(lc.pars$ages, names(e0m)))
result[[sex]]$nx <- lc.pars$nx
if(keep.lt) {
result[[sex]]$sr <- matrix(LCres$Sr, nrow=resnage,
dimnames=list(age.groups, names(e0m)))
result[[sex]]$Lx <- matrix(LCres$LLm, nrow=resnage,
dimnames=list(age.groups, names(e0m)))
result[[sex]]$lx <- matrix(LCres$lx, nrow=resnage,
dimnames=list(age.groups, names(e0m)))
} else {
result[[sex]]$sr <- NULL
result[[sex]]$Lx <- NULL
result[[sex]]$lx <- NULL
}
}
return(result)
}
PPgp/MortCast documentation built on June 12, 2022, 12:53 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .kranges ultimate.bx rotate.leecarter lileecarter.estimate .finish.bx bx.estimate leecarter.estimate
Documented in leecarter.estimate lileecarter.estimate rotate.leecarter ultimate.bx
#' @title Lee-Carter Estimation
#' @description Estimate Lee-Carter parameters (Lee and Carter 1992).
#' @details
#' The function estimates parameters of \eqn{\log(m_x(t)) = a_x + b_x k(t) + \epsilon_x(t)} (Lee and Carter 1992).
#' The argument \code{ax.index} determines which time periods to use to
#' estimate the \eqn{a_x} parameter, while \code{ax.smooth} controls if
#' the resulting \eqn{a_x} should be smoothened over ages (see Sevcikova et al. 2016 for details).
#'
#' @param mx A matrix of age-specific mortality rates where rows correspond to age groups
#' and columns correspond to time periods. Rownames define the starting ages of the age groups.
#' @param ax.index A vector of column indices of \code{mx} to be used to estimate the \eqn{a_x} parameter.
#' By default all time periods are used.
#' @param ax.smooth Logical allowing to smooth the \eqn{a_x} over ages.
#' @param ax.smooth.df Degree of freedom for smoothing if \code{ax.smooth} is \code{TRUE}.
#' Default is half the length of \eqn{a_x}.
#' @param bx.postprocess Logical determining if numerical anomalies in \eqn{b_x} should be dealt with.
#' @param nx Size of age groups. By default ages are determined by rownames of \code{mx}. This argument is only used if
#' \code{mx} has no rownames. If \code{nx} is 5, the age groups are interpreted as 0, 1, 5, 10, \dots. For \code{nx} equals 1,
#' the age groups are interpreted as 0, 1, 2, 3, \dots.
#' @return List with elements \code{ax}, \code{bx} and \code{kt} corresponding to the estimated parameters.
#' @export
#' @seealso \code{\link{mortcast}}, \code{\link{lileecarter.estimate}}
#' @references
#' Lee, R. D. and Carter, L. (1992). Modeling and forecasting the time series of
#' US mortality. Journal of the American Statistical Association, 87, 659-671.
#'
#' Sevcikova H., Li N., Kantorova V., Gerland P., Raftery A.E. (2016).
#' Age-Specific Mortality and Fertility Rates for Probabilistic Population Projections.
#' In: Schoen R. (eds) Dynamic Demographic Analysis. The Springer Series on Demographic Methods
#' and Population Analysis, vol 39. Springer, Cham
#'
#' @examples
#' data(mxM, package = "wpp2017")
#' mx <- subset(mxM, name == "Netherlands")[,4:16]
#' rownames(mx) <- c(0,1, seq(5, 100, by=5))
#' lc.ax.avg <- leecarter.estimate(mx)
#' lc.ax.last <- leecarter.estimate(mx, ax.index=ncol(mx))
#' plot(lc.ax.avg$ax, type="l")
#' lines(lc.ax.last$ax, col="blue")
#'
leecarter.estimate <- function(mx, ax.index = NULL, ax.smooth = FALSE,
ax.smooth.df = NULL, bx.postprocess = TRUE, nx = 5) {
if(length(dim(mx)) < 2 || ncol(mx) < 2) stop("mx must be a matrix (age x time) of at least two columns.")
lmx <- log(mx)
if(length(dim(lmx))==0)
lmx <- as.matrix(lmx)
if(is.null(rownames(lmx))) # default ages
rownames(lmx) <- if(nx > 1) c(0,1, seq(nx, by=nx, length=nrow(lmx)-2)) else 0:(nrow(lmx)-1)
nest <- ncol(lmx)
if(is.null(ax.index)) ax.index <- 1:nest
ax <- apply(lmx[,ax.index, drop=FALSE], 1, sum, na.rm=TRUE) / length(ax.index)
if(ax.smooth) {
if (is.null(ax.smooth.df)) ax.smooth.df <- ceiling(length(ax)/2)
ax.sm <- smooth.spline(ax, df = ax.smooth.df)$y
ax[-1] <- ax.sm[-1] # keep value of the first age group
}
kt <- apply(lmx, 2, sum) - sum(ax)
bx <- bx.estimate(lmx, ax, kt, postprocess = bx.postprocess)
return(list(ax = ax, bx = bx, kt = kt, nx = nx))
}
bx.estimate <- function(lmx, ax, kt, postprocess=TRUE) {
x2 <- sum(kt*kt)
x1 <- rep(NA, nrow(lmx))
for (i in 1:nrow(lmx))
x1[i] <- sum((lmx[i,]-ax[i])*kt)
bx <- x1/x2
if(postprocess)
bx <- .finish.bx(bx)
names(bx) <- rownames(lmx)
return(bx)
}
.finish.bx <- function(bx) {
negbx <- which(bx <= 0)
lnegbx <- length(negbx)
if(lnegbx > 0 && negbx[1] == 1) {
bx[1] <- 0
negbx <- if(lnegbx > 1) negbx[2:lnegbx] else c()
lnegbx <- length(negbx)
}
while(lnegbx > 0) {
bx[negbx] <- 0.5 * bx[negbx-1]
negbx <- which(bx < 0)
lnegbx <- length(negbx)
}
lbx <- length(bx)
for (i in 1:(lbx-1)) {
if (bx[lbx+1 - i] == 0) bx[lbx+1 - i] <- bx[lbx - i]
}
bx <- bx/sum(bx) # must sum to 1
return(bx)
}
#' @title Coherent Lee-Carter Estimation
#' @description Estimate coherent Lee-Carter parameters (Li and Lee 2005).
#' @details The coherent Lee-Carter parameters for male and female mortality rates
#' share the same \eqn{b_x} which is the average of the age-specific
#' \eqn{b_x} parameters.
#'
#' The function in addition computes the ultimate \eqn{b^u_x} as defined in
#' Li et al. (2013) based on the coherent \eqn{b_x}.
#' @param mxM A matrix of male age-specific mortality rates where rows correspond to age groups
#' and columns correspond to time periods. For 5-year age groups, the first and second rows corresponds to
#' 0-1 and 1-5 age groups, respectively. Rownames define the starting ages of the respective groups.
#' @param mxF A matrix of female mortality rates of the same shape as \code{mxM}.
#' @param nx Size of age groups. Should be either 5 or 1.
#' @param ... Additional arguments passed to \code{\link{leecarter.estimate}}.
#' @return List containing elements \code{bx} (coherent \eqn{b_x} parameter),
#' \code{ultimate.bx} (ultimate \eqn{b^u_x} parameter), \code{ages} (age groups), \code{nx} (age group interval), and
#' lists \code{female} and \code{male}, each with the Lee-Carter parameters.
#' @export
#'
#' @references
#' Li, N. and Lee, R. D. (2005). Coherent mortality forecasts for a group of populations:
#' An extension of the Lee-Carter method. Demography, 42, 575-594.
#'
#' Li, N., Lee, R. D. and Gerland, P. (2013). Extending the Lee-Carter method to model the rotation
#' of age patterns of mortality decline for long-term projections. Demography, 50, 2037-2051.
#'
#' @examples
#' data(mxM, mxF, package = "wpp2017")
#' country <- "Germany"
#' mxm <- subset(mxM, name == country)[,4:16]
#' mxf <- subset(mxF, name == country)[,4:16]
#' rownames(mxm) <- rownames(mxf) <- c(0,1, seq(5, 100, by=5))
#' lc <- lileecarter.estimate(mxm, mxf)
#' plot(lc$bx, type="l")
#' lines(lc$ultimate.bx, lty=2)
#'
lileecarter.estimate <- function(mxM, mxF, nx = 5, ...) {
lc.male <- leecarter.estimate(mxM, nx = nx, ...)
lc.female <- leecarter.estimate(mxF, nx = nx, ...)
if(length(lc.female$ax) != length(lc.male$ax) || length(lc.female$kt) != length(lc.male$kt))
stop("Mismatch in dimensions of male and female mortality. Check the mxM and mxF arguments.")
bx <- (lc.male$bx + lc.female$bx)/2
ages <- as.integer(names(lc.female$bx))
# correct nx if wrong input
nx <- ages[4] - ages[3]
return(list(bx = bx, ultimate.bx = ultimate.bx(bx), ages = ages, nx = nx,
female=list(ax=lc.female$ax, bx=bx, kt=lc.female$kt, sex.bx = lc.female$bx),
male=list(ax=lc.male$ax, bx=bx, kt=lc.male$kt, sex.bx = lc.male$bx)
))
}
#' @title Rotated Lee-Carter
#' @description Rotate the Lee-Carter parameter \eqn{b_x} over time to reach an ultimate \eqn{b^u_x},
#' as described in Li et al. (2013).
#'
#' @param bx A vector of the Lee-Carter \eqn{b_x} parameter
#' (from e.g. \code{\link{lileecarter.estimate}} or \code{\link{leecarter.estimate}}).
#' @param ultimate.bx A vector of the ultimate \eqn{b^u_x} parameter as defined in Li, Lee, Gerland (2013)
#' (obtained using \code{\link{lileecarter.estimate}} or \code{\link{ultimate.bx}}).
#' @param e0 A time series of life expectancies.
#' @param e0l Level of life expectancy at which the rotation starts.
#' @param e0u Level of life expectancy at which the rotation finishes.
#' @param p Exponent of the smooth function.
#' @return Function \code{rotate.leecarter} returns a matrix of rotated \eqn{B_x(t)} where rows correspond to age groups and columns
#' correspond to time periods (given by the vector \code{e0}).
#' @export
#'
#' @references
#' Li, N., Lee, R. D. and Gerland, P. (2013). Extending the Lee-Carter method to model the rotation
#' of age patterns of mortality decline for long-term projections. Demography, 50, 2037-2051.
#'
#' @examples
#' data(mxF, mxM, e0Fproj, e0Mproj, package = "wpp2017")
#' country <- "Japan"
#' mxm <- subset(mxM, name == country)[,4:16]
#' mxf <- subset(mxF, name == country)[,4:16]
#' e0f <- as.numeric(subset(e0Fproj, name == country)[-(1:2)])
#' e0m <- as.numeric(subset(e0Mproj, name == country)[-(1:2)])
#' rownames(mxm) <- rownames(mxf) <- c(0,1, seq(5, 100, by=5))
#' lc <- lileecarter.estimate(mxm, mxf)
#' rotlc <- rotate.leecarter(lc$bx, lc$ultimate.bx, (e0f + e0m)/2)
#' plot(lc$bx, type="l")
#' lines(lc$ultimate.bx, col="red")
#' for(i in 1:ncol(rotlc)) lines(rotlc[,i], col="grey")
#'
rotate.leecarter <- function(bx, ultimate.bx, e0, e0l = 80, e0u = 102, p = 0.5) {
npred <- length(e0)
wt <- (e0 - e0l)/(e0u-e0l)
wst <- (0.5*(1+(sin(pi/2*(2*wt-1)))))^p
Bxt <- matrix(NA, nrow=length(bx), ncol=npred)
for(t in 1:npred) {
Bxt[,t] <- switch(cut(e0[t], c(0, e0l, e0u, 9999), labels=FALSE, right=FALSE),
bx,
(1-wst[t])*bx + wst[t]*ultimate.bx,
ultimate.bx)
}
dimnames(Bxt) <- list(names(bx), names(e0))
return(Bxt)
}
#' @rdname rotate.leecarter
#' @return Function \code{ultimate.bx} returns a vector of the ultimate \eqn{b^u_x}.
#' @export
ultimate.bx <- function(bx) {
# ultimate bx (Li, Lee, Gerland 2013)
bux <- bx
lbx <- length(bx)
idx15 <- which(names(bx) == 15)
idx60 <- which(names(bx) == 60)
bux[1:idx60] <- mean(bux[idx15:idx60])
bux[(idx60+1):lbx] <- bux[(idx60+1):lbx] * (bux[idx60]/bux[idx60+1]) # adjust so that b(65-69)=b(60-64)
return(bux/sum(bux)) # must sum to 1
}
.kranges <- function(Bxt, axm, axf) {
lmin <- -12
lmax <- 0
machine.max <- log(.Machine$double.xmax)
machine.min <- log(.Machine$double.xmin)
npred <- ncol(Bxt)
kranges <- list(male=list(), female=list())
ax <- list(male=axm, female=axf)
bx.lim <- rbind(apply(Bxt, 2, function(x) min(x[x>0])),
apply(Bxt, 2, function(x) max(x[x>0])))
for(sex in c("male", "female")) {
this.lmin <- floor(min(lmin, min(ax[[sex]])-0.1))
kranges[[sex]]$kl <- pmin((lmax - apply(ax[[sex]], 2, max))/bx.lim[2,], machine.max)
kranges[[sex]]$ku <- pmax((this.lmin - apply(ax[[sex]], 2, min))/bx.lim[1,], machine.min)
}
return(kranges)
}
#' @title Coherent Rotated Lee-Carter Prediction
#' @description Predict age-specific mortality rates using the coherent rotated Lee-Carter method.
#' @details This function implements Steps 6-9 of Algorithm 2 in Sevcikova et al. (2016).
#' It uses the abridged or unabridged life table function to find the level of mortality that coresponds to the given
#' life expectancy. Thus, it can be used for both, mortality for 5- or 1-year age groups.
#' @param e0m A time series of future male life expectancy.
#' @param e0f A time series of future female life expectancy.
#' @param lc.pars A list of coherent Lee-Carter parameters with elements \code{bx}, \code{ultimate.bx},
#' \code{ages}, \code{nx},
#' \code{female} and \code{male} as returned by \code{\link{lileecarter.estimate}}.
#' The \code{female} and \code{male} objects are again lists that should contain a vector
#' \code{ax} and optionally a matrix \code{axt} if the \eqn{a_x} parameter
#' needs to be defined as time dependent. In such a case, rows are age groups and columns are
#' time periods corresponding to the length of the \code{e0f} and \code{e0m} vectors.
#' @param rotate If \code{TRUE} the rotation method of \eqn{b_x} is used as described in Li et al. (2013).
#' @param keep.lt Logical. If \code{TRUE} additional life table columns are kept in the
#' resulting object.
#' @param constrain.all.ages By default the method constrains the male mortality to be above female
#' mortality for old ages if the male life expectancy is below the female life expectancy. Setting
#' this argument to \code{TRUE} causes this constraint to be applied to all ages.
#' @param \dots Additional life table arguments.
#' @return List with elements \code{female} and \code{male}, each of which contains a matrix \code{mx}
#' with the predicted mortality rates. If \code{keep.lt} is \code{TRUE}, it also
#' contains matrices \code{sr} (survival rates), and life table quantities \code{Lx} and \code{lx}.
#' @export
#'
#' @seealso \code{\link{rotate.leecarter}}, \code{\link{leecarter.estimate}}, \code{\link{lileecarter.estimate}},
#' \code{\link{mortcast.blend}}
#' @references
#' Li, N., Lee, R. D. and Gerland, P. (2013). Extending the Lee-Carter method to model the rotation
#' of age patterns of mortality decline for long-term projections. Demography, 50, 2037-2051.
#'
#' Sevcikova H., Li N., Kantorova V., Gerland P., Raftery A.E. (2016).
#' Age-Specific Mortality and Fertility Rates for Probabilistic Population Projections.
#' In: Schoen R. (eds) Dynamic Demographic Analysis. The Springer Series on Demographic Methods
#' and Population Analysis, vol 39. Springer, Cham
#'
#' @examples
#' # estimate parameters from historical mortality data (5-year age groups)
#' data(mxM, mxF, e0Fproj, e0Mproj, package = "wpp2017")
#' country <- "Brazil"
#' mxm <- subset(mxM, name == country)[,4:16]
#' mxf <- subset(mxF, name == country)[,4:16]
#' rownames(mxm) <- rownames(mxf) <- c(0,1, seq(5, 100, by=5))
#' lc <- lileecarter.estimate(mxm, mxf)
#'
#' # project into future for given levels of life expectancy
#' e0f <- as.numeric(subset(e0Fproj, name == country)[-(1:2)])
#' e0m <- as.numeric(subset(e0Mproj, name == country)[-(1:2)])
#' pred <- mortcast(e0m, e0f, lc)
#'
#' # plot first projection in black and the remaining ones in grey
#' plot(lc$ages, pred$female$mx[,1], type="b", log="y", ylim=range(pred$female$mx),
#' ylab="female mx", xlab="Age", main=paste(country, "(5-year age groups)"), cex=0.5)
#' for(i in 2:ncol(pred$female$mx)) lines(lc$ages, pred$female$mx[,i], col="grey")
#'
#' # similarly for 1-year age groups
#' # derive toy 1-year mx using model life tables at given level of e0
#' mxm1y <- mlt(seq(65, 71, length = 4), sex = "male", nx = 1)
#' mxf1y <- mlt(seq(73, 78, length = 4), sex = "female", nx = 1)
#'
#' # estimate parameters
#' lc1y <- lileecarter.estimate(mxm1y, mxf1y, nx = 1)
#'
#' # project into the future
#' pred1y <- mortcast(e0m, e0f, lc1y)
#'
#' # plot first projection in black and the remaining ones in grey
#' plot(lc1y$ages, pred1y$female$mx[,1], type="b", log="y", ylim=range(pred1y$female$mx),
#' ylab="female mx", xlab="Age", main="1-year age groups", cex=0.5)
#' for(i in 2:ncol(pred1y$female$mx)) lines(lc1y$ages, pred1y$female$mx[,i], col="grey")
#'
mortcast <- function (e0m, e0f, lc.pars, rotate = TRUE, keep.lt = FALSE,
constrain.all.ages = FALSE, ...) {
get.a0rule <- function(a0rule = c("ak", "cd"), ...)
list(ak = 1, cd = 2)[[match.arg(a0rule)]]
# if e0 is a data.frame, convert to vector (it would not drop dimension without as.matrix)
if(length(dim(e0m)) > 0) e0m <- drop(as.matrix(e0m))
if(length(dim(e0f)) > 0) e0f <- drop(as.matrix(e0f))
# prepare for computation
e0 <- list(female=e0f, male=e0m)
npred <- length(e0f)
nage <- length(lc.pars$ages)
if(lc.pars$nx > 1) {
resnage <- nage-1 # number of age groups of the resulting matrices
age.groups <- lc.pars$ages[-2] # group 0-1 collapsed into 0-5
} else {
resnage <- nage # all ages
age.groups <- lc.pars$ages
}
a0cat <- get.a0rule(...)
zeromatsr <- matrix(0, nrow=resnage, ncol=npred)
zeromatmx <- matrix(0, nrow=nage, ncol=npred)
# in an abridged case, lx, Lx and sr will be returned for 5-year intervals
ressex <- list(mx=zeromatmx, lx=zeromatsr, sr=zeromatsr, Lx=zeromatsr)
result <- list(female = ressex, male = ressex)
# rotate bx if needed
if(rotate)
Bxt <- rotate.leecarter(lc.pars$bx, lc.pars$ultimate.bx, (e0f + e0m)/2)
else
Bxt <- matrix(lc.pars$bx, nrow=nage, ncol=npred)
# allow for time-dependent ax
for(sex in c("female", "male"))
if(is.null(lc.pars[[sex]]$axt))
lc.pars[[sex]]$axt <- matrix(lc.pars[[sex]]$ax, nrow=nrow(Bxt), ncol=npred)
# compute ranges for k(t)
kranges <- .kranges(Bxt, lc.pars$male$axt, lc.pars$female$axt)
#Get the projected kt from e0, and make projection of Mx
for (sex in c("female", "male")) { # iterate over female and male (order matters because of the constrain)
LCres <- .C("LC", as.integer(npred), as.integer(c(female=2, male=1)[sex]), as.integer(nage), as.integer(lc.pars$nx),
as.numeric(lc.pars[[sex]]$axt), as.numeric(Bxt), as.numeric(e0[[sex]]),
Kl=as.numeric(kranges[[sex]]$kl), Ku=as.numeric(kranges[[sex]]$ku),
# 1 for constraining old ages only; 2 for constraining all ages
constrain=as.integer((sex == "male") * ((sex == "male") + (constrain.all.ages == TRUE))),
FMx=as.numeric(result$female$mx), FEop=as.numeric(e0$female), a0rule = as.integer(a0cat),
LLm = as.numeric(result[[sex]]$Lx), Sr=as.numeric(result[[sex]]$sr),
lx=as.numeric(result[[sex]]$lx), Mx=as.numeric(result[[sex]]$mx))
result[[sex]]$mx <- matrix(LCres$Mx, nrow=nage,
dimnames=list(lc.pars$ages, names(e0m)))
result[[sex]]$nx <- lc.pars$nx
if(keep.lt) {
result[[sex]]$sr <- matrix(LCres$Sr, nrow=resnage,
dimnames=list(age.groups, names(e0m)))
result[[sex]]$Lx <- matrix(LCres$LLm, nrow=resnage,
dimnames=list(age.groups, names(e0m)))
result[[sex]]$lx <- matrix(LCres$lx, nrow=resnage,
dimnames=list(age.groups, names(e0m)))
} else {
result[[sex]]$sr <- NULL
result[[sex]]$Lx <- NULL
result[[sex]]$lx <- NULL
}
}
return(result)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.