Goals

Data analysis has interconnected steps.

Each update…

…can invalidate other work.

Do you hunt for all the changes yourself?

• Messy and prone to human error.
• Not reproducible.

Do you rerun everything from scratch?

• Takes too long.
• Too frustrating.

Pipeline toolkits for large computation

• Lots of tools: github.com/pditommaso/awesome-pipeline.
• Few designed for R.

Small example analysis project: GDP

• Which quantity is a better predictor of GDP per capita, life expectancy or population?
  • Gapminder data: observations on multiple countries from 1952 to 2007.
  • Bayesian regeression with rstanarm: straightforward inference and interpretation.
• But Bayesian methods can be computationally expensive!

Let’s use the drake R package.

The drake plan: steps of the workflow.

library(drake)
plan <- drake_plan(
  life_exp = fit_model(covariate = "life_exp", dataset = dataset),
  log_pop = fit_model(covariate = "log_pop", dataset = dataset),
  dataset = prepare_gapminder(gapminder::gapminder),
  plots = gather_plots(life_exp, log_pop),
  report = render(
    knitr_in("report.Rmd"),
    output_file = file_out("report.html")
  )
)

The plan is just a data frame.

plan
## # A tibble: 5 x 2
##   target   command                                                         
##   <chr>    <chr>                                                           
## 1 life_exp "fit_model(covariate = \"life_exp\", dataset = dataset)"        
## 2 log_pop  "fit_model(covariate = \"log_pop\", dataset = dataset)"         
## 3 dataset  prepare_gapminder(gapminder::gapminder)                         
## 4 plots    gather_plots(life_exp, log_pop)                                 
## 5 report   "render(knitr_in(\"report.Rmd\"), output_file = file_out(\"repo…

Support for the plan

# Load your packages: library(drake), library(gapminder), etc.
source("packages.R")

# Define your functions: prepare_gapminder(), fit_model(), etc.
source("functions.R") 

Run your workflow.

make(plan)
## target dataset
## target log_pop
## target life_exp
## target plots
## target report

Output file report.html

Get targets from the cache.

readd(dataset)
## # A tibble: 1,704 x 3
##    life_exp log_gdp_percap log_pop
##       <dbl>          <dbl>   <dbl>
##  1    -2.37          -1.21   0.113
##  2    -2.26          -1.17   0.170
##  3    -2.13          -1.14   0.236
##  4    -1.97          -1.15   0.309
##  5    -1.81          -1.25   0.387
##  6    -1.63          -1.20   0.467
##  7    -1.52          -1.03   0.377
##  8    -1.44          -1.14   0.423
##  9    -1.38          -1.36   0.525
## 10    -1.37          -1.37   0.717
## # … with 1,694 more rows

Find things to improve.

loadd(plots)
plot(plots)

Go back and change a function.

gather_plots
## function (...) 
## {
##     grobs <- lapply(list(...), plot_model)
##     arrangeGrob(grobs = grobs)
## }

# Update in functions.R
# and source("functions.R") again.
gather_plots <- function(...) {
  grobs <- lapply(list(...), plot_model)
  arrangeGrob(grobs = grobs, ncol = 2)
}

Which targets need an update?

config <- drake_config(plan)
outdated(config)
## [1] "plots"  "report"

vis_drake_graph(config, targets_only = TRUE)

vis_drake_graph()

Run only the parts that need to change.

make(plan)
## target plots
## target report

Life expectancy has a stronger association.

loadd(plots)
plot(plots)

Reproducibilty is about confidence and trust.

• Tangible evidence that your results match the underlying code and data:

make(plan)
## All targets are already up to date.

outdated(config)
## character(0)

Scale up to many targets.

• Larger example: Produc dataset.
• Predict gross state product using all possible 3-covariate models.

library(Ecdat)

# We want to predict "gsp" based on the other variables.
predictors <- setdiff(colnames(Produc), "gsp")

# We will try all combinations of three covariates.
combos <- combn(predictors, 3)

# We have a grid of covariates.
x1 <- rlang::syms(combos[1, ])
x2 <- rlang::syms(combos[2, ])
x3 <- rlang::syms(combos[3, ])

