|
| 1 | +library(macpan2); library(tidyverse) |
| 2 | +sir = mp_tmb_library("starter_models", "sir", package = "macpan2") |
| 3 | +data = read_csv("https://raw.githubusercontent.com/canmod/macpan-workshop/refs/heads/main/data/early_covid_on_reports.csv") |
| 4 | +cleaned_data = (data |
| 5 | + |> rename(time = date, matrix = var) |
| 6 | + |> select(-province) |
| 7 | +) |
| 8 | +sir_covid = mp_tmb_insert(sir |
| 9 | + |
| 10 | + ## Modify the the part of the model that |
| 11 | + ## is executed once "during" each iteration |
| 12 | + ## of the simulation loop. |
| 13 | + , phase = "during" |
| 14 | + |
| 15 | + ## Insert the new expressions at the |
| 16 | + ## end of each iteration. |
| 17 | + ## (Inf ~ infinity ~ end) |
| 18 | + , at = Inf |
| 19 | + |
| 20 | + ## Specify the model for under-reporting as |
| 21 | + ## a simple fraction of the number of new |
| 22 | + ## cases every time-step. |
| 23 | + , expressions = list(report ~ report_prob * infection) |
| 24 | + |
| 25 | + ## Update defaults to something a little more |
| 26 | + ## reasonable for early covid in Ontario. |
| 27 | + , default = list( |
| 28 | + gamma = 1/14 ## 2-week recovery |
| 29 | + , beta = 2.5/14 ## R0 = 2.5 |
| 30 | + , report_prob = 0.1 ## 10% of cases get reported |
| 31 | + , N = 1.4e7 ## ontario population |
| 32 | + , I = 50 ## start with undetected infections |
| 33 | + ) |
| 34 | +) |
| 35 | +cal = mp_tmb_calibrator(sir_covid, cleaned_data |
| 36 | + , traj = "report" |
| 37 | + , par = c("beta", "gamma") |
| 38 | +) |
| 39 | +mp_optimize(cal) |
| 40 | +fitted_data = mp_trajectory(cal) |
| 41 | +(ggplot() |
| 42 | + + geom_line(aes(time, value), data = fitted_data) |
| 43 | + + geom_point(aes(time, value), data = cleaned_data) |
| 44 | +) |
0 commit comments