Nothing
#' @title Kannisto Method
#' @description Extrapolate given mortality rates using the original Kannisto method.
#' @details
#' The function first estimates the original Kannisto parameters
#' by passing mortality rates for age groups \code{est.ages} into
#' the \code{\link{kannisto.estimate}} function.
#' The estimated parameters are then passed to the projection function
#' \code{\link{kannisto.predict}} to extrapolate into ages \code{proj.ages}.
#' Lastly, the input mortality object is extended by results for the extrapolated ages.
#' If \code{proj.ages} contains age groups that are included in \code{mx},
#' values for those age groups are overwritten.
#' @param mx A vector or matrix of mortality rates. If it is a matrix,
#' rows correspond to age groups with rownames identifying the corresponding age as integers.
#' By default five-years age groups are assigned to rows if rownames are not given.
#' @param est.ages A vector of integers identifying the ages to be used
#' for estimation. It should be a subset of rownames of \code{mx}. Change the defaults if
#' 1-year age groups are used (see Example).
#' @param proj.ages A vector of integers identifying the age groups for which mortality rates
#' are to be projected. Change the defaults if 1-year age groups are used (see Example).
#' @return A vector or matrix containing the input mortality object \code{mx}
#' extended by the extrapolated age groups.
#' @export
#' @references
#' Thatcher, A. R., Kannisto, V. and Vaupel, J. W. (1998). The Force of Mortality at Ages 80 to 120,
#' volume 5 of Odense Monographs on Population Aging Series. Odense, Denmark: Odense University Press.
#'
#' @seealso \code{\link{kannisto.estimate}}, \code{\link{kannisto.predict}}, \code{\link{cokannisto}}
#' @examples
#' data(mxM, package = "wpp2017")
#' mx <- subset(mxM, name == "Burkina Faso")[,-(1:3)]
#' rownames(mx) <- c(0,1, seq(5, 100, by=5))
#' mxnew <- kannisto(mx)
#' ages <- as.integer(rownames(mxnew))
#' plot(ages, mxnew[,"2095-2100"], type="l", log="y",
#' xlab="age", ylab="mx", col="red")
#' lines(ages, mxnew[,"2010-2015"])
#'
#' # Kannisto for 1-year age groups
#' # derive toy 1-year mx using model life tables at e0 of 70
#' mx1y <- mlt(70, sex = "male", nx = 1)
#' # Pretend we only observed mx for ages 0:100.
#' # Use 90-99 for estimation and extend mx from 100 to 140
#' mx1ynew <- kannisto(mx1y[1:100, , drop = FALSE], est.ages = 90:99, proj.ages = 100:140)
#' # Plot the new mx for old ages
#' plot(80:140, mx1ynew[81:141], type = "l", xlab="age", ylab="mx", col="red")
#' # Check how it compares to the original mx that was not used in the estimation
#' lines(100:130, mx1y[101:nrow(mx1y)])
#'
kannisto <- function(mx, est.ages = seq(80, 95, by=5),
proj.ages = seq(100, 130, by=5)) {
ages <- if(length(dim(mx)) == 0) names(mx) else rownames(mx)
Mxe <- as.matrix(mx)
if(is.null(ages)) # default ages
ages <- as.character(c(0,1, seq(5, by=5, length=nrow(Mxe)-2)))
ages.num <- as.integer(ages)
rownames(Mxe) <- ages
est.ages.char <- as.character(est.ages)
if(any(!est.ages.char %in% rownames(Mxe)))
stop("est.ages are not included in mx. Check the rownames of mx.")
est.data <- Mxe[est.ages.char,, drop = FALSE]
kann.pars <- apply(est.data, 2, kannisto.estimate, ages = est.ages)
kanncoefs <- lapply(kann.pars, function(x) x$coefficients)
res <- lapply(kanncoefs, kannisto.predict, ages=proj.ages)
mres <- sapply(res, cbind)
all.ages <- as.character(sort(unique(c(as.integer(ages), proj.ages))))
proj.ages.char <- as.character(proj.ages)
resMx <- matrix(NA, nrow=length(all.ages), ncol=ncol(Mxe),
dimnames=list(all.ages, colnames(Mxe)))
rownames(resMx) <- all.ages
resMx[proj.ages.char,] <- mres
orig.ages <- rownames(Mxe)[rownames(Mxe) %in% rownames(resMx) & !rownames(Mxe) %in% proj.ages.char]
resMx[orig.ages,] <- Mxe[orig.ages,]
if(length(dim(mx)) == 0) { # convert to vector
resMxv <- resMx[,1]
names(resMxv) <- rownames(resMx)
resMx <- resMxv
}
return(resMx)
}
#' @title Coherent Kannisto Method
#' @description Extrapolate given mortality rates into higher ages
#' using the Coherent Kannisto method as described in Sevcikova et al. (2016).
#' @details
#' The function first estimates the coherent Kannisto parameters
#' by passing mortality rates for age groups \code{est.ages} into
#' the \code{\link{cokannisto.estimate}} function.
#' The estimated parameters are then passed to the projection function
#' \code{\link{kannisto.predict}} to extrapolate into ages \code{proj.ages}.
#' Lastly, the input mortality objects are extended by results for the extrapolated ages.
#' If \code{proj.ages} contains age groups that are included in \code{mxM} and \code{mxF},
#' values for those age groups are overwritten.
#' @param mxM A vector or matrix of male mortality rates. If it is a matrix,
#' rows correspond to age groups with rownames identifying the corresponding age as integers.
#' By default five-years age groups are assigned to rows if rownames are not given.
#' @param mxF A vector or matrix of female mortality rates. Its length or dimension
#' should be the same \code{mxM}.
#' @param est.ages A vector of integers identifying the ages to be used
#' for estimation. It should be a subset of rownames of \code{mxM}. Change the defaults if
#' 1-year age groups are used (see Example in \code{\link{kannisto}}).
#' @param proj.ages A vector of integers identifying the age groups for which mortality rates
#' are to be projected. Change the defaults if
#' 1-year age groups are used (see Example in \code{\link{kannisto}}).
#' @return A list of two vectors or matrices (for male and female) containing the input motality
#' objects extended by the extrapolated age groups.
#' @export
#'
#' @seealso \code{\link{cokannisto.estimate}}, \code{\link{kannisto.predict}}
#' @references
#' 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, mxF, package = "wpp2017")
#' country <- "South Africa"
#' mxm <- subset(mxM, name == country)[,-(1:3)]
#' mxf <- subset(mxF, name == country)[,-(1:3)]
#' rownames(mxm) <- rownames(mxf) <- c(0,1, seq(5, 100, by=5))
#' mxnew <- cokannisto(mxm, mxf)
#' ages <- as.integer(rownames(mxnew$male))
#' plot(ages, mxnew$male[,"2095-2100"], type="l", log="y",
#' xlab="age", ylab="mx", col="blue", main=country)
#' lines(ages, mxnew$female[,"2095-2100"], col="red")
#' lines(ages, mxnew$male[,"2010-2015"], lty=2, col="blue")
#' lines(ages, mxnew$female[,"2010-2015"], lty=2, col="red")
#' legend("bottomright", legend=c("male 2010-2015", "female 2010-2015",
#' "male 2095-2100", "female 2095-2100"), bty="n",
#' col=rep(c("blue", "red"),2), lty=c(2,2,1,1))
#'
cokannisto <- function(mxM, mxF,
est.ages = seq(80, 95, by=5),
proj.ages = seq(100, 130, by=5)) {
ages <- if(length(dim(mxM)) == 0) names(mxM) else rownames(mxM)
MxeM <- as.matrix(mxM)
MxeF <- as.matrix(mxF)
if(is.null(ages)) # default ages
ages <- as.character(c(0,1, seq(5, by=5, length=nrow(MxeM)-2)))
ages.num <- as.integer(ages)
rownames(MxeM) <- rownames(MxeF) <- ages
est.ages.char <- as.character(est.ages)
if(any(!est.ages.char %in% rownames(MxeM)))
stop("est.ages are not included in mxM.")
est.data <- rbind(MxeM[est.ages.char,,drop=FALSE], MxeF[est.ages.char,,drop=FALSE])
nest <- length(est.ages)
kann.pars <- apply(est.data, 2,
function(x) cokannisto.estimate(x[1:nest], x[(nest+1):length(x)],
ages = est.ages, fitted = FALSE))
male.pars <- lapply(kann.pars, function(x) x[["male"]]$coefficients)
resM <- lapply(male.pars, kannisto.predict, ages=proj.ages)
matresM <- sapply(resM, cbind)
female.pars <- lapply(kann.pars, function(x) x[["female"]]$coefficients)
resF <- lapply(female.pars, kannisto.predict, ages=proj.ages)
matresF <- sapply(resF, cbind)
all.ages <- as.character(sort(unique(c(as.integer(ages), proj.ages))))
proj.ages.char <- as.character(proj.ages)
resMxM <- matrix(NA, nrow=length(all.ages), ncol=ncol(MxeM),
dimnames=list(all.ages, colnames(MxeM)))
resMxF <- matrix(NA, nrow=length(all.ages), ncol=ncol(MxeF),
dimnames=list(all.ages, colnames(MxeF)))
resMxM[proj.ages.char,] <- matresM
resMxF[proj.ages.char,] <- matresF
orig.ages <- rownames(MxeM)[rownames(MxeM) %in% rownames(resMxM) & !rownames(MxeM) %in% proj.ages.char]
resMxM[orig.ages,] <- MxeM[orig.ages,]
resMxF[orig.ages,] <- MxeF[orig.ages,]
if(length(dim(mxM)) == 0) { # convert to vector
resMxMv <- resMxM[,1]
resMxFv <- resMxF[,1]
names(resMxMv) <- names(resMxFv) <- rownames(resMxM)
resMxM <- resMxMv
resMxF <- resMxFv
}
return(list(male=resMxM, female=resMxF))
}
#' @title Kannisto Estimation
#' @description Estimate the Kannisto parameters (Thatcher et al. 1998).
#' @details Given the Kannisto equation \eqn{logit(m_x) = \log(c) + dx}{logit(mx) = log(c) + dx},
#' the function estimates the \eqn{c} and \eqn{d} parameters using
#' values of \code{ages} as the covariate \eqn{x}.
#' @param mx A vector of mortality rates.
#' @param ages A vector of ages corresponding to \code{mx}. These can be indices of age groups or raw ages.
#' @return List with the following components:
#' \describe{
#' \item{\code{coefficients}:}{named vector with Kannisto coefficients \eqn{c} and \eqn{d}.}
#' \item{\code{fitted.values}:}{the fitted values}
#' \item{\code{residuals}:}{input rates minus the fitted values}
#'}
#' @export
#' @references
#' Thatcher, A. R., Kannisto, V. and Vaupel, J. W. (1998). The Force of Mortality at Ages 80 to 120,
#' volume 5 of Odense Monographs on Population Aging Series. Odense, Denmark: Odense University Press.
#'
#' @seealso \code{\link{kannisto.predict}}, \code{\link{kannisto}}, \code{\link{cokannisto.estimate}}
#'
#' @examples
#' data(mxM, package = "wpp2017")
#' mx <- subset(mxM, name == "Canada")[,"2010-2015"]
#' kannisto.estimate(mx[18:21], ages = 18:21)
#'
kannisto.estimate <- function(mx, ages){
y <- log(mx) - log(1-mx)
x <- ages
coefs <- coefficients(lm(y ~ x))
cf <- c(c=exp(coefs[[1]]), d=coefs[['x']])
fitted <- kannisto.predict(cf, ages)
return(list(coefficients = cf,
fitted.values = fitted,
residuals = mx - fitted))
}
#' @title Coherent Kannisto Estimation
#' @description Estimate the coherent Kannisto parameters as described in Sevcikova et al. (2016).
#' @details Given the Kannisto equation \eqn{logit(m_x) = \log(c) + dx}{logit(mx) = log(c) + dx},
#' the Coherent Kannisto method estimates the \eqn{d} parameter jointly for male and female
#' data, in order to prevent mortality crossovers in higher ages.
#' @param mxM A vector of male mortality rates.
#' @param mxF A vector of female mortality rates.
#' @param ages A vector of ages corresponding to \code{mxM} and \code{mxF}.
#' @param fitted Logical. If \code{TRUE} the fitted values and residuals are returned.
#' @return List of two lists, one for male and one for female. Each of the two lists contains the following components:
#' \describe{
#' \item{\code{coefficients}:}{named vector with the coherent Kannisto coefficients \eqn{c} and \eqn{d}. The \eqn{d} values are the same in both lists.}
#' \item{\code{fitted.values}:}{the fitted values (not included if \code{fitted} is \code{FALSE})}
#' \item{\code{residuals}:}{input rates minus the fitted values (not included if \code{fitted} is \code{FALSE})}
#'}
#' @export
#'
#' @seealso \code{\link{cokannisto}}, \code{\link{kannisto.predict}}, \code{\link{kannisto}}
#'
#' @references
#' 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, mxF, package = "wpp2017")
#' country <- "Brazil"
#' mxm <- subset(mxM, name == country)[,"2010-2015"]
#' mxf <- subset(mxF, name == country)[,"2010-2015"]
#' cokannisto.estimate(mxm[18:21], mxf[18:21], ages = 18:21)
#'
cokannisto.estimate <- function(mxM, mxF, ages, fitted = TRUE){
mx <- c(mxM, mxF)
y <- log(mx) - log(1-mx)
x <- c(ages, ages)
g <- c(rep(1, length(mxM)), rep(0, length(mxF)))
coefs <- coefficients(lm(y ~ g + x))
female.coefs <- c(c=exp(coefs[[1]]), d=coefs[['x']])
male.coefs <- c(c=exp(coefs[[1]] + coefs[['g']]), d=coefs[['x']])
result <- list(female = list(coefficients = female.coefs),
male = list(coefficients = male.coefs))
if(fitted) {
fittedF <- kannisto.predict(female.coefs, ages)
result$female$fitted.values <- fittedF
result$female$residuals <- mxF - fittedF
fittedM <- kannisto.predict(male.coefs, ages)
result$male$fitted.values <- fittedM
result$male$residuals <- mxM - fittedM
}
return(result)
}
#' @title Kannisto Prediction
#' @description Given estimated Kannisto parameters (coherent or original),
#' it predicts mortality rates for given ages.
#' @details Given parameters \eqn{c} and \eqn{d} in \code{pars},
#' the function uses the Kannisto equation \eqn{logit(m_x) = \log(c) + dx}{logit(mx) = log(c) + dx},
#' to predict mortality rates for age groups \eqn{x} given by \code{ages}.
#' @param pars A named vector with Kanisto coefficients \eqn{c} and \eqn{d}
#' (e.g. result of \code{\link{kannisto.estimate}} or \code{\link{cokannisto.estimate}}).
#' @param ages A vector of ages to make prediction for. These can be indices of age groups or raw ages,
#' but on the same scale as used in the estimation.
#' @return Vector of predicted mortality rates.
#' @export
#' @seealso \code{\link{cokannisto}}, \code{\link{kannisto.estimate}}, \code{\link{cokannisto.estimate}}
#' @references
#' Thatcher, A. R., Kannisto, V. and Vaupel, J. W. (1998). The Force of Mortality at Ages 80 to 120,
#' volume 5 of Odense Monographs on Population Aging Series. Odense, Denmark: Odense University Press.
#'
#' @examples
#' data(mxM, mxF, package = "wpp2017")
#' mxm <- subset(mxM, name == "Germany")[,"2010-2015"]
#' ages <- c(0, 1, seq(5, 130, by=5))
#'
#' # using original Kannisto parameters
#' pars <- kannisto.estimate(mxm[18:21], ages = ages[18:21])
#' mxm.pred <- kannisto.predict(pars$coefficients, ages = ages[22:28])
#' plot(ages, c(mxm[1:21], mxm.pred), type="l", log="y",
#' xlab="age", ylab="mx")
#'
#' # Coherent Kannisto
#' mxf <- subset(mxF, name == "Germany")[,"2010-2015"]
#' copars <- cokannisto.estimate(
#' mxm[18:21], mxf[18:21], ages = ages[18:21])
#' cmxm.pred <- kannisto.predict(copars[["male"]]$coefficients, ages = ages[22:28])
#' cmxf.pred <- kannisto.predict(copars[["female"]]$coefficients, ages = ages[22:28])
#' plot(ages, c(mxm[1:21], cmxm.pred), type="l", log="y",
#' xlab="age", ylab="mx", col="blue")
#' lines(ages, c(mxf[1:21], cmxf.pred), col="red")
#'
kannisto.predict <- function(pars, ages){
numer <- pars["c"] * exp(pars["d"] * ages)
return(numer/(1 + numer))
}
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.