Reproducible computation at scale in R

class: center, middle, inverse, title-slide

# Reproducible computation at scale in R
### Will Landau

---

## Large statistical computation

* [Bayesian data analysis](https://mc-stan.org/)
* [Bayesian network meta-analysis](https://bookdown.org/MathiasHarrer/Doing_Meta_Analysis_in_R/bayesian-network-meta-analysis.html)
* [Graph-based multiple comparison procedures](https://github.com/kornl/gMCP)
* [Subgroup identification](https://cran.r-project.org/web/packages/TSDT/index.html)
* [Predictive modeling](http://appliedpredictivemodeling.com/computing)
* [Deep neural networks](https://keras.rstudio.com/)
* [PK/PD modeling](https://github.com/nlmixrdevelopment/nlmixr)
* Clinical trial simulation
* Target identification

???

Thank you all for coming, and thank you to nyhackr and Lander Analytics for the opportunity to speak today.

I come from the life sciences, and we develop ambitious computational workflows for Statistics and data science. There's a lot of Bayesian analysis, machine learning, simulation, and prediction. Other domains have similar workloads, and we need to think about both efficiency and reproducibility right from the start.

---

## Common features

1. Heavy use of the [R language](https://www.r-project.org/).
2. Long runtimes.
3. Multiple sub-tasks.
4. Frequent changes to code and data.

<img src = "./images/sisyphus.svg" align="left" style="border: none; box-shadow: none; height: 325px; text-align: center;">
<br>
<a href="https://openclipart.org/detail/275842/sisyphus-overcoming-silhouette">https://openclipart.org/detail/275842/sisyphus-overcoming-silhouette</a>

???

Many of these projects require long runtimes. Methods like Markov chain Monte Carlo and deep neural nets are computationally expensive. It could take hours or even days just to fit a single model. That's fine if you're only going to run the project once, or at regularly scheduled times. But if the code is still under development, it's easy to get trapped in a vicious Sisyphean cycle.

---

## Interconnected tasks
<center>
<img src = "./images/workflow.png" align="middle" style="border: none; box-shadow: none; text-align: center;">
</center>

???

A large workflow has a large number of moving parts. We have datasets that we preprocess or simulate, analyses of those datasets, and summaries of the analyses.

---

## Changes

???

If you change any one of these parts - whether it's a bugfix, a tweak to a model, or some new data -

---

## Consequences

???

Then everything that depends on it is no longer valid, and you need to rerun the computation to bring the results back up to date. This is seriously frustrating when you're in development and you're still making a constant stream of changes to code and data in real time. If every change means you need to rerun the project, there's no way the results can keep up...

---

## Pipeline tools and workflow managers

- Tons exist already: [github.com/pditommaso/awesome-pipeline](https://github.com/pditommaso/awesome-pipeline).
- Most are language-agnostic or designed for Python or the shell.

???

...unless you use a pipeline tool. There are pipeline tools for production which resemble Apache Airflow, and there are pipeline tools for development which resemble GNU Make. Today, I'm going to focus on Make-like tools because those are the ones I think are designed for this part of the process. It's an action-packed space, and there are a lot of great options. But unfortunately, there's not a whole lot for R.

---

## What distinguishes `targets`?

* Fundamentally designed for R.
* Supports a clean, modular, function-oriented programming style.
* Abstracts files as R objects and automatically manages data.
* Surpasses the permanent limitations of its predecessor, [`drake`](https://github.com/ropensci/drake): <https://wlandau.github.io/targets/articles/need.html>

???

That's where targets comes in. targets is a Make-like pipeline tool that is fundamentally designed for R. You can call it from an R session, it supports a clean, idiomatic, function-oriented style of programming, and it helps you store and retrieve your results. Most importantly, it gets you out of the Sisyphean loop of long computation, enhances reproducibility, and takes the frustration out of data science.

---

## Example targets workflow: validate a Bayesian model

$$
`\begin{aligned}
y_i &\stackrel{\text{iid}}{\sim} \text{Normal}(\alpha + \beta \cdot x_i, \ \sigma^2) \\
\alpha &\sim \text{Normal}(0, 1) \\
\beta &\sim \text{Normal}(0, 1) \\
\sigma &\sim \text{HalfCauchy}(0, 1) \\
\end{aligned}`
$$

1. Simulate 1000 datasets from the model.
1. Fit the model to each simulated dataset.
1. Verify that in 90% of the simulation reps, the `$\beta$` parameter values used for data generation are covered in 90% credible intervals. (Should repeat for other parameters too.)

---

## Stan model

```r
data {
  int <lower = 1> n;
  vector[n] x;
  vector[n] y;
}
parameters {
  real alpha;
  real beta;
  real <lower = 0> sigma;
}
model {
  y ~ normal(alpha + x * beta, sigma);
  alpha ~ normal(0, 1);
  beta ~ normal(0, 1);
  sigma ~ cauchy(0, 1);
}
```

---

## Function for data generation

```r
simulate_data_continuous <- function() {
  alpha <- rnorm(1, 0, 1)
  beta <- rnorm(1, 0, 1)
  sigma <- rhcauchy(1, 1)
  x <- rnorm(100, 1, 1) # continuous covariate
  y <- rnorm(100, alpha + x * beta, sigma)
  sim <- basename(tempfile(pattern = "sim"))
  tibble(x = x, y = y, beta_true = beta, sim = sim)
}
```

---

## Function to fit the model

```r
fit_model <- function(data, model_file) {
  stan_data <- list(x = data$x, y = data$y, n = nrow(data))
  truth <- data$beta_true[1]
  model <- cmdstan_model(model_file)
  fit <- model$sample(data = stan_data, refresh = 0)
  fit$summary() %>%
    filter(variable == "beta") %>%
    mutate(beta_true = truth, cover_beta = q5 < truth & truth < q95)
}
```

---

## Utility functions

```r
compile_model <- function(model_file) {
  quiet(cmdstan_model(model_file))
  model_file
}

map_sims <- function(data, model_file) {
  data %>%
    group_by(sim) %>%
    group_modify(~quiet(fit_model(.x, model_file = model_file))) %>%
    ungroup()
}

quiet <- function(code) {
  sink(nullfile())
  on.exit(sink())
  suppressMessages(code)
}
```

---

## File structure

```r
├── run.sh
├── run.R
*├── _targets.R
├── _targets/
├── sge.tmpl
├── R
*│   ├── functions.R
*│   └── utils.R
├── stan
│   └── model.stan
└── report.Rmd
```

---

## Setup in _targets.R

```r
# Beginning of _targets.R
library(targets)
library(tarchetypes)
source("R/functions.R")
source("R/utils.R")
options(tidyverse.quiet = TRUE)
options(clustermq.scheduler = "sge", clustermq.template = "sge.tmpl")
tar_option_set(
  packages = c("cmdstanr", "extraDistr", "rmarkdown", "tidyverse")
)
```

---

## Pipeline in _targets.R

```r
# End of _targets.R
tar_pipeline(
  tar_target(
    model_file,
    quiet(compile_model("stan/model.stan")),
    format = "file",
    deployment = "main"
  ),
  tar_target(index_batch, seq_len(25)),
  tar_target(index_sim, seq_len(40)),
  tar_target(
    data_continuous,
    map_dfr(index_sim, ~simulate_data_continuous()),
    pattern = map(index_batch),
    format = "fst_tbl"
  )
)
```

---

## Inspect the graph

```r
tar_visnetwork()
```

---

## Run the targets

```r
tar_make()
#> ● run target batch_index
#> ● run target model_file
#> ● run branch data_continuous_f385e0ef
#> ● run branch data_continuous_0217a462
#> ● run branch data_continuous_a06bc09e
#> ● run branch data_continuous_c1e661c3
#> ...
```

---

## Now all targets are up to date.

```r
tar_make()
#> ✓ skip target batch_index
#> ✓ skip target model_file
#> ✓ skip branch data_continuous_f385e0ef
#> ✓ skip branch data_continuous_0217a462
#> ✓ skip branch data_continuous_a06bc09e
#> ✓ skip branch data_continuous_c1e661c3
#> ...
#> ✓ Already up to date.
```

---

## Now all targets are up to date.

```r
tar_visnetwork()
```

---

## Inspect the data

```r
tar_read(data_continuous, branches = 1) # see also tar_load()
#> # A tibble: 200 x 4
#>         x       y beta_true rep            
#>     <dbl>   <dbl>     <dbl> <chr>          
#>  1  0.446  1.37      -0.870 rep3841223f956f
#>  2  2.01  -0.361     -0.870 rep3841223f956f
#>  3  1.05  -0.183     -0.870 rep3841223f956f
#>  4  0.734 -0.891     -0.870 rep3841223f956f
#>  5  0.567  1.18      -0.870 rep3841223f956f
#>  6  1.30   0.306     -0.870 rep3841223f956f
#>  7  0.649  0.0242    -0.870 rep3841223f956f
#>  8  0.678  0.793     -0.870 rep3841223f956f
#>  9  1.05   0.688     -0.870 rep3841223f956f
#> 10 -0.437  2.25      -0.870 rep3841223f956f
#> # … with 190 more rows
```

---

## Add the models to the pipeline

```r
# _targets.R
# setup...
tar_pipeline(
  # existing targets...
  tar_target(
    fit_continuous,
    map_sims(data_continuous, model_file = model_file),
    pattern = map(data_continuous),
    format = "fst_tbl"
  )
)
```

---

## Check the dependency graph

```r
tar_visnetwork()
```

---

## Run the models on a cluster

---

## Inspect the model results

```r
tar_read(fit_continuous)
#> # A tibble: 8 x 13
#>   rep   variable    mean  median      sd     mad     q5     q95
#>   <chr> <chr>      <dbl>   <dbl>   <dbl>   <dbl>  <dbl>   <dbl>
#> 1 rep3… beta     -0.765  -0.765  0.0868  0.0856  -0.905 -0.619 
#> 2 rep3… beta      1.03    1.03   0.0840  0.0852   0.892  1.17  
#> 3 rep3… beta     -0.560  -0.560  0.0309  0.0306  -0.612 -0.511 
#> 4 rep3… beta      0.859   0.859  0.0302  0.0306   0.810  0.908 
#> 5 rep3… beta     -0.0980 -0.0962 0.111   0.109   -0.280  0.0792
#> 6 rep3… beta      0.462   0.463  0.0567  0.0565   0.369  0.555 
#> 7 rep3… beta      0.316   0.318  0.0848  0.0840   0.178  0.454 
#> 8 rep3… beta      0.499   0.499  0.00727 0.00709  0.486  0.511 
#> # … with 5 more variables: rhat <dbl>, ess_bulk <dbl>,
#> #   ess_tail <dbl>, beta_true <dbl>, cover_beta <lgl>
```

---

## Inspect the graph

```r
tar_visnetwork(label = c("time", "size", "branches"))
```

---

## Add discrete data

```r
# functions.R
simulate_data_discrete <- function() {
  alpha <- rnorm(1, 0, 1)
  beta <- rnorm(1, 0, 1)
  sigma <- rhcauchy(1, 1)
  x <- rbinom(100, 1, 0.5) # discrete covariate
  y <- rnorm(100, alpha + x * beta, sigma)
  sim <- basename(tempfile(pattern = "sim"))
  tibble(x = x, y = y, beta_true = beta, sim = sim)
}
```

---

## Add an R Markdown report

---

## Add new targets

```r
# _targets.R
# setup...
tar_pipeline(
  # existing targets...
  tar_target(
    data_discrete,
    map_dfr(index_sim, ~simulate_data_discrete()),
    pattern = map(index_batch),
    format = "fst_tbl"
  ),
  tar_target(
    fit_discrete,
    map_sims(data_discrete, model_file = model_file),
    pattern = map(data_discrete),
    format = "fst_tbl"
  ),
  tar_render(report, "report.Rmd") # from tarchetypes
)
```

---

## New targets

```r
tar_visnetwork(label = c("time", "size", "branches"))
```

---

## Only the new targets run

```r
tar_make_clustermq(workers = 25)
#> ✓ skip target batch_index
#> ✓ skip target model_file
#> ✓ skip branch data_continuous_f385e0ef
#> ✓ skip branch data_continuous_0217a462
#> ...
#> ✓ skip branch fit_continuous_5f87e2ff
#> ✓ skip branch fit_continuous_d5dfa376
#> ...
#> ● run branch data_discrete_06241955
#> ● run branch data_discrete_88664909
#> ...
#> ● run branch fit_discrete_e31322f6
#> ● run branch fit_discrete_a2fab7f4
#> ...
#> ● run target report
```

---

## Change a function

```r
# functions.R
simulate_data_discrete <- function(...) {
  alpha <- rnorm(1, 0, 1)
  beta <- rnorm(1, 0, 1)
  sigma <- rhcauchy(1, 1)
* x <- rbinom(100, 1, 0.3) # previously 0.5
  y <- rnorm(100, alpha + x * beta, sigma)
  sim <- basename(tempfile(pattern = "sim"))
  tibble(x = x, y = y, beta_true = beta, sim = sim)
}
```

---

## Downstream targets are automatically invalidated

```r
tar_visnetwork()
```

---

## Only the affected targets rerun

---

## Inspect the report

---

## Inspect the report

---

## Evidence of reproducibility

```r
tar_visnetwork()
```

---

## Resources

* Get [`targets`](https://github.com/wlandau/targets):

```r
install.packages("remotes")
remotes::install_github("wlandau/targets")
```

* Code: <https://github.com/wlandau/targets-stan>
* These slides: <https://wlandau.github.io/nyhackr2020>
* Tutorial materials: <https://github.com/wlandau/targets-tutorial>
* Development repository: <https://github.com/wlandau/targets>
* Full user manual: <https://wlandau.github.io/targets-manual/>
* Reference website: <https://wlandau.github.io/targets/>

???

There are several resources to learn about targets. There's a reference website, an online user manual, and a repository with the example code from today.

---

## Live coding: continuously deploy a Shiny app

```r
tar_pipeline(
  ...,
  tar_target(
    results_file,
    export_results(continuous = fit_continuous, discrete = fit_discrete),
    format = "file"
  ),
  tar_target(app_source, "app.R", format = "file"),
  tar_target(deploy, deploy_app(app_source, results_file), deployment = "main")
)
```

---

## What about `drake`?

* `drake` is still an excellent choice for pipeline management, but it has permanent user-side limitations.
* `targets` was created to overcome these limitations and create a smoother user experience.
    1. Stronger guardrails by design.
    1. A friendlier, lighter, more transparent data management system.
    1. Show which *functions* are up to date.
    1. More flexible dynamic branching.
    1. Improved parallel efficiency.
    1. Designed for custom user-side [metaprogramming](https://wlandau.github.io/targets-manual/branching.html#metaprogramming) and target archetypes: <https://wlandau.github.io/tarchetypes/>.
* The statement of need describes the details: <https://wlandau.github.io/targets/articles/need.html>.

???

But what about drake? drake already does all these things, and it's still an excellent choice for Make-like pipeline management. But it does have permanent user-side limitations. We've been developing, improving, expanding, and refining drake for several years, and we've reached the point where the most important problems to tackle are exactly the problems we cannot solve in this tool. It's just too big and too set in its ways, and its architecture was originally designed around assumptions that no longer hold up. To overcome these permanent limitations, we need a new tool that borrows from drake's journey and advances the user experience beyond what drake is capable of, and that new tool is targets. targets has stronger guardrails, lighter data management, greater transparency around data and the process of watching for changes, more flexible dynamic branching, better parallel efficiency, and design that lets us build on top of it more easily. The targets package website has a statement of need describing the changes in more detail.

---

## Guardrails in `targets`

* The only way to use `targets` is the correct way.
* Main guardrails:
    1. Always run in a fresh R process (unless you deliberately configure `targets` for debugging).
    2. Require a `_targets.R` configuration file in the project root.
    3. Require the `_targets/` data store to always be in the project root.

???

Let's start with guardrails. The only way to use targets is the correct way. Unless you deliberately opt out for debugging purposes, targets always does its work in a fresh clean reproducible R process to avoid unpredictably invalidating targets. It's also paternalistic about your working directory and data storage. It removes flexibility where it doesn't belong, and it keeps users from getting themselves into trouble.

---

## `drake`'s cache

```
.drake/
├── config/
├── data/
├───── 17bfcef645301416.rds
├───── 21935c86f12692e2.rds
├───── 37caf5df2892cfc4.rds
├───── ...
├── drake/
├── keys/
├───── memoize/
├───── meta/
├───── objects/
├───── ...
└── scratch/
```

???

Now for data management. drake's cache is a large file system. It has hundreds of tiny bookkeeping files, it gets heavier and heavier over time unless you deliberately run garbage collection, and it's not human-readable. If you commit this thing to Git, or try to share with your colleagues, or try to diagnose cryptic errors that come from corrupted data, you're going to have a hard time.

---

## The data store in `targets`

```
_targets/
├── meta/
├───── meta
├───── progress
├── objects/
├───── target_name_1
├───── target_name_2
├───── target_name_3
└───── ...
```

???

targets simplifies data storage. It only stores what it needs to. There's one file per target in the objects/ folder, one data.table with metadata, and one data.table with runtime progress information. The data store is much lighter and much more portable, and the file names all make sense. And it survives corruptions far more easily. Of all these files, only the metadata file is sacred. Anything else can break and the project can recover just by rerunning one or two targets.

---

## Seamless AWS S3 integration

<https://wlandau.github.io/targets-manual/cloud.html>

```r
# _targets.R
tar_option_set(resources = list(bucket = "my-bucket-name"))
tar_pipeline(
  tar_target(dataset, get_large_dataset(), format = "aws_fst_tbl"),
  tar_target(analysis, analyze_dataset(dataset), format = "aws_qs")
)
```

```r
# R session
tar_make()
tar_read(dataset)
```

???

The simpler data store paved the way for seamless, Metaflow-like integration with Amazon S3. There's a little config, but after that, `targets` will automatically upload the return values to an S3 bucket. And retrieving cloud data feels exactly the same as reading local data.

---

## Show which functions are out of date

![](./images/graph_imports.png)

???

drake tells you which targets are out of date, but it cannot point to the precise functions and global objects that changed since the last run of the pipeline. targets provides this information, which significantly decreases frustration and increases reproducibility. In this example dependency graph, we can see that the create_plot() function changed, which is one of the reasons targets hist and report are no longer up to date.

---

## Dynamic branching with `dplyr`

```r
library(dplyr)
library(targets)
data.frame(
  x = seq_len(6),
  id = rep(letters[seq_len(3)], each = 2)
) %>%
  group_by(id) %>%
* tar_group()
#> # A tibble: 6 x 3
#> # Groups:   id [3]
#>       x id    tar_group
#>   <int> <chr>     <int>
#> 1     1 a             1
#> 2     2 a             1
#> 3     3 b             2
#> 4     4 b             2
#> 5     5 c             3
#> 6     6 c             3
```

???

Now for dynamic branching. drake began as a static tool, whereas targets is fundamentally dynamic. In fact, targets finally achieves one of the most frequently requested features for drake: dynamic branching over arbitrary subsets of data frames. First, we define a grouped data frame with `dplyr::group_by()`, `targets::tar_group()`, and `iteration = "group"`.

---

## Define a target with groups.

```r
tar_target(
  data,
  data.frame(
    x = seq_len(6),
    id = rep(letters[seq_len(3)], each = 2)
  ) %>%
    group_by(id) %>%
*   tar_group(),
* iteration = "group"
)
```

???

Then, every downstream target that branches over it will automatically create one branch for each group of rows.

---

## Inefficient dynamic branching in drake

![](./images/dynamic_drake.png)

???

The efficiency of dynamic branching is also better. Dynamic sub-targets in drake are strictly nested inside their parents, so all the sub-targets in a group must finish before any downstream sub-target can begin. But targets knows how to advance forward even if not all the upstream branches are finished. This behavior saves runtime and conserves computing resources and gets work done faster..

---

## Efficient dynamic branching in targets

![](./images/dynamic_targets.png)

???

---

## Composable dynamic branching patterns

```r
# _targets.R
tar_pipeline(
  tar_target(w, seq_len(2)),
  tar_target(x, head(letters, 3)),
  tar_target(y, head(LETTERS, 3)),
  tar_target(
    z,
    data.frame(w = w, x = x, y = y),
*   pattern = cross(w, map(x, y))
  )
)
```

---

## Composable dynamic branching patterns

* `pattern = cross(w, map(x, y))` is equivalent to `crossing(w, nesting(x, y))` from `tidyr`.

```r
tar_make(reporter = "silent")

tar_read(z)
#>   w x y
#> 1 1 a A
#> 2 1 b B
#> 3 1 c C
#> 4 2 a A
#> 5 2 b B
#> 6 2 c C
```

* Idea from [@glep](https://github.com/glep) and [@djbirke](https://github.com/djbirke) on GitHub.

---

## New pattern types (easy to add more)

* `map()`: iterate over one or more targets in sequence.
* `cross()`: iterate over combinations of slices of targets.
* `head()`: restrict branching to the first few elements.
* `tail()`: restrict branching to the last few elements.
* `sample()`: restrict branching to a random subset of elements.

---

## Dynamic branching emulation

```r
tar_pattern(cross(sample(x, 2), map(y, z)), x = 25, y = 3, z = 3)
#> # A tibble: 6 x 3
#>   x     y     z    
#>   <chr> <chr> <chr>
#> 1 x_14  y_1   z_1  
#> 2 x_14  y_2   z_2  
#> 3 x_14  y_3   z_3  
#> 4 x_4   y_1   z_1  
#> 5 x_4   y_2   z_2  
#> 6 x_4   y_3   z_3
```

---

## Metaprogramming

* `tar_target_raw()` avoids non-standard evaluation and supports third-party metaprogramming.
* The following are equivalent ways to define a target.

```r
# For most users:
tar_target(data, simulate_data(), pattern = map(index))

# For developers who metaprogram reusable pipeline archetypes:
tar_target_raw(
  "data",
  quote(simulate_data()),
  pattern = quote(map(index))
)
```

???

Lastly, targets is easier to extend and build on. Unlike drake, targets can avoid non-standard evaluation and domain-specific languages. The tar_target_raw() function lets you define targets programmatically and opens up a lot of flexibility. The online manual walks you through how to to use tar_target_raw() for custom static branching.

---

## Target archetypes

* The `tarchetypes` package has helpers for commonly used targets: <https://wlandau.github.io/tarchetypes/>

Function | Target archetype
---|---
`tar_render()` | Render a dependency-aware R Markdown report.
`tar_knit()` | Run a dependency-aware `knitr` report.
`tar_change()` | Always run a target when a custom object changes.
`tar_force()` | Always run a target when a custom condition is true.
`tar_skip()` | Never run a target when a custom condition is true.
`tar_plan()` | Simplified `drake`-like syntax for `targets` pipelines.

???

In addition, there's an external package called `tarchetypes` that uses these metaprogramming capabilities to write shorthand for commonly used targets and pipelines. These archetypes abstract away configuration details and make it easier to write concise readable pipelines.