From 67a883009ef30e753e9586a53254c10fac10f1d4 Mon Sep 17 00:00:00 2001 From: Hana Sevcikova Date: Thu, 20 Jun 2024 17:49:49 -0700 Subject: [PATCH] argument show.median added --- ChangeLog | 4 +- R/plot_functions.R | 102 +++++++++++++++++++++++++++------------ man/plot-trajectories.Rd | 4 +- 3 files changed, 74 insertions(+), 36 deletions(-) diff --git a/ChangeLog b/ChangeLog index ed16bae..7ccef41 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,10 +1,10 @@ -7.4-2.90xx (2024-06-05) +7.4-2.90xx (2024-06-20) ----- 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), +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) diff --git a/R/plot_functions.R b/R/plot_functions.R index 43fe2ee..a155f7a 100644 --- a/R/plot_functions.R +++ b/R/plot_functions.R @@ -565,7 +565,7 @@ 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, - traj.index = NULL, show.mean = 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'), @@ -666,79 +666,117 @@ 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" + } } - lines(x2, tfr.median, type='l', col=col[3], lwd=lwd[3]) - lty <- 1 + 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 + } + # 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) + max.lty <- max(lty) + } + main.legend <- "" + 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]) - lty <- c(max(lty)+1, 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(show.mean){ - # plot mean - tfr.mean <- get.mean.from.prediction(tfr.pred, country$index, country$code) + 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(max(lty)+1, lty) + lty <- c(lty, max.lty+1) + max.lty <- max(lty) } - if (half.child.variant) { - lty <- c(lty, max(lty)+1) + 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]) } diff --git a/man/plot-trajectories.Rd b/man/plot-trajectories.Rd index e3732ad..ff14bf1 100644 --- a/man/plot-trajectories.Rd +++ b/man/plot-trajectories.Rd @@ -14,7 +14,7 @@ tfr.trajectories.plot(tfr.pred, country, pi = c(80, 95), half.child.variant = TRUE, nr.traj = NULL, adjusted.only = TRUE, typical.trajectory = FALSE, mark.estimation.points = FALSE, - traj.index = NULL, show.mean = 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'), @@ -39,7 +39,7 @@ tfr.trajectories.table(tfr.pred, country, pi = c(80, 95), \item{typical.trajectory}{Logical. If \code{TRUE} one trajectory is shown that has the smallest distance to the median.} \item{mark.estimation.points}{Logical. If \code{TRUE}, points that were not used in the estimation (phase I) are shown in lighter color than points in phase II and III.} \item{traj.index}{Vector of trajectory indices to show. If not given, the trajectories are selected using equidistant spacing.} - \item{show.mean}{Logical indicating if the mean of the distribution should be shown.} + \item{show.mean, show.median}{Logical indicating if the mean or/and the median of the distribution should be shown.} \item{xlim, ylim, type, xlab, ylab}{Graphical parameters passed to the \code{plot} function.} \item{lwd, col}{Vector of six elements giving the line width and color for: 1. observed data, 2. imputed missing data, 3. median, 4. quantiles, 5. half-child variant, 6. trajectories.} \item{show.legend}{Logical controlling whether the legend should be drawn.}