Nothing
DLcurve.plot.all <- function (mcmc.list = NULL, sim.dir = NULL,
output.dir = file.path(getwd(), "DLcurves"),
output.type="png",
burnin = NULL, verbose = FALSE, ...) {
if(!file.exists(output.dir)) dir.create(output.dir, recursive=TRUE)
if(is.null(mcmc.list)) mcmc.list <- get.tfr.mcmc(sim.dir=sim.dir, verbose=verbose, burnin=burnin)
mcmc.list <- get.mcmc.list(mcmc.list)
meta <- mcmc.list[[1]]$meta
.do.plot.all.country.loop(as.character(country.names(meta)), meta, output.dir, DLcurve.plot, output.type=output.type,
file.prefix='DLplot', plot.type='DL graph', verbose=verbose, mcmc.list = mcmc.list, burnin = burnin, ...)
}
stop.if.country.not.DL <- function(country.obj, meta) {
if (!is.element(country.obj$index, meta$id_DL))
stop('Country ', country.obj$name, ' not estimated because no decline observed.')
}
tfr.world.dlcurves <- function(x, mcmc.list, burnin=NULL, countryUc=NULL, ...) {
# Get the hierarchical DL curves with U_c for a given country (countryUc)
# If countryUc is not given, take the middle point of the U_c prior.
if(inherits(mcmc.list, 'bayesTFR.prediction')) {
if(!is.null(burnin) && burnin != mcmc.list$burnin)
warning('Prediction was generated with different burnin. Burnin set to ', mcmc.list$burnin)
burnin <- 0 # because burnin was already cut of the traces
}
if(is.null(burnin)) burnin <- 0
mcmc.list <- get.mcmc.list(mcmc.list)
country <- if(!is.null(countryUc)) get.country.object(countryUc, mcmc.list[[1]]$meta)$code else countryUc
return(tfr.get.dlcurves(x, mcmc.list, country.code=NULL, country.index=NULL, burnin=burnin, countryUc=country, ...))
}
tfr.country.dlcurves <- function(x, mcmc.list, country, burnin=NULL, ...) {
# Get country-specific DL curves.
# It's a wrapper around tfr.get.dlcurves for easier usage.
if(inherits(mcmc.list, 'bayesTFR.prediction')) {
if(!is.null(burnin) && burnin != mcmc.list$burnin)
warning('Prediction was generated with different burnin. Burnin set to ', mcmc.list$burnin)
burnin <- 0 # because burnin was already cut of the traces
}
if(is.null(burnin)) burnin <- 0
mcmc.list <- get.mcmc.list(mcmc.list)
country.obj <- get.country.object(country, mcmc.list[[1]]$meta)
if(is.null(country.obj$code))
stop("Country ", country, " not found.")
return(tfr.get.dlcurves(x, mcmc.list, country.code=country.obj$code, country.index=country.obj$index, burnin=burnin, ...))
}
tfr.get.dlcurves <- function(x, mcmc.list, country.code, country.index, burnin=0, nr.curves=2000, predictive.distr=FALSE,
return.sigma=FALSE, countryUc=NULL) {
# If country.code is null, get world distribution. In such a case, countryUc gives the country code
# that should be used for retrieving U_c. If not given, the middle point between U_c prior is taken.
.get.sig <- function(i, sigma, S, a, b) return(sigma + (x[i] - S)*ifelse(x[i] > S, -a, b))
.get.sig.distr <- function(i, traces)
return(apply(traces, 1, function(y)
return(pmax(.get.sig(i, y['sigma0'], y['S_sd'], y['a_sd'], y['b_sd']), mcmc$meta$sigma0.min))))
dlc <- sigma.all <- c()
cspec <- TRUE
Uvalue <- NULL
if(length(mcmc.list) == 0) stop("No MCMCs available.")
if(!is.null(country.code) && !is.element(country.index, mcmc.list[[1]]$meta$id_Tistau)) {
U.var <- paste0("U_c", country.code)
d.var <- paste0("d_c", country.code)
Triangle_c4.var <- paste0("Triangle_c4_c", country.code)
gamma.vars <- paste0("gamma_", 1:3, "_c", country.code)
} else {
U.var <- "U"
d.var <- "d"
Triangle_c4.var <- "Triangle_c4"
gamma.vars <- paste0("gamma_", 1:3)
alpha.vars <- paste0('alpha_',1:3)
delta.vars <- paste0('delta_',1:3)
if(!is.null(country.code))
Uvalue = get.observed.tfr(country.index, mcmc.list[[1]]$meta,
'tfr_matrix_all')[mcmc.list[[1]]$meta$tau_c[country.index]]
else {
if(!is.null(countryUc)) U.var <- paste0("U_c", countryUc)
else Uvalue <- mcmc.list[[1]]$meta$U.c.low.base + (mcmc.list[[1]]$meta$U.up - mcmc.list[[1]]$meta$U.c.low.base)/2
}
cspec <- FALSE
}
# Compute the quantiles on a sample of at least 2000.
nr.curves.from.mc <- if (!is.null(nr.curves)) ceiling(max(nr.curves, 2000)/length(mcmc.list))
else NULL
for (mcmc in mcmc.list) {
th.burnin <- get.thinned.burnin(mcmc,burnin)
thincurves.mc <- get.thinning.index(nr.points=nr.curves.from.mc,
all.points=mcmc$length - th.burnin)
if(cspec) # country specific
traces <- data.frame(load.tfr.parameter.traces.cs(mcmc, country.code,
burnin=th.burnin,
thinning.index=thincurves.mc$index))
else { #Tistau country. Use hierarchical distr.
traces <- data.frame(load.tfr.parameter.traces(mcmc,
burnin=th.burnin,
thinning.index=thincurves.mc$index))
ltraces <- nrow(traces)
for (i in 1:3){
gamma <- rnorm(ltraces, mean = traces[,alpha.vars[i]],
sd = traces[,delta.vars[i]])
traces <- cbind(traces, gamma)
colnames(traces)[ncol(traces)] <- gamma.vars[i]
}
Triangle4_tr_s <- rnorm(ltraces, mean = traces[,'Triangle4'], sd = traces[, 'delta4'])
traces <- cbind(traces,
Triangle_c4=( mcmc$meta$Triangle_c4.up*exp(Triangle4_tr_s) + mcmc$meta$Triangle_c4.low)/(1+exp(Triangle4_tr_s)))
d_tr_s <- rnorm(ltraces, mean = traces[,'chi'], sd = traces[,'psi'])
if(is.null(Uvalue)) {
traces <- cbind(traces, data.frame(load.tfr.parameter.traces.cs(mcmc, countryUc, par.names = "U",
burnin=th.burnin,
thinning.index=thincurves.mc$index)))
} else traces <- cbind(traces, U=rep(Uvalue, ltraces))
traces <- cbind(traces, d=(mcmc$meta$d.up*(exp(d_tr_s) + mcmc$meta$d.low)/(1+exp(d_tr_s))))
}
theta <- (traces[, U.var] - traces[, Triangle_c4.var] ) *
exp(traces[, gamma.vars, drop=FALSE])/apply(exp(traces[,gamma.vars, drop=FALSE]), 1, sum)
theta <- cbind(theta, traces[, Triangle_c4.var], traces[, d.var])
dl <- t(apply(theta, 1, DLcurve, tfr = x, p1 = mcmc$meta$dl.p1, p2 = mcmc$meta$dl.p2,
annual = get.item(mcmc$meta, "annual.simulation", FALSE)))
if(length(x) == 1) dl <- t(dl)
if(predictive.distr || return.sigma) {
wp.traces <- load.tfr.parameter.traces(mcmc,
burnin=th.burnin,
thinning.index=thincurves.mc$index,
par.names=c('sigma0', 'S_sd', 'a_sd', 'b_sd'))
sigma_eps <- NULL
for(j in 1:length(x))
sigma_eps <- cbind(sigma_eps, .get.sig.distr(j, wp.traces))
if(predictive.distr) {
errors <- t(apply(sigma_eps, 1, function(sig) rnorm(dim(dl)[2],0,sig)))
if(length(x) == 1 && all(dim(errors) == rev(dim(dl)))) errors <- t(errors)
dlc <- rbind(dlc, dl+errors)
} else {
dlc <- rbind(dlc, dl)
sigma.all <- rbind(sigma.all, sigma_eps)
}
} else dlc <- rbind(dlc, dl)
}
return (if(!return.sigma) dlc else list(dl=dlc, sigma=sigma.all))
}
.match.colors.with.default <- function(col, default.col) {
ldcol <- length(default.col)
lcol <- length(col)
if(lcol >= ldcol) return(col)
col <- rep(col, ldcol)
col[(lcol+1):ldcol] <- default.col[(lcol+1):ldcol]
return(col)
}
DLcurve.plot <- function (mcmc.list, country, burnin = NULL, pi = 80, tfr.max = 10,
nr.curves = NULL, predictive.distr=FALSE, ylim = NULL, xlab = "TFR (reversed)", ylab = "TFR decrement",
main = NULL, show.legend=TRUE, col=c('black', 'red', "#00000020"), ...
)
{
if(inherits(mcmc.list, 'bayesTFR.prediction')) {
if(!is.null(burnin) && burnin != mcmc.list$burnin)
warning('Prediction was generated with different burnin. Burnin set to ', mcmc.list$burnin)
burnin <- 0 # because burnin was already cut of the traces
}
col <- .match.colors.with.default(col, c('black', 'red', "#00000020"))
if(is.null(burnin)) burnin <- 0
mcmc.list <- get.mcmc.list(mcmc.list)
meta <- mcmc.list[[1]]$meta
country.obj <- get.country.object(country, meta)
if(is.null(country.obj$code)) stop("Country ", country, " not found.")
country <- country.obj
if (!is.null(mcmc.list[[1]]$uncertainty) && mcmc.list[[1]]$uncertainty)
{
mcmc.list.tmp <- list(meta = meta, mcmc.list = mcmc.list)
obs.data <- unlist(array(get.tfr.estimation(mcmc.list.tmp, country.obj$code, probs = 0.5)$tfr_quantile[,1]))
}
else
obs.data <- get.observed.tfr(country$index, meta, 'tfr_matrix_observed', 'tfr_matrix_all')[1:meta$T_end_c[country$index]]
#stop.if.country.not.DL(country, meta)
tfr_plot <- seq(0, tfr.max, 0.1)
dlc <- tfr.get.dlcurves(tfr_plot, mcmc.list, country$code, country$index, burnin, nr.curves,
predictive.distr=predictive.distr)
miny <- min(dlc)
maxy <- max(dlc)
decr <- -diff(obs.data)
dl.obs.idx <- if(max(meta$tau_c[country$index],1) >= meta$lambda_c[country$index]) c()
else seq(max(meta$tau_c[country$index],1), meta$lambda_c[country$index]-1)
p3.obs.idx <- if(meta$lambda_c[country$index]>length(decr)) c() else seq(meta$lambda_c[country$index], length(decr))
maxy <- max(maxy, decr, na.rm=TRUE)
miny <- min(miny, decr, na.rm=TRUE)
thincurves <- get.thinning.index(nr.curves, dim(dlc)[1])
ltype <- "l"
if (thincurves$nr.points == 0) {
ltype <- "n"
thincurves$index <- 1
}
if (is.null(main)) main <- country$name
if (is.null(ylim)) ylim <- c(miny, maxy)
# plot trajectories
plot(dlc[thincurves$index[1], ] ~ tfr_plot, col = col[3],
type = ltype,
xlim = c(max(tfr_plot), min(tfr_plot)),
ylim = ylim, ylab = ylab, xlab = xlab, main = main, ...
)
if (thincurves$nr.points > 1) {
for (i in 2:thincurves$nr.points) {
lines(dlc[thincurves$index[i], ] ~ tfr_plot, col = col[3])
}
}
# plot quantiles
dl50 <- apply(dlc, 2, quantile, 0.5)
lines(dl50 ~ tfr_plot, col = col[2], lwd = 2)
lty <- 2:(length(pi) + 1)
for (i in 1:length(pi)) {
al <- (1 - pi[i]/100)/2
dlpi <- apply(dlc, 2, quantile, c(al, 1 - al))
lines(dlpi[1, ] ~ tfr_plot, col = col[2], lty = lty[i],
lwd = 2)
lines(dlpi[2, ] ~ tfr_plot, col = col[2], lty = lty[i],
lwd = 2)
}
# plot observed data
obs.legend <- list(legend=c(), pch=c(), bg=c(), col=c())
if(length(dl.obs.idx) > 0 && any(!is.na(decr[dl.obs.idx]))) {
points(decr[dl.obs.idx] ~ obs.data[-length(obs.data)][dl.obs.idx], pch = 19, col=col[1])
obs.legend$legend <- c(obs.legend$legend, 'Phase II data')
obs.legend$pch <- c(obs.legend$pch, 19)
obs.legend$col <- c(obs.legend$col, col[1])
#obs.legend$bg <- c(obs.legend$bg, col[1])
}
endpI <- if(length(dl.obs.idx)==0) max(meta$tau_c[country$index],1) else dl.obs.idx[1]-1
if((endpI > 1) && any(!is.na(decr[1:endpI]))) {# draw phase I points
points(decr[1:endpI] ~ obs.data[-length(obs.data)][1:endpI], pch=0, col=col[1])
obs.legend$legend <- c(obs.legend$legend, 'Phase I data')
obs.legend$pch <- c(obs.legend$pch, 0)
obs.legend$col <- c(obs.legend$col, col[1])
#obs.legend$bg <- c(obs.legend$bg, col[1])
}
if(length(p3.obs.idx)>0) {
points(decr[p3.obs.idx] ~ obs.data[-length(obs.data)][p3.obs.idx], pch = 2, col=col[1]) # draw phase III points
obs.legend$legend <- c(obs.legend$legend, 'Phase III data')
obs.legend$pch <- c(obs.legend$pch, 2)
obs.legend$col <- c(obs.legend$col, col[1])
#obs.legend$bg <- c(obs.legend$bg, 'grey')
}
if(show.legend)
legend("topright", legend = c("median", paste("PI", pi), obs.legend$legend),
lty = c(1, lty, rep(0,length(obs.legend$pch))), bty = "n",
col = c(rep(col[2], 1+length(pi)), obs.legend$col),
pch=c(rep(-1,length(pi)+1), obs.legend$pch),
#bg=c(rep(-1,length(pi)+1), obs.legend$bg),
)
}
.get.trajectories.table <- function(tfr.pred, country, obs.data, pi, pred.median, cqp, half.child.variant=FALSE,
uncertainty=FALSE, adjusted = TRUE) {
l <- tfr.pred$nr.projections
obs.data <- obs.data[!is.na(obs.data)]
x1 <- as.integer(names(obs.data))
year.step <- ifelse(get.item(tfr.pred$mcmc.set$meta, "annual.simulation", FALSE), 1, 5)
x2 <- seq(max(x1)+year.step, by=year.step, length=l)
if (!uncertainty)
tfr <- as.matrix(obs.data, ncol=1)
else
{
tmp <- get.tfr.estimation(mcmc.list = tfr.pred$mcmc.set, country = country$code,
probs = c(0.5, sort(c((100-pi)/200, 1-(100-pi)/200))), adjust = adjusted)
tfr <- as.matrix(tmp$tfr_quantile)[,1:(1+2*length(pi))]
}
rownames(tfr) <- x1
pred.table <- matrix(NA, ncol=2*length(pi)+1, nrow=l)
pred.table[,1] <- pred.median[2:(l+1)]
colnames(pred.table) <- c('median', rep(NA,ncol(pred.table)-1))
if (uncertainty) colnames(tfr) <- c('median', rep(NA,ncol(pred.table)-1))
idx <- 2
for (i in 1:length(pi)) {
al <- (1-pi[i]/100)/2
if (!is.null(cqp[[i]])) {
pred.table[,idx:(idx+1)] <- t(cqp[[i]][,2:(l+1)])
} else{
pred.table[,idx:(idx+1)] <- matrix(NA, nrow=l, ncol=2)
}
colnames(pred.table)[idx:(idx+1)] <- c(al, 1-al)
idx <- idx+2
}
rownames(pred.table) <- x2
cn <- colnames(pred.table)[2:ncol(pred.table)]
pred.table[,2:ncol(pred.table)] <- pred.table[,cn[order(cn)]]
colnames(pred.table)[2:ncol(pred.table)] <- cn[order(cn)]
if (uncertainty) colnames(tfr)[2:ncol(tfr)] <- cn[order(cn)]
if(half.child.variant) {
up.low <- get.half.child.variant(median=c(0, pred.table[,1]))
pred.table <- cbind(pred.table, t(up.low[,2:ncol(up.low)]))
colnames(pred.table)[(ncol(pred.table)-1):ncol(pred.table)] <- c('-0.5child', '+0.5child')
if (uncertainty)
{
tfr <- cbind(tfr, matrix(NA, nrow = nrow(tfr), ncol = ncol(pred.table) - ncol(tfr)))
colnames(tfr)[(ncol(tfr)-1):ncol(tfr)] <- c('-0.5child', '+0.5child')
}
}
return(rbind(cbind(tfr, matrix(NA, nrow=nrow(tfr), ncol=ncol(pred.table)-ncol(tfr))), pred.table))
}
tfr.trajectories.table <- function(tfr.pred, country, pi=c(80, 95), half.child.variant=TRUE, adjusted = TRUE) {
if (missing(country)) {
stop('Argument "country" must be given.')
}
country.obj <- get.country.object(country, tfr.pred$mcmc.set$meta)
if(is.null(country.obj$code)) stop("Country ", country, " not found.")
country <- country.obj
uncertainty <- FALSE
if ((length(tfr.pred$mcmc.set$mcmc.list)>0 && !is.null(tfr.pred$mcmc.set$mcmc.list[[1]]$uncertainty) &&
tfr.pred$mcmc.set$mcmc.list[[1]]$uncertainty) || (country$index %in% tfr.pred$mcmc.set$meta$extra))
uncertainty <- TRUE
obs.data <- get.data.for.country.imputed(tfr.pred, country$index)
if(!is.null(tfr.pred$present.year.index)) obs.data <- obs.data[1:min(length(obs.data), tfr.pred$present.year.index.all)]
pred.median <- get.median.from.prediction(tfr.pred, country$index, country$code, adjusted = adjusted)
trajectories <- get.trajectories(tfr.pred, country$code, adjusted = adjusted)
cqp <- list()
for (i in 1:length(pi))
cqp[[i]] <- get.traj.quantiles(tfr.pred, country$index, country$code, trajectories$trajectories, pi[i],
est.uncertainty = uncertainty, adjusted = adjusted)
return(.get.trajectories.table(tfr.pred, country, obs.data, pi, pred.median, cqp, half.child.variant, uncertainty, adjusted = adjusted))
}
get.typical.trajectory.index <- function(trajectories) {
med <- apply(trajectories, 1, median, na.rm=TRUE)
sumerrors <- apply(abs(trajectories - med), 2, sum)
sorterrors <- order(sumerrors)
return(sorterrors[round(length(sorterrors)/2, 0)])
}
get.trajectories <- function(tfr.pred, country, nr.traj=NULL, adjusted=TRUE, base.name='traj', typical.trajectory=FALSE) {
traj.file <- file.path(tfr.pred$output.directory, paste(base.name, '_country', country, '.rda', sep=''))
if (!file.exists(traj.file)) return(list(trajectories=NULL))
load(traj.file)
if(typical.trajectory) {
traj.idx <- get.typical.trajectory.index(trajectories)
} else {
thintraj <- get.thinning.index(nr.traj, dim(trajectories)[2])
if (thintraj$nr.points == 0) return(list(trajectories=NULL))
traj.idx <- thintraj$index
}
if(!is.null(trajectories)) {
if(adjusted) {
shift <- get.tfr.shift(country, tfr.pred)
if(!is.null(shift)) trajectories <- trajectories + shift
}
rownames(trajectories) <- get.prediction.years(tfr.pred$mcmc.set$meta, dim(trajectories)[1],
present.year.index=tfr.pred$present.year.index)
}
return(list(trajectories=trajectories, index=traj.idx))
}
get.quantile.from.prediction <- function(tfr.pred, quantile, country.index, country.code=NULL, adjusted=TRUE,
est.uncertainty = FALSE) {
quant.values <- tfr.pred$quantiles[country.index, as.character(quantile),]
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=0.5, adjust = adjusted)
unc.last.time <- which(tfr.est$tfr_quantile$year == dimnames(tfr.pred$quantiles)[[3]][1])
quant.values[1] <- unlist(tfr.est$tfr_quantile[unc.last.time, 1])
}
if (!adjusted) return(quant.values)
shift <- get.tfr.shift(country.code, tfr.pred)
if(!is.null(shift)) quant.values <- quant.values + shift
return(quant.values)
}
get.median.from.prediction <- function(tfr.pred, country.index, country.code=NULL, adjusted=TRUE, ...) {
return(get.quantile.from.prediction(tfr.pred, quantile=0.5, country.index=country.index,
country.code=country.code, adjusted=adjusted, ...))
}
get.traj.quantiles <- function(tfr.pred, country.index, country.code, trajectories=NULL, pi=80,
adjusted=TRUE, est.uncertainty = FALSE) {
al <- (1-pi/100)/2
quantile.values <- as.numeric(dimnames(tfr.pred$quantiles)[[2]])
alidx<-round(quantile.values,6)==round(al,6)
cqp <- NULL
if (any(alidx)) { # pre-saved quantiles
alidx2 <- round(quantile.values,6)==round(1-al,6)
cqp <- rbind(tfr.pred$quantiles[country.index, alidx,],
tfr.pred$quantiles[country.index, alidx2,])
} else { # non-standard quantiles
reload <- FALSE
if (is.null(trajectories)) {
if(tfr.pred$nr.traj > 0) reload <- TRUE
} else {
if (dim(trajectories)[2] < 2000 && tfr.pred$nr.traj > dim(trajectories)[2]) reload <- TRUE
}
if(reload) {
#load 2000 trajectories maximum for computing quantiles
traj.reload <- get.trajectories(tfr.pred, country.code, 2000)
trajectories <- traj.reload$trajectories
}
if (!is.null(trajectories)) {
cqp <- apply(trajectories, 1,
quantile, c(al, 1-al), na.rm = TRUE)
}
}
if (est.uncertainty && has.est.uncertainty(tfr.pred$mcmc.set))
{ # replace quantiles from the first value (present year) with estimated uncertainty
tfr.est <- get.tfr.estimation(mcmc.list=tfr.pred$mcmc.set, country = country.code, probs=c(al, 1-al), adjust = adjusted)
unc.last.time <- which(tfr.est$tfr_quantile$year == dimnames(tfr.pred$quantiles)[[3]][1])
cqp[,1] <- unlist(tfr.est$tfr_quantile[unc.last.time, 1:2])
}
if(!adjusted) return(cqp)
shift <- get.tfr.shift(country.code, tfr.pred)
if(!is.null(shift)) cqp <- cqp + matrix(shift, nrow=nrow(cqp), ncol=ncol(cqp), byrow=TRUE)
return(cqp)
}
tfr.trajectories.plot.all <- function(tfr.pred,
output.dir=file.path(getwd(), 'TFRtrajectories'),
output.type="png", main=NULL, verbose=FALSE, ...) {
# plots TFR trajectories for all countries
.do.plot.all(tfr.pred$mcmc.set$meta, output.dir, tfr.trajectories.plot, output.type=output.type,
verbose=verbose, tfr.pred=tfr.pred, ...)
}
.do.plot.all.country.loop <- function(country.names, meta, output.dir, func, output.type="png",
file.prefix='TFRplot', plot.type='TFR graph', country.table=NULL,
main=NULL, verbose=FALSE, ...) {
if(!file.exists(output.dir)) dir.create(output.dir, recursive=TRUE)
postfix <- output.type
if(output.type=='postscript') postfix <- 'ps'
main.arg <- main
for (country in country.names) {
country.obj <- if(!is.null(meta)) get.country.object(country, meta)
else get.country.object(country, country.table=country.table)
if(verbose)
cat('Creating', plot.type, 'for', country, '(', country.obj$code, ')\n')
if(!is.null(main) && grepl('XXX', main, fixed=TRUE))
main.arg <- gsub('XXX', as.character(country.obj$name), main, fixed=TRUE)
do.call(output.type, list(file.path(output.dir,
paste(file.prefix,'_c', country.obj$code, '.', postfix, sep=''))))
do.call(func, list(country=country.obj$code, main=main.arg, ...))
dev.off()
}
if(verbose)
cat('\nPlots stored into', output.dir, '\n')
}
.do.plot.all <- function(meta, ...) {
# processes plotting function func for all countries
.do.plot.all.country.loop(country.names(meta), meta, ...)
}
get.half.child.variant <- function(median, increment=c(0, 0.25, 0.4, 0.5)) {
l <- length(median)
lincr <- length(increment)
upper <- lower <- c()
for (i in 1:l) {
upper <- c(upper, median[i]+increment[min(i,lincr)])
lower <- c(lower, median[i]-increment[min(i,lincr)])
}
return(rbind(lower, upper))
}
tfr.estimation.plot <- function(mcmc.list = NULL, country = NULL, sim.dir = NULL,
burnin = 0, thin = 1, pis = c(80, 95), plot.raw = TRUE,
grouping = 'source', save.image=TRUE, plot.dir = 'Estimation.plot',
adjust = TRUE, country.code = deprecated(), ISO.code = deprecated())
{
if (lifecycle::is_present(country.code)) {
lifecycle::deprecate_warn("7.1-0", "tfr.estimation.plot(country.code)", "tfr.estimation.plot(country)")
country <- country.code
}
if (lifecycle::is_present(ISO.code)) {
lifecycle::deprecate_warn("7.1-0", "tfr.estimation.plot(ISO.code)", "tfr.estimation.plot(country)")
country <- ISO.code
}
if (is.null(mcmc.list))
mcmc.list <- get.tfr.mcmc(sim.dir)
if (is.null(mcmc.list)) {
warning('MCMC does not exist.')
return(NULL)
}
if(inherits(mcmc.list, 'bayesTFR.prediction')) {
if(burnin != mcmc.list$burnin && burnin != 0)
warning('Prediction was generated with different burnin. Burnin set to ', mcmc.list$burnin)
burnin <- 0 # because burnin was already cut of the traces
mcmc.list <- mcmc.list$mcmc.set
}
meta <- mcmc.list$meta
if (is.null(mcmc.list$mcmc.list[[1]]$uncertainty) || !mcmc.list$mcmc.list[[1]]$uncertainty)
{
stop("MCMC does not consider uncertainty of past TFR.")
}
tfr.object <- get.tfr.estimation(mcmc.list=mcmc.list, country=country, sim.dir=sim.dir,
burnin=burnin, thin=thin, probs=sort(c((1-pis/100)/2, 0.5, pis/100 + (1-pis/100)/2)),
adjust = adjust)
country.obj <- get.country.object(country, meta)
quantile_tbl <- tfr.object$tfr_quantile
names(quantile_tbl)[1:(1 + 2 * length(pis))] <- paste0("Q", sort(c((100-pis)/2, 50, pis + (100-pis)/2)))
names.col <- paste0("Q", sort(c((100-pis)/2, 50, pis + (100-pis)/2)))
requireNamespace('ggplot2')
q <- ggplot2::ggplot(data=quantile_tbl) + ggplot2::xlab("year") + ggplot2::ylab("TFR")
q <- q + ggplot2::geom_ribbon(ggplot2::aes_string(x="year", ymin=names.col[1], ymax=names.col[length(names.col)]), alpha=0.2, fill='red') +
ggplot2::geom_line(ggplot2::aes_string(x="year", y="Q50"), size = 0.8, color="red") +
ggplot2::geom_point(ggplot2::aes_string(x="year", y="Q50"), size = 1, color="red") +
ggplot2::ggtitle(country.obj$name)
if (length(pis) > 1)
q <- q + ggplot2::geom_ribbon(ggplot2::aes_string(x="year", ymin=names.col[2], ymax=names.col[length(names.col)-1]), alpha=0.3, fill='red')
if (plot.raw)
{
if (country.obj$index %in% meta$extra) raw.data <- meta$raw_data_extra[[country.obj$index]]
else raw.data <- meta$raw_data.original[meta$raw_data.original$country_code == country.obj$code, ]
if(! grouping %in% colnames(raw.data) && !is.null(meta$covariates) && meta$covariates[1] %in% colnames(raw.data))
grouping <- meta$covariates[1]
if(grouping %in% colnames(raw.data)) {
ngroups <- t(unique(subset(raw.data, select=grouping)))
q <- q + ggplot2::geom_point(mapping = ggplot2::aes_string(x="year", y="tfr", color=grouping, shape=grouping),
data=raw.data, size=2.5) + ggplot2::scale_shape_manual(values=rep(15:18, len=length(ngroups)))
} else {
q <- q + ggplot2::geom_point(mapping = ggplot2::aes_string(x="year", y="tfr"), data=raw.data, size=2.5)
warning("No grouping of raw data used because column ", grouping, " not available. Use argument 'grouping' to group raw data.")
}
}
wpp.data <- get.observed.tfr(country.obj$index, meta, "tfr_matrix_all")
wpp.data <- data.frame(year=quantile_tbl$year, tfr = as.numeric(wpp.data))
#wpp.data <- wpp.data[wpp.data$year %% 5 == 3,]
q <- q + ggplot2::geom_line(data = wpp.data, ggplot2::aes_string(x="year", y="tfr"), size = 0.8) + ggplot2::theme_bw() +
ggplot2::geom_point(data = wpp.data, ggplot2::aes_string(x="year", y="tfr"), size = 0.7)
# re-draw the median
q <- q + ggplot2::geom_line(ggplot2::aes_string(x="year", y="Q50"), size = 0.8, color="red") +
ggplot2::geom_point(ggplot2::aes_string(x="year", y="Q50"), size = 1, color="red")
if (save.image)
{
if(!dir.exists(plot.dir)) dir.create(plot.dir)
pdf(file = file.path(plot.dir, paste0('tfr_country_', country.obj$code, ".pdf")), width=10, height=5)
print (q)
dev.off()
}
return(q)
}
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,
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'),
show.legend=TRUE, add=FALSE, uncertainty=FALSE, col_unc="purple", ...
) {
# lwd/col is a vector of 6 line widths/colors for:
# 1. observed data, 2. imputed missing data, 3. median, 4. quantiles, 5. +- 0.5 child, 6. trajectories
if (missing(country)) {
stop('Argument "country" must be given.')
}
if(length(lwd) < 6) {
lwd <- rep(lwd, 6)
lwd[6] <- 1
}
col <- .match.colors.with.default(col, c('black', 'green', 'red', 'red', 'blue', '#00000020'))
country.obj <- get.country.object(country, tfr.pred$mcmc.set$meta)
if(is.null(country.obj$code)) stop("Country ", country, " not found.")
if (uncertainty)
{
tfr.object <- get.tfr.estimation(mcmc.list=tfr.pred$mcmc.set, country = country.obj$code,
probs=sort(c((1-pi/100)/2, 0.5, pi/100 + (1-pi/100)/2)))
}
country <- country.obj
tfr_observed <- get.observed.tfr(country$index, tfr.pred$mcmc.set$meta, 'tfr_matrix_observed', 'tfr_matrix_all')
T_end_c <- tfr.pred$mcmc.set$meta$T_end_c
if(!is.null(tfr.pred$present.year.index.all)) T_end_c <- pmin(T_end_c, tfr.pred$present.year.index.all)
tfr_matrix_reconstructed <- get.tfr.reconstructed(tfr.pred$tfr_matrix_reconstructed, tfr.pred$mcmc.set$meta)
suppl.T <- length(tfr_observed) - tfr.pred$mcmc.set$meta$T_end
y1.part1 <- tfr_observed[1:T_end_c[country$index]]
y1.is.not.na <- which(!is.na(y1.part1))
y1.part1 <- y1.part1[y1.is.not.na]
lpart1 <- length(y1.part1)
y1.part2 <- NULL
lpart2 <- min(tfr.pred$mcmc.set$meta$T_end, tfr.pred$present.year.index) - T_end_c[country$index] + suppl.T
if (lpart2 > 0) {
p2idx <- (T_end_c[country$index]+1-suppl.T):nrow(tfr_matrix_reconstructed)
y1.part2 <- tfr_matrix_reconstructed[p2idx,country$index]
names(y1.part2) <- rownames(tfr_matrix_reconstructed)[p2idx]
}
x1 <- as.integer(c(names(y1.part1), names(y1.part2)))
x2 <- as.numeric(dimnames(tfr.pred$quantiles)[[3]])
trajectories <- get.trajectories(tfr.pred, country$code, nr.traj, typical.trajectory=typical.trajectory)
# plot historical data: observed
if (!add) {
if(is.null(xlim)) xlim <- c(min(x1,x2), max(x1,x2))
if(is.null(ylim))
{
ylim <- c(0, max(trajectories$trajectories, y1.part1, y1.part2, na.rm=TRUE))
if (uncertainty)
{
ylim[2] <- max(ylim[2], max(tfr.object$tfr_quantile[,-ncol(tfr.object$tfr_quantile), with = FALSE]))
}
}
if(is.null(main)) main <- country$name
plot(xlim, ylim, type='n', xlim=xlim, ylim=ylim, ylab=ylab, xlab=xlab, main=main,
panel.first = grid())
}
points.x <- x1[1:lpart1]
points.y <- y1.part1
if (mark.estimation.points) {
tfr.est <- get.observed.tfr(country$index, tfr.pred$mcmc.set$meta, 'tfr_matrix')[1:T_end_c[country$index]][y1.is.not.na]
elim.idx <- c()
# Phase I
end.na <- which(!is.na(tfr.est))
end.na <- if(length(end.na)==0) length(tfr.est) else end.na[1]
if(end.na > 1) {
na.idx <- 1:end.na
points(points.x[na.idx], points.y[na.idx], type=type, lwd=lwd[1], col=rgb(t(col2rgb(col[1])/255), alpha=0.1), ...)
elim.idx <- c(elim.idx, na.idx[-end.na])
}
# Phase III
p3.idx <- if(tfr.pred$mcmc.set$meta$lambda_c[country$index]>=length(tfr.est)) c() else seq(tfr.pred$mcmc.set$meta$lambda_c[country$index], length(tfr.est))
if(length(p3.idx) > 0) {
points(points.x[p3.idx], points.y[p3.idx], type=type, lwd=lwd[1], col=rgb(t(col2rgb(col[1])/255), alpha=0.3), pch = 4, ...)
elim.idx <- c(elim.idx, p3.idx)
}
if(length(elim.idx) > 0) {
points.x <- points.x[-elim.idx]
points.y <- points.y[-elim.idx]
}
}
if (!uncertainty || mark.estimation.points)
points(points.x, points.y, type=type, lwd=lwd[1], col=col[1], ...)
if(lpart2 > 0) { # imputed values
lines(x1[(lpart1+1): length(x1)], y1.part2, pch=2, type='b', col=col[2], lwd=lwd[2])
lines(x1[lpart1:(lpart1+1)], c(y1.part1[lpart1], y1.part2[1]), col=col[2], lwd=lwd[2]) # connection between the two parts
}
# plot trajectories
if(!is.null(trajectories$trajectories)) {
for (i in 1:length(trajectories$index)) {
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) {
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]
}
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))
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])
}
}
}
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
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)
}
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])
}
}
}
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')
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 (half.child.variant) {
lty <- c(lty, max(lty)+1)
llty <- length(lty)
up.low <- get.half.child.variant(median=tfr.median)
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')
cols <- c(cols, col[5])
lwds <- c(lwds, lwd[5])
}
if(show.legend) {
pch <- rep(-1, length(legend))
if (!uncertainty)
{
legend <- c(legend, 'observed TFR')
cols <- c(cols, col[1])
lty <- c(lty, 1)
pch <- c(pch, 1)
lwds <- c(lwds, lwd[1])
}
if(!uncertainty && (lpart2 > 0)) {
legend <- c(legend, 'imputed TFR')
cols <- c(cols, col[2])
lty <- c(lty, 1)
pch <- c(pch, 2)
lwds <- c(lwds, lwd[2])
}
if (uncertainty)
{
legend <- c(legend, if(adjusted.only) 'est. with uncertainty' else 'adj. estimates')
lty <- c(lty, 1)
pch <- c(pch, -1)
cols <- c(cols, col_unc)
lwds <- c(lwds, lwd[1])
if(!adjusted.only) {
legend <- c(legend, 'BHM estimates')
lty <- c(lty, unadj.lty)
cols <- c(cols, col_unc)
pch <- c(pch, -1)
lwds <- c(lwds, lwd[1])
}
}
legend('bottomleft', legend=legend, lty=lty, bty='n', col=cols, pch=pch, lwd=lwds)
}
}
extract.plot.args <- function(...) {
# split '...' into plot arguments and the rest
all.plot.args <- names(formals(plot.default))
args <- list(...)
which.plot.args <- pmatch(names(args), all.plot.args)
is.fun.arg <- is.na(which.plot.args)
return(list(plot.args=args[!is.fun.arg], other.args=args[is.fun.arg]))
}
do.plot.tfr.partraces <- function(mcmc.list, func, par.names, main.postfix='', chain.ids=NULL,
nr.points=NULL, dev.ncol=5, ...) {
mcmc.list <- get.mcmc.list(mcmc.list)
if (is.null(chain.ids)) {
nr.chains <- length(mcmc.list)
chain.ids <- 1:nr.chains
} else {
nr.chains <- length(chain.ids)
}
pars <- list()
iter <- rep(NA, nr.chains)
mclen <- rep(0, nr.chains)
# split '...' into function arguments and plot arguments
split.args <- extract.plot.args(...)
fun.args <- split.args$other.args
plot.args <- split.args$plot.args
thin <- fun.args$thin
fun.args$thin <- NULL
if(is.null(fun.args$burnin)) fun.args$burnin <- 0
orig.burnin <- fun.args$burnin
i <- 1
for(chain in chain.ids) {
mcmc <- mcmc.list[[chain]]
if (!is.null(thin) || mcmc$thin > 1) {
consolidated.burn.thin <- burn.and.thin(mcmc, orig.burnin,
if (is.null(thin)) mcmc$thin else thin)
fun.args$burnin <- consolidated.burn.thin$burnin
if (!is.null(consolidated.burn.thin$index)) fun.args$thinning.index <- consolidated.burn.thin$index
else {
if(!is.null(thin)) {
thin <- max(thin, mcmc$thin)
fun.args$thinning.index <- unique(round(seq(1,mcmc$length-consolidated.burn.thin$burnin,
by=thin/mcmc$thin)))
} else {
fun.args$thinning.index <- seq(1,mcmc$length-consolidated.burn.thin$burnin)
mclen[i] <- length(fun.args$thinning.index)
}
}
}
pars[[i]] <- eval(do.call(func, c(list(mcmc, par.names=par.names), fun.args)))
pars[[i]] <- filter.traces(pars[[i]], par.names)
iter[i] <- mcmc$finished.iter
if (i==1) {
par.names.l <- length(colnames(pars[[1]]))
maxy <- rep(NA, par.names.l)
miny <- rep(NA, par.names.l)
}
ipara<-1
for (para in colnames(pars[[i]])) {
maxy[ipara] <- max(maxy[ipara], pars[[i]][,para], na.rm=TRUE)
miny[ipara] <- min(miny[ipara], pars[[i]][,para], na.rm=TRUE)
ipara <- ipara+1
}
i <- i+1
}
if (par.names.l < dev.ncol) {
ncols <- par.names.l
nrows <- 1
} else {
ncols <- dev.ncol
nrows <- ceiling(par.names.l/dev.ncol)
}
par.cur <- par(mfrow=c(nrows,ncols))
col <- 1:nr.chains
maxx<-max(iter)
if(is.null(plot.args$xlim)) plot.args$xlim <- c(1+orig.burnin,maxx)
if(is.null(plot.args$xlab)) plot.args$xlab <- 'iterations'
if(is.null(plot.args$ylab)) plot.args$ylab <- ''
ylim <- plot.args$ylim
ipara <- 1
for (para in colnames(pars[[1]])) {
#mx <- length(pars[[1]][,para])
if (is.null(thin)) {
maxmclen <- max(mclen)
xindex <- if(maxmclen > 0) seq(1+orig.burnin, maxx, length=maxmclen)
else (1+orig.burnin):maxx
} else xindex <- seq(1+orig.burnin, maxx, by=thin)
thinpoints <- get.thinning.index(nr.points, length(xindex))
if (thinpoints$nr.points > 0) {
plot.args$ylim <- if(is.null(ylim)) c(miny[ipara], maxy[ipara]) else ylim
do.call('plot', c(list(xindex[thinpoints$index],
pars[[1]][thinpoints$index, para],
main=paste(para, main.postfix),
col=col[1], type='l'), plot.args))
if (nr.chains > 1) {
for (i in 2:nr.chains) {
lines(xindex[thinpoints$index],
pars[[i]][thinpoints$index,para], col=col[i], type='l')
}
}
}
ipara <- ipara+1
}
par(par.cur)
#stop('')
}
tfr.partraces.plot <- function(mcmc.list=NULL, sim.dir=file.path(getwd(), 'bayesTFR.output'),
chain.ids=NULL, par.names=tfr.parameter.names(trans=TRUE),
nr.points=NULL, dev.ncol=5, low.memory=TRUE, ...) {
if (is.null(mcmc.list))
mcmc.list <- get.tfr.mcmc(sim.dir, low.memory=low.memory)
do.plot.tfr.partraces(mcmc.list, 'load.tfr.parameter.traces', chain.ids=chain.ids,
nr.points=nr.points, par.names=par.names, dev.ncol=dev.ncol, ...)
}
tfr.partraces.cs.plot <- function(country, mcmc.list=NULL, sim.dir=file.path(getwd(), 'bayesTFR.output'),
chain.ids=NULL, par.names=tfr.parameter.names.cs(trans=TRUE),
nr.points=NULL, dev.ncol=3, low.memory=TRUE, ...) {
if (is.null(mcmc.list))
mcmc.list <- get.tfr.mcmc(sim.dir, low.memory=low.memory)
mcmc.list <- get.mcmc.list(mcmc.list)
country.obj <- get.country.object(country, mcmc.list[[1]]$meta)
if (is.null(country.obj$name))
stop('Country ', country, ' not found.')
stop.if.country.not.DL(country.obj, mcmc.list[[1]]$meta)
do.plot.tfr.partraces(mcmc.list, 'load.tfr.parameter.traces.cs',
main.postfix=paste('(',country.obj$name,')', sep=''), chain.ids=chain.ids, nr.points=nr.points,
country=country.obj$code, par.names=par.names, dev.ncol=dev.ncol, ...)
}
tfr3.partraces.plot <- function(mcmc.list=NULL, sim.dir=file.path(getwd(), 'bayesTFR.output'),
chain.ids=NULL, par.names=tfr3.parameter.names(),
nr.points=NULL, dev.ncol=3, low.memory=TRUE, ...) {
if (is.null(mcmc.list))
mcmc.list <- get.tfr3.mcmc(sim.dir, low.memory=low.memory)
else if(inherits(mcmc.list, 'bayesTFR.prediction'))
stop('Function not available for bayesTFR.prediction objects.')
tfr.partraces.plot(mcmc.list, sim.dir=NULL, chain.ids=chain.ids, par.names=par.names,
nr.points=nr.points, dev.ncol=dev.ncol, ...)
}
tfr3.partraces.cs.plot <- function(country, mcmc.list=NULL, sim.dir=file.path(getwd(), 'bayesTFR.output'),
chain.ids=NULL, par.names=tfr3.parameter.names.cs(),
nr.points=NULL, dev.ncol=2, low.memory=TRUE, ...) {
if (is.null(mcmc.list))
mcmc.list <- get.tfr3.mcmc(sim.dir, low.memory=low.memory)
else if(inherits(mcmc.list, 'bayesTFR.prediction'))
stop('Function not available for bayesTFR.prediction objects.')
mcmc.list <- get.mcmc.list(mcmc.list)
country.obj <- get.country.object(country, mcmc.list[[1]]$meta)
if (is.null(country.obj$name))
stop('Country ', country, ' not found.')
do.plot.tfr.partraces(mcmc.list, 'load.tfr.parameter.traces.cs',
main.postfix=paste('(',country.obj$name,')', sep=''), chain.ids=chain.ids, nr.points=nr.points,
country=country.obj$code, par.names=par.names, dev.ncol=dev.ncol, ...)
}
do.plot.tfr.pardensity <- function(mcmc.list, func, par.names, par.names.ext, main.postfix='',
func.args=NULL, chain.ids=NULL, burnin=NULL, dev.ncol=5, ...) {
if(inherits(mcmc.list, 'bayesTFR.prediction')) {
if(!is.null(burnin) && burnin != mcmc.list$burnin)
warning('Prediction was generated with different burnin. Burnin set to ', mcmc.list$burnin,
'.\n Use a bayesTFR.mcmc.set object as the first argument, if the original traces should be used.')
burnin <- 0 # because burnin was already cut of the traces
if (!is.null(chain.ids) && max(chain.ids) > 1) {
warning('Thinned traces from all chains used for plotting density.\n For selecting individual chains, use a bayesTFR.mcmc.set object as the first argument.')
chain.ids <- NULL
}
}
if(is.null(burnin)) burnin <- 0
mcmc.list <- get.mcmc.list(mcmc.list)
if (!is.null(chain.ids)) mcmc.list <- mcmc.list[chain.ids]
par.names.l <- length(par.names.ext)
if (par.names.l < dev.ncol) {
ncols <- par.names.l
nrows <- 1
} else {
ncols <- dev.ncol
nrows <- ceiling(par.names.l/dev.ncol)
}
args <- extract.plot.args(...)
par.cur <- par(mfrow=c(nrows,ncols))
tfr_flag <- FALSE
for (para in par.names) {
if (!tfr_flag && length(grep('tfr_*', para) > 0))
{
para <- 'tfr'
tfr_flag <- TRUE
}
values <- eval(do.call(func, c(list(mcmc.list, par.names=para, burnin=burnin), func.args)))
values <- filter.traces(values, par.names)
for (par.name in colnames(values)) {
dens <- do.call('density', c(list(values[,par.name]), args$other.args))
do.call('plot', c(list(dens, main=paste(par.name, main.postfix)), args$plot.args))
}
}
par(par.cur)
}
tfr.pardensity.plot <- function(mcmc.list=NULL, sim.dir=file.path(getwd(), 'bayesTFR.output'),
chain.ids=NULL, par.names=tfr.parameter.names(trans=TRUE),
burnin=NULL, dev.ncol=5, low.memory=TRUE, ...) {
if (is.null(mcmc.list))
mcmc.list <- get.tfr.mcmc(sim.dir, low.memory=low.memory)
par.names.ext <- get.full.par.names(par.names, tfr.parameter.names.extended())
if(length(par.names.ext) <= 0)
stop('Parameter names are not valid parameters.\nUse function tfr.parameter.names(...) or valid parameter names.')
do.plot.tfr.pardensity(mcmc.list, 'get.tfr.parameter.traces', chain.ids=chain.ids, par.names=par.names,
par.names.ext=par.names.ext,
burnin=burnin, dev.ncol=dev.ncol, ...)
}
tfr.pardensity.cs.plot <- function(country, mcmc.list=NULL, sim.dir=file.path(getwd(), 'bayesTFR.output'),
chain.ids=NULL, par.names=tfr.parameter.names.cs(trans=TRUE),
burnin=NULL, dev.ncol=3, low.memory=TRUE, ...) {
if (is.null(mcmc.list))
mcmc.list <- get.tfr.mcmc(sim.dir, low.memory=low.memory)
mcmc.l <- get.mcmc.list(mcmc.list)
country.obj <- get.country.object(country, mcmc.l[[1]]$meta)
if (is.null(country.obj$name))
stop('Country ', country, ' not found.')
stop.if.country.not.DL(country.obj, mcmc.l[[1]]$meta)
par.names.ext <- get.full.par.names.cs(par.names,
tfr.parameter.names.cs.extended(country.obj$code))
if(length(par.names.ext) <= 0 && length(grep('tfr_*', par.names)) <= 0)
stop('Parameter names are not valid country-specific parameters.\nUse function tfr.parameter.names.cs(...) or valid parameter names.')
else if (length(par.names.ext) <= 0)
par.names.ext <- paste0('tfr_c', country.obj$code)
do.plot.tfr.pardensity(mcmc.list, 'get.tfr.parameter.traces.cs', chain.ids=chain.ids, par.names=par.names,
par.names.ext=par.names.ext,
main.postfix=paste('(',country.obj$name,')', sep=''),
func.args=list(country.obj=country.obj),
burnin=burnin, dev.ncol=dev.ncol, ...)
}
tfr3.pardensity.plot <- function(mcmc.list=NULL, sim.dir=file.path(getwd(), 'bayesTFR.output'),
chain.ids=NULL, par.names=tfr3.parameter.names(),
burnin=NULL, dev.ncol=3, low.memory=TRUE, ...) {
if (is.null(mcmc.list))
mcmc.list <- get.tfr3.mcmc(sim.dir, low.memory=low.memory)
else if(inherits(mcmc.list, 'bayesTFR.prediction'))
stop('Function not available for bayesTFR.prediction objects.')
do.plot.tfr.pardensity(mcmc.list, 'get.tfr.parameter.traces', chain.ids=chain.ids, par.names=par.names,
par.names.ext=par.names,
burnin=burnin, dev.ncol=dev.ncol, ...)
}
tfr3.pardensity.cs.plot <- function(country, mcmc.list=NULL, sim.dir=file.path(getwd(), 'bayesTFR.output'),
chain.ids=NULL, par.names=tfr3.parameter.names.cs(),
burnin=NULL, dev.ncol=2, low.memory=TRUE, ...) {
if (is.null(mcmc.list))
mcmc.list <- get.tfr3.mcmc(sim.dir, low.memory=low.memory)
else if(inherits(mcmc.list, 'bayesTFR.prediction'))
stop('Function not available for bayesTFR.prediction objects.')
mcmc.l <- get.mcmc.list(mcmc.list)
country.obj <- get.country.object(country, mcmc.l[[1]]$meta)
if (is.null(country.obj$name))
stop('Country ', country, ' not found.')
par.names.ext <- get.full.par.names.cs(par.names, paste(par.names, '_c', country.obj$code, sep=''))
if(length(par.names.ext) <= 0)
stop('Parameter names are not valid country-specific parameters.\nUse function tfr3.parameter.names.cs(...) or valid parameter names.')
do.plot.tfr.pardensity(mcmc.list, 'get.tfr.parameter.traces.cs', chain.ids=chain.ids, par.names=par.names,
par.names.ext=par.names.ext,
main.postfix=paste('(',country.obj$name,')', sep=''),
func.args=list(country.obj=country.obj),
burnin=burnin, dev.ncol=dev.ncol, ...)
}
".get.gamma.pars" <- function(pred, ...) UseMethod (".get.gamma.pars")
.get.gamma.pars.bayesTFR.prediction <- function(pred, ...) {
# estimated by
# library(MASS)
# data <- pred$tfr_matrix_reconstructed[12,]
# gd <- fitdistr(data-min(data)+0.05, densfun='gamma')
# min(data) is 0.95
return(list(gamma.pars=list(shape=1.87, rate=0.94), gamma.shift=0.95-0.05, min.value=0.7,
max.value=NULL))
}
get.tfr.map.parameters <- function(pred, tfr.range=NULL, nr.cats=50, same.scale=TRUE,
quantile=0.5, ...) {
map.pars <- list(pred=pred, quantile=quantile, ...)
if (same.scale) {
gp <- .get.gamma.pars(pred)
data <- pred$quantiles[,as.character(quantile),1]
q <- if(is.null(tfr.range)) c(min(pmax(data,gp$gamma.shift)), max(data))-gp$gamma.shift
else c(max(gp$gamma.shift, tfr.range[1]-gp$gamma.shift), max(gp$gamma.shift, tfr.range[2]-gp$gamma.shift))
min.max.q <- pgamma(q, shape=gp$gamma.pars[['shape']], rate=gp$gamma.pars[['rate']])
quantiles <- seq(min.max.q[1], min.max.q[2], length=nr.cats)
quant.values <- c(gp$min.value,
qgamma(quantiles, shape=gp$gamma.pars[['shape']], rate=gp$gamma.pars[['rate']])+gp$gamma.shift)
if(!is.null(gp$max.value) && gp$max.value > max(quant.values)) quant.values <- c(quant.values, gp$max.value)
if(!is.null(tfr.range)) {
if(tfr.range[1] < quant.values[1])
quant.values <- c(tfr.range[1], quant.values)
if(tfr.range[1] > quant.values[1])
quant.values <- quant.values[quant.values >= tfr.range[1]]
last <- length(quant.values)
if(tfr.range[2] > quant.values[last])
quant.values <- c(quant.values, tfr.range[2])
if(tfr.range[2] < quant.values[last])
quant.values <- quant.values[quant.values <= tfr.range[2]]
}
col <- rainbow(500, start=0, end=0.67)[seq(500, 1, length=length(quant.values)-1)]
map.pars$catMethod <- quant.values
} else {
col <- rainbow(500, start=0, end=0.67)[seq(500, 1, length=nr.cats)]
map.pars$numCats <- nr.cats
}
map.pars$colourPalette <- col
return(map.pars)
}
bdem.map.all <- function(pred, output.dir, type='tfr', output.type='png', range=NULL, nr.cats=50, same.scale=TRUE,
quantile=0.5, file.prefix='TFRwrldmap_', ...) {
if(!file.exists(output.dir)) dir.create(output.dir, recursive=TRUE)
all.years <- get.all.prediction.years(pred)
output.args <- list()
postfix <- output.type
if(output.type=='postscript') postfix <- 'ps'
filename.arg <- 'filename'
if(output.type=='postscript' || output.type=='pdf') {filename.arg<-'file'}
else{output.args[['width']] <- 1000}
map.pars <- get.tfr.map.parameters(pred, tfr.range=range, nr.cats=nr.cats, same.scale=same.scale,
quantile=quantile, ...)
for (year in all.years) {
output.args[[filename.arg]] <- file.path(output.dir,
paste(file.prefix, year, '.', postfix, sep=''))
do.call(paste(type, '.map', sep=''), c(list(year=year, device=output.type,
device.args=output.args), map.pars))
dev.off()
}
cat('Maps written into', output.dir, '\n')
}
tfr.map.all <- function(pred, output.dir, output.type='png', tfr.range=NULL, nr.cats=50, same.scale=TRUE,
quantile=0.5, file.prefix='TFRwrldmap_', ...) {
bdem.map.all(pred=pred, output.dir=output.dir, type='tfr', output.type=output.type, range=tfr.range,
nr.cats=nr.cats, same.scale=same.scale, quantile=quantile, file.prefix=file.prefix, ...)
}
".map.main.default" <- function(pred, ...) UseMethod (".map.main.default")
.map.main.default.bayesTFR.prediction <- function(pred, ...) return('TFR: quantile')
par.names.for.worldmap <- function(pred, ...) UseMethod ("par.names.for.worldmap")
par.names.for.worldmap.bayesTFR.prediction <- function(pred, ...) {
return(c('lambda', tfr.parameter.names.cs.extended()))
}
"get.data.for.worldmap" <- function(pred, ...) UseMethod ("get.data.for.worldmap")
get.data.for.worldmap.bayesTFR.prediction <- function(pred, quantile=0.5, year=NULL, par.name=NULL,
adjusted=FALSE, projection.index=1, pi=NULL, ...) {
meta <- pred$mcmc.set$meta
quantiles <- quantile
if (!is.null(pi)) {
qlower <- (1-pi/100)/2
quantiles <- c(quantile, qlower, 1-qlower)
}
if(!is.null(par.name)) { # data are parameter values
if (!is.element(par.name, par.names.for.worldmap(pred)))
stop('Illegal par.name. Allowed values:',
paste(par.names.for.worldmap(pred), collapse=', '))
data <- c()
if (par.name == 'lambda') {
tfr <- get.data.imputed(pred)
tfr.years <- get.estimation.years(meta)
all.years <- c(tfr.years, get.prediction.years(meta, pred$nr.projections+1, pred$present.year.index)[-1])
nr.data <- pred$nr.projections+dim(tfr)[1]
for (country in 1:get.nr.countries(meta)) {
country.obj <- get.country.object(country, meta, index=TRUE)
tfr.and.pred.median <- c(tfr[,country],
get.quantile.from.prediction(pred, quantile, country.obj$index, country.obj$code,
adjusted=adjusted)[-1])
lambda <- all.years[find.lambda.for.one.country(tfr.and.pred.median, nr.data)]
data <- c(data, lambda)
}
codes <- meta$regions$country_code
} else {
con <- textConnection("sout", "w", local=TRUE) # redirect output (to get rid of coda's messages)
for (country in get.countries.index(meta)) {
country.obj <- get.country.object(country, meta, index=TRUE)
sink(con, type='message')
s <- summary(coda.list.mcmc(pred$mcmc.set, country=country.obj$code,
par.names=NULL, par.names.cs=par.name, thin=1, burnin=0), quantiles = quantiles)
sink(type='message')
data <- rbind(data, s$quantiles)
}
close(con)
codes <- meta$regions$country_code[get.countries.index(meta)]
}
projection.index <- 1
projection <- TRUE
period <- paste('Parameter', par.name, 'for')
} else { # data are TFRs
projection <- TRUE
if(!is.null(year)) {
ind.proj <- get.predORest.year.index(pred, year)
if(! 'index' %in% names(ind.proj))
stop('Projection year ', year, ' not found.')
projection.index <- ind.proj['index']
projection <- ind.proj['is.projection']
}
if(projection) {
if(!all(is.element(as.character(quantiles), dimnames(pred$quantiles)[[2]])))
stop('Some of the quantiles ', paste(quantiles, collapse=', '), ' not found.\nAvailable: ',
paste(dimnames(pred$quantiles)[[2]], collapse=', '),
'\nCheck arguments "quantile" and "pi".')
data <- pred$quantiles[, as.character(quantiles), projection.index]
if(adjusted) data <- data + get.tfr.shift.all(pred, projection.index)
period <- get.prediction.periods(meta, projection.index,
present.year.index=pred$present.year.index)[projection.index]
} else {
data <- get.data.imputed(pred)[projection.index, ]
period <- get.tfr.periods(meta)[projection.index]
}
codes <- meta$regions$country_code
}
if(adjusted) period <- paste(period, 'adjusted')
rownames(data) <- NULL
low<-NULL
up<-NULL
res <- data
if(!is.null(dim(data))) {
res <- data[,1]
if(ncol(data) > 1) {
low <- data[,2]
up <- data[,3]
}
}
return(list(period=period, data=res, country.codes=codes, lower=low, upper=up))
}
tfr.map <- function(pred, quantile=0.5, year=NULL, par.name=NULL, adjusted=FALSE,
projection.index=1, device='dev.new', main=NULL,
resolution=c("coarse","low","less islands","li","high"),
device.args=NULL, data.args=NULL, ...
) {
resolution <- match.arg(resolution)
#if(resolution=='high') require(rworldxtra)
data.period <- do.call(get.data.for.worldmap, c(list(pred, quantile, year=year,
par.name=par.name, adjusted=adjusted, projection.index=projection.index), data.args))
#data.period.base <- do.call(get.data.for.worldmap, c(list(pred, quantile, year=2013,
# par.name=par.name, adjusted=adjusted, projection.index=projection.index), data.args))
#data <- (data.period$data - data.period.base$data)/1000
data <- data.period$data
period <- data.period$period
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'), ]
#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
sPDF$UN[sPDF$ISO3=="CYN"] <- 196
## Kosovo -> assign to Serbia
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 (SOM) -> fixed in rworlmap version 1.3.6
#sPDF$UN[sPDF$ISO3=="SOL"] <- 706
#mtfr <- joinCountryData2Map(tfr, joinCode='UN', nameJoinColumn='un')
# join sPDF with tfr
mtfr <- rep(NA, length(sPDF$UN))
valididx <- which(is.element(sPDF$UN, tfr$un))
mtfr[valididx] <- tfr$tfr[sapply(sPDF$UN[valididx], function(x,y) which(y==x), tfr$un)]
sPDF$tfr <- mtfr
if(is.null(main)) {
main <- paste(period, .map.main.default(pred, data.period), quantile)
}
if (device != 'dev.cur')
do.call(rworldmap::mapDevice, c(list(device=device), device.args))
mapParams<-rworldmap::mapCountryData(sPDF, nameColumnToPlot='tfr', addLegend=FALSE, mapTitle=main, ...
)
# Default for legendIntervals changed in rworlmap 1.3.6 from "page" to "data". Therefore need to pass it explicitly here.
do.call(rworldmap::addMapLegend, c(mapParams, legendWidth=0.5, legendMar=2, legendLabels='all', legendIntervals = "page"))
#do.call(addMapLegend, c(mapParams, legendWidth=0.5, legendMar=2, legendLabels='all', sigFigs=2, legendShrink=0.8, tcl=-0.3, digits=1))
}
tfr.ggmap <- function(pred, quantile=0.5, year=NULL, par.name=NULL, adjusted=FALSE,
projection.index=1, main=NULL, data.args=NULL, viridis.option = "B",
nr.cats = 10, same.scale = FALSE, plot = TRUE, file.name = NULL, plot.size = 4, ...
) {
# function for quantile transformation
make_quantile_trans <- function(x, format = scales::label_number()) {
name <- paste0("quantiles_of_", deparse1(substitute(x)))
xs <- sort(x)
N <- length(xs)
transform <- function(x) findInterval(x, xs)/N # find the last element that is smaller
inverse <- function(q) xs[1+floor(q*(N-1))]
scales::trans_new(
name = name,
transform = transform,
inverse = inverse,
breaks = function(x, n = 5) inverse(scales::extended_breaks()(transform(x), n)),
minor_breaks = function(x, n = 5) inverse(scales::regular_minor_breaks()(transform(x), n)),
format = format,
domain = xs[c(1, N)]
)
}
for(pkg in c("ggplot2", "sf", "spData", "scales"))
requireNamespace(pkg)
data.period <- do.call(get.data.for.worldmap, c(list(pred, quantile, year=year,
par.name=par.name, adjusted=adjusted, projection.index=projection.index), data.args))
data <- data.period$data
period <- data.period$period
tfr <- data.frame(cbind(un=data.period$country.codes, tfr=data))
e <- new.env(parent = emptyenv())
data("iso3166", envir=e)
data("world", package = "spData", envir=e)
tfr <- merge(tfr, e$iso3166[, c("uncode", "charcode")], by.x = "un", by.y = "uncode")
#e$world <- e$world[- which(!e$world$iso_a2 %in% tfr$charcode),]
# align with UN countries
e$world <- e$world[- which(e$world$name_long == "Antarctica"), c("iso_a2", "name_long", "geom")] # remove Antarctica
e$world$iso_a2[e$world$name_long == "Somaliland"] <- "SO" # set Somaliland as Somalia
e$world$iso_a2[e$world$name_long == "Northern Cyprus"] <- "CY" # set Northern Cyprus as Cyprus
tfr <- tfr[tfr$charcode %in% e$world$iso_a2,]
e$world <- merge(e$world, tfr, by.x = "iso_a2", by.y = "charcode", all = TRUE)
if(same.scale & is.null(par.name)) {
all.data <- pred$quantiles[,as.character(quantile),]
all.data <- all.data[data.period$country.codes %in% e$world$un]
} else all.data <- e$world$tfr
if(is.null(main))
main <- paste(period, .map.main.default(pred, data.period), quantile)
# g <- ggplot(world) + geom_sf(aes(fill = tfr), colour = "grey", lwd = 0.1) +
# #scale_fill_viridis_c(option = "A", direction = -1, breaks = scales::breaks_pretty()) +
# theme(legend.position="bottom") + # coord_sf(default_crs = "+proj=robin +ellps=WGS84", label_axes = "--EN")
# scale_fill_distiller(palette = "YlOrRd", direction = 1, breaks = round(quantile(world$tfr, probs = seq(0, 1, length = 10)), 2),
# #breaks = scales::breaks_pretty(n = 10),
# trans = make_quantile_trans(world$tfr)) # trans = "reverse",
# g <- g + theme(axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank())
# g2 <- g + guides(fill = guide_colourbar(title = "", barwidth = unit(0.5, "npc", data = g), barheight = 0.5,
# direction = "horizontal"))
#g5 <- g2 + ggtitle("(3) Palette: gradient from yellow to red") +
# scale_fill_gradient(low = "yellow", high = "red", breaks = round(quantile(world$tfr, probs = seq(0, 1, length = 10)), 2), trans = make_quantile_trans(world$tfr))
world.rob <- sf::st_transform(e$world, "+proj=robin +ellps=WGS84")
grobin <- ggplot2::ggplot(world.rob) + ggplot2::geom_sf(ggplot2::aes_string(fill = "tfr"), colour = "grey", lwd = 0.1, ...) +
ggplot2::coord_sf(datum = NA) +
ggplot2::scale_fill_viridis_c(option = viridis.option, direction = -1, na.value= "white",
breaks = round(quantile(all.data, probs = seq(0, 1, length = nr.cats), na.rm = TRUE), 2),
trans = make_quantile_trans(all.data),
limits = range(all.data)) +
ggplot2::theme(legend.position="bottom", axis.text.x = ggplot2::element_blank(), axis.text.y = ggplot2::element_blank(),
axis.ticks = ggplot2::element_blank())
grobin <- grobin + ggplot2::guides(fill = ggplot2::guide_colourbar(title = "", barwidth = ggplot2::unit(0.5, "npc", data = grobin),
barheight = 0.5, direction = "horizontal")) +
ggplot2::ggtitle(main)
if(plot == TRUE){
plot_ratio <- 2.360463 # derived via library(tmaptools); get_asp_ratio(world.rob)
if(is.null(file.name)){
grDevices::dev.new(width = plot_ratio*plot.size, height = plot.size)
print(grobin)
} else {
ggplot2::ggsave(file.name, grobin, width = plot_ratio*plot.size, height = plot.size)
}
}
invisible(grobin)
}
tfr.map.gvis <- function(pred, year=NULL, quantile=0.5, pi=80, par.name=NULL,
adjusted=FALSE, ...)
bdem.map.gvis(pred, year=year,
quantile=quantile, pi=pi, par.name=par.name, adjusted=adjusted, ...)
"bdem.map.gvis" <- function(pred, ...) UseMethod ("bdem.map.gvis")
bdem.map.gvis.bayesTFR.prediction <- function(pred, year=NULL, quantile=0.5, pi=80,
par.name=NULL, html.file=NULL, adjusted=FALSE, ...) {
.do.gvis.bdem.map('TFR', 'BHM for Total Fertility Rate<br>', pred, year=year,
quantile=quantile, pi=pi, par.name=par.name, adjusted=adjusted, ...)
}
.do.gvis.bdem.map <- function(what, title, pred, year=NULL, quantile=0.5, pi=80,
par.name=NULL, adjusted=FALSE, ...) {
data.period <- get.data.for.worldmap(pred, quantile, year=year,
par.name=par.name, projection.index=1, adjusted=adjusted, pi=pi, ...)
mapdata <- data.period$data
period <- data.period$period
lower <- data.period$lower
upper <- data.period$upper
un <- data.period$country.codes
countries.table <- get.countries.table(pred)
if(!is.null(par.name)) what <- par.name
e <- new.env(parent = emptyenv())
data('iso3166', envir=e)
unmatch <- match(un, e$iso3166$uncode)
unidx <- which(!is.na(unmatch))
ct.idx <- sapply(un[unidx], function(x, y) which(y==x), countries.table$code)
country.names <- countries.table$name[ct.idx]
#remove problematic characters
country.names <- iconv(country.names, "latin1", "ASCII", "?")
data <- data.frame(un=un[unidx], name=country.names,
iso=e$iso3166$charcode[unmatch][unidx])
data[[what]] <- mapdata[unidx]
if(!is.null(lower)) { # confidence intervals defined
lower.name <- paste('lower_', pi, sep='')
upper.name <- paste('upper_', pi, sep='')
data[[lower.name]] <- round(lower[unidx], 2)
data[[upper.name]] <- round(upper[unidx], 2)
data$pi <- paste(e$iso3166$charcode[unmatch][unidx], ': ', pi, '% PI (', data[[lower.name]], ', ',
data[[upper.name]], ')', sep='')
hovervar <- 'pi'
} else { # no confidence intervals
lower.name <- 'lower'
upper.name <- 'upper'
data[[lower.name]] <- data[[upper.name]] <- rep(NA, length(unidx))
hovervar <- ''
}
col <- c('#0000CC', '#00CCFF', '#33FF66', '#FFFF66', '#FF9900', '#FF3300')
geo <- googleVis::gvisGeoChart(data, locationvar="iso", colorvar=what, hovervar=hovervar,
options=list(height=500, width=900, dataMode='regions',
colors=paste('[', paste(shQuote(col), collapse=', '), ']')))
#geo$html$caption <- paste(title, 'in', period,'<br>\n')
geo$html$caption <- paste(title, period, .map.main.default(pred, data.period), quantile)
bdem.data <- data[,c('un', 'iso', 'name', what, lower.name, upper.name)]
gvis.table <- googleVis::gvisTable(bdem.data,
options=list(width=600, height=600, page='disable', pageSize=198))
page <- list(type="MapTable",
chartid=format(Sys.time(), "BdemMap-%Y-%m-%d-%H-%M-%S"),
html=list(Header=geo$html$header,
Chart1=geo$html$chart,
Caption1=geo$html$caption,
Chart2=gvis.table$html$chart,
Caption2=gvis.table$html$caption,
Footer=gvis.table$html$footer)
)
class(page) <- list("gvis", class(page))
plot(page)
invisible(page)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.