This vignette is adapted from the homepage of the SimEngine website.

```
library(SimEngine)
#> Loading required package: magrittr
#> Welcome to SimEngine! Full package documentation can be found at:
#> https://avi-kenny.github.io/SimEngine
```

**SimEngine** is an open-source R package for structuring, maintaining, running, and debugging statistical simulations on both local and cluster-based computing environments. Emphasis is placed on thorough documentation, and scalability.

The goal of many statistical simulations is to test how a new statistical method performs against existing methods. Most statistical simulations include three basic phases: (1) generate some data, (2) run one or more methods using the generated data, and (3) compare the performance of the methods.

To briefly illustrate how these phases are implemented using **SimEngine**, we will use the example of estimating the average treatment effect of a drug in the context of a randomized controlled trial (RCT).

The simulation object (an R object of class *sim_obj*) will contain all data, functions, and results related to your simulation.

```
library(SimEngine)
<- new_sim() sim
```

Most simulations will involve one or more functions that create a dataset designed to mimic some real-world data structure. Here, we write a function that simulates data from an RCT in which we compare a continuous outcome (e.g. blood pressure) between a treatment group and a control group. We generate the data by looping through a set of patients, assigning them randomly to one of the two groups, and generating their outcome according to a simple model.

```
# Code up the dataset-generating function
<- function (num_patients) {
create_rct_data <- data.frame(
df "patient_id" = integer(),
"group" = character(),
"outcome" = double(),
stringsAsFactors = FALSE
)for (i in 1:num_patients) {
<- ifelse(sample(c(0,1), size=1)==1, "treatment", "control")
group <- ifelse(group=="treatment", -7, 0)
treatment_effect <- rnorm(n=1, mean=130, sd=5) + treatment_effect
outcome <- list(i, group, outcome)
df[i,]
}return (df)
}
# Test the function
create_rct_data(5)
#> patient_id group outcome
#> 1 1 treatment 126.4852
#> 2 2 treatment 125.2483
#> 3 3 treatment 120.4302
#> 4 4 control 138.7341
#> 5 5 treatment 128.4791
```

With **SimEngine**, any functions that you declare (or load via `source()`

) are automatically added to your simulation object when the simulation runs. In this example, we test two different estimators of the average treatment effect. For simplicity, we code this as a single function and use the `type`

argument to specify which estimator we want to use, but you could also write two separate functions. The first estimator uses the known probability of being assigned to the treatment group (0.5), whereas the second estimator uses an estimate of this probability based on the observed data. Don’t worry too much about the mathematical details; the important thing is that both methods attempt to take in the dataset generated by the `create_rct_data()`

function and return an estimate of the treatment effect, which in this case is *-7*.

```
# Code up the estimators
<- function(df, type) {
est_tx_effect <- nrow(df)
n <- sum(df$outcome * (df$group=="treatment"))
sum_t <- sum(df$outcome * (df$group=="control"))
sum_c if (type=="est1") {
<- 0.5
true_prob return ( sum_t/(n*true_prob) - sum_c/(n*(1-true_prob)) )
else if (type=="est2") {
} <- sum(df$group=="treatment") / n
est_prob return ( sum_t/(n*est_prob) - sum_c/(n*(1-est_prob)) )
}
}
# Test out the estimators
<- create_rct_data(10000)
df est_tx_effect(df, "est1")
#> [1] -9.259742
est_tx_effect(df, "est2")
#> [1] -6.930036
```

Often, we want to run the same simulation multiple times (with each run referred to as a “simulation replicate”), but with certain things changed. In this example, perhaps we want to vary the number of patients and the method used to estimate the average treatment effect. We refer to the things that vary as “simulation levels”. By default, **SimEngine** will run our simulation 10 times for each level combination. Below, since there are two methods and three values of num_patients, we have six level combinations and so **SimEngine** will run a total of 60 simulation replicates. Note that we make extensive use of the pipe operators (`%>%`

and `%<>%`

) from the **magrittr** package; if you have never used pipes, check out the magrittr documentation.

```
%<>% set_levels(
sim estimator = c("est1", "est2"),
num_patients = c(50, 200, 1000)
)
```

The simulation script is a function that runs a single simulation replicate and returns the results. Within a script, you can reference the current simulation level values using the variable *L*. For example, when the first simulation replicate is running, `L$estimator`

will equal “est1” and `L$num_patients`

will equal 50. In the last simulation replicate, `L$estimator`

will equal “est2” and `L$num_patients`

will equal 1,000. Your script will automatically have access to any functions that you created earlier.

```
%<>% set_script(function() {
sim <- create_rct_data(L$num_patients)
df <- est_tx_effect(df, L$estimator)
est return (list("est"=est))
})
```

Your script should always return a named list, although your list can be complex and contain dataframes, multiple levels of nesting, etc. Note that you can also code your estimators as separate functions and call them from within the script using `use_method`

.

This controls options related to your entire simulation, such as the number of simulation replicates to run for each level combination and how to parallelize your code. This is also where you should specify any packages your simulation needs (instead of using `library()`

or `require()`

). This is discussed in detail on the `set_config`

page. We set `num_sim`

to 10, and so **SimEngine** will run a total of 60 simulation replicates (10 for each level combination).

```
%<>% set_config(
sim num_sim = 10,
n_cores = 2,
parallel = "outer",
packages = c("ggplot2", "stringr")
)#>
#> Attaching package: 'ggplot2'
#> The following object is masked from 'package:SimEngine':
#>
#> vars
```

All 600 replicates are run at once and results are stored in the simulation object.

```
%<>% run()
sim #> Done. No errors or warnings detected.
```

Once the simulations have finished, use the `summarize()`

function to calculate common summary statistics, such as bias, variance, MSE, and coverage.

```
%>% summarize(
sim bias = list(truth=-7, estimate="est"),
mse = list(truth=-7, estimate="est")
)#> level_id estimator num_patients bias_est MSE_est
#> 1 1 est1 50 -16.98962339 1388.1446784
#> 2 2 est2 50 -0.57840518 4.0898524
#> 3 3 est1 200 6.22177351 653.7529757
#> 4 4 est2 200 -0.43567402 0.4212225
#> 5 5 est1 1000 -0.26793128 44.1567650
#> 6 6 est2 1000 -0.06323944 0.1279339
```

In this example, we see that the MSE of estimator 1 is much higher than that of estimator 2 and that MSE decreases with increasing sample size for both estimators, as expected.

You can also directly access the results for individual simulation replicates.

```
head(sim$results)
#> sim_uid level_id sim_id estimator num_patients runtime est
#> 1 1 1 1 est1 50 0.010332823 -47.299229
#> 2 2 1 2 est1 50 0.006906986 -35.046499
#> 3 3 1 3 est1 50 0.004701138 -68.928312
#> 4 4 1 4 est1 50 0.013962984 -35.526029
#> 5 5 1 5 est1 50 0.004229069 -16.109661
#> 6 6 1 6 est1 50 0.003988981 2.783245
```