Skip to contents

Fit a Model

Source: R/bage_mod-methods.R
fit.bage_mod.Rd

Derive the posterior distribution for a model.

Usage

# S3 method for class 'bage_mod'
fit(
  object,
  method = c("standard", "inner-outer"),
  vars_inner = NULL,
  optimizer = c("multi", "nlminb", "BFGS", "CG"),
  quiet = TRUE,
  start_oldpar = FALSE,
  ...
)

Arguments

object

A bage_mod object, created with mod_pois(), mod_binom(), or mod_norm().

method

Estimation method. Current choices are "standard" (the default) and "inner-outer". See below for details.

vars_inner

Names of variables to use for inner model when method is "inner-outer". If NULL(the default)vars_inner` is the age, sex/gender, and time variables.

optimizer

Which optimizer to use. Current choices are "multi", "nlminb", "BFGS", and "CG". Default is "multi". See below for details.

quiet

Whether to suppress warnings and progress messages from the optimizer. Default is TRUE.

start_oldpar

Whether the optimizer should start at previous estimates. Used only when fit() is being called on a fitted model. Default is FALSE.

...

Not currently used.

Value

A bage_mod object

Estimation methods

When method is "standard" (the default), all parameters, other than the lowest-level rates, probabilities, or means are jointly estimated within TMB.

When method is "inner-outer", estimation is carried out in multiple steps, which, in large models, can sometimes reduce computation times. In Step 1, the data is aggregated across all dimensions other than those specified in var_inner, and a model for the inner variables is fitted to the data. In Step 2, the data is aggregated across the remaining variables, and a model for the outer variables is fitted to the data. In Step 3, values for dispersion are calculated. Parameter estimates from steps 1, 2, and 3 are then combined. "inner-outer" methods are still experimental, and may change in future.

Optimizer

The choices for the optimizer argument are:

  • "multi" Try "nlminb", and if that fails, restart from the parameter values where "nlminb" stopped, using "BFGS". The default.

  • "nlminb" stats::nlminb()

  • "BFGS" stats::optim() using method "BFGS".

  • "GC" stats::optim() using method "CG" (conjugate gradient).

See also

Examples

## specify model
mod <- mod_pois(injuries ~ age + sex + year,
                data = nzl_injuries,
                exposure = popn)

## examine unfitted model
mod
#> 
#>     ------ Unfitted Poisson model ------
#> 
#>    injuries ~ age + sex + year
#> 
#>   exposure = popn
#> 
#>         term  prior along n_par n_par_free
#>  (Intercept) NFix()     -     1          1
#>          age   RW()   age    12         12
#>          sex NFix()     -     2          2
#>         year   RW()  year    19         19
#> 
#>  disp: mean = 1
#> 
#>  n_draw var_time var_age var_sexgender
#>    1000     year     age           sex
#> 

## fit model
mod <- fit(mod)
#> Building log-posterior function...
#> Finding maximum...
#> Drawing values for hyper-parameters...

## examine fitted model
mod
#> 
#>     ------ Fitted Poisson model ------
#> 
#>    injuries ~ age + sex + year
#> 
#>   exposure = popn
#> 
#>         term  prior along n_par n_par_free std_dev
#>  (Intercept) NFix()     -     1          1       -
#>          age   RW()   age    12         12    0.76
#>          sex NFix()     -     2          2    0.72
#>         year   RW()  year    19         19    0.08
#> 
#>  disp: mean = 1
#> 
#>  n_draw var_time var_age var_sexgender optimizer
#>    1000     year     age           sex    nlminb
#> 
#>  time_total time_max time_draw iter converged                    message
#>        0.25     0.11      0.11   13      TRUE   relative convergence (4)
#> 

## extract rates
aug <- augment(mod)
#> Drawing `.expected`...
#> Drawing `.fitted`...
aug
#> # A tibble: 912 × 9
#>    age   sex    ethnicity  year injuries  popn .observed
#>    <fct> <chr>  <chr>     <int>    <int> <int>     <dbl>
#>  1 0-4   Female Maori      2000       12 35830 0.000335 
#>  2 5-9   Female Maori      2000        6 35120 0.000171 
#>  3 10-14 Female Maori      2000        3 32830 0.0000914
#>  4 15-19 Female Maori      2000        6 27130 0.000221 
#>  5 20-24 Female Maori      2000        6 24380 0.000246 
#>  6 25-29 Female Maori      2000        6 24160 0.000248 
#>  7 30-34 Female Maori      2000       12 22560 0.000532 
#>  8 35-39 Female Maori      2000        3 22230 0.000135 
#>  9 40-44 Female Maori      2000        6 18130 0.000331 
#> 10 45-49 Female Maori      2000        6 13770 0.000436 
#> # ℹ 902 more rows
#> # ℹ 2 more variables: .fitted <rdbl<1000>>, .expected <rdbl<1000>>

## extract hyper-parameters
comp <- components(mod)
comp
#> # A tibble: 37 × 4
#>    term        component level                   .fitted
#>    <chr>       <chr>     <chr>              <rdbl<1000>>
#>  1 (Intercept) effect    (Intercept)  -2.5 (-4.2, -0.76)
#>  2 age         effect    0-4          -2.4 (-4.1, -0.77)
#>  3 age         effect    5-9           -3.7 (-5.4, -2.1)
#>  4 age         effect    10-14           -3.3 (-5, -1.6)
#>  5 age         effect    15-19        -1.6 (-3.2, 0.093)
#>  6 age         effect    20-24         -1.4 (-3.1, 0.21)
#>  7 age         effect    25-29        -1.6 (-3.2, 0.092)
#>  8 age         effect    30-34       -1.7 (-3.3, -0.027)
#>  9 age         effect    35-39       -1.7 (-3.3, -0.025)
#> 10 age         effect    40-44       -1.7 (-3.3, -0.025)
#> # ℹ 27 more rows