merged with master

.Rbuildignore

@@ -2,3 +2,6 @@
 ^README\.md$
 ^.github$
 ^data-raw
+^tests/last.dump.rda
+^tests/.clustersize
+^tests/.proc

ChangeLog

@@ -1,3 +1,23 @@
+5.3-0/1 (11/05/2024)
+-----
+Added dataset include_2024.
+
+Default for par.names and par.names.cs in e0.raftery.diag changed to all parameters.
+
+Fixed bug in e0.joint.plot for annual prediction object.
+
+Fixed bug in setting time index when imputation is present in e0.predict.subnat().
+
+Added argument use.wpp.data to run.e0.mcmc(). 
+
+Added the mean variant into summary output.
+
+e0.trajectories.plot() got arguments to plot the means (show.mean), the medians (show.median),
+as well as selected trajectories (traj.index).
+
+Added option of keeping multiple predictions directories in one simulation directory
+(argument "subdir" added to various functions).
+
 5.2-0 (09/15/2023)
 -----
 Annual subnational projections are now possible, via the argument "annual"

DESCRIPTION

@@ -1,10 +1,20 @@
 Package: bayesLife
 Type: Package
 Title: Bayesian Projection of Life Expectancy
-Version: 5.2-0
-Date: 2023-09-15
-Author: Hana Sevcikova, Adrian Raftery, Jennifer Chunn 
-Maintainer: Hana Sevcikova <[email protected]>
+Version: 5.3-1
+Date: 2024-12-05
+Authors@R: c(person(given = "Hana",
+                    family = "Sevcikova",
+                    role = c("cre", "aut"),
+                    email = "[email protected]"),
+             person(given = "Adrian",
+                      family = "Raftery",
+                      role = "aut",
+                      email = "[email protected]"),
+             person(given = "Jennifer",
+                      family = "Chunn",
+                      role = "aut")
+            )
 Description: Making probabilistic projections of life expectancy for all countries of the world, using a Bayesian hierarchical model <doi:10.1007/s13524-012-0193-x>. Subnational projections are also supported.
 Depends: 
     bayesTFR (>= 7.3-0),

NAMESPACE

