Skip to content

Commit

Permalink
Merge branch 'master' into cran
Browse files Browse the repository at this point in the history
  • Loading branch information
@hanase
hanase committed Nov 8, 2024
2 parents 6df97a8 + 9c7116a commit a118438
Show file tree
Hide file tree
Showing 27 changed files with 622 additions and 160 deletions.
24 changes: 24 additions & 0 deletions ChangeLog
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,27 @@
7.4-4 (2024-11-01)
-----
Two countries added to small countries in include_2024.

7.4-3 (2024-10-31)
-----
Added option of keeping multiple predictions directories in one simulation directory
(argument "subdir" added to various functions).

Added option of shifting projection to wpp by matching means.

Added dataset include_2024.

Fixed bug related to a combination of ar.phase2 and missing data (thanks to Shunqi Cheng).

Added the mean variant into summary outputs.

tfr.trajectories.plot() got arguments to plot the means (show.mean), medians (show.median),
as well as selected trajectories (traj.index).

7.4-1/2 (2023-10-17)
-----
Removed rgdal from suggested packages.

7.4-0 (2023-09-14)
-----
Annual subnational projections are now possible, controlled by the argument "annual"
Expand Down
31 changes: 26 additions & 5 deletions DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -1,12 +1,33 @@
Package: bayesTFR
Version: 7.4-0
Date: 2023-09-14
Version: 7.4-4
Date: 2024-11-01
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]>
Authors@R: c(person(given = "Hana",
family = "Sevcikova",
role = c("cre", "aut"),
email = "[email protected]"),
person(given = "Leontine",
family = "Alkema",
role = "aut"),
person(given = "Peiran",
family = "Liu",
role = "aut"),
person(given = "Adrian",
family = "Raftery",
role = "aut",
email = "[email protected]"),
person(given = "Bailey",
family = "Fosdick",
role = "aut"),
person(given = "Patrick",
family = "Gerland",
role = "aut",
email = "[email protected]"
)
)
Depends: R (>= 3.5.0)
Imports: mvtnorm, MASS, coda, graphics, grDevices, stats, utils, wpp2019, data.table, lifecycle
Suggests: rworldmap, snowFT, googleVis, rgdal, sp, wpp2017, wpp2015, wpp2012, wpp2010, ggplot2, sf, spData, scales
Suggests: rworldmap, snowFT, googleVis, sp, wpp2017, wpp2015, wpp2012, wpp2010, ggplot2, sf, spData, scales
Description: Making probabilistic projections of total fertility rate for all countries of the world, using a Bayesian hierarchical model <doi:10.1007/s13524-011-0040-5> <doi:10.18637/jss.v106.i08>. Subnational probabilistic projections are also supported <doi:10.4054/DemRes.2018.38.60>.
License: GPL-3 | file LICENSE
URL: https://bayespop.csss.washington.edu
Expand Down
2 changes: 2 additions & 0 deletions NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,7 @@ export(
get.thinned.tfr.mcmc,
get.tfr.prediction,
has.tfr.prediction,
available.tfr.predictions,
get.tfr.convergence.all,
get.tfr.convergence,
get.tfr3.convergence.all,
Expand Down Expand Up @@ -102,6 +103,7 @@ export(
tfr.country.dlcurves,
tfr.dl.coverage,
get.median.from.prediction,
get.mean.from.prediction,
bdem.parameter.traces,
get.traj.ascii.header,
get.nr.countries,
Expand Down
27 changes: 21 additions & 6 deletions R/get_outputs.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -222,14 +222,23 @@ create.thinned.tfr.mcmc <- function(mcmc.set, thin=1, burnin=0, output.dir=NULL,
return(thin.mcmc.set)
}

has.tfr.prediction <- function(mcmc=NULL, sim.dir=NULL) {
has.tfr.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.tfr.prediction <- function(mcmc=NULL, sim.dir=NULL, mcmc.dir=NULL) {
available.tfr.predictions <- function(mcmc=NULL, sim.dir=NULL, full.names = FALSE){
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.')
all.dirs <- list.dirs(sim.dir, recursive = FALSE, full.names = TRUE)
pred.dirs <- names(unlist(sapply(all.dirs, function(x) list.files(x, pattern = "^prediction.rda$"))))
if(length(pred.dirs) == 0 || full.names == TRUE) return(pred.dirs)
return(sapply(pred.dirs, basename, USE.NAMES = FALSE))
}

get.tfr.prediction <- function(mcmc=NULL, sim.dir=NULL, mcmc.dir=NULL, subdir = "predictions") {
############
# Returns an object of class bayesTFR.prediction
# Set mcmc.dir to NA, if the prediction object should not have a pointer
Expand All @@ -238,11 +247,15 @@ get.tfr.prediction <- function(mcmc=NULL, sim.dir=NULL, mcmc.dir=NULL) {
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.tfr.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)
bayesTFR.prediction$output.directory <- normalizePath(output.dir)
Expand Down Expand Up @@ -365,7 +378,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
Expand Down
6 changes: 3 additions & 3 deletions R/mcmc_ini.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -95,9 +95,9 @@ find.lambda.for.one.country <- function(tfr, T_end, annual = FALSE) {
if(annual) # convert lambda from 5-year scale to annual scale
{
lambda <- min(which(year.bin == lambda) + 2, Tendorig)
if (length(year.bin) - lambda < 5) { # if in the last time period, set it to the end of the period
lambda <- length(year.bin)
while(is.na(tfrorig[lambda])) lambda <- lambda - 1 # move it before the last NA if any
if (Tendorig - lambda < 5) { # if in the last (observed) time period, set it to the end of the period
lambda <- Tendorig
#while(is.na(tfrorig[lambda])) lambda <- lambda - 1 # move it before the last NA if any - not needed anymore
}
}

Expand Down
130 changes: 97 additions & 33 deletions R/plot_functions.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down Expand Up @@ -551,7 +564,8 @@ tfr.estimation.plot <- function(mcmc.list = NULL, country = NULL, sim.dir = NULL
tfr.trajectories.plot <- function(tfr.pred, country, pi=c(80, 95),
half.child.variant=TRUE, nr.traj=NULL,
adjusted.only = TRUE, typical.trajectory=FALSE,
mark.estimation.points=FALSE,
mark.estimation.points=FALSE,
traj.index = NULL, show.mean = FALSE, show.median = TRUE,
xlim=NULL, ylim=NULL, type='b',
xlab='Year', ylab='TFR', main=NULL, lwd=c(2,2,2,2,2,1),
col=c('black', 'green', 'red', 'red', 'blue', '#00000020'),
Expand Down Expand Up @@ -595,7 +609,11 @@ tfr.trajectories.plot <- function(tfr.pred, country, pi=c(80, 95),
}
x1 <- as.integer(c(names(y1.part1), names(y1.part2)))
x2 <- as.numeric(dimnames(tfr.pred$quantiles)[[3]])

if(!is.null(traj.index)) nr.traj <- length(traj.index)
trajectories <- get.trajectories(tfr.pred, country$code, nr.traj, typical.trajectory=typical.trajectory)
if(!is.null(traj.index) && !is.null(trajectories$trajectories)) trajectories$index <- traj.index

# plot historical data: observed
if (!add) {
if(is.null(xlim)) xlim <- c(min(x1,x2), max(x1,x2))
Expand Down Expand Up @@ -648,69 +666,116 @@ tfr.trajectories.plot <- function(tfr.pred, country, pi=c(80, 95),
lines(x2, trajectories$trajectories[,trajectories$index[i]], type='l', col=col[6], lwd=lwd[6])
}
}
# plot median
tfr.median <- get.median.from.prediction(tfr.pred, country$index, country$code)
if(uncertainty) {
if(uncertainty)
unc.last.time <- which(tfr.object$tfr_quantile$year == x2[1])
tfr.median[1] <- unlist(tfr.object$tfr_quantile[unc.last.time,])[length(pi)+1]

# extract median & mean
tfr.median <- tfr.mean <- tfr.main.proj <- NULL
if(show.median)
tfr.median <- get.median.from.prediction(tfr.pred, country$index, country$code)
if(show.mean)
tfr.mean <- get.mean.from.prediction(tfr.pred, country$index, country$code)

# set the main projection (solid line)
main.proj.name <- ""
if(!is.null(tfr.median)){
tfr.main.proj <- tfr.median
main.proj.name <- "median"
} else {
if(!is.null(tfr.mean)){
tfr.main.proj <- tfr.mean
main.proj.name <- "mean"
}
}
lty <- c()
if(!is.null(tfr.main.proj)){
if(uncertainty) # replace last observed with estimated median
tfr.main.proj[1] <- unlist(tfr.object$tfr_quantile[unc.last.time,])[length(pi)+1]
# draw main projection
lines(x2, tfr.main.proj, type='l', col=col[3], lwd=lwd[3])
lty <- 1
}
lines(x2, tfr.median, type='l', col=col[3], lwd=lwd[3])
lty <- 1

# plot given CIs
if(length(pi) > 0){
lty <- c(lty, 2:(length(pi)+1))
pi.lty <- 2:(length(pi)+1)
lty <- c(lty, pi.lty)
for (i in 1:length(pi)) {
cqp <- get.traj.quantiles(tfr.pred, country$index, country$code, trajectories$trajectories, pi[i],
est.uncertainty = uncertainty)
if (!is.null(cqp)) {
lines(x2, cqp[1,], type='l', col=col[4], lty=lty[i+1], lwd=lwd[4])
lines(x2, cqp[2,], type='l', col=col[4], lty=lty[i+1], lwd=lwd[4])
lines(x2, cqp[1,], type='l', col=col[4], lty=pi.lty[i], lwd=lwd[4])
lines(x2, cqp[2,], type='l', col=col[4], lty=pi.lty[i], lwd=lwd[4])
}
}
}
max.lty <- if(length(lty) == 0) 1 else max(lty)

if (uncertainty)
{
col_median <- length(pi)+1
lines(tfr.object$tfr_quantile$year, as.data.frame(tfr.object$tfr_quantile)[, col_median], type='l', col=col_unc, lwd=lwd[3])
if(!adjusted.only) { # plot unadjusted estimation median
unadj.lty <- max(lty)+1
unadj.lty <- max.lty+1
tfr.object.unadj <- get.tfr.estimation(mcmc.list=tfr.pred$mcmc.set, country = country.obj$code,
probs=0.5, adjust = FALSE)
lines(tfr.object.unadj$tfr_quantile$year, as.data.frame(tfr.object.unadj$tfr_quantile)$V1, type='l', col=col_unc, lwd=lwd[3], lty = max(lty)+1)
lines(tfr.object.unadj$tfr_quantile$year, as.data.frame(tfr.object.unadj$tfr_quantile)$V1, type='l', col=col_unc, lwd=lwd[3], lty = unadj.lty)
}
if(length(pi) > 0) {
for (i in 1:length(pi)) {
lines(tfr.object$tfr_quantile$year, as.data.frame(tfr.object$tfr_quantile)[, length(pi)+1-i], type='l', col=col_unc, lty=lty[i+1], lwd=lwd[4])
lines(tfr.object$tfr_quantile$year, as.data.frame(tfr.object$tfr_quantile)[, length(pi)+1+i], type='l', col=col_unc, lty=lty[i+1], lwd=lwd[4])
lines(tfr.object$tfr_quantile$year, as.data.frame(tfr.object$tfr_quantile)[, length(pi)+1-i], type='l', col=col_unc, lty=pi.lty[i], lwd=lwd[4])
lines(tfr.object$tfr_quantile$year, as.data.frame(tfr.object$tfr_quantile)[, length(pi)+1+i], type='l', col=col_unc, lty=pi.lty[i], lwd=lwd[4])
}
}
}
legend <- c()
cols <- c()
lwds <- c()
if(!adjusted.only) { # plot unadjusted median
bhm.median <- get.median.from.prediction(tfr.pred, country$index, country$code, adjusted=FALSE)
lines(x2, bhm.median, type='l', col=col[3], lwd=lwd[3], lty=max(lty)+1)
legend <- c(legend, 'BHM median')
if(!adjusted.only) { # plot unadjusted median & mean
if(main.proj.name == "mean"){
bhm.main <- get.mean.from.prediction(tfr.pred, country$index, country$code, adjusted=FALSE)
bhm.main.name <- 'BHM mean'
} else {
bhm.main <- get.median.from.prediction(tfr.pred, country$index, country$code, adjusted=FALSE)
bhm.main.name <- 'BHM median'
}
lines(x2, bhm.main, type='l', col=col[3], lwd=lwd[3], lty=max.lty+1)
legend <- c(legend, bhm.main.name)
cols <- c(cols, col[3])
lwds <- c(lwds, lwd[3])
lty <- c(max(lty)+1, lty)
lty <- c(max.lty+1, lty)
max.lty <- max(lty)
}
median.legend <- if(adjusted.only) 'median' else 'adj. median'
legend <- c(legend, median.legend, if(length(pi) > 0) paste0(pi, '% PI') else c())
cols <- c(cols, col[3], rep(col[4], length(pi)))
lwds <- c(lwds, lwd[3], rep(lwd[4], length(pi)))
if (half.child.variant) {
lty <- c(lty, max(lty)+1)
if(main.proj.name != ""){
main.legend <- if(adjusted.only) main.proj.name else paste('adj.', main.proj.name)
legend <- c(legend, main.legend)
cols <- c(cols, col[3])
lwds <- c(lwds, lwd[3])
}
legend <- c(legend, if(length(pi) > 0) paste0(pi, '% PI') else c())
cols <- c(cols, rep(col[4], length(pi)))
lwds <- c(lwds, rep(lwd[4], length(pi)))
if(show.median && show.mean){
# plot mean in addition to median
lines(x2, tfr.mean, type='l', col=col[3], lwd=1, lty=max(lty)+1)
legend <- c(legend, 'mean')
cols <- c(cols, col[3])
lwds <- c(lwds, 1)
lty <- c(lty, max.lty+1)
max.lty <- max(lty)
}
if (half.child.variant && !is.null(tfr.main.proj)) {
lty <- c(lty, max.lty+1)
max.lty <- max(lty)
llty <- length(lty)
up.low <- get.half.child.variant(median=tfr.median)
up.low <- get.half.child.variant(median=tfr.main.proj)
if(uncertainty) {
up.low <- up.low[,-1]
x2t <- x2[-1]
} else x2t <- x2
lines(x2t, up.low[1,], type='l', col=col[5], lty=lty[llty], lwd=lwd[5])
lines(x2t, up.low[2,], type='l', col=col[5], lty=lty[llty], lwd=lwd[5])
legend <- c(legend, '+/- 0.5 child')
legend <- c(legend, paste(main.proj.name, '+/- 0.5 child'))
cols <- c(cols, col[5])
lwds <- c(lwds, lwd[5])
}
Expand Down Expand Up @@ -1210,11 +1275,10 @@ tfr.map <- function(pred, quantile=0.5, year=NULL, par.name=NULL, adjusted=FALSE
tfr <- data.frame(cbind(un=data.period$country.codes, tfr=data))
map <- rworldmap::getMap(resolution=resolution)
#first get countries excluding Antarctica which crashes spTransform (says the help page for joinCountryData2Map)
sPDF <- map[-which(map$ADMIN=='Antarctica'), ]
if(requireNamespace("rgdal", quietly=TRUE)) {
#transform map to the Robinson projection
sPDF <- sp::spTransform(sPDF, CRSobj=sp::CRS("+proj=robin +ellps=WGS84"))
}
sPDF <- map[-which(map$ADMIN=='Antarctica'), ]
#transform map to the Robinson projection
sPDF <- sp::spTransform(sPDF, CRSobj = sp::CRS("+proj=robin +ellps=WGS84"))

## recode missing UN codes and UN member states
sPDF$UN <- sPDF$ISO_N3
## N. Cyprus -> assign to Cyprus
Expand All @@ -1223,7 +1287,7 @@ tfr.map <- function(pred, quantile=0.5, year=NULL, par.name=NULL, adjusted=FALSE
sPDF$UN[sPDF$ISO3=="KOS"] <- 688
## W. Sahara -> no UN numerical code assigned in Natural Earth map (its ISO3 changed in rworlmap 1.3.6)
sPDF$UN[sPDF$ISO3=="ESH"] <- 732
## Somaliland -> assign to Somalia -> fixed in rworlmap version 1.3.6
## Somaliland -> assign to Somalia (SOM) -> fixed in rworlmap version 1.3.6
#sPDF$UN[sPDF$ISO3=="SOL"] <- 706

#mtfr <- joinCountryData2Map(tfr, joinCode='UN', nameJoinColumn='un')
Expand Down
Loading

0 comments on commit a118438

Please sign in to comment.