starter model

canmod · Oct 9, 2024 · 2d0d7f8 · 2d0d7f8
1 parent 88f9635
commit 2d0d7f8
Show file tree

Hide file tree

Showing 3 changed files with 241 additions and 59 deletions.
diff --git a/inst/starter_models/sir_mosquito/README.Rmd b/inst/starter_models/sir_mosquito/README.Rmd
@@ -22,45 +22,13 @@ knitr::opts_chunk$set(
 
 This [model](https://homepages.warwick.ac.uk/~masfz/ModelingInfectiousDiseases/Chapter4/Program_4.4/index.html) from [@keeling2011modeling] describes the transmission dynamics between humans and mosquitoes, focusing on a vector-borne disease.
 
-The code in this article uses the following packages.
+# Packages Used
 
 ```{r packages, message=FALSE, warning=FALSE}
 library(macpan2)
 library(ggplot2)
 ```
 
-# States
-
-| **Variable** | **Description**                                  |
-|--------------|--------------------------------------------------|
-| $S_H$          | Number of susceptible humans                     |
-| $I_H$          | Number of infected humans                        |
-| $S_M$          | Number of susceptible mosquitoes                 |
-| $I_M$          | Number of infected mosquitoes                    |
-
-# Parameters
-
-| **Variable**  | **Description**                                  |
-|---------------|--------------------------------------------------|
-| $r$             | Rate at which humans are bitten |
-| $\beta_{ij}$             | Transmission probability (following a bite) to species $i$ from species $j$.|
-| $\gamma_i$             | Recovery rate for host species $i$.                          |
-| $\nu_i$             | Birth rate for host species $i$.                                 |
-| $\mu_i$             | Death rate for species $i$                    |
-
-# Dynamics
-
-The model is described by the following differential equations:
-
-$$
-\begin{align*}
-\frac{dS_H}{dt} &= \nu_H - r \beta_{HM} S_H I_M - \mu_H S_H \\
-\frac{dI_H}{dt} &= r \beta_{HM} S_H I_M - \gamma_H I_H - \mu_H I_H \\
-\frac{dS_M}{dt} &= \nu_M - r \beta_{MH} S_M I_H - \mu_M S_M \\
-\frac{dI_M}{dt} &= r \beta_{MH} S_M I_H - \gamma_M I_M - \mu_M I_M \\
-\end{align*}
-$$
-
 # Model Specification
 
 This model has been specified in the `sir_mosquito` directory [here](https://github.com/canmod/macpan2/blob/main/inst/starter_models/sir_mosquito/tmb.R) and is accessible from the `macpan2` model library (see [Example Models](https://canmod.github.io/macpan2/articles/example_models.html) for details). 
@@ -108,6 +76,39 @@ layout = mp_layout_paths(spec)
 ```
 
 
+# States
+
+| **Variable** | **Description**                                  |
+|--------------|--------------------------------------------------|
+| $S_H$          | Number of susceptible humans                     |
+| $I_H$          | Number of infected humans                        |
+| $S_M$          | Number of susceptible mosquitoes                 |
+| $I_M$          | Number of infected mosquitoes                    |
+
+# Parameters
+
+| **Variable**  | **Description**                                  |
+|---------------|--------------------------------------------------|
+| $r$             | Rate at which humans are bitten |
+| $\beta_{ij}$             | Transmission probability (following a bite) to species $i$ from species $j$.|
+| $\gamma_i$             | Recovery rate for host species $i$.                          |
+| $\nu_i$             | Birth rate for host species $i$.                                 |
+| $\mu_i$             | Death rate for species $i$                    |
+
+# Dynamics
+
+The model is described by the following differential equations:
+
+$$
+\begin{align*}
+\frac{dS_H}{dt} &= \nu_H - r \beta_{HM} S_H I_M - \mu_H S_H \\
+\frac{dI_H}{dt} &= r \beta_{HM} S_H I_M - \gamma_H I_H - \mu_H I_H \\
+\frac{dS_M}{dt} &= \nu_M - r \beta_{MH} S_M I_H - \mu_M S_M \\
+\frac{dI_M}{dt} &= r \beta_{MH} S_M I_H - \gamma_M I_M - \mu_M I_M \\
+\end{align*}
+$$
+
+
 # Simulation
 
 Simulation of this and other models depends on the kind of [state update](https://canmod.github.io/macpan2/reference/mp_euler) that you use. Here we use RK4, which is a standard ODE solver.

diff --git a/inst/starter_models/sir_mosquito/README.md b/inst/starter_models/sir_mosquito/README.md
@@ -2,11 +2,12 @@ Mosquito-Vector SIR
 ================
 Steve Walker
 
+-   <a href="#packages-used" id="toc-packages-used">Packages Used</a>
+-   <a href="#model-specification" id="toc-model-specification">Model
+    Specification</a>
 -   <a href="#states" id="toc-states">States</a>
 -   <a href="#parameters" id="toc-parameters">Parameters</a>
 -   <a href="#dynamics" id="toc-dynamics">Dynamics</a>
--   <a href="#model-specification" id="toc-model-specification">Model
-    Specification</a>
 -   <a href="#simulation" id="toc-simulation">Simulation</a>
 -   <a href="#references" id="toc-references">References</a>
 
@@ -16,13 +17,36 @@ from ([Keeling and Rohani 2011](#ref-keeling2011modeling)) describes the
 transmission dynamics between humans and mosquitoes, focusing on a
 vector-borne disease.
 
-The code in this article uses the following packages.
+# Packages Used
 
 ``` r
 library(macpan2)
 library(ggplot2)
 ```
 
+# Model Specification
+
+This model has been specified in the `sir_mosquito` directory
+[here](https://github.com/canmod/macpan2/blob/main/inst/starter_models/sir_mosquito/tmb.R)
+and is accessible from the `macpan2` model library (see [Example
+Models](https://canmod.github.io/macpan2/articles/example_models.html)
+for details).
+
+``` r
+spec = mp_tmb_library("starter_models", "sir_mosquito", package = "macpan2")
+```
+
+This specification can be used to draw the following flow diagram using
+code found in the [source for this
+article](https://github.com/canmod/macpan2/blob/main/inst/starter_models/sir_mosquito/README.Rmd).
+The dashed lines show the state dependence of infection, which is the
+only way that the human and mosquito sub-models interact (e.g.,
+mosquitos never become humans, they just infect them). The per-capita
+infection rate for humans depends on the number of infectious mosquitos,
+and vice versa.
+
+![](./figures/diagram-1.png)<!-- -->
+
 # States
 
 | **Variable** | **Description**                  |
@@ -55,29 +79,6 @@ $$
 \end{align*}
 $$
 
-# Model Specification
-
-This model has been specified in the `sir_mosquito` directory
-[here](https://github.com/canmod/macpan2/blob/main/inst/starter_models/sir_mosquito/tmb.R)
-and is accessible from the `macpan2` model library (see [Example
-Models](https://canmod.github.io/macpan2/articles/example_models.html)
-for details).
-
-``` r
-spec = mp_tmb_library("starter_models", "sir_mosquito", package = "macpan2")
-```
-
-This specification can be used to draw the following flow diagram using
-code found in the [source for this
-article](https://github.com/canmod/macpan2/blob/main/inst/starter_models/sir_mosquito/README.Rmd).
-The dashed lines show the state dependence of infection, which is the
-only way that the human and mosquito sub-models interact (e.g.,
-mosquitos never become humans, they just infect them). The per-capita
-infection rate for humans depends on the number of infectious mosquitos,
-and vice versa.
-
-![](./figures/diagram-1.png)<!-- -->
-
 # Simulation
 
 Simulation of this and other models depends on the kind of [state

diff --git a/inst/starter_models/sir_waning/README.Rmd b/inst/starter_models/sir_waning/README.Rmd
@@ -92,4 +92,184 @@ $$
 \end{align*}
 $$
 
+# Calibration
+
+## Simulate fake data
+
+We modify the specification so that it is different from the default library model that we will calibrate. We will use the [`mp_rk4`](https://canmod.github.io/macpan2/reference/mp_euler.html) ODE solver.
+
+# simulator object
+sir_waning = mp_simulator(  
+    model = spec
+  , time_steps = time_steps
+  , outputs = c("I","S", "waning_immunity")
+)
+
+## -------------------------
+## specify objective function
+## -------------------------
+
+# negative log likelihood
+obj_fn = ~ -sum(dpois(I_obs, rbind_time(I, I_obs_times)))
+
+# update simulator to create new variables 
+# I_obs and I_obs_times and initialize
+sir_waning$update$matrices(
+    I_obs = empty_matrix
+  , I_obs_times = empty_matrix
+)
+
+# update simulator to include this function
+sir_waning$replace$obj_fn(obj_fn)
+
+## -------------------------
+## parameterize model
+## -------------------------
+
+# choose which parameter(s) to estimate - log(beta) and phi
+sir_waning$update$transformations(Log("beta"))
+sir_waning$update$transformations(Log("phi"))
+
+sir_waning$replace$params(c(log(spec$default$beta),log(spec$default$phi)),
+                          c("log_beta","log_phi")
+                          )
+
+sir_waning
+
+
+## -------------------------
+## simulate fake data
+## -------------------------
+
+# beta value to simulate data with
+true_beta = 0.3
+
+# phi value to simulate data with
+true_phi = 0.09
+
+## simulate observed data using true parameters
+observed_data = sir_waning$report(c(log(true_beta),log(true_phi)))
+
+## compute incidence for observed data
+I_obs = rpois(time_steps, subset(observed_data, matrix == "I", select = c(value)) %>% pull())
+I_obs_times = subset(observed_data, matrix == "I", select = c(time)) %>% pull()
+
+if (interactive()) {
+  plot(I_obs, type = "l", las = 1)
+}
+
+
+## -------------------------
+## update simulator with fake data to fit to
+## -------------------------
+
+sir_waning$update$matrices(
+    I_obs = I_obs
+  , I_obs_times = I_obs_times
+)
+
+## -------------------------
+## plot likelihood surface (curve)
+## -------------------------
+
+# plot surface as contours
+if (interactive()) {
+  log_betas = seq(from = log(0.1), to = log(1), length = 100)
+  log_phis = seq(from = log(1e-3), to = log(0.2), length = 100)
+  x_y = expand.grid(log_betas, log_phis) %>% setNames(c("log_betas","log_phis"))
+
+  ll = apply(
+      x_y
+    , 1
+    , function(z) {sir_waning$objective(z["log_betas"], z["log_phis"])}
+  ) 
+
+  dat_for_plot <- cbind(x_y, ll)
+
+  ggplot(dat_for_plot, aes(log_betas, log_phis, z=ll)) +
+    geom_contour_filled()+
+    ## add true parameter values to compare
+    geom_vline(xintercept = log(true_beta), col='red')+
+    geom_hline(yintercept = log(true_phi), col='red')
+
+
+}
+
+## -------------------------
+## fit parameters
+## -------------------------
+
+## optimize and check convergence
+## warning message, but converges
+sir_waning$optimize$nlminb()
+
+## plot observed vs predicted
+if (interactive()) {
+
+  ## estimates are close to true values
+  print(sir_waning$current$params_frame())
+  print(paste0("exp(default beta) ",exp(sir_waning$current$params_frame()$default[1])))
+  print(paste0("exp(current beta) ",exp(sir_waning$current$params_frame()$current[1])))
+  print(paste0("exp(default phi) ",exp(sir_waning$current$params_frame()$default[2])))
+  print(paste0("exp(current phi) ",exp(sir_waning$current$params_frame()$current[2])))
+
+
+  data_to_plot <- (cbind(as.numeric(I_obs),1:time_steps)
+                   %>% data.frame()
+                   %>% setNames(c("value","time"))
+                   %>% mutate(type="observed")
+  ) %>% union(sir_waning$report() 
+              %>% filter(matrix=="I") 
+              %>% select(time,value)
+              %>% mutate(type="predicted")
+  )
+
+  ggplot(data_to_plot, aes(x=time, y=value, col=type))+
+    geom_line()+
+    theme_bw()+
+    ylab("I")
+
+}
+
+## -------------------------
+## exploring
+## -------------------------
+
+## plot S
+if (interactive()) {
+  plot(sir_waning$report() %>% filter(matrix=="S") %>% select(time,value),
+       type="l", las = 1, ylab='S')
+}
+
+## plot waning immunity (phi*R)
+if (interactive()) {
+  plot(sir_waning$report() %>% filter(matrix=="waning_immunity") %>% select(time,value), type = "l", las = 1, ylab='Waning Immunity')
+}
+
+## -------------------------
+## exercise random effects a bit
+## -------------------------
+
+## In real life we would add an informative or partially informative
+## prior for the random effect, but here we just use uniform for log phi.
+
+matrix_version = "1.6-5"
+if (packageVersion("Matrix") >= matrix_version) {
+  sir_waning$reset$params()
+  sir_waning$replace$params(log(spec$default$beta), "log_beta")
+  sir_waning$replace$random(log(spec$default$phi), "log_phi")
+  sir_waning$optimize$nlminb()
+
+  ## They both seem to be reasonably close to the true values
+  print(sir_waning$sdreport())
+  print(exp(sir_waning$sdreport()$par.fixed))
+  print(exp(sir_waning$sdreport()$par.random))
+  print(sir_waning$current$params_frame())
+  print(sir_waning$current$random_frame())
+  if (require(tmbstan) & require(broom.mixed)) {
+    mp_tmbstan_coef(sir_waning, conf.int = TRUE, drop.pars = character())
+  }
+}
+
+
 # References