[help] Workflow for parameter recovery simulation studies?

[venpopov]

Help

• I understand and agree to https://books.ropensci.org/targets/help.html.

Description

My use case is similar to the example with sensitivity analysis for bayesian models. Generally, I would have some data generating function that takes parameters and a model fitting function that takes data and reports parameter estimates (doesn't matter bayesian or not, both workflows). In parameter recovery studies, I would typically:

• generate a grid of parameters (could be hundreds/thousands of parameter combinations)
• generate a dataset for each parameter set
• estimate the model for each dataset,
• repeat this N times for each parameter set

I am currently trying to implement this with tar_map_rep like so (toy example):

library(targets)
library(tarchetypes)
tar_option_set(
 packages = c("bmm", "dplyr", "purrr")
)

# toy example
pars <- expand.grid(
 n = c(20, 50, 100),
 kappa = seq(1, 16, 5),
 p_mem = seq(0.1, 1, 0.5)
)

transform_for_mixtur <- function(error) {
 data.frame(response = error, target = 0, id = 1)
}

recover_mixture2p <- function(n, kappa, p_mem) {
 suppressMessages(
   bmm::rmixture2p(n, kappa, p_mem) |>
     transform_for_mixtur() |>
     mixtur::fit_mixtur(model = "2_component", unit = "radians") |>
     select(kappa, p_t) |>
     rename(kappa_est = kappa, p_mem_est = p_t) |>
     mutate(n = n, kappa = kappa, p_mem = p_mem)
 )
}

list(
 tar_map_rep(
   name = results,
   command = recover_mixture2p(n, kappa, p_mem),
   values = pars,
   batches = 2,
   reps = 3
 )
)

This seems to work well, but even with this toy example, the DAG is unusable:

The realistic example would have:

• a grid of 10000+ parameter combinations
• 1000 replicates (e.g. batches = 20, reps = 50)
• be run in parallel

I have two questions:

• Is the way I am doing it the best way for my stated purposes?
• If so, is there a way to make the DAG more useful? I read through the help pages and I don't see a way to make it less cluttered in a case like this with tar_map_rep

[wlandau]

Some thoughts:

• For so many parameter combinations and replications, I would recommend using dynamic branching with batching to keep the total number of targets down in the low tens of thousands.
• On top of that, you might use static branching for different kinds of models or different simulation scenarios.
• Since there are so many replications, I recommend being careful about the size of the output data for each iteration. Rather than returning all the MCMC samples, I recommend returning compact summary statistics, convergence diagnostics, specific posterior probabilities of interest, etc.
• With each replication, I also recommend returning metadata about the kind of scenario or model or prior etc. that the iteration ran.
• targets >= 1.9.0 is much more efficient than version 1.8.0 for these kinds of large pipelines. Unfortunately version 1.9.0 is taking time to reach CRAN, but for now you can install it with remotes::install_github("ropensci/targets").

[venpopov]

With the pure dynamic branching approach, is this how you would do it?

list(
  tar_target(index_batch, seq_len(40)),
  tar_target(
    name = pars, 
    command = expand.grid(
       n = c(20, 50, 100),
       kappa = seq(1, 16, 5),
       p_mem = seq(0.1, 1, 0.5)
    )
  ),
  tar_target(
    data_continuous,
    purrr::pmap(pars, recover_mixture2p),
    pattern = map(index_batch)
  )
)

I was trying to model it after the example in the batching page. But if I understand this correctly, it would only have a branch for each batch, but within that do all the parameters in a single branch? And how would you include replications, in addition to the batches? I'm having a bit of trouble figuring out how the different pieces fit together.

Everything else we're on the same page, I already do it exactly like that ;)

[wlandau]

You might consider gridding out all your replications up front and then batching them with tarchetypes::tar_group_count(). Here's an sketch that groups all the rows of the pars grid into 100 batches.

list(
  tarchetypes::tar_group_count(
    name = pars, 
    count = 100,
    command = expand.grid(
       n = c(20, 50, 100),
       kappa = seq(1, 16, 5),
       p_mem = seq(0.1, 1, 0.5),
       replication = seq_len(100)
    )
  ),
  tar_target(
    data_continuous,
    purrr::pmap(pars, recover_mixture2p),
    pattern = map(pars)
  )
)

Alternatively, you could assign an explicit batch index column with tar_group_by(). In this case, each element of batch_index (of whatever you call it and pass to tar_group_by()) will mark each batch for anything that branches over it downstream.

list(
  tarchetypes::tar_group_by(
    name = pars, 
    batch_index,
    command = expand.grid(
       n = c(20, 50, 100),
       kappa = seq(1, 16, 5),
       p_mem = seq(0.1, 1, 0.5),
       batch_index = seq_len(100)
    )
  ),
  tar_target(
    data_continuous,
    purrr::pmap(pars, recover_mixture2p),
    pattern = map(pars)
  )
)

[venpopov]

Thank you for the suggestions! These would be useful to start experimenting with to find a workflow that fits this type of work most fruitfully.