@@ -36,6 +36,7 @@ export(
 	create.thinned.e0.mcmc,
 	get.e0.prediction,
 	has.e0.prediction,
+	available.e0.predictions,
 	get.e0.jmale.prediction,
 	has.e0.jmale.prediction,
 	get.rege0.prediction,

R/diagnostics.R

@@ -1,7 +1,7 @@
 e0.raftery.diag <- function(mcmc = NULL, 
 							 sim.dir = file.path(getwd(), 'bayesLife.output'),
 							 burnin = 0, country = NULL,
-							 par.names = NULL, par.names.cs = NULL,
+							 par.names = NA, par.names.cs = NA,
 							 country.sampling.prop = 1,
 							 verbose = TRUE, ...) {
     mcmc.set <- if (is.null(mcmc)) get.e0.mcmc(sim.dir = sim.dir, low.memory = TRUE) else mcmc

R/get_outputs.R

@@ -12,6 +12,7 @@ get.e0.mcmc <- function(sim.dir = file.path(getwd(), 'bayesLife.output'),
 	load(file = mcmc.file.path)
 	bayesLife.mcmc.meta$output.dir <- normalizePath(sim.dir)
 	if(is.null(bayesLife.mcmc.meta$annual.simulation)) bayesLife.mcmc.meta$annual.simulation <- FALSE
+	if(is.null(bayesLife.mcmc.meta$use.wpp.data)) bayesLife.mcmc.meta$use.wpp.data <- TRUE
 	if(is.null(bayesLife.mcmc.meta$mcmc.options)) {
 	    bayesLife.mcmc.meta <- .convert.meta.from.legacy.form(bayesLife.mcmc.meta)
 	}
@@ -92,14 +93,19 @@ e0.mcmc <- function(mcmc.set, chain.id=1) return (mcmc.set$mcmc.list[[chain.id]]
 e0.mcmc.list <- function(mcmc.set, chain.ids=NULL) 
 	return(bayesTFR::tfr.mcmc.list(mcmc.set=mcmc.set, chain.ids=chain.ids))
 
-has.e0.prediction <- function(mcmc=NULL, sim.dir=NULL) {
+has.e0.prediction <- function(mcmc=NULL, sim.dir=NULL, subdir = "predictions") {
 	if (!is.null(mcmc)) sim.dir <- if(is.character(mcmc)) mcmc else mcmc$meta$output.dir
 	if (is.null(sim.dir)) stop('Either mcmc or directory must be given.')
-	if(file.exists(file.path(sim.dir, 'predictions', 'prediction.rda'))) return(TRUE)
+	if(file.exists(file.path(sim.dir, subdir, 'prediction.rda'))) return(TRUE)
 	return(FALSE)
 }
 
-get.e0.prediction <- function(mcmc=NULL, sim.dir=NULL, joint.male=FALSE, mcmc.dir=NULL) {
+available.e0.predictions <- function(mcmc=NULL, sim.dir=NULL, full.names = FALSE){
+    return(bayesTFR::available.tfr.predictions(mcmc=mcmc, sim.dir = sim.dir, full.names = full.names))
+}
+
+get.e0.prediction <- function(mcmc=NULL, sim.dir=NULL, joint.male=FALSE, mcmc.dir=NULL, 
+                              subdir = "predictions") {
 	############
 	# Returns an object of class bayesLife.prediction
 	# Set mcmc.dir to NA, if the prediction object should not have a pointer 
@@ -108,11 +114,15 @@ get.e0.prediction <- function(mcmc=NULL, sim.dir=NULL, joint.male=FALSE, mcmc.di
 	if (!is.null(mcmc)) 
 		sim.dir <- if(is.character(mcmc)) mcmc else mcmc$meta$output.dir
 	if (is.null(sim.dir)) stop('Either mcmc or directory must be given.')
-	output.dir <- file.path(sim.dir, 'predictions')
+	output.dir <- file.path(sim.dir, subdir)
 	pred.file <- file.path(output.dir, 'prediction.rda')
 	if(!file.exists(pred.file)) {
 		warning('File ', pred.file, ' does not exist.')
-		return(NULL)
+	    if(length((alt.preds <- available.e0.predictions(sim.dir = sim.dir))) > 0){
+	        output.dir <- file.path(sim.dir, alt.preds[1])
+	        pred.file <- file.path(output.dir, 'prediction.rda')
+	        warning('Extracting predictions from ', alt.preds[1])
+	    } else return(NULL)
 	}
 	load(file=pred.file)
 	bayesLife.prediction$output.directory <- output.dir

R/plot_functions.R

@@ -106,8 +106,13 @@ e0.joint.plot <- function(e0.pred, country, pi=95, years, nr.points=500,
 	if(!has.e0.jmale.prediction(e0.pred)) 
 		stop('A male prediction does not exist for the given prediction object. Run e0.jmale.predict.')
 	start.year <- as.integer(dimnames(e0.pred$quantiles)[[3]][1])
-	years.obs <- years[years <= start.year+2]
-	years.pred <- years[years > start.year+2]
+	if(e0.pred$mcmc.set$meta$annual.simulation){
+	    years.obs <- years[years <= start.year]
+	    years.pred <- years[years > start.year]
+	} else {
+	    years.obs <- years[years <= start.year+2]
+	    years.pred <- years[years > start.year+2]
+	}
 	years.idx <- unlist(lapply(years.pred, bayesTFR:::get.prediction.year.index, pred=e0.pred))
 	years.idx <- years.idx[years.idx > 1]
 	years.obs.idx <- unlist(lapply(years.obs, bayesTFR:::get.estimation.year.index, meta=e0.pred$mcmc.set$meta))
@@ -195,6 +200,7 @@ e0.trajectories.plot.all <- function(e0.pred,
 
 e0.trajectories.plot <- function(e0.pred, country, pi=c(80, 95), both.sexes=FALSE,
 								  nr.traj=NULL, adjusted.only = TRUE, typical.trajectory=FALSE,
+								  traj.index = NULL, show.mean = FALSE, show.median = TRUE,
 								  xlim=NULL, ylim=NULL, type='b', 
 								  xlab='Year', ylab='Life expectancy at birth', main=NULL, 
 								  lwd=c(2,2,2,2,1), col=c('black', 'green', 'red', 'red', '#00000020'),
@@ -320,24 +326,44 @@ e0.trajectories.plot <- function(e0.pred, country, pi=c(80, 95), both.sexes=FALS
 		e0pred <- pred[[ipred]]
 		this.col <- plotcols[[ipred]]
 		meta <- e0pred$mcmc.set$meta
+		if(!is.null(traj.index)) nr.traj <- length(traj.index)
+		e0.median <- e0.mean <- e0.main.proj <- NULL
 		if(do.average) {
 			trajectories <- get.e0.trajectories.object(pred, country$code, nr.traj=nr.traj, typical.trajectory=typical.trajectory, pi=pi)
-			e0.median <- trajectories$median
+			if(show.median)
+			    e0.median <- trajectories$median
+			if(show.mean)
+			    e0.mean <- apply(trajectories$trajectories, 1, mean, na.rm=TRUE)
 		} else {
 			trajectories <- get.e0.trajectories.object(e0pred, country$code, nr.traj=nr.traj, typical.trajectory=typical.trajectory)
-			e0.median <- bayesTFR::get.median.from.prediction(e0pred, country$index, country$code)
+			if(show.median)
+			    e0.median <- bayesTFR::get.median.from.prediction(e0pred, country$index, country$code)
+			if(show.mean)
+			    e0.mean <- bayesTFR::get.mean.from.prediction(e0pred, country$index, country$code)
+		}
+		if(!is.null(traj.index) && !is.null(trajectories$trajectories)) trajectories$index <- traj.index
+		# set the main projection (solid line)
+		main.proj.name <- ""
+		if(!is.null(e0.median)){
+		    e0.main.proj <- e0.median
+		    main.proj.name <- "median"
+		} else {
+		    if(!is.null(e0.mean)){
+		        e0.main.proj <- e0.mean
+		        main.proj.name <- "mean"
+		    }
 		}
 		cqp <- list()
 		if(ipred > 1) add <- TRUE
 		if(!add)
 			ylim.loc <- c(min(if (!is.null(trajectories$trajectories))
 							trajectories$trajectories[,trajectories$index]
 					  	else NULL, 
-					  	ylim.loc[1], e0.median, na.rm=TRUE), 
+					  	ylim.loc[1], e0.main.proj, na.rm=TRUE), 
 				  	max(if (!is.null(trajectories$trajectories))
 							trajectories$trajectories[,trajectories$index]
 					  else NULL, 
-					  ylim.loc[2], e0.median, na.rm=TRUE))
+					  ylim.loc[2], e0.main.proj, na.rm=TRUE))
 		if(length(pi) > 0) {
 			for (i in 1:length(pi)) {
 				if(do.average) cqp[[i]] <- trajectories$quantiles[[i]]
@@ -381,11 +407,16 @@ e0.trajectories.plot <- function(e0.pred, country, pi=c(80, 95), both.sexes=FALS
 					col=this.col[5], lwd=lwd[5])
 			}
 		}
-		# plot median	
-		lines(plot.data[[ipred]]$pred.x, e0.median, type='l', col=this.col[3], lwd=lwd[3])
-		legend <- if(adjusted.only) 'median' else 'adj. median'
-		if(do.both.sexes) legend <- paste(lowerize(get.sex.label(meta)), legend)
-		lty <- 1
+		legend <- lty <- lwds <- cols <- c()
+		# plot main projection
+		if(!is.null(e0.main.proj)){
+		    lines(plot.data[[ipred]]$pred.x, e0.main.proj, type='l', col=this.col[3], lwd=lwd[3])
+		    legend <- if(adjusted.only) main.proj.name else paste('adj.', main.proj.name)
+		    if(do.both.sexes) legend <- paste(lowerize(get.sex.label(meta)), legend)
+		    lty <- 1
+		    lwds <- lwd[3]
+		    cols <- this.col[3]
+		}
 		# plot given CIs
 		if(length(pi) > 0) {
 			tlty <- 2:(length(pi)+1)
@@ -401,18 +432,33 @@ e0.trajectories.plot <- function(e0.pred, country, pi=c(80, 95), both.sexes=FALS
 		legend <- c(legend, paste('observed', if(do.both.sexes) paste(lowerize(get.sex.label(meta)), 'e0') else 'e0'))
 		lty <- c(lty, 1)
 		pchs <- c(rep(-1, length(legend)-1), pch[1])
-		lwds <- c(lwd[3], rep(lwd[4], length(pi)), lwd[1])
-		cols <- c(this.col[3], rep(this.col[4], length(pi)), this.col[1])
-		if(!adjusted.only) { # plot unadjusted median
-			bhm.median <- bayesTFR::get.median.from.prediction(e0pred, country$index, country$code, adjusted=FALSE)
-			lines(plot.data[[ipred]]$pred.x, bhm.median, type='l', col=this.col[3], lwd=lwd[3], lty=max(lty)+1)
-			bhm.leg <- 'BHM median'
-			if(do.both.sexes) bhm.leg <- paste(lowerize(get.sex.label(meta)), bhm.leg)
-			legend <- c(legend, bhm.leg)
+		lwds <- c(lwds, rep(lwd[4], length(pi)), lwd[1])
+		cols <- c(cols, rep(this.col[4], length(pi)), this.col[1])
+		if(!adjusted.only) { # plot unadjusted median / mean
+		    if(main.proj.name == "mean"){
+		        bhm.main <- bayesTFR::get.mean.from.prediction(e0pred, country$index, country$code, adjusted=FALSE)
+		        bhm.main.name <- 'BHM mean'
+		    } else {
+			    bhm.main <- bayesTFR::get.median.from.prediction(e0pred, country$index, country$code, adjusted=FALSE)
+			    bhm.main.name <- 'BHM median'
+		    }
+			lines(plot.data[[ipred]]$pred.x, bhm.main, type='l', col=this.col[3], lwd=lwd[3], lty=max(lty)+1)
+			if(do.both.sexes) bhm.main.name <- paste(lowerize(get.sex.label(meta)), bhm.main.name)
+			legend <- c(legend, bhm.main.name)
 			cols <- c(cols, this.col[3])
 			lwds <- c(lwds, lwd[3])
 			lty <- c(lty, max(lty)+1)
 		}
+		if(show.median && show.mean){
+		    # plot mean in addition to median
+		    lines(plot.data[[ipred]]$pred.x, e0.mean, type='l', col=this.col[3], lwd=1, lty=max(lty)+1)
+		    mean.leg <- 'mean'
+		    if(do.both.sexes) mean.leg <- paste(lowerize(get.sex.label(meta)), mean.leg)
+		    legend <- c(legend, mean.leg)
+		    cols <- c(cols, this.col[3])
+		    lwds <- c(lwds, 1)
+		    lty <- c(lty, max(lty)+1)
+		}
 		if(lpart2 > 0) {
 			legend <- c(legend, paste('imputed', if(do.both.sexes) paste(lowerize(get.sex.label(meta)), 'e0') else 'e0'))
 			cols <- c(cols, this.col[2])

R/project_subnat.R

@@ -1,6 +1,6 @@
 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, annual = NULL,
+                               end.year=2100, start.year=NULL, subdir = "predictions", 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 
@@ -30,7 +30,7 @@ e0.predict.subnat <- function(countries, my.e0.file, sim.dir=file.path(getwd(),
         return(compute.alpha.ar1(...))
 
     method <- match.arg(method)
-    wpred <- get.e0.prediction(sim.dir) # contains national projections
+    wpred <- get.e0.prediction(sim.dir, subdir = subdir) # contains national projections
     wdata <- wpred$e0.matrix.reconstructed
     wmeta <- wpred$mcmc.set$meta
     if(!is.null(seed)) set.seed(seed)
@@ -153,7 +153,7 @@ e0.predict.subnat <- function(countries, my.e0.file, sim.dir=file.path(getwd(),
         PIs_cqp <- array(NA, c(nr.reg, length(quantiles.to.keep), nrow(wtrajs)),
                      dimnames=list(meta$regions$country_code, dimnames(wpred$quantiles)[[2]], dimnames(wtrajs)[[1]]))
         mean_sd <- array(NA, c(nr.reg, 2, nrow(wtrajs)))
-        #meta$Tc.index <- .get.Tcindex(meta$e0.matrix, cnames = meta$regions$country_name)
+        meta$Tc.index <- .get.Tcindex(meta$e0.matrix, cnames = meta$regions$country_name, stop.if.less.than2 = FALSE) # allow just one data point
         country.char <- as.character(country.obj$code)
         e0reconstructed <- meta$e0.matrix
 
@@ -171,7 +171,7 @@ e0.predict.subnat <- function(countries, my.e0.file, sim.dir=file.path(getwd(),
                 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
+                meta$Tc.index[[region]] <- meta$Tc.index[[region]][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)