# 7. Simulation

# 1 Introduction

Simulation studies can be used to assess the performance of a model. The basic idea is to generate some parameter values, use these parameters to generate some data, use the data to try to infer the original parameter values, and then see how close the inferred parameter values are to the actual parameter values.

Figure 1.1: Simulation study of a model

Function `report_sim()` automates the process of doing a simulation study. We are still experimenting with `report_sim()`, and the interface may change. Suggestions are welcome, ideally through raising an issue here.

# 2 Estimation model matches data generation model

The most straightforward type of simulation is when the `estimation model' used to do the inference matches the`data generation model’ used to create the data. Even when the estimation model matches the data generation model, the inferred values for the parameters will not exactly reproduce the actual values, since data is drawn at random, and provides a noisy signal about parameters it was generated from. However, if the experiment is repeated many times, with a different randomly-drawn dataset each time, the errors should more or less average out at zero, 50% credible intervals should contain the true values close to 50% of the time, and 95% credible intervals should contain the true values close to 95% of the time.

To illustrate, we use investigate the performance of a model of divorce rates in New Zealand. We reduce the number of ages and time periods to speed up the calculations.

    library(bage)
    library(dplyr, warn.conflicts = FALSE)
    library(poputils)

    divorces_small <- nzl_divorces |>
      filter(age_upper(age) < 40,
             time >= 2016) |>
      droplevels()

    mod <- mod_pois(divorces ~ age * sex + time,
                    data = divorces_small,
                        exposure = population)
    mod
    #>
    #>     ------ Unfitted Poisson model ------
    #>
    #>    divorces ~ age * sex + time
    #>
    #>   exposure = population
    #>
    #>         term  prior along n_par n_par_free
    #>  (Intercept) NFix()     -     1          1
    #>          age   RW()   age     4          4
    #>          sex NFix()     -     2          2
    #>         time   RW()  time     6          6
    #>      age:sex   RW()   age     8          8
    #>
    #>  disp: mean = 1
    #>
    #>  n_draw var_time var_age var_sexgender
    #>    1000     time     age           sex

To do the simulation study, we pass the model to `report_sim()`. If only one model is supplied, `report_sim()` assumes that that model should be used as the estimation model and as the data generation model. By default `report_sim()` repeats the experiment 100 times, generating a different dataset each time.

    set.seed(0)
    res <- report_sim(mod_est = mod)
    res
    #> $components
    #> # A tibble: 9 × 7
    #>   term        component  .error .cover_50 .cover_95 .length_50 .length_95
    #>   <chr>       <chr>       <dbl>     <dbl>     <dbl>      <dbl>      <dbl>
    #> 1 (Intercept) effect     0.0408     0.48      0.93       1.18       3.41
    #> 2 age         effect    -0.106      0.485     0.94       1.41       4.14
    #> 3 age         hyper      0.0407     0.41      0.83       0.584      2.25
    #> 4 sex         effect    -0.126      0.475     0.945      1.13       3.30
    #> 5 time        effect    -0.0602     0.482     0.97       1.23       3.64
    #> 6 time        hyper     -0.0120     0.37      0.89       0.412      1.20
    #> 7 age:sex     effect     0.163      0.445     0.924      1.42       4.08
    #> 8 age:sex     hyper     -0.0466     0.49      0.88       0.479      1.65
    #> 9 disp        disp      -0.0689     0.52      0.96       0.244      0.733
    #>
    #> $augment
    #> # A tibble: 2 × 7
    #>   .var      .observed     .error .cover_50 .cover_95 .length_50 .length_95
    #>   <chr>         <dbl>      <dbl>     <dbl>     <dbl>      <dbl>      <dbl>
    #> 1 .fitted        115.  -0.000202     0.506     0.946    0.00857     0.0246
    #> 2 .expected      115. -49.3          0.477     0.929   93.5       352.

The output from `report_sim()` is a list of two data frames. The first data frame contains results for parameters associated with the `components()` function: main effects and interactions, associated hyper-parameters, and dispersion. The second data frame contains results for parameters associated with the `augment()` function: the lowest-level rates parameters.

