set.seed(seed <- 1) library("bayesPop") # ============ Section 3.2 ============ # setting iter to 1000 and nr.traj to 100 takes about 1/2 hour for all three steps # to complete if run sequentially. Step 1 and 2 can be run at the same time. iter <- 1000 # number of MCMC iterations nr.traj <- 100 # number of population trajectories # 1. TFR estimation and prediction sim.dir.tfr <- file.path(getwd(), "TFRprojections") run.tfr.mcmc(iter = iter, nr.chains = 1, thin = 1, output.dir = sim.dir.tfr, seed = seed) run.tfr3.mcmc(sim.dir = sim.dir.tfr, iter = iter, nr.chains = 1, thin = 1, seed = seed) tfr.predict(sim.dir = sim.dir.tfr, nr.traj = nr.traj, burnin = iter/2, burnin3 = iter/2, seed = seed) # 2. e0 estimation and prediction sim.dir.e0 <- file.path(getwd(), "e0_projections") run.e0.mcmc(sex = "Female", iter = iter, nr.chains = 1, thin = 1, output.dir = sim.dir.e0, seed = seed) e0.predict(sim.dir = sim.dir.e0, nr.traj = nr.traj, burnin = iter/2, seed = seed) # Will generate a few warnings about NaNs produced which can be ignored. Note # that the above function does not generate projections for HIV/AIDS countries. # 3. Population prediction sim.dir.pop <- file.path(getwd(), "pop_projections") pop.pred <- pop.predict(output.dir = sim.dir.pop, inputs = list(tfr.sim.dir = sim.dir.tfr, e0F.sim.dir = sim.dir.e0, e0M.sim.dir = "joint_"), keep.vital.events = TRUE, verbose = TRUE) # Will generate warnings for HIV/AIDS countries that e0 not available. # retrieve the prediction object pop.pred <- get.pop.prediction(sim.dir.pop) # ============ Section 3.3 ============ country <- "Netherlands" summary(pop.pred, country) # Figure 3 left pdf("traj_Netherlands.pdf", height = 5, width = 5) par(mgp = c(1.7, .5, 0), mar = c(1.6, 2.6, 1.5, .1)) pop.trajectories.plot(pop.pred, country = country, sum.over.ages = TRUE) dev.off() # Figure 3 right pdf("traj_Netherlands_male0_14.pdf", height = 5, width = 5) par(mgp = c(1.7, .5, 0), mar = c(1.6, 2.6, 1.5, .1)) pop.trajectories.plot(pop.pred, country = country, sex = "male", age = 1:3, sum.over.ages = TRUE) dev.off() # Figure 4 left pdf("trajbyage_Netherlands_2100.pdf", height = 5, width = 5) par(mgp = c(1.7, .5, 0), mar = c(3.8, 2.6, 1.5, .1)) pop.byage.plot(pop.pred, country = country, year = 2100) dev.off() # Figure 4 right pdf("trajbyage_Netherlands_1960_2100.pdf", height = 5, width = 5) par(mgp = c(1.7, .5, 0), mar = c(3.8, 2.6, 1.5, .1)) pop.byage.plot(pop.pred, country = country, year = 2060, pi = 80, nr.traj = 50) pop.byage.plot(pop.pred, country = country, year = 1960, add = TRUE, col = "blue", show.legend = FALSE) dev.off() # Cohort projections pop.cohorts.plot(pop.pred, country = country) cohort.data <- cohorts(pop.pred, country = country) head(names(cohort.data)) # Figure 5 pdf("cohorts_Netherlands.pdf", height = 3, width = 8) pop.cohorts.plot(pop.pred, cohort.data = cohort.data, cohorts = c(1980, 2000, 2020)) dev.off() # Age categories which(pop.pred$ages < 15) pop.pred$ages[1:3] # ============ Section 4 ========== # Classic pyramid, Figure 6 left pdf("pyr_Netherlands.pdf", height = 5, width = 5) pop.pyramid(pop.pred, country, c(2100, 2015), age = 1:23) dev.off() # Trajectories pyramid, Figure 6 right pdf("pyrtraj_Netherlands.pdf", height = 5, width = 5) pop.trajectories.pyramid(pop.pred, country, c(2100, 2025, 1950), age = 1:23, pi = 95, nr.traj = 0, proportion = TRUE) dev.off() # Generate pyramids for all countries and two sets of years pop.pyramidAll(pop.pred, year = list(c(2100, 2015), c(2050, 2015)), age = 1:23, output.dir = "mypyramids") # Pyramids for user-defined data data <- read.table(file.path(find.package("bayesPop"), "ex-data", "popestimates_WAKing.txt"), header = TRUE, row.names = 1) head(data) WA <- data[, c("WA.male", "WA.female")] colnames(WA) <- c("M", "F") Ki <- data[, c("King.male", "King.female")] colnames(Ki) <- c("M", "F") pyr <- get.bPop.pyramid(list(WA, Ki), LRcolnames = c("M", "F"), legend = c("Washington", "King County")) # Figure 7 left pdf("pyr_WAKing.pdf", height = 5, width = 6) plot(pyr, main = "Population in 2011", pyr2.par = list(height = 0.7, col = "violet", border = "violet")) dev.off() # Figure 7 right pyr.prop <- get.bPop.pyramid(list(WA, Ki), is.proportion = NA, LRcolnames = c("M", "F"), legend = c("Washington", "King County")) pdf("pyr_WAKing_prop.pdf", height = 5, width = 6) pop.pyramid(pyr.prop, main = "Population in 2011 (proportions)", pyr1.par = list(col = "lightgreen", border = "lightgreen", density = 30), pyr2.par = list(col = "darkred", border = "darkred", density = 50, height = 0.3)) dev.off() # ============ Section 5 ========== # Checking dimension of basic component B <- get.pop("BNL{}", pop.pred, observed = FALSE) P <- get.pop("PNL{4:10}", pop.pred, observed = FALSE) all(dim(B) == dim(P)) # Figure 8 left expr <- "pop.apply(PDE_F{4:10}, gmedian, cats=seq(15, by=5, length=8))" pdf("traj_expr_median_childbearingage.pdf", height = 5, width = 5) par(mgp = c(1.7, .5, 0), mar = c(1.6, 1.6, 1.5, .1)) pop.trajectories.plot(pop.pred, nr.traj = 20, expression = expr, main = "Median age of women in childbearing ages") dev.off() # Figure 8 right expr <- mac.expression(country = "DE") expr pdf("plot_mac.pdf", height = 5, width = 5) par(mgp = c(1.7, .5, 0), mar = c(1.6, 1.6, 1.5, .5)) pop.trajectories.plot(pop.pred, nr.traj = 20, expression = expr, main = "Mean age of childbearing") dev.off() # Figure 9 left pdf("traj_expr_births_per_woman.pdf", height = 5, width = 5) par(mgp = c(1.7, .5, 0), mar = c(3, 1.6, 1.5, .1)) pop.byage.plot(pop.pred, expression = "BDE{} / PDE_F{4:10}", year = 2030, nr.traj = 20, main = "Births by mother's age per woman - 2030") dev.off() # Figure 9 right pdf("plot_cohort_births.pdf", height = 5, width = 5) par(mgp = c(1.7, .5, 0), mar = c(3, 1.6, 1.5, .1)) pop.cohorts.plot(pop.pred, expression = "BDE{}", cohorts = 2015, legend.pos = "topleft", main = "Births for female cohort 2015-2020") dev.off() # map (will take some time when called for the first time) pop.map(pop.pred, expression = "MXXX[-1]", year = 2050) # ============ Section 6 ========== # generate aggregations for World, Africa, Europe, Northern # America, Asia, Latin Am. pop.aggr <- pop.aggregate(pop.pred, input.type = "country", regions = c(900, 903, 908, 905, 935, 904), verbose = TRUE) pop.aggr$countries pop.aggr <- get.pop.aggregation(sim.dir.pop, name = "country") # Figure 10 pdf("traj_aggr.pdf", width = 10, height = 6) par(mfrow = c(2,2), mgp = c(1.7, .5, 0), mar = c(1.6, 2.6, 1.5, .1)) for (country in c(903, 908, 905, 935)) pop.trajectories.plot(pop.aggr, country, sum.over.ages = TRUE) dev.off() # Figure 11 left pdf("pyr_world.pdf", width = 6, height = 5) pop.pyramid(pop.aggr, 900, year = c(2100, 2015), proportion = TRUE) dev.off() # Figure 11 right pdf("trajbyage_aggr.pdf", width = 6, height = 5) par(mgp = c(1.7, .5, 0), mar = c(3.8, 1.6, 1.5, .1)) pop.byage.plot(pop.aggr, expression = "P908{} / P900{}", year = 2100, main = "Proportion of Europeans to world population", pi = 80, ylim = c(0, 0.35), show.legend = FALSE) pop.byage.plot(pop.aggr, expression = "P908{} / P900{}", year = 2015, add = TRUE, show.legend = FALSE, col = "blue") legend("topleft", legend = c(2100, 2015), col = c("red", "blue"), lty = 1) # combining non-agregated regions with aggregated ones dev.off() pop.trajectories.plot(pop.pred, expression = "PIND / P900", main = "Proportion of population of India to the world population")