calibration sir example

inst/starter_models/lotka_volterra/README.md

@@ -0,0 +1,19 @@
+---
+title: "Lotka-Volterra"
+index_entry: "simple two-species competition model"
+author: Jennifer Freeman
+---
+
+- $X$ - number of individuals in species $X$  
+- $Y$ - number of individuals in species $Y$ 
+- $r_i$ - growth rate of species $I$
+- $a_{ij}$ - intra/inter-specific density dependence, ``effect of species $j$ on species $i$'' (Hastings, 1997)
+
+$$
+\begin{align*}
+\frac{dX}{dt} &= r_x X (1 - a_{xx}X - a_{xy}Y) \\
+\frac{dY}{dt} &= r_y Y (1 - a_{yy}Y - a_{yx}X)
+\end{align*}
+$$
+
+Hastings, A. (1997). Competition. In: Population Biology. Springer, New York, NY. https://doi.org/10.1007/978-1-4757-2731-9_7

inst/starter_models/lotka_volterra/model.R

@@ -0,0 +1,39 @@
+library(macpan2)
+
+## expression lists
+###################################################
+
+## free form computations for convenience
+## none
+
+## absolute flow rates (per time only)
+## reference 
+flow_rates = list(
+  ## growth rate of species X and Y
+    growth_x ~ rx * X
+  , growth_y ~ ry * Y
+  ## intraspecific effect
+  ## species X on X
+  , intraspecific_x ~ growth_x * axx * X
+  ## species Y on Y
+  , intraspecific_y ~ growth_y * ayy * Y
+  ## interspecific effect
+  ## species X on Y
+  , interspecific_xy ~ growth_x * ayx * Y
+  ## species Y on X
+  , interspecific_yx ~ growth_y * axy * X
+)
+
+## state updates
+state_updates = list(
+    X ~ X + growth_x - intraspecific_x - interspecific_xy
+  , Y ~ Y + growth_y - intraspecific_y - interspecific_yx
+)
+
+## simple unstructured scalar expression
+expr_list =  ExprList(
+  during = c(
+      flow_rates
+    , state_updates
+  )
+)

inst/starter_models/macpan_base/model.R

@@ -0,0 +1,52 @@
+library(macpan2)
+
+## expression lists
+###################################################
+
+## free form computations for convenience
+computations = list(
+  N ~ sum(S, E, Ia, Ip, Im, Is, R, H, ICUs, ICUd, H2, D)
+)
+
+## absolute flow rates (per time only)
+flow_rates = list(
+    SE ~ S * (beta0 / N) * (Ia * Ca + Ip * Cp + Im * Cm * (1 - iso_m) + Is * Cs *(1 - iso_s))
+  , EIa ~ E * alpha * sigma
+  , EIp ~ E * (1 - alpha)* sigma
+  , IaR ~ Ia * gamma_a
+  , IpIm ~ Ip * mu * gamma_p
+  , ImR ~ Im * gamma_m
+  , IpIs ~ Ip * (1 - mu) * gamma_p
+  , IsICUs ~ Is * (1 - nonhosp_mort) * (1 - phi1) * (1 - phi2) * gamma_s
+  , ICUsH2 ~ ICUs * psi1
+  , H2R ~ H2 * psi3
+  , IsH ~ Is * (1 - nonhosp_mort) * phi1 * gamma_s
+  , IsICUd ~ Is * (1 - nonhosp_mort) * (1 - phi1) * phi2 * gamma_s
+  , ICUdD ~ ICUd * psi2
+  , HR ~ H * rho
+)
+
+## state updates
+state_updates = list(
+    S ~ S - SE
+  , E ~ E + SE - EIa - EIp
+  , Ia ~ Ia + EIa - IaR
+  , Ip ~ Ip + EIp - IpIm - IpIs
+  , Im ~ Im + IpIm
+  , Is ~ Is + IpIs - IsICUs - IsH - IsICUd
+  , R ~ R + IaR + H2R + HR
+  , H ~ H + IsH - HR
+  , ICUs ~ ICUs + IsICUs - ICUsH2
+  , ICUd ~ ICUd + IsICUd - ICUdD
+  , H2 ~ H2 + ICUsH2 - H2R
+  , D ~ D + ICUdD
+)
+
+## simple unstructured scalar expression
+expr_list =  ExprList(
+    before = computations
+  , during = c(
+      flow_rates
+    , state_updates
+  )
+)

inst/starter_models/seir/model.R

@@ -0,0 +1,33 @@
+library(macpan2)
+
+## expression lists
+###################################################
+
+## free form computations for convenience
+computations = list(
+  N ~ sum(S, E, I, R)
+)
+
+## absolute flow rates (per time only)
+flow_rates = list(
+    infection ~ S * I * beta / N
+  , exposure ~ alpha * E
+  , recovery ~ gamma * I
+)
+
+## state updates
+state_updates = list(
+    S ~ S - exposure
+  , E ~ E + exposure - infection
+  , I ~ I + infection - recovery
+  , R ~ R + recovery
+)
+
+## simple unstructured scalar expression
+expr_list =  ExprList(
+    before = computations
+  , during = c(
+      flow_rates
+    , state_updates
+  )
+)

inst/starter_models/sir/example.R

@@ -15,15 +15,21 @@ tmb_simulator  ## note the new expression before the simulation loop
 time_steps = 100L
 true_beta = 0.4
 
+## set time_steps value
 tmb_simulator$replace$time_steps(time_steps)
