added mean variant to summaries

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.

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 ([email protected]), Leontine Alkema ([email protected]), Peiran Liu ([email protected]), Adrian Raftery ([email protected]), Bailey Fosdick ([email protected]), Patrick Gerland ([email protected])
 Maintainer: Hana Sevcikova <[email protected]>

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

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) {

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),

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"

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

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}