Skip to contents

Simulation Study of a Model

Source: R/sim.R
report_sim.Rd

Use simulated data to assess the performance of an estimation model.

Usage

report_sim(
  mod_est,
  mod_sim = NULL,
  method = c("standard", "inner-outer"),
  vars_inner = NULL,
  n_sim = 100,
  point_est_fun = c("median", "mean"),
  widths = c(0.5, 0.95),
  report_type = c("short", "long", "full"),
  n_core = 1
)

Arguments

mod_est

The model whose performance is being assessed. An object of class bage_mod.

mod_sim

The model used to generate the simulated data. If no value is supplied, mod_est is used.

method

Estimation method used for mod_est. See fit().

vars_inner

Variables used in inner model with "inner-outer"estimation method. See fit().

n_sim

Number of sets of simulated data to use. Default is 100.

point_est_fun

Name of the function to use to calculate point estimates. The options are "mean" and "median". The default is "mean".

widths

Widths of credible intervals. A vector of values in the interval (0, 1]. Default is c(0.5, 0.95).

report_type

Amount of detail in return value. Options are "short" and "long". Default is "short".

n_core

Number of cores to use for parallel processing. If n_core is 1 (the default), no parallel processing is done.

Value

A named list with a tibble called "components" and a tibble called "augment".

Warning

The interface for report_sim() is still under development and may change in future.

See also

Examples

## results random, so set seed
set.seed(0)

## make data - outcome variable (deaths here)
## needs to be present, but is not used
data <- data.frame(region = c("A", "B", "C", "D", "E"),
                   population = c(100, 200, 300, 400, 500),
                   deaths = NA)

## simulation with estimation model same as
## data-generating model
mod_est <- mod_pois(deaths ~ region,
                    data = data,
                    exposure = population) |>
  set_prior(`(Intercept)` ~ Known(0))
report_sim(mod_est = mod_est,
           n_sim = 10) ## in practice should use larger value
#> $components
#> # A tibble: 4 × 7
#>   term        component  .error .cover_50 .cover_95 .length_50 .length_95
#>   <chr>       <chr>       <dbl>     <dbl>     <dbl>      <dbl>      <dbl>
#> 1 (Intercept) effect     0           1         1         0           0   
#> 2 region      effect    -0.0348      0.62      0.98      0.995       2.92
#> 3 region      hyper      0.0360      0.4       0.9       0.582       1.94
#> 4 disp        disp       0.0496      0.5       0.9       0.765       2.84
#> 
#> $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        4.41 -0.0102      0.6       0.96      0.115      0.340
#> 2 .expected      4.41 -0.545       0.62      0.98      2.92      15.4  
#> 

## simulation with estimation model different
## from data-generating model
mod_sim <- mod_est |>
  set_prior(region ~ N(s = 2))
report_sim(mod_est = mod_est,
           mod_sim = mod_sim,
           n_sim = 10)
#> $components
#> # A tibble: 4 × 7
#>   term        component  .error .cover_50 .cover_95 .length_50 .length_95
#>   <chr>       <chr>       <dbl>     <dbl>     <dbl>      <dbl>      <dbl>
#> 1 (Intercept) effect     0            1        1         0           0   
#> 2 region      effect     0.0944       0.6      0.98      1.16        3.35
#> 3 region      hyper     -0.479        0.2      0.8       0.644       2.11
#> 4 disp        disp       0.109        0.5      0.8       0.977       4.46
#> 
#> $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        157.  0.200      0.58      0.96      0.382       1.10
#> 2 .expected      157. -1.27       0.6       0.98    146.        685.  
#>