+
+## feed log(true_beta) to the simulator because we have
+## already specified log-transformation of this parameter
 observed_data = tmb_simulator$report(log(true_beta))
+
+## .mats_to_return is set to "I", so observed_data$value is
+## the prevalence (density of I) over time
 observed_data$value = rpois(time_steps, observed_data$value)
 
 if (interactive()) {
   plot(observed_data$value, type = "l", las = 1)
 }
 
-
 ## -------------------------
 ## update simulator with fake data to fit to
 ## -------------------------
@@ -35,7 +41,7 @@ tmb_simulator$update$matrices(
 
 
 ## -------------------------
-## plot likelihood surface
+## plot likelihood surface (curve)
 ## -------------------------
 
 if (interactive()) {
@@ -51,11 +57,12 @@ if (interactive()) {
 
 
 ## -------------------------
-## optimize the model
+## fit parameters
 ## -------------------------
 
 tmb_simulator$optimize$nlminb()
 
+## plot observed vs predicted value
 if (interactive()) {
   print(tmb_simulator$current$params_frame())
   plot(observed_data$value, type = "l", las = 1)

inst/starter_models/sir_demo/model.R

@@ -0,0 +1,36 @@
+library(macpan2)
+
+## expression lists
+###################################################
+
+## free form computations for convenience
+computations = list(
+  N ~ sum(S, I, R)
+)
+
+## absolute flow rates (per time only)
+flow_rates = list(
+    birth ~ birth_rate * N
+  , infection ~ S * I * beta / N
+  , recovery ~ gamma * I
+)
+
+## state updates
+state_updates = list(
+    S ~ S - infection + birth - death_rate * S
+  , I ~ I + infection - recovery - death_rate * I
+  , R ~ R + recovery - death_rate * R
+)
+
+## simple unstructured scalar expression
+## in the general case,
+## birth rate and death rate can differ, N might not be constant
+## nothing can be computed before simulation loop
+expr_list =  ExprList(
+    during = c(
+        computations
+      , flow_rates
+      , state_updates
+  )
+)
+

inst/starter_models/sir_waning/model.R

@@ -0,0 +1,32 @@
+library(macpan2)
+
+## expression lists
+###################################################
+
+## free form computations for convenience
+computations = list(
+  N ~ sum(S, I, R)
+)
+
+## absolute flow rates (per time only)
+flow_rates = list(
+    infection ~ S * I * beta / N
+  , recovery ~ gamma * I
+  , waning_immunity ~ wane * R
+)
+
+## state updates
+state_updates = list(
+    S ~ S - infection + waning_immunity
+  , I ~ I + infection - recovery
+  , R ~ R + recovery - waning_immunity
+)
+
+## simple unstructured scalar expression
+expr_list =  ExprList(
+    before = computations
+  , during = c(
+      flow_rates
+    , state_updates
+  )
+)

inst/starter_models/ww/model.R

@@ -0,0 +1,59 @@
+library(macpan2)
+
+## expression lists
+###################################################
+
+## free form computations for convenience
+computations = list(
+  N ~ sum(S, E, Ia, Ip, Im, Is, R, H, ICUs, ICUd, H2, D, W, A)
+)
+
+## absolute flow rates (per time only)
+flow_rates = list(
+    SE ~ S * (beta0 / N) * (Ia * Ca + Ip * Cp + Im * Cm * (1 - iso_m) + Is * Cs *(1 - iso_s))
+  , EIa ~ E * alpha * sigma
+  , EIp ~ E * (1 - alpha)* sigma
+  , IaR ~ Ia * gamma_a
+  , IpIm ~ Ip * mu * gamma_p
+  , ImR ~ Im * gamma_m
+  , IpIs ~ Ip * (1 - mu) * gamma_p
+  , IsICUs ~ Is * (1 - nonhosp_mort) * (1 - phi1) * (1 - phi2) * gamma_s
+  , ICUsH2 ~ ICUs * psi1
+  , H2R ~ H2 * psi3
+  , IsH ~ Is * (1 - nonhosp_mort) * phi1 * gamma_s
+  , IsICUd ~ Is * (1 - nonhosp_mort) * (1 - phi1) * phi2 * gamma_s
+  , ICUdD ~ ICUd * psi2
+  , HR ~ H * rho
+  , IaW ~ Ia * nu # or * or N?
+  , IpW ~ Ip * nu 
+  , ImW ~ Im * nu
+  , IsW ~ Is * nu
+  , WA ~ W * xi
+)
+
+## state updates
+state_updates = list(
+    S ~ S - SE
+  , E ~ E + SE - EIa - EIp
+  , Ia ~ Ia + EIa - IaR
+  , Ip ~ Ip + EIp - IpIm - IpIs
+  , Im ~ Im + IpIm
+  , Is ~ Is + IpIs - IsICUs - IsH - IsICUd
+  , R ~ R + IaR + H2R + HR
+  , H ~ H + IsH - HR
+  , ICUs ~ ICUs + IsICUs - ICUsH2
+  , ICUd ~ ICUd + IsICUd - ICUdD
+  , H2 ~ H2 + ICUsH2 - H2R
+  , D ~ D + ICUdD
+  , W ~ W + IaW + IpW + ImW + IsW - WA
+  , A ~ A + WA
+)
+
+## simple unstructured scalar expression
+expr_list =  ExprList(
+    before = computations
+  , during = c(
+      flow_rates
+    , state_updates
+  )
+)