Skip to content

Commit

Permalink
argument show.median added
Browse files Browse the repository at this point in the history
  • Loading branch information
hanase committed Jun 21, 2024
1 parent 2cf3d4b commit 67a8830
Show file tree
Hide file tree
Showing 3 changed files with 74 additions and 36 deletions.
4 changes: 2 additions & 2 deletions ChangeLog
Original file line number Diff line number Diff line change
@@ -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)
Expand Down
102 changes: 70 additions & 32 deletions R/plot_functions.R
Original file line number Diff line number Diff line change
Expand Up @@ -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'),
Expand Down Expand Up @@ -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])
}
Expand Down
4 changes: 2 additions & 2 deletions man/plot-trajectories.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -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'),
Expand All @@ -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.}
Expand Down

0 comments on commit 67a8830

Please sign in to comment.