From be90a0222ff725270f9c3ec726d4d9aaa26fbb15 Mon Sep 17 00:00:00 2001 From: Hana Sevcikova Date: Fri, 11 Aug 2023 15:20:10 -0700 Subject: [PATCH 1/3] annual subnational projections --- R/project_subnat.R | 128 +++++++++++++++++++++++++++++++++------------ R/run_mcmc.R | 7 ++- R/wpp_data.R | 2 +- 3 files changed, 101 insertions(+), 36 deletions(-) diff --git a/R/project_subnat.R b/R/project_subnat.R index ac79e0e..b7f310b 100644 --- a/R/project_subnat.R +++ b/R/project_subnat.R @@ -1,7 +1,7 @@ e0.predict.subnat <- function(countries, my.e0.file, sim.dir=file.path(getwd(), 'bayesLife.output'), method = c("ar1", "shift", "scale"), predict.jmale = FALSE, my.e0M.file = NULL, - end.year=2100, start.year=NULL, output.dir = NULL, nr.traj=NULL, seed = NULL, - ar.pars = c(rho = 0.95, U = 82.5, a = 0.0482, b = -0.0154), + end.year=2100, start.year=NULL, output.dir = NULL, annual = NULL, + nr.traj=NULL, seed = NULL, ar.pars = NULL, save.as.ascii = 0, verbose = TRUE, jmale.estimates = NULL, ...) { # Run subnational projections, using the Scale AR(1) model applied to a national bayesLife simulation # sim.dir is the world-national simulation. Set output.dir to store results somewhere else. @@ -29,20 +29,39 @@ e0.predict.subnat <- function(countries, my.e0.file, sim.dir=file.path(getwd(), compute.alpha.shift <- function(...) return(compute.alpha.ar1(...)) - joint.male.est <- subnat.gap.estimates() - method <- match.arg(method) wpred <- get.e0.prediction(sim.dir) # contains national projections wdata <- wpred$e0.matrix.reconstructed wmeta <- wpred$mcmc.set$meta - year.step <- if(wmeta$annual.simulation) 1 else 5 if(!is.null(seed)) set.seed(seed) if (is.null(output.dir)) output.dir <- wmeta$output.dir quantiles.to.keep <- as.numeric(dimnames(wpred$quantiles)[[2]]) - e <- new.env() + + wannual <- wmeta$annual.simulation + if(is.null(annual)) annual <- wannual + year.step <- if(annual) 1 else 5 + wyear.step <- if(wannual) 1 else 5 + + joint.male.est <- subnat.gap.estimates(annual) + + ar.pars.default <- c(rho = 0.95, U = 82.5, a = 0.0482, b = -0.0154) + if(annual) { # Raftery's conversion from 5-year AR(1) parameters to 1-year parameters + ar.pars.default["a"] <- ar.pars.default["a"] * (1-ar.pars.default["rho"]^(2/5))/(1-ar.pars.default["rho"]^2) + ar.pars.default["b"] <- ar.pars.default["b"] * (1-ar.pars.default["rho"]^(2/5))/(1-ar.pars.default["rho"]^2) + ar.pars.default["rho"] <- ar.pars.default["rho"]^(1/5) + } + # Make sure all AR(1) parameters are available. If not, take the default values. + if(is.null(ar.pars)) ar.pars <- ar.pars.default + for(par in names(ar.pars.default)) if(! par %in% names(ar.pars)) ar.pars[par] <- ar.pars.default[par] + result <- list() orig.nr.traj <- nr.traj sex <- list(F = 'female', M = 'male') + + # set start year and present year of the national simulation + wstarty.global <- if(is.null(start.year)) as.integer(dimnames(wdata)[[1]][wpred$present.year.index]) + wyear.step else start.year + wpresenty.global <- wstarty.global - wyear.step + for (country in countries) { country.obj <- get.country.object(country, wmeta) if(is.null(country.obj$code)) { @@ -50,40 +69,71 @@ e0.predict.subnat <- function(countries, my.e0.file, sim.dir=file.path(getwd(), next } if(verbose) cat('\n', country.obj$name, ': predicting e0 for', sex[[wmeta$sex]]) - starty <- if(is.null(start.year)) as.integer(dimnames(wdata)[[1]][wpred$present.year.index]) + year.step else start.year + + we0 <- wdata[, country.obj$index] + we0obsy <- as.integer(names(we0))[-length(we0)] + + wstarty <- wstarty.global + wpresenty <- wpresenty.global + + starty <- wstarty # set subnational start year to the national start year (might get overwritten below) + presenty <- starty - year.step # subnational present year (might get overwritten below) + + if(!wannual){ # national prediction is 5-year + # snap start year and present year to be in the middle of the corresponding 5-years interval + seqy <- seq(min(we0obsy)-3, wpred$end.year, by=5) + wstarty <- seqy[cut(wstarty, seqy, labels=FALSE)] + 3 + wpresenty <- seqy[cut(wstarty, seqy, labels=FALSE)-1] + 3 + if(!annual){ # subnational prediction is 5-year + starty <- wstarty + presenty <- wpresenty + } + } else { # national prediction is annual + if(!annual){ # subnational prediction is 5-year + # adjust start year and present year to be in the middle of the corresponding 5-years interval + seqy <- seq(1700, wpred$end.year, by=5) + starty <- seqy[cut(wstarty, seqy, labels=FALSE)] + 3 + presenty <- starty - 5 + } + } + + # get various meta values for the subnational simulation meta <- e0.mcmc.meta.ini.subnat(wmeta, country = country.obj$code, my.e0.file = my.e0.file, - start.year = 1900, present.year = starty - year.step, verbose = verbose) - class(meta) <- "bayesLife.mcmc.meta" + start.year = 1900, present.year = presenty, annual = annual, verbose = verbose) + + if(is.null(meta$e0.matrix)) + stop("No data available. It might be a mismatch between the available observed years and start.year. ", + "The data file should contain year ", if(annual) presenty else paste(presenty - 3, presenty + 2, sep = "-"), + ". Or change the argument start.year.") + this.output.dir <- file.path(output.dir, paste0('subnat_', method), paste0('c', country.obj$code)) outdir <- file.path(this.output.dir, 'predictions') meta$output.dir <- this.output.dir + # extract national trajectories and thin them if needed wtrajs <- get.e0.trajectories(wpred, country.obj$code) nr.traj <- orig.nr.traj if(is.null(nr.traj)) nr.traj <- ncol(wtrajs) thinning.index <- round(seq(1, ncol(wtrajs), length = nr.traj)) wtrajs <- wtrajs[as.integer(rownames(wtrajs)) <= end.year, thinning.index] nr.traj <- ncol(wtrajs) - wyears <- as.integer(rownames(wtrajs)) - wend.year <- max(wyears) - we0 <- wdata[, country.obj$index] - we0obsy <- as.integer(names(we0))[-length(we0)] - if(!wmeta$annual.simulation){ - seqy <- seq(min(we0obsy) - 3, max(wyears) + 2, by = 5) - midy <- seqy + 3 - } else{ - seqy <- midy <-(min(we0obsy)-1):max(wyears) - } - presenty <- midy[cut(starty, seqy, labels=FALSE)-1] - if(any(wyears < presenty)) # remove time periods from national trajectories before present year - wtrajs <- wtrajs[-which(wyears < presenty),] - if(presenty < min(wyears)) { # add observed data to national trajectories if present year is not there - adde0 <- we0[-length(we0)][we0obsy >= presenty] - adddata <- matrix(adde0, ncol=nr.traj, nrow=length(adde0), - dimnames=list(names(adde0), colnames(wtrajs))) - wtrajs <- rbind(adddata, wtrajs) + + # attach observed data to the trajectories, so that imputation can be done using all data + adde0 <- matrix(we0, ncol = nr.traj, nrow = length(we0), + dimnames = list(names(we0), colnames(wtrajs))) + wtrajs.all <- rbind(adde0[! rownames(adde0) %in% rownames(wtrajs),], wtrajs) + + if(annual && !wannual) { # interpolate national trajectories + wtrajs.all <- bayesTFR:::interpolate.trajectories(wtrajs.all) + } else { + if(!annual && wannual) # average annual national trajectories + wtrajs.all <- bayesTFR:::average.trajectories(wtrajs.all) } + # remove everything from the national trajectories before present year + yrs <- as.integer(rownames(wtrajs.all)) + wtrajs <- wtrajs.all[-which(yrs < presenty),] + if(!presenty %in% rownames(meta$e0.matrix)) { # add NAs to e0 matrix addyears <- sort(seq(presenty, max(as.integer(rownames(meta$e0.matrix)) + year.step), by=-year.step)) adddata <- matrix(NA, nrow=length(addyears), ncol=ncol(meta$e0.matrix), @@ -97,6 +147,7 @@ e0.predict.subnat <- function(countries, my.e0.file, sim.dir=file.path(getwd(), bayesLife.mcmc.meta <- meta store.bayesLife.meta.object(bayesLife.mcmc.meta, this.output.dir) + # prepare for subnational simulation this.nr.project <- nrow(wtrajs) - 1 nr.reg <- get.nr.countries(meta) PIs_cqp <- array(NA, c(nr.reg, length(quantiles.to.keep), nrow(wtrajs)), @@ -105,24 +156,26 @@ e0.predict.subnat <- function(countries, my.e0.file, sim.dir=file.path(getwd(), #meta$Tc.index <- .get.Tcindex(meta$e0.matrix, cnames = meta$regions$country_name) country.char <- as.character(country.obj$code) e0reconstructed <- meta$e0.matrix + + # iterate over subnational units for(region in 1:nr.reg) { reg.obj <- get.country.object(region, meta, index=TRUE) regcode.char <- as.character(reg.obj$code) rege0 <- get.observed.e0(region, meta, 'e0.matrix') rege0.last <- rep(rege0[length(rege0)], nr.traj) c.first <- do.call(paste0("compute.alpha.", method), - list(rege0.last, we0[names(we0) %in% names(rege0.last)])) # set of initial scales + list(rege0.last, wtrajs.all[rownames(wtrajs.all) %in% names(rege0.last)])) # set of initial scales #stop("") if(is.na(rege0.last[1])) { # impute where NA's at the end for(i in length(rege0):1) if(!is.na(rege0[i])) break - widx <- which(names(we0) %in% names(rege0[i])) - c.first <- rep(do.call(paste0("compute.alpha.", method), list(rege0[i], we0[widx])), + widx <- which(rownames(wtrajs.all) %in% names(rege0[i])) + c.first <- rep(do.call(paste0("compute.alpha.", method), list(rege0[i], wtrajs.all[widx,])), nr.traj) # set of initial scales meta$Tc.index[region] <- i imptraj <- matrix(NA, nrow = length(rege0) - i, ncol = nr.traj) # trajectory matrix for imputation for(tr in 1:nr.traj) { # iterate over trajectories imp.time <- i:(length(rege0)-1) - natval <- we0[widx + (imp.time - i + 1)] + natval <- wtrajs.all[widx + (imp.time - i + 1), tr] impval <- do.call(paste0("generate.trajectory.", method), list(natval, c.first[tr])) imptraj[imp.time - i + 1, tr] <- impval @@ -135,6 +188,7 @@ e0.predict.subnat <- function(countries, my.e0.file, sim.dir=file.path(getwd(), rege0.last <- imptraj[nrow(imptraj),] } # end of imputation + # initiate matrix for holding the projections e0.pred <- matrix(NA, nrow = this.nr.project + 1, ncol = nr.traj, dimnames = list(rownames(wtrajs), NULL)) e0.pred[1,] <- rege0.last @@ -143,7 +197,9 @@ e0.predict.subnat <- function(countries, my.e0.file, sim.dir=file.path(getwd(), e0.pred[-1, s] <- do.call(paste0("generate.trajectory.", method), list(wtrajs[-1,s], c.first[s])) } + # save resulting trajectories trajectories <- e0.pred + rownames(trajectories) <- rownames(wtrajs) save(trajectories, file = file.path(outdir, paste0('traj_country', reg.obj$code, '.rda'))) # compute quantiles @@ -161,7 +217,7 @@ e0.predict.subnat <- function(countries, my.e0.file, sim.dir=file.path(getwd(), output.directory = normalizePath(outdir), mcmc.set=list(meta=meta, mcmc.list=list()), nr.projections=this.nr.project, - burnin=NA, end.year=wend.year, start.year=starty, + burnin=NA, end.year = max(as.integer(rownames(wtrajs))), start.year=starty, method = method, ar.pars = ar.pars, present.year.index=present.year.index, present.year.index.all=present.year.index, @@ -229,8 +285,14 @@ e0.jmale.predict.subnat <- function(e0.pred, estimates = NULL, gap.lim = c(0,18) invisible(bayesLife.prediction) } -subnat.gap.estimates <- function() +subnat.gap.estimates <- function(annual = FALSE){ + if(annual) + return(list(eq1 = list(coefficients = c(-0.2572068905, e0.1953 = -0.0002524553, Gprev = 0.9914796728, e0 = 0.0049163005, e0d75 = -0.0193942510), + sigma = 0.1723969, dof = 2.311855), + eq2 = list(coefficients = c(Gprev = 1), sigma = 0.3429213, dof = NULL))) + # 5-year estimates return(list(eq1 = list(coefficients = c(-1.076940553, e0.1953 = 0.002739663, Gprev = 0.942276337, e0 = 0.019435863, e0d75 = -0.099589171), sigma = 0.3902612, dof = 4.058482), eq2 = list(coefficients = c(Gprev = 1), sigma = 0.4296343, dof = NULL)) ) +} \ No newline at end of file diff --git a/R/run_mcmc.R b/R/run_mcmc.R index 5956cd8..c58e761 100644 --- a/R/run_mcmc.R +++ b/R/run_mcmc.R @@ -711,11 +711,14 @@ e0.mcmc.ini.extra <- function(mcmc, countries, index.replace = NULL) { } e0.mcmc.meta.ini.subnat <- function(meta, country, start.year = 1950, present.year = 2020, - my.e0.file = NULL, annual.simulation = FALSE, verbose = FALSE) { + my.e0.file = NULL, annual = NULL, verbose = FALSE) { meta$start.year <- start.year meta$present.year <- present.year + if(is.null(annual)) annual <- meta$annual.simulation + meta$annual.simulation <- annual + data <- get.wpp.e0.subnat(country, start.year = start.year, present.year = present.year, - my.e0.file = my.e0.file, annual = annual.simulation) + my.e0.file = my.e0.file, annual = annual) data$nr.countries <- ncol(data$e0.matrix) data$Tc.index <- .get.Tcindex(data$e0.matrix, cnames=data$regions$country_name, stop.if.less.than2 = FALSE) data$subnat <- TRUE diff --git a/R/wpp_data.R b/R/wpp_data.R index 6354be2..5cef579 100644 --- a/R/wpp_data.R +++ b/R/wpp_data.R @@ -203,7 +203,7 @@ get.wpp.e0.subnat.joint <- function(country, meta, my.e0.file) { LEXmatrix.regions <- bayesTFR:::get.observed.time.matrix.and.regions( data_countries, loc_data, start.year = meta$start.year, - present.year = meta$present.year, annual = meta$annual, + present.year = meta$present.year, annual = meta$annual.simulation, datacolnames=c(country.code='reg_code', country.name='name', reg.name='reg_name', reg.code='NA', area.name='country', area.code='country_code')) From dcd51343ca49b78e7fd2c9de5561476dcb6142a8 Mon Sep 17 00:00:00 2001 From: Hana Sevcikova Date: Fri, 15 Sep 2023 12:34:28 -0700 Subject: [PATCH 2/3] added data check --- R/wpp_data.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/wpp_data.R b/R/wpp_data.R index 5cef579..da27e38 100644 --- a/R/wpp_data.R +++ b/R/wpp_data.R @@ -171,7 +171,8 @@ get.wpp.e0.subnat <- function(country, start.year=1950, present.year=2010, my.e0 data_countries <- data[include | locations$prediction.only,] nr_countries_estimation <- sum(include) - #stop("") + if(nrow(data_countries) == 0) + stop("No data for country ", country, " available.") LEXmatrix.regions <- bayesTFR:::get.observed.time.matrix.and.regions( data_countries, loc_data, start.year = start.year, present.year = present.year, annual = annual, From bb2caac3dbed03e54a4d355bd8311afa5c6b3837 Mon Sep 17 00:00:00 2001 From: Hana Sevcikova Date: Fri, 15 Sep 2023 14:18:20 -0700 Subject: [PATCH 3/3] edited docs for subnational projection --- ChangeLog | 5 +++++ DESCRIPTION | 4 ++-- man/bayesLife-package.Rd | 4 ++-- man/e0.predict.subnat.Rd | 25 +++++++++++++++++-------- tests/test_functions.R | 38 +++++++++++++++++++++++++++++++++++++- 5 files changed, 63 insertions(+), 13 deletions(-) diff --git a/ChangeLog b/ChangeLog index 062044d..0ebb032 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,8 @@ +5.2-0 (09/15/2023) +----- +Annual subnational projections are now possible, via the argument "annual" +in e0.predict.subnat(). + 5.1-1 (02/04/2023) ----- Allow use of annual supplemental data. diff --git a/DESCRIPTION b/DESCRIPTION index d8c7b91..fb288ce 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: bayesLife Type: Package Title: Bayesian Projection of Life Expectancy -Version: 5.1-1 -Date: 2023-02-04 +Version: 5.2-0 +Date: 2023-09-15 Author: Hana Sevcikova, Adrian Raftery, Jennifer Chunn Maintainer: Hana Sevcikova Description: Making probabilistic projections of life expectancy for all countries of the world, using a Bayesian hierarchical model . Subnational projections are also supported. diff --git a/man/bayesLife-package.Rd b/man/bayesLife-package.Rd index bf88242..25c0db1 100644 --- a/man/bayesLife-package.Rd +++ b/man/bayesLife-package.Rd @@ -35,7 +35,7 @@ Existing simulation results can be accessed using the \link{get.e0.mcmc} functio For a table with countries included in the mcmc or prediction object, the function \link[bayesTFR]{get.countries.table} can be used in the same way as in \pkg{bayesTFR}. -Historical data are taken from one of the packages \pkg{wpp2019} (default), \pkg{wpp2017}, \pkg{wpp2015}, \pkg{wpp2012} or \pkg{wpp2010}, depending on users settings. +Historical data are taken from one of the packages \pkg{wpp2019} (default), \pkg{wpp2017}, \pkg{wpp2015}, \pkg{wpp2012} or \pkg{wpp2010}, depending on users settings. For more recent data, package \pkg{wpp2022} can be installed from Github (@PPgp). } @@ -57,7 +57,7 @@ A. E. Raftery, N. Li, H. Sevcikova, P. Gerland, G. K. Heilig (2012). Bayesian p A. E. Raftery, N. Lalic, P. Gerland (2014). Joint Probabilistic Projection of Female and Male Life Expectancy. Demographic Research, 30:795-822. -H. Sevcikova, A. E. Raftery (2021). Probabilistic Projection of Subnational Life Expectancy. Journal of Official Statistics, in press. +H. Sevcikova, A. E. Raftery (2021). Probabilistic Projection of Subnational Life Expectancy. Journal of Official Statistics, , Vol. 37, no. 3, 591-610. %Lalic, N. (2011). Master's thesis at the Department of Statistics, University of Washington. } diff --git a/man/e0.predict.subnat.Rd b/man/e0.predict.subnat.Rd index fb006c2..8056cb2 100644 --- a/man/e0.predict.subnat.Rd +++ b/man/e0.predict.subnat.Rd @@ -15,15 +15,15 @@ e0.predict.subnat(countries, my.e0.file, method = c("ar1", "shift", "scale"), predict.jmale = FALSE, my.e0M.file = NULL, end.year = 2100, start.year = NULL, output.dir = NULL, - nr.traj = NULL, seed = NULL, - ar.pars = c(rho = 0.95, U = 82.5, a = 0.0482, b = -0.0154), - save.as.ascii = 0, verbose = TRUE, jmale.estimates = NULL, \dots) + annual = NULL, nr.traj = NULL, seed = NULL, + ar.pars = NULL, save.as.ascii = 0, verbose = TRUE, + jmale.estimates = NULL, \dots) e0.jmale.predict.subnat(e0.pred, estimates = NULL, gap.lim = c(0,18), max.e0.eq1.pred = 86, my.e0.file = NULL, save.as.ascii = 0, verbose = TRUE) -subnat.gap.estimates() +subnat.gap.estimates(annual = FALSE) } \arguments{ @@ -36,9 +36,15 @@ subnat.gap.estimates() \item{end.year}{End year of the projections.} \item{start.year}{Start year of the projections. By default, projections start at the same time point as the national projections.} \item{output.dir}{Directory into which the resulting prediction objects and the trajectories are stored. See below for details.} + \item{annual}{Logical indicating if the subnational projection should be on an annual scale or a 5-year scale. By default, +the scale is matched to the national simulation given by \code{sim.dir}. If given, the scale must match to the scale of the subnational data provided in \code{my.e0.file}. +If the subnational and national scales are not the same, +the national trajectories are either interpolated (if \code{annual = TRUE} and the national simulation is not annual) or averaged to 5-year values +(if \code{annual = FALSE} and the national simulation is annual).} \item{nr.traj}{Number of trajectories to be generated. If \code{NULL}, the number of trajectories in the national projections is used.} \item{seed}{Seed of the random number generator. If \code{NULL} no seed is set. It can be used to generate reproducible projections.} -\item{ar.pars}{List containing the parameter estimates for the AR(1) method (i.e. if \code{method = "ar1"}, default). It must have elements called \code{rho}, \code{U}, \code{a} and \code{b}. See the reference paper for details on the estimation.} +\item{ar.pars}{Named vector containing the parameter estimates for the AR(1) method (i.e. if \code{method = "ar1"}, default). If given, it must have elements called \code{rho}, \code{U}, \code{a} and \code{b}. See the reference paper for details on the estimation. By default for a 5-year simulation, \code{c(rho = 0.95, U = 82.5, a = 0.0482, b = -0.0154)} is used. +For an annual simulation these default parameters are scaled to \code{c(rho = 0.9898, U = 82.5, a = 0.01, b = -0.0032)}, see details below.} \item{save.as.ascii}{Either a number determining how many trajectories should be converted into an ASCII file, or \dQuote{all} in which case all trajectories are converted. By default no conversion is performed.} \item{verbose}{Logical switching log messages on and off.} \item{jmale.estimates, estimates}{List with estimates for the female-male gap model. The default values are retrieved using the function \code{subnat.gap.estimates()}.} @@ -50,11 +56,14 @@ subnat.gap.estimates() \details{ The \code{e0.predict.subnat} function implements the methodology described in Sevcikova and Raftery (2021). Given a set of national bayesLife projections, it applies one of the methods (AR(1), Shift or Scale) to each national trajectory and each subregion of given countries which yields subnational e0 projections. -The file on subnational data passed into \code{my.e0.file} and \code{my.e0M.file} has to have a column \dQuote{country_code} with numerical values corresponding to countries given in the argument \code{countries}, and column \dQuote{reg_code} giving the numerical identifier of each subregion. Column \dQuote{name} should be used for subregion name, and column \dQuote{country_name} for country name. An optional column \dQuote{include_code} can be used to eliminate entries from processing. Entries with values of 1 or 2 will be included, all others will be ignored. Column \dQuote{last.observed} can be used to define which time period contains the last observed data point (given as integer, e.g. year in the middle of the time period). Remaining columns define the time periods, e.g. \dQuote{2000-2005}, \dQuote{2005-2010}. The package contains an example of such dataset, see Example below. +The file on subnational data passed into \code{my.e0.file} and \code{my.e0M.file} has to have a column \dQuote{country_code} with numerical values corresponding to countries given in the argument \code{countries}, and column \dQuote{reg_code} giving the numerical identifier of each subregion. Column \dQuote{name} should be used for subregion name, and column \dQuote{country_name} for country name. An optional column \dQuote{include_code} can be used to eliminate entries from processing. Entries with values of 1 or 2 will be included, all others will be ignored. Column \dQuote{last.observed} can be used to define which time period contains the last observed data point (given as integer, e.g. year in the middle of the time period). Remaining columns define the time periods, e.g. \dQuote{2000-2005}, \dQuote{2005-2010} for a 5-year simulation, or \dQuote{2020}, \dQuote{2021} for an annual simulation. The package contains an example of such dataset, see Example below. + +The default AR(1) parameters for the \dQuote{ar1} method were designed for a 5-year simulation, see Sevcikova & Raftery (2021) for more detail. These are \eqn{\rho = 0.95, U = 82.5, a = 0.0482, b = -0.0154}. If an annual AR(1) process is desired, we use the following conversion for the autoregressive parameter \eqn{\rho} and the \eqn{a} and \eqn{b} parameters: +\eqn{\rho^* = \rho^{(1/5)}, a^* = a * \sqrt{((1 - \rho^{(2/5)})/(1 - \rho^2))}, b^* = b * \sqrt{((1 - \rho^{(2/5)})/(1 - \rho^2))}}. The \eqn{U} parameter stays the same for both processes. Thus, the annual parameters are \code{c(rho = 0.9898, U = 82.5, a = 0.01, b = -0.0032)}. Note that if the \code{ar.pars} argument is specified by the user, it is assumed that the parameters have been scaled appropriately and thus, no conversion takes place. Argument \code{output.dir} gives a location on disk where results of the function should be stored. If it is \code{NULL} (default), results are stored in the same directory as the national projections. In both cases a subdirectory called \dQuote{subnat_\code{method}} is created in which each country has its own subfolder with the country code in its name. Each such subfolder contains the same type of outputs as in the national case generated using \code{\link{e0.predict}}, most importantly a directory \dQuote{predictions} with trajectories for each region. -If the argument \code{predict.jmale} is \code{TRUE}, the \code{e0.predict.subnat} invokes the \code{e0.jmale.predict.subnat} function for each country. However, one can call the \code{e0.jmale.predict.subnat} function explicitly. It applies the female-male gap model to regions of one country. See \code{\link{e0.jmale.predict}} for more detail on the model. The default covariates of the model are not estimated on the fly. They were estimated externally using subnational data for about 30 countries and can be viewed using \code{subnat.gap.estimates()}. +If the argument \code{predict.jmale} is \code{TRUE}, the \code{e0.predict.subnat} invokes the \code{e0.jmale.predict.subnat} function for each country. However, one can call the \code{e0.jmale.predict.subnat} function explicitly. It applies the female-male gap model to regions of one country. See \code{\link{e0.jmale.predict}} for more detail on the model. The default covariates of the model are not estimated on the fly. They were estimated externally using subnational data for about 30 countries and can be viewed using \code{subnat.gap.estimates()}, either for estimates derived from 5-year data (default) or annual data (\code{annual = TRUE}). } \value{ @@ -68,7 +77,7 @@ Even though the resulting object contains subnational results, the names of its } \references{ -H. Sevcikova, A. E. Raftery (2021). Probabilistic Projection of Subnational Life Expectancy. Journal of Official Statistics, in press. +H. Sevcikova, A. E. Raftery (2021). Probabilistic Projection of Subnational Life Expectancy. Journal of Official Statistics, Vol. 37, no. 3, 591-610. } \author{ diff --git a/tests/test_functions.R b/tests/test_functions.R index f91f767..0eb0936 100644 --- a/tests/test_functions.R +++ b/tests/test_functions.R @@ -720,8 +720,44 @@ test.subnational.predictions <- function() { trajs <- get.e0.trajectories(preds[["36"]], "Victoria") stopifnot(all(dim(trajs) == c(7, 30))) - test.ok(test.name) unlink(sim.dir, recursive=TRUE) + + # Annual subnational projections + my.sub.file.annual.F <- tempfile() + my.sub.file.annual.M <- tempfile() + # female data + datF <- data.table::fread(my.sub.file) # load the data and save only the last time period, pretending this is a single year data point + datF <- datF[, .(country_code, country_name, reg_code, name, `2013` = `2010-2015`)] + data.table::fwrite(datF, file = my.sub.file.annual.F, sep = "\t") # save it + # male data + datM <- data.table::copy(datF)[, `2013` := `2013`*0.95] + data.table::fwrite(datM, file = my.sub.file.annual.M, sep = "\t") + + sim.dir <- tempfile() + preds <- e0.predict.subnat(c(36, 124, 360), my.e0.file=my.sub.file.annual.F, sim.dir=nat.dir, + output.dir=sim.dir, start.year = 2016, # will impute 2 years + end.year=2030, annual = TRUE, + predict.jmale = TRUE, my.e0M.file = my.sub.file.annual.M) + # Retrieve trajectories + trajs <- get.e0.trajectories(preds[["124"]], "Alberta") + stopifnot(all(dim(trajs) == c(16, 30))) # from 2015 to 2030 + stopifnot(all(c(2015, 2030) %in% as.integer(rownames(trajs)))) + + pred <- get.rege0.prediction(sim.dir, 360) + spred <- summary(pred, "Papua") + stopifnot(max(spred$projection.years) == 2030) + + predM <- get.e0.jmale.prediction(pred) + t <- tfr.trajectories.table(predM, "Bali") + stopifnot(max(as.integer(rownames(t))) == 2030) + stopifnot(all(c(2015, 2024, 2029) %in% as.integer(rownames(t)))) + + unlink(sim.dir, recursive=TRUE) + unlink(my.sub.file.annual.F) + unlink(my.sub.file.annual.M) + + test.ok(test.name) + } test.run.annual.simulation <- function(wpp.year = 2019) {