diff --git a/ChangeLog b/ChangeLog index e7f8e12f..683d78e4 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,7 +1,9 @@ -7.4-2.9001 (2024-02-19) +7.4-2.90xx (2024-06-05) ----- Fixed bug related to a combination of ar.phase2 and missing data (thanks to Shunqi Cheng). +Added the mean variant into summary outputs. + 7.4-1/2 (2023-10-17) ----- Removed rgdal from suggested packages. diff --git a/DESCRIPTION b/DESCRIPTION index 1d9becfe..ea93a2cc 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: bayesTFR -Version: 7.4-2.9001 -Date: 2024-02-19 +Version: 7.4-2.9002 +Date: 2024-06-05 Title: Bayesian Fertility Projection Author: Hana Sevcikova (hanas@uw.edu), Leontine Alkema (alkema@nus.edu.sg), Peiran Liu (prliu@uw.edu), Adrian Raftery (raftery@uw.edu), Bailey Fosdick (bfosdick@uw.edu), Patrick Gerland (gerland@un.org) Maintainer: Hana Sevcikova diff --git a/R/get_outputs.R b/R/get_outputs.R index 10fd7ac6..1c998703 100644 --- a/R/get_outputs.R +++ b/R/get_outputs.R @@ -365,7 +365,9 @@ get.tfr.estimation <- function(mcmc.list=NULL, country=NULL, sim.dir=NULL, output <- list(tfr_table=tfr_table, country.obj = country.obj) if (!is.null(probs)) { - tfr_quantile <- apply(tfr_table, 2, quantile, probs = probs) + if(length(probs) == 1 && probs == "mean"){ + tfr_quantile <- apply(tfr_table, 2, mean) + } else tfr_quantile <- apply(tfr_table, 2, quantile, probs = probs) if (is.null(dim(tfr_quantile))) tfr_quantile <- matrix(tfr_quantile, nrow = 1) tfr_quantile <- data.table::data.table(t(tfr_quantile)) start.year <- mcmc.list$meta$start.year diff --git a/R/plot_functions.R b/R/plot_functions.R index f6e74785..5ca80e6a 100644 --- a/R/plot_functions.R +++ b/R/plot_functions.R @@ -369,6 +369,19 @@ get.median.from.prediction <- function(tfr.pred, country.index, country.code=NUL return(get.quantile.from.prediction(tfr.pred, quantile=0.5, country.index=country.index, country.code=country.code, adjusted=adjusted, ...)) } +get.mean.from.prediction <- function(tfr.pred, country.index, country.code=NULL, adjusted=TRUE, + est.uncertainty = FALSE) { + mean.values <- tfr.pred$traj.mean.sd[country.index, 1,] + if(est.uncertainty && has.est.uncertainty(tfr.pred$mcmc.set)){ # get the right value for present year + tfr.est <- get.tfr.estimation(mcmc.list=tfr.pred$mcmc.set, country = country.code, probs="mean", adjust = adjusted) + unc.last.time <- which(tfr.est$tfr_quantile$year == dimnames(tfr.pred$quantiles)[[3]][1]) + mean.values[1] <- unlist(tfr.est$tfr_quantile[unc.last.time, 1]) + } + if (!adjusted) return(mean.values) + shift <- get.tfr.shift(country.code, tfr.pred) + if(!is.null(shift)) mean.values <- mean.values + shift + return(mean.values) +} get.traj.quantiles <- function(tfr.pred, country.index, country.code, trajectories=NULL, pi=80, adjusted=TRUE, est.uncertainty = FALSE) { diff --git a/R/predict_tfr.R b/R/predict_tfr.R index 36881cac..51ab44c7 100644 --- a/R/predict_tfr.R +++ b/R/predict_tfr.R @@ -1118,12 +1118,12 @@ get.projection.summary.header.bayesTFR.prediction <- function(pred, ...) "get.UN.variant.names" <- function(pred, ...) UseMethod("get.UN.variant.names") get.UN.variant.names.bayesTFR.prediction <- function(pred, ...) - return(c('BHM median', 'BHM80 lower', 'BHM80 upper', 'BHM95 lower', 'BHM95 upper', 'Low', + return(c('BHM median', 'BHM80 lower', 'BHM80 upper', 'BHM95 lower', 'BHM95 upper', 'BHM mean', 'Low', 'High', 'Constant fertility')) "get.friendly.variant.names" <- function(pred, ...) UseMethod("get.friendly.variant.names") get.friendly.variant.names.bayesTFR.prediction <- function(pred, ...) - return(c('median', 'lower 80', 'upper 80', 'lower 95', 'upper 95', '-0.5child', '+0.5child', 'constant')) + return(c('median', 'lower 80', 'upper 80', 'lower 95', 'upper 95', 'mean', '-0.5child', '+0.5child', 'constant')) get.wpp.revision.number <- function(pred) { wpps <- c(2008, 2010, 2012, seq(2015, by = 2, length = 3), seq(2022, by = 2, length = 10)) @@ -1173,7 +1173,7 @@ do.write.projection.summary <- function(pred, output.dir, revision=NULL, indicat UN.variant.names <- get.UN.variant.names(pred) friendly.variant.names <- get.friendly.variant.names(pred) # the code is dependent on the following order of the variants (it's presumed): - # median, lower 80, upper 80, lower 95, upper 95, -1/2child, +1/2 child, constant + # median, lower 80, upper 80, lower 95, upper 95, mean, -1/2child, +1/2 child, constant nr.var <- length(UN.variant.names) result1 <- result2 <- NULL @@ -1184,7 +1184,10 @@ do.write.projection.summary <- function(pred, output.dir, revision=NULL, indicat if(est.uncertainty){ # add uncertainty est.object <- get.tfr.estimation(mcmc.list = pred$mcmc.set, country = country.obj$code, probs = c(0.5, 0.1, 0.9, 0.025, 0.975), adjust = adjusted) + est.object.mean <- get.tfr.estimation(mcmc.list = pred$mcmc.set, country = country.obj$code, + probs = "mean", adjust = adjusted) this.tfr.unc <- as.matrix(est.object$tfr_quantile)[1:ltfr, -ncol(est.object$tfr_quantile)] + this.tfr.unc <- cbind(this.tfr.unc, mean = unlist(est.object.mean$tfr_quantile[1:ltfr, 1])) #this.tfr <- unlist(est.object$tfr_quantile[,1])[1:ltfr] } this.result1 <- cbind( @@ -1194,15 +1197,18 @@ do.write.projection.summary <- function(pred, output.dir, revision=NULL, indicat # prediction median <- get.median.from.prediction(pred, country.obj$index, country.obj$code, adjusted=adjusted, est.uncertainty = est.uncertainty) + mean <- get.mean.from.prediction(pred, country.obj$index, country.obj$code, adjusted=adjusted, + est.uncertainty = est.uncertainty) proj.result <- rbind(median, get.traj.quantiles(pred, country.obj$index, country.obj$code, pi=80, adjusted=adjusted, est.uncertainty = est.uncertainty), get.traj.quantiles(pred, country.obj$index, country.obj$code, pi=95, adjusted=adjusted, - est.uncertainty = est.uncertainty)) + est.uncertainty = est.uncertainty), + mean) if(any(friendly.variant.names == '-0.5child')) proj.result <- rbind(proj.result, get.half.child.variant(median)) - #browser() + proj.result <- round(rbind(proj.result, rep(median[1], nr.proj)), precision) # constant variant if(est.uncertainty){ # user-friendly output contains observed years as well @@ -1219,7 +1225,6 @@ do.write.projection.summary <- function(pred, output.dir, revision=NULL, indicat result1 <- rbind(result1, this.result1) for(ivar in 1:nr.var) { this.obs.tfr <- if(est.uncertainty) obs.tfr[ivar,] else this.tfr - #browser() result2 <- rbind(result2, cbind(revision=rep(revision, nr.proj.all), variant=rep(e$UN_variants[e$UN_variants$Vshort==UN.variant.names[ivar],'VarID'], nr.proj.all), country=rep(country.obj$code, nr.proj.all), diff --git a/data-raw/UN_variants.txt b/data-raw/UN_variants.txt index 26602326..7bc15975 100644 --- a/data-raw/UN_variants.txt +++ b/data-raw/UN_variants.txt @@ -22,3 +22,4 @@ 13 204 "BHM80 lower" "BHM80 lower" "Standard" 13 206 "BHM95 upper" "BHM95 upper" "Standard" 13 207 "BHM95 lower" "BHM95 lower" "Standard" +13 208 "BHM mean" "BHM mean" "Standard" diff --git a/data/UN_variants.rda b/data/UN_variants.rda index af9da415..3399da64 100644 Binary files a/data/UN_variants.rda and b/data/UN_variants.rda differ diff --git a/man/get.tfr.estimation.Rd b/man/get.tfr.estimation.Rd index cad29870..760c5aa9 100644 --- a/man/get.tfr.estimation.Rd +++ b/man/get.tfr.estimation.Rd @@ -17,17 +17,17 @@ get.tfr.estimation(mcmc.list = NULL, country = NULL, \item{sim.dir}{Directory with the MCMC simulation results. Only used if \code{mcmc.list} is \code{NULL}.} \item{burnin}{Burn-in for getting trajectories and quantiles. A positive burn-in \eqn{x} will remove first \eqn{x} iterations from each chain.} \item{thin}{Thin for getting trajectories and quantiles. Thinning level \eqn{x} greater than 1 will store one iteration per \eqn{x} samples.} - \item{probs}{A vector of numbers between \eqn{[0,1]} specifying which estimation quantiles should be outputted. If it is set to \code{NULL} no quantiles are returned.} + \item{probs}{A vector of numbers between \eqn{[0,1]} specifying which estimation quantiles should be outputted. If it is set to \code{NULL} no quantiles are returned. It can be set to the word \dQuote{mean}, in which case the estimation mean is outputted.} \item{adjust}{Logical indicating whether the adjusted median and trajectories should be returned.} \item{country.code, ISO.code}{Deprecated arguments. Use argument \code{country} instead.} } \details{ -This function is used to obtain the TFR estimation trajectories as well as corresponding quantiles if the \code{mcmc.list} has been obtained while taking account for uncertainty about the past, i.e. \code{uncertainty=TRUE} in \code{\link{run.tfr.mcmc}}. Quantiles are included in the results if \code{probs} is not \code{NULL}. +This function is used to obtain the TFR estimation trajectories as well as corresponding quantiles if the \code{mcmc.list} has been obtained while taking account for uncertainty about the past, i.e. \code{uncertainty=TRUE} in \code{\link{run.tfr.mcmc}}. Quantiles or the mean are included in the results if \code{probs} is not \code{NULL}. } \value{ \item{tfr_table}{Table storing the trajectories. It is a matrix with rows equal to number of trajectories, and column equal to number of time periods.} \item{country.obj}{A list storing information about the country which the trajectories and quantiles correspond to. It corresponds to the output of \code{\link{get.country.object}}.} - \item{tfr_quantile}{Optional. A data.table object, storing the quantiles of estimates for each time period as specified by the \code{probs} argument. The time periods are contained in the \code{year} column.} + \item{tfr_quantile}{Optional. A data.table object, storing the quantiles or the mean of estimates for each time period as specified by the \code{probs} argument. The time periods are contained in the \code{year} column.} } \author{ Peiran Liu, Hana Sevcikova diff --git a/man/write.projection.summary.Rd b/man/write.projection.summary.Rd index 8580981e..4900bb17 100644 --- a/man/write.projection.summary.Rd +++ b/man/write.projection.summary.Rd @@ -1,12 +1,12 @@ \name{write.projection.summary} -\Rdversion{1.1} +%\Rdversion{1.1} \alias{write.projection.summary} %- Also NEED an '\alias' for EACH other topic documented here. \title{ Writing Projection Summary Files } \description{ -The function creates two files containing projection summaries, such as the median, the lower and upper bound of the 80 and 90\% probability intervals, respectively, the +/- 0.5 child variant and the constant variant. One file is in a user-friendly format, whereas the other is in a UN-specific format with internal coding of the time and the variants. In addition, a file containing some of the model parameters is created. +The function creates two files containing projection summaries, such as the median, mean, the lower and upper bound of the 80 and 90\% probability intervals, respectively, the +/- 0.5 child variant and the constant variant. One file is in a user-friendly format, whereas the other is in a UN-specific format with internal coding of the time and the variants. In addition, a file containing some of the model parameters is created. } \usage{ write.projection.summary(dir = file.path(getwd(), "bayesTFR.output"), @@ -24,7 +24,7 @@ write.projection.summary(dir = file.path(getwd(), "bayesTFR.output"), The first file that the function creates is called \file{projection_summary_user_friendly.csv} (or \file{projection_summary_user_friendly_adjusted.csv} if \code{adjusted=TRUE}), it is a comma-separated table with the following columns: \describe{\item{country_name}{country name} \item{country_code}{country code} - \item{variant}{name of the variant, such as \dQuote{median}, \dQuote{lower 80}, \dQuote{upper 80}, \dQuote{lower 95}, \dQuote{upper 95}, \dQuote{-0.5child}, \dQuote{+0.5child}, \dQuote{constant}} + \item{variant}{name of the variant, such as \dQuote{median}, \dQuote{lower 80}, \dQuote{upper 80}, \dQuote{lower 95}, \dQuote{upper 95}, \dQuote{mean}, \dQuote{-0.5child}, \dQuote{+0.5child}, \dQuote{constant}} \item{\emph{period1:}}{e.g. \dQuote{2005-2010}: TFR for the first time period. If \code{est.uncertainty} is \code{TRUE}, the first time period is the first observed time period. Otherwise it is the last observed one.} \item{\emph{period2:}}{e.g. \dQuote{2010-2015}: TFR for the second time period} \item{\dots}{further columns with TFR projections}