Simulation Study of a Model

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

When report_type is "short", a tibble with the following columns:

• component. Part of model. See Details. "par" is the rate, probability, or mean parameter from the likelihood.

• vals_sim. Simulated value for parameter, averaged across all simulations and cells.

• error_point_est. Point estimate minus simulation-true value, averaged across all simulations and cells.

• cover. Actual proportion of simulation-true values that fall within each type of interval, averaged across all simulations and cells.

When report_type is "long", a tibble with the following columns:

• component. Part of model. See components(). "par" is the rate, probability, or mean parameter from the likelihood.

• term. Category within component.

• level. Category within term.

• vals_sim. Simulated values for parameter, stored in an rvec.

• error_point_est. Point estimates minus simulation-true values, stored in an rvec.

• cover. Actual proportions of simulation-true values falling within each type of interval, stored in an rvec.

See also

• mod_pois(), mod_binom(), mod_norm() Specify a model

• components(), augment() Draw from joint prior or posterior distribution of model

• replicate_data() Generate replicate data for a model

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.0825       0.6       1        1.03        3.00
#> 3 region      hyper      0.0219       0.3       0.9      0.606       1.96
#> 4 disp        disp      -0.336        0.6       1        0.882       3.16
#> 
#> $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        1.63 -0.00249      0.54      0.94     0.0779      0.229
#> 2 .expected      1.63 -0.258        0.6       1        1.26        5.16 
#> 

## 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.276        0.5      0.96      1.31        3.81
#> 3 region      hyper     -0.500        0.2      0.8       0.666       2.08
#> 4 disp        disp      -0.0676       0.5      0.8       1.09        3.83
#> 
#> $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        12.3 -0.00339      0.54      0.96      0.113      0.331
#> 2 .expected      12.3 -6.28         0.5       0.96      6.27      27.9  
#>