Scale up to many targets.

• Experimental interface coming to drake >= 7.0.0.

plan <- drake_plan(
  model = target(
    fit_gsp_model(gsp ~ x1 + x2 + x3, data = Ecdat::Produc),
    transform = map(x1 = !!x1, x2 = !!x2, x3 = !!x3)
  ),
  rmspe_i = target(
    get_rmspe(model, Ecdat::Produc),
    transform = map(model)
  ),
  rmspe = target(
    do.call(rbind, rmspe_i),
    transform = combine()
  ),
  plot = ggsave(file_out("rmspe.pdf"), plot_rmspe(rmspe))
)

Scale up to many targets.

## # A tibble: 170 x 2
##    target              command                                             
##    <chr>               <chr>                                               
##  1 model_state_year_p… fit_gsp_model(gsp ~ state + year + pcap, data = Ecd…
##  2 model_state_year_h… fit_gsp_model(gsp ~ state + year + hwy, data = Ecda…
##  3 model_state_year_w… fit_gsp_model(gsp ~ state + year + water, data = Ec…
##  4 model_state_year_u… fit_gsp_model(gsp ~ state + year + util, data = Ecd…
##  5 model_state_year_pc fit_gsp_model(gsp ~ state + year + pc, data = Ecdat…
##  6 model_state_year_e… fit_gsp_model(gsp ~ state + year + emp, data = Ecda…
##  7 model_state_year_u… fit_gsp_model(gsp ~ state + year + unemp, data = Ec…
##  8 model_state_pcap_h… fit_gsp_model(gsp ~ state + pcap + hwy, data = Ecda…
##  9 model_state_pcap_w… fit_gsp_model(gsp ~ state + pcap + water, data = Ec…
## 10 model_state_pcap_u… fit_gsp_model(gsp ~ state + pcap + util, data = Ecd…
## # … with 160 more rows

Scale up to many targets.

Persistent parallel workers

# If running locally:
options(clustermq.scheduler = "multicore")

# If on a SLURM cluster (other clusters also supported):
drake_hpc_template_file("slurm_clustermq.tmpl") # Modify by hand.
options(
  clustermq.scheduler = "slurm",
  clustermq.template = "slurm_clustermq.tmpl"
)

# Either way, make() is the same.
make(plan, jobs = 2, parallelism = "clustermq")

Persistent parallel workers

Transient parallel workers

# If running locally:
future::plan(future::multisession)

# If on a SLURM cluster (other clusters also supported):
library(future.batchtools)
drake_hpc_template_file("slurm_batchtools.tmpl") # Modify by hand.
future::plan(batchtools_slurm, template = "slurm_batchtools.tmpl")

# Either way, make() is the same.
make(plan, jobs = 2, parallelism = "future")

Distributed computing: transient workers

Get drake and get help.

install.packages("drake") # CRAN release

library(remotes)
install_github("ropensci/drake") # Development version

library(drake)
drake_example("gapminder") # Download today's code.
drake_example("gsp")

• These slides (source)
• Quick reference guide
• Technical manual and tutorials
• File an issue.
• Contribute code.
• Discuss at rOpenSci.

Thanks

• Maëlle Salmon
• Ben Marwick
• Julia Lowndes
• Peter Slaughter
• Jenny Bryan
• Rich FitzJohn
• Stefanie Butland

• Jarad Niemi
• Kirill Müller
• Kendon Bell
• Alex Axthelm
• Jasper Clarkberg
• Henrik Bengtsson
• Eric Nantz
• Ian Watson
• Daniel Falster
• Chan-Yub Park
• Tiernan Martin

References

• drake
• rOpenSci
• Gapminder and the gapminder package
• Produc dataset dataset from the Ecdat package
• List of pipeline toolkits
• Images (desk, Sisyphus) taken from from openclipart.org under CC0 1.0, unlimited commercial use.
• The embedded parallel computing animations were created using PowToon and are used only as explicitly permitted in the Terms and Conditions of Use. The maintainers of and contributors to drake make no copyright claim to any of the constituent graphicical elements that comprise the animations.