Climate Change and Population Modeling in R

October 8, 2017

(This article was first published on R Views, and kindly contributed to R-bloggers)

A recent paper in nature climate change: Less than 2°C warming by 2100 unlikely (Raftery et al. 2017), concludes that the goal of the Paris Agreement is unlikely to be met.

Although the conclusion is disheartening, the paper advances the science of climate modeling by developing a joint Bayesian hierarchical model for Gross Domestic Product per capita and carbon intensity. This ensemble of models, in turn, depends on the availability of probabilistic population projections developed by the BayesPop Project at the University of Washington and available on CRAN.

The linkage between CO2 emissions and population modeling comes through the use of a version of the Kaya identity, which “expresses future emissions levels in a country as a product of three components: population, GDP per capita, and carbon intensity (CO2 emissions per unit of GDP).”

The 3700 lines of R code that implement the models underlying the paper are available here on GitHub. Unless you work in this aspect of climate modeling, mastering this code seems likely to be a formidable task. A good first step, however, might be to work through the population models developed by BayesPop. These models comprise a series of R packages:

  • bayesTFR implements a Bayesian hierarchical model to make projections of total fertility for all of the countries in the world.
  • bayesLife implements a Bayesian Hierarchical model to make life expectancy projections for all the countries in the world is described here.
  • bayesPop uses bayesTFR and bayesLife to generate population projections for all of the world’s countries.
  • bayesDem provides graphical user interface for bayesTFR, bayesLife and bayesPop.
  • wppExployer implements a shiny app for exploring the data in on the World Population Prospects This data is available in R packages wpp2017, wpp2015, wpp2012, and wpp2010,

The wppExployer shiny app is also available online.

To give you a feel for how these models work, we use the bayesTFR package to estimate and plot fertility rates for the United States. (Note that, for brevity, the informative console output of the simulations in the code below have been suppressed.)

# This command produces output data, such as in the directory ex-data
sim.dir <- tempfile()
# Phase II MCMCs
m <- run.tfr.mcmc(nr.chains=1, iter=60, output.dir=sim.dir, seed=1, verbose=TRUE)
## Loading required package: wpp2017
## Warning: package 'wpp2017' was built under R version 3.4.2
# Phase III MCMCs (not included in the package)
m3 <- run.tfr3.mcmc(sim.dir=sim.dir, nr.chains=2, iter=100, thin=1, seed=1, verbose=TRUE)
# Prediction
pred <- tfr.predict(m, burnin=30, burnin3=50, verbose=TRUE)
summary(pred, country='United States of America')
unlink(sim.dir, recursive=TRUE)
tfr.trajectories.plot(pred,country="United States of America")

The transparent, reproducible work being done by Raftery and his collaborators, and the other members of the BayesPop Project, is a great gift to the R community. Anyone with an appetite for R code and the patience to learn something of Bayesian hierarchical models can gain some insight into one of the most important scientific problems of our day.

To leave a comment for the author, please follow the link and comment on their blog: R Views. offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Search R-bloggers


Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)