As can be seen in the results, the errors do not average out at exactly zero, 50% credible intervals do not contain the true value exactly 50% of the time, and 95% credible intervals do not contain the true value exactly 95% of the time. However, increasing the number of simulations from the default value of 100 to, say, 1000 will reduce the average size of the errors closer to zero, and bring the actual coverage rates closer to their advertised values. When larger values of `n_sim` are used, it can be helpful to use parallel processing to speed up calculations, which is done through the `n_core` argument.

One slightly odd feature of the results is that the mean for `.expected` is very large. This reflects the fact that the data generation model draws some extreme values. We are developing a set of more informative priors that should avoid this behavior in future versions of bage.

# 3 Estimation model different from data generation model

In actual applications, no estimation model ever perfectly describes the true data generating process. It can therefore be helpful to see how robust a given model is to misspecification, that is, to cases where the estimation model differs from the data generation model.

With `report_sim()`, this can be done by using one model for the `mod_est` argument, and a different model for the `mod_sim` argument.

Consider, for instance, a case where the time effect is generated from an `AR1()` prior, while the estimation model continues to use the default value of a `RW()` prior,

    mod_ar1 <- mod |>
      set_prior(time ~ AR1())
    mod_ar1
    #>
    #>     ------ Unfitted Poisson model ------
    #>
    #>    divorces ~ age * sex + time
    #>
    #>   exposure = population
    #>
    #>         term  prior along n_par n_par_free
    #>  (Intercept) NFix()     -     1          1
    #>          age   RW()   age     4          4
    #>          sex NFix()     -     2          2
    #>         time  AR1()  time     6          6
    #>      age:sex   RW()   age     8          8
    #>
    #>  disp: mean = 1
    #>
    #>  n_draw var_time var_age var_sexgender
    #>    1000     time     age           sex

We set the `mod_sim` argument to `mod_ar1` and generate the report.

    res_ar1 <- report_sim(mod_est = mod, mod_sim = mod_ar1)
    res_ar1
    #> $components
    #> # A tibble: 8 × 7
    #>   term        component  .error .cover_50 .cover_95 .length_50 .length_95
    #>   <chr>       <chr>       <dbl>     <dbl>     <dbl>      <dbl>      <dbl>
    #> 1 (Intercept) effect    -0.121      0.53      0.92       1.17       3.40
    #> 2 age         effect     0.115      0.495     0.96       1.43       4.17
    #> 3 age         hyper      0.0631     0.43      0.81       0.592      2.28
    #> 4 sex         effect    -0.109      0.495     0.955      1.12       3.30
    #> 5 time        effect    -0.0705     0.605     0.96       1.20       3.58
    #> 6 age:sex     effect     0.117      0.471     0.956      1.44       4.12
    #> 7 age:sex     hyper     -0.0773     0.52      0.86       0.482      1.64
    #> 8 disp        disp      -0.0664     0.48      0.93       0.224      0.673
    #>
    #> $augment
    #> # A tibble: 2 × 7
    #>   .var      .observed    .error .cover_50 .cover_95 .length_50 .length_95
    #>   <chr>         <dbl>     <dbl>     <dbl>     <dbl>      <dbl>      <dbl>
    #> 1 .fitted        90.6 -0.000130     0.510     0.945    0.00819     0.0235
    #> 2 .expected      90.6  0.710        0.492     0.937   29.2        93.3

In this case, although actual coverage for the hyper-parameters (in the `components` part of the results) now diverges from the advertised coverage, coverage for the low-level rates (in the `augment` part of the results) is still close to advertised coverage.

# 4 The relationship between `report_sim()` and `replicate_data()`

Functions `report_sim()` and `replicate_data()` overlap, in that both use simulated data to provide insights into model performance. Their aims are, however, different. Typically, `report_sim()` is used before fitting a model, to assess its performance across a random selection of possible datasets, while `replicate_data()` is used after fitting a model, to assess its performance on the dataset to hand.