This function is a fairly flexible way of running simulation studies in trialr, and beyond. It essentially uses delegates to perform this pattern:

for i in 1:N:
  data = get_data_func()
  fit = fit_model_func(data)
  if summarise_func is null:
    sims[i] = fit
  else
    sims[i] = summarise_func(data, fit)
  end
loop
return sims
trialr_simulate(
  N,
  get_data_func,
  fit_model_func,
  summarise_func = NULL,
  num_logs = 10,
  num_saves = NULL,
  save_func = NULL
)

Arguments

N

integer, number of simulated iterations to run.

get_data_func

Function that takes no parameters and returns a sampled dataset to be analysed. I.e. the call signature is f().

fit_model_func

Function that accepts the output of get_data_func as the sole parameter and fits the model or performs the analysis, returning an object of arbitrary type.

summarise_func

Optional. If provided, this function should accept the ouputs of get_data_func and fit_model_func as parameters 1 & 2 and perform some post-fit processing or simplification. The result of this call is the output from iteration i. If omitted, the fit object from fit_model_func is simply used as the output from iteration i.

num_logs

Number of log messages to receive about progress. NULL to suppress logging. E.g. if N=100 and num_logs=10, you will get log messages when i=10, 20, 30, etc.

num_saves

Number of interimittent saves to attempt. NULL to suppress saving E.g. if N=100 and num_saves=10, the save_func delegate will be called after iteration i=10, 20, 30, etc.

save_func

Optional. Function that takes the interim list of simulated objects as the sole parameter and saves them somehow. This, combined with num_saves, allows periodic saving of in-progress results to avoid complete data loss if the simulation study fails for some reason.

Value

list of length N. The items in the list are as returned by summarise_func or fit_model_func.

Examples

get_data_func <- function() { group_sizes <- rbinom(n = 5, size = 50, prob = c(0.1, 0.3, 0.3, 0.2, 0.1)) group_responses <- rbinom(n = 5, size = group_sizes, prob = c(0.2, 0.5, 0.2, 0.2, 0.2)) list( group_responses = group_responses, group_sizes = group_sizes, mu_mean = gtools::logit(0.1), mu_sd = 1, tau_alpha = 2, tau_beta = 20 ) } fit_model_func <- function(data) { data <- append(data, list(refresh = 0)) do.call(stan_hierarchical_response_thall, args = data) } summarise_func <- function(data, fit) { # Probability that estimate response rate exceeds 30% unname(colMeans(as.data.frame(fit, 'prob_response') > 0.3)) } if (FALSE) { sims <- trialr_simulate(N = 20, get_data_func, fit_model_func, summarise_func) # Posterior probabilities that the response rate in each cohort exceeds 30%: do.call(rbind, sims) # Cohorts are in columns; simulated iterations are in rows. }