From 043501e277fffd82077ac6955c7de818c0d16079 Mon Sep 17 00:00:00 2001 From: Ivan Williams Date: Mon, 23 Sep 2024 16:01:04 -0300 Subject: [PATCH] include HMD_old_logquad --- NAMESPACE | 2 + R/log_quad_augm.R | 173 +++++++++++++++++++++++++++++++++++++++ man/HMD_old_logquad.Rd | 83 +++++++++++++++++++ man/logquad_augmented.Rd | 35 ++++++++ 4 files changed, 293 insertions(+) create mode 100644 R/log_quad_augm.R create mode 100644 man/HMD_old_logquad.Rd create mode 100644 man/logquad_augmented.Rd diff --git a/NAMESPACE b/NAMESPACE index cd4760c5..0f09fb05 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,6 +3,7 @@ export("%>%") export(ADM) export(AGEN) +export(HMD_old_logquad) export(ID) export(IRD) export(OPAG) @@ -71,6 +72,7 @@ export(is_age_redundant) export(is_age_sequential) export(is_single) export(loess_smth1) +export(logquad_augmented) export(lt_a_closeout) export(lt_a_pas) export(lt_a_un) diff --git a/R/log_quad_augm.R b/R/log_quad_augm.R new file mode 100644 index 00000000..ab96f471 --- /dev/null +++ b/R/log_quad_augm.R @@ -0,0 +1,173 @@ + +# author: IW -------------------------------------------------------------- + +#' HMD pattern for adult ages. +#' @description Adjust rates in oldest ages using HMD pattern, based on log-quad method. +#' @details One possible scenario when mortality data on last ages is not reliable, is to use a mortality pattern with some known index in previous ages. +#' This function gives a HMD pattern based on 5q0 and 45q15, using log-quad model. Additionally, a value on mortality between 60 and 75 can be included to make a better adjustment in level. +#' @param nMx numeric. Vector of mortality rates in abridged age classes. +#' @param Age integer. Single ages (abridged not allowed). +#' @param Sex character. Either male \code{"m"}, female \code{"f"}, or both \code{"b"}. +#' @param q0_5 numeric. Probability of death from born to age 5. By default implicit values in `nMx` should be entered. +#' @param q15_45 numeric. Probability of death from age 15 to age 60. By default implicit values in `nMx` should be entered. +#' @param fitted_logquad Optional, defaults to \code{NULL}. An object of class +#' \code{wilmoth}. If full HMD is not enough, one +#' can fit a Log-Quadratic (\url{https://github.com/mpascariu/MortalityEstimate}) model +#' based on any other collection of life tables; +#' @param q60_15 numeric. Probability of death from age 60 to age 75. When external information on those ages level is available, +#' can be included to increase parameter `ax` from log-quad model in last ages (Li, 2003). +#' @param Age_transition integer. Form which age should transition to HMD pattern starts. +#' @param window_transition integer. Number of ages to the left and to the right of `Age_transition` to do a log-linear transition in rates. +#' @param plot_comparison Show or not a plot with the result. +#' @param ... Other arguments to be passed on to the \code{lt_single} function. +#' @export +#' @return life table as in \code{lt_single} function. +#' @examples +#' # Mortality rates from UN Chilean with e0=70. Wat would be the rates based on HMD pattern? +#' # In this case making a transition of 10 years at age 80, and returning an OAG=100. +#' \dontrun{ +#' lt <- DemoToolsData::modelLTx1 +#' lt <- lt[lt$family == "Chilean" & lt$sex == "female" & lt$e0 == 70,] +#' chilean70_adjHMD <- HMD_old_logquad(nMx = lt$mx1, +#' Age = lt$age, +#' Sex = "f", +#' q0_5 = 1 - lt$lx1[lt$age==5]/lt$lx1[lt$age==0], +#' q15_45 = 1 - lt$lx1[lt$age==60]/lt$lx1[lt$age==15], +#' Age_transition = 80, +#' window_transition = 10, +#' plot_comparison = TRUE, +#' OAnew = 100) +#' # We know (as an example) that q60_15 is .5 higher than what HMD pattern would be. +#' chilean70_adjHMD_augm <- HMD_old_logquad(nMx = lt$mx1, +#' Age = lt$age, +#' Sex = "f", +#' q0_5 = 1 - lt$lx1[lt$age==5]/lt$lx1[lt$age==0], +#' q15_45 = 1 - lt$lx1[lt$age==60]/lt$lx1[lt$age==15], +#' q60_15 = (1 - lt$lx1[lt$age==75]/lt$lx1[lt$age==60]) * 1.5, +#' Age_transition = 80, window_transition = 10, +#' OAnew = 100, plot_comparison = TRUE) +#' } + +HMD_old_logquad <- function(nMx, Age = NULL, + Sex = "b", + q0_5 = NULL, q15_45 = NULL, + q60_15 = NULL, + Age_transition = 80, + window_transition = 3, + plot_comparison = FALSE, + fitted_logquad = NULL, + ...){ + + # check if age is not complete + if(!is_single(Age)) stop("Need rates by single age.") + if(is.null(Age)) Age <- 0:length(nMx) + + # check if an optional fitted_logquad is specified + if(is.null(fitted_logquad)){ + if(Sex == "b"){ + fitted_logquad <- DemoTools::fitted_logquad_b + } + if(Sex == "f"){ + fitted_logquad <- DemoTools::fitted_logquad_f + } + if(Sex == "m"){ + fitted_logquad <- DemoTools::fitted_logquad_m + } + } + + # load the log-quad model already fitted in the package (different from MortalityEstimate) and estimate with input data + logquad_model <- lt_model_lq(Sex = Sex, q0_5 = q0_5, q15_45 = q15_45, fitted_logquad = fitted_logquad) + mx_logquad_5q0_45q15 <- logquad_model$lt + + # augmented method if was asked + if(!is.null(q60_15)){ + mx_logquad_5q0_45q15 <- tryCatch({ + logquad_augmented(coeffs = fitted_logquad$coefficients, + k = logquad_model$values$k, + Age = logquad_model$lt$Age, + q0_5 = q0_5, q60_15 = q60_15, Sex = Sex)}, + error = function(e) { + warning("Augmented log-quad was not possible. Revise q60_15. Returned basic log-quad.") + mx_logquad_5q0_45q15}) + } + + # make it single + mx_logquad_5q0_45q15 <- lt_abridged2single(nMx = mx_logquad_5q0_45q15$nMx, + Age = mx_logquad_5q0_45q15$Age, + lx = mx_logquad_5q0_45q15$lx, + Sex = Sex, + ...) + + # smooth transition with a given length window + Ages <- 0:max(mx_logquad_5q0_45q15$Age) + Age_smooth <- (Age_transition-window_transition):(Age_transition+window_transition) + Age_not_smooth <- Ages[!Ages %in% Age_smooth] + nMx_to_smooth_transition <- data.frame(Age = Age_not_smooth, + nMx = c(nMx[Agemax(Age_smooth)])) + nMx_interpolated <- exp(stats::approx(x = nMx_to_smooth_transition$Age, + y = log(nMx_to_smooth_transition$nMx), + xout = Age_smooth)$y) + smooth_transtition_nMx <- dplyr::bind_rows( + nMx_to_smooth_transition, + data.frame(Age = Age_smooth, nMx = nMx_interpolated)) %>% + dplyr::arrange(Age) + + # plot diagnostic + if(plot_comparison){ + df <- + rbind( + data.frame(Age = Age, nMx = nMx, Type = "Input"), + data.frame(Age = smooth_transtition_nMx$Age, nMx = smooth_transtition_nMx$nMx, Type = "Adjusted"), + data.frame(Age = Age_smooth, nMx = nMx_interpolated, Type = "Transition")) %>% + ggplot2::ggplot(ggplot2::aes(x = Age, y = nMx, color = Type)) + + ggplot2::geom_line() + + ggplot2::geom_vline(xintercept = Age_transition, linetype = "dashed", color = "grey") + + ggplot2::scale_y_log10() + + ggplot2::theme_bw() + } + + # rebuild lt and finish, with input OAnew if was not defined as additional input argument + lt_extra_arguments <- list(...) + if(!("OAnew" %in% names(lt_extra_arguments))){OAnew <- max(Age)} else {OAnew <- lt_extra_arguments$OAnew} + mx_logquad_5q0_45q15 <- lt_single_mx(nMx = smooth_transtition_nMx$nMx, + Age = smooth_transtition_nMx$Age, + Sex = Sex, + OAnew = OAnew) + return(mx_logquad_5q0_45q15) +} + +#' Augmented logquad +#' @description Adjust rates in oldest ages that comes from a HMD model, using an external estimate of 15q60 (Li, 2014). As an example see\code{\link[DemoTools]{HMD_old_logquad}}. +#' @details Parameter \code{a(x)} is augmented based on an external estimate of 15q60. +#' @param coeffs data.frame. Columns \code{a(x)}, \code{b(x)}, \code{c(x)} and \code{v(x)} from fitted logquad model. See \code{fitted_logquad_b}. +#' @param k numeric. Adult mortality related value from log-quad estimatation based on one or two input parameters. See \code{lt_model_lq}. +#' @param Sex character. Either male \code{"m"}, female \code{"f"}, or both \code{"b"}. +#' @param Age integer. Abridged lower bound ages. Same length than rows in `coeffs`. +#' @param q0_5 numeric. Probability of death from born to age 5. +#' @param q60_15 numeric. Probability of death from age 60 to age 75. +#' @param ... Other arguments to be passed on to the \code{lt_abridged} function. +#' @export +#' @return life table as in \code{lt_abridged} function. +#' @references See [Li (2014)](https://www.un.org/development/desa/pd/content/estimating-life-tables-developing-countries). + +logquad_augmented <- function(coeffs, k, q0_5, q60_15, Sex = "b", Age, ...){ + + # arrange parameters and get rates from the model + params <- c(1, log(q0_5), log(q0_5)^2, k) + mx_logquad <- exp(rowSums(as.matrix(coeffs) %*% diag(params))) + + # adjust ax with the ratio between input value and implicit value from model + q60_15_hat <- q60_15 + age_q60_15 <- which(Age == 60) + q60_15 <- 1 - exp(-5*(sum(mx_logquad[c(age_q60_15, age_q60_15+1, age_q60_15+2)]))) + coeffs_hat <- coeffs + coeffs_hat$ax[age_q60_15:nrow(coeffs_hat)] <- coeffs_hat$ax[age_q60_15:nrow(coeffs_hat)] + log(log(1-q60_15_hat)/log(1-q60_15)) + + # new rates with changed ax + mx_logquad_augm <- exp(rowSums(as.matrix(coeffs_hat) %*% diag(params))) + lt_logquad_augm <- lt_abridged(nMx = mx_logquad_augm, Age = Age, Sex = Sex, ...) + + # output + return(lt_logquad_augm) +} \ No newline at end of file diff --git a/man/HMD_old_logquad.Rd b/man/HMD_old_logquad.Rd new file mode 100644 index 00000000..dd6f75f8 --- /dev/null +++ b/man/HMD_old_logquad.Rd @@ -0,0 +1,83 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/log_quad_augm.R +\name{HMD_old_logquad} +\alias{HMD_old_logquad} +\title{HMD pattern for adult ages.} +\usage{ +HMD_old_logquad( + nMx, + Age = NULL, + Sex = "b", + q0_5 = NULL, + q15_45 = NULL, + q60_15 = NULL, + Age_transition = 80, + window_transition = 3, + plot_comparison = FALSE, + fitted_logquad = NULL, + ... +) +} +\arguments{ +\item{nMx}{numeric. Vector of mortality rates in abridged age classes.} + +\item{Age}{integer. Single ages (abridged not allowed).} + +\item{Sex}{character. Either male \code{"m"}, female \code{"f"}, or both \code{"b"}.} + +\item{q0_5}{numeric. Probability of death from born to age 5. By default implicit values in \code{nMx} should be entered.} + +\item{q15_45}{numeric. Probability of death from age 15 to age 60. By default implicit values in \code{nMx} should be entered.} + +\item{q60_15}{numeric. Probability of death from age 60 to age 75. When external information on those ages level is available, +can be included to increase parameter \code{ax} from log-quad model in last ages (Li, 2003).} + +\item{Age_transition}{integer. Form which age should transition to HMD pattern starts.} + +\item{window_transition}{integer. Number of ages to the left and to the right of \code{Age_transition} to do a log-linear transition in rates.} + +\item{plot_comparison}{Show or not a plot with the result.} + +\item{fitted_logquad}{Optional, defaults to \code{NULL}. An object of class +\code{wilmoth}. If full HMD is not enough, one +can fit a Log-Quadratic (\url{https://github.com/mpascariu/MortalityEstimate}) model +based on any other collection of life tables;} + +\item{...}{Other arguments to be passed on to the \code{lt_single} function.} +} +\value{ +life table as in \code{lt_single} function. +} +\description{ +Adjust rates in oldest ages using HMD pattern, based on log-quad method. +} +\details{ +One possible scenario when mortality data on last ages is not reliable, is to use a mortality pattern with some known index in previous ages. +This function gives a HMD pattern based on 5q0 and 45q15, using log-quad model. Additionally, a value on mortality between 60 and 75 can be included to make a better adjustment in level. +} +\examples{ +# Mortality rates from UN Chilean with e0=70. Wat would be the rates based on HMD pattern? +# In this case making a transition of 10 years at age 80, and returning an OAG=100. +\dontrun{ +lt <- DemoToolsData::modelLTx1 +lt <- lt[lt$family == "Chilean" & lt$sex == "female" & lt$e0 == 70,] +chilean70_adjHMD <- HMD_old_logquad(nMx = lt$mx1, + Age = lt$age, + Sex = "f", + q0_5 = 1 - lt$lx1[lt$age==5]/lt$lx1[lt$age==0], + q15_45 = 1 - lt$lx1[lt$age==60]/lt$lx1[lt$age==15], + Age_transition = 80, + window_transition = 10, + plot_comparison = TRUE, + OAnew = 100) +# We know (as an example) that q60_15 is .5 higher than what HMD pattern would be. +chilean70_adjHMD_augm <- HMD_old_logquad(nMx = lt$mx1, + Age = lt$age, + Sex = "f", + q0_5 = 1 - lt$lx1[lt$age==5]/lt$lx1[lt$age==0], + q15_45 = 1 - lt$lx1[lt$age==60]/lt$lx1[lt$age==15], + q60_15 = (1 - lt$lx1[lt$age==75]/lt$lx1[lt$age==60]) * 1.5, + Age_transition = 80, window_transition = 10, + OAnew = 100, plot_comparison = TRUE) +} +} diff --git a/man/logquad_augmented.Rd b/man/logquad_augmented.Rd new file mode 100644 index 00000000..af9923a4 --- /dev/null +++ b/man/logquad_augmented.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/log_quad_augm.R +\name{logquad_augmented} +\alias{logquad_augmented} +\title{Augmented logquad} +\usage{ +logquad_augmented(coeffs, k, q0_5, q60_15, Sex = "b", Age, ...) +} +\arguments{ +\item{coeffs}{data.frame. Columns \code{a(x)}, \code{b(x)}, \code{c(x)} and \code{v(x)} from fitted logquad model. See \code{fitted_logquad_b}.} + +\item{k}{numeric. Adult mortality related value from log-quad estimatation based on one or two input parameters. See \code{lt_model_lq}.} + +\item{q0_5}{numeric. Probability of death from born to age 5.} + +\item{q60_15}{numeric. Probability of death from age 60 to age 75.} + +\item{Sex}{character. Either male \code{"m"}, female \code{"f"}, or both \code{"b"}.} + +\item{Age}{integer. Abridged lower bound ages. Same length than rows in \code{coeffs}.} + +\item{...}{Other arguments to be passed on to the \code{lt_abridged} function.} +} +\value{ +life table as in \code{lt_abridged} function. +} +\description{ +Adjust rates in oldest ages that comes from a HMD model, using an external estimate of 15q60 (Li, 2014). As an example see\code{\link[DemoTools]{HMD_old_logquad}}. +} +\details{ +Parameter \code{a(x)} is augmented based on an external estimate of 15q60. +} +\references{ +See \href{https://www.un.org/development/desa/pd/content/estimating-life-tables-developing-countries}{Li (2014)}. +}