Skip to content

Commit 2dd897e

Browse files
Merge pull request #251 from canmod/bmb_calibration
calibration vignette tweaks
2 parents 7a22a5a + f5a93af commit 2dd897e

File tree

1 file changed

+15
-12
lines changed

1 file changed

+15
-12
lines changed

vignettes/calibration.Rmd

Lines changed: 15 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,7 @@ options(macpan2_verbose = FALSE)
3232

3333
Before reading this article on calibrating models to data, please first look at the [quickstart guide](https://canmod.github.io/macpan2/articles/quickstart) and the article on the [model library](https://canmod.github.io/macpan2/articles/example_models).
3434

35-
## Hello World
35+
## Hello, World
3636

3737
We'll do the first thing you should always do when trying out a new fitting procedure: simulate clean, nice data from the model and see if you can recover something close to the true parameters.
3838

@@ -54,7 +54,8 @@ sir_simulator = mp_simulator(sir_spec
5454
, outputs = c("S", "I", "R")
5555
, default = list(N = 300, R = 100, beta = 0.25, gamma = 0.1)
5656
)
57-
sir_results = mp_trajectory(sir_simulator)
57+
sir_results = mp_trajectory(sir_simulator) |>
58+
mutate(across(matrix, ~factor(., levels = c("S", "I", "R"))))
5859
(sir_results
5960
|> ggplot(aes(time, value, colour = matrix))
6061
+ geom_line()
@@ -87,14 +88,14 @@ The next step is to produce an object that can be calibrated through optimizatio
8788
```{r calibrator}
8889
sir_calibrator = mp_tmb_calibrator(sir_spec
8990
, data = sir_prevalence
90-
, traj = "I"
91+
, traj = list(I = mp_poisson()),
9192
, par = c("beta", "R")
9293
, default = list(N = 300)
9394
)
9495
print(sir_calibrator)
9596
```
9697

97-
Note that the calibrator has a few new expressions that deal with comparisons with data. In particular it is the objective function that we will optimize. But before that we can do a sanity check to make sure that the default values give a reasonable-looking trajectory.
98+
The calibrator has a few new expressions that deal with comparisons with data; in particular, it defines the objective function that we will optimize. But before that we can do a sanity check to make sure that the default values give a reasonable-looking trajectory.
9899

99100
```{r plot1}
100101
(sir_calibrator
@@ -107,32 +108,34 @@ Note that the calibrator has a few new expressions that deal with comparisons wi
107108

108109
### Step 2: do the fit
109110

110-
Doing the fit is straightforward.
111+
Doing the fit is straightforward; this calls a nonlinear optimizer built into base R (`nlminb` by default), starting from the default values specified in the calibrator.
111112

112113
```{r sir_fit, results = "hide"}
113114
mp_optimize(sir_calibrator)
114115
```
115116

116-
Note that the `mp_optimize` function has modified the `sir_calibrator` object, which now contains the new fitted parameter values and the results of the optimization.
117+
The `mp_optimize` function has **modified** the `sir_calibrator` object in place; it now contains the new fitted parameter values and the results of the optimization.
117118

118119
### Step 3: check the fit
119120

120-
We can print the results of the optimizer (`nlminb` in this case) using the `mp_optimizer_output` function. **Always check the value of the convergence code** (if it's not 0, then something *may* have gone wrong ...).
121+
We can print the results of the optimizer (`nlminb` in this case) using the `mp_optimizer_output` function. **Always check the value of the convergence code** (if it's not 0, then something may have gone wrong ...).
121122

122123
```{r check_fit}
123124
mp_optimizer_output(sir_calibrator)
124125
mp_optimize(sir_calibrator)
125126
mp_optimizer_output(sir_calibrator, what="all")
126127
```
127128

128-
As mentioned above, the best-fit parameters are stored internally, and we can get information about
129-
them using the `mp_tmb_coef` function. (Note that if you get a message about the `broom.mixed` package, please install it. `mp_tmb_coef` is a wrapper for this function).
129+
As mentioned above, the best-fit parameters are stored internally; we can get information about
130+
them using the `mp_tmb_coef` function. (If you get a message about the `broom.mixed` package, please install it. `mp_tmb_coef` is a wrapper for `broom.mixed::tidy()`).
131+
130132
```{r params}
131133
sir_estimates = mp_tmb_coef(sir_calibrator, conf.int = TRUE)
132-
print(sir_estimates)
134+
print(sir_estimates, digits = 3)
133135
```
134136

135-
These correspond pretty well to the known true values of the simulation model.
137+
This parameter corresponds pretty well to the known true values we used to simulate.
138+
136139
```{r default_reminder}
137140
mp_default(sir_simulator) |> filter(matrix %in% sir_estimates$mat)
138141
```
@@ -173,4 +176,4 @@ $$
173176

174177
Given these assumptions we choose $\mathbf b$ to maximize the resulting likelihood function, and use functionality from the `TMB` package (and sometimes the `tmbstan`/`rstan` packages) to do statistical inference on the fitted parameters and trajectories.
175178

176-
We recognize that this statistical model will often be overly restrictive. The `macpan2` package has a developer interface that is much more flexible, allowing for more detailed control over `TMB`, `tmbstan`, and `rstan`. This interface allows for arbitrary likelihood functions, prior distributions, parameter transformations, flexible parameter time-variation models, random effects and more. See [here](https://canmod.github.io/macpan2/articles/calibration_advanced.html) and [here](https://canmod.github.io/macpan2/articles/time_varying_parameters.html) for more information, although because these guides describe a developer interface the instructions may be unclear to some or many readers. Our plan is to continue adding interface layers, such as the interface described in this vignette, so that more of `macpan2` can be exposed to users.
179+
This statistical model will often be too simple. The `macpan2` package has an extremely flexible developer interface that allows for more detailed control over `TMB`, `tmbstan`, and `rstan`. This interface allows for arbitrary likelihood functions, prior distributions, parameter transformations, flexible parameter time-variation models, random effects and more. See [here](https://canmod.github.io/macpan2/articles/calibration_advanced.html) and [here](https://canmod.github.io/macpan2/articles/time_varying_parameters.html) for more information, although because these guides describe a developer interface the instructions may be unclear to some or many readers. Our plan is to continue adding interface layers, such as the interface described in this vignette, so that more of `macpan2` can be exposed to users.

0 commit comments

Comments
 (0)