Code
+
@@ -863,9 +864,9 @@
-
Code
+
@@ -889,10 +890,10 @@
-
Code
+
Code
@@ -900,15 +901,15 @@
-
Code
+
Code
@@ -917,37 +918,37 @@
-
+
Code
viewof reset = Inputs.button([
- ["Reset all sliders", () => {
- set(viewof S0, reset_S)
- set(viewof I0, reset_I)
- set(viewof beta_input, reset_beta)
- }]
-])
-viewof S0 = Inputs.range(
- [500, 15000],
- {value: reset_S, step: 1, label: md`${tex`S(0)`}`}
-)
-
-viewof I0 = Inputs.range(
- [0.001, 50],
- {value: reset_I, step: 0.001, label: md`${tex`I(0)`}`}
-)
-
-viewof beta_input = Inputs.range(
- [1, 100],
- {value: reset_beta, step: 0.001, label: md`${tex`\beta (\times 10^{-3})`}`}
-)
-
-viewof site_select = Inputs.select(
- ["Site 1", "Site 2", "Site 3"],
- {label: "Select a site:"}
-)
-
-// convert to daily time scale as easier to manipulate
-beta = (beta_input / 365) * (10 ** (-3))
-
-md`${tex`R_0 = ${R0_str}`}`Code
viewof reset = Inputs.button([
+ ["Reset all sliders", () => {
+ set(viewof S0, reset_S)
+ set(viewof I0, reset_I)
+ set(viewof beta_input, reset_beta)
+ }]
+])
+viewof S0 = Inputs.range(
+ [500, 15000],
+ {value: reset_S, step: 1, label: md`${tex`S(0)`}`}
+)
+
+viewof I0 = Inputs.range(
+ [0.001, 50],
+ {value: reset_I, step: 0.001, label: md`${tex`I(0)`}`}
+)
+
+viewof beta_input = Inputs.range(
+ [1, 100],
+ {value: reset_beta, step: 0.001, label: md`${tex`\beta (\times 10^{-3})`}`}
+)
+
+viewof site_select = Inputs.select(
+ ["Site 1", "Site 2", "Site 3"],
+ {label: "Select a site:"}
+)
+
+// convert to daily time scale as easier to manipulate
+beta = (beta_input / 365) * (10 ** (-3))
+
+md`${tex`R_0 = ${R0_str}`}`
@@ -1000,11 +1001,11 @@
-
Code
+
Code
@@ -1035,8 +1036,8 @@
-
Code
+
@@ -1053,8 +1054,8 @@
-
Code
+
@@ -1071,13 +1072,13 @@
-
Code
+
Code
@@ -1085,21 +1086,21 @@
-
Code
+
Code
@@ -1107,7 +1108,7 @@
-
Code
+
@@ -1115,7 +1116,7 @@
-
Code
+
@@ -1123,30 +1124,30 @@
-
diff --git a/_book/r-session-03_files/figure-html/unnamed-chunk-30-1.png b/_book/r-session-03_files/figure-html/unnamed-chunk-30-1.png
index dbcc19f..3848a80 100644
Binary files a/_book/r-session-03_files/figure-html/unnamed-chunk-30-1.png and b/_book/r-session-03_files/figure-html/unnamed-chunk-30-1.png differ
diff --git a/_book/r-session-03_files/figure-html/unnamed-chunk-31-1.png b/_book/r-session-03_files/figure-html/unnamed-chunk-31-1.png
index 2f71477..0f0976b 100644
Binary files a/_book/r-session-03_files/figure-html/unnamed-chunk-31-1.png and b/_book/r-session-03_files/figure-html/unnamed-chunk-31-1.png differ
diff --git a/_book/r-session-03_files/figure-html/unnamed-chunk-34-1.png b/_book/r-session-03_files/figure-html/unnamed-chunk-34-1.png
index df4ff47..2e230d4 100644
Binary files a/_book/r-session-03_files/figure-html/unnamed-chunk-34-1.png and b/_book/r-session-03_files/figure-html/unnamed-chunk-34-1.png differ
diff --git a/_book/r-session-03_files/figure-html/unnamed-chunk-35-1.png b/_book/r-session-03_files/figure-html/unnamed-chunk-35-1.png
index 899cbab..93133ae 100644
Binary files a/_book/r-session-03_files/figure-html/unnamed-chunk-35-1.png and b/_book/r-session-03_files/figure-html/unnamed-chunk-35-1.png differ
diff --git a/_book/r-session-03_files/figure-html/unnamed-chunk-38-1.png b/_book/r-session-03_files/figure-html/unnamed-chunk-38-1.png
index 3c276b2..b98e979 100644
Binary files a/_book/r-session-03_files/figure-html/unnamed-chunk-38-1.png and b/_book/r-session-03_files/figure-html/unnamed-chunk-38-1.png differ
diff --git a/_book/r-session-03_files/figure-html/unnamed-chunk-4-1.png b/_book/r-session-03_files/figure-html/unnamed-chunk-4-1.png
index 7f19d74..d6856ea 100644
Binary files a/_book/r-session-03_files/figure-html/unnamed-chunk-4-1.png and b/_book/r-session-03_files/figure-html/unnamed-chunk-4-1.png differ
diff --git a/_book/r-session-03_files/figure-html/unnamed-chunk-43-1.png b/_book/r-session-03_files/figure-html/unnamed-chunk-43-1.png
index cec9d20..45a0d9c 100644
Binary files a/_book/r-session-03_files/figure-html/unnamed-chunk-43-1.png and b/_book/r-session-03_files/figure-html/unnamed-chunk-43-1.png differ
diff --git a/_book/r-session-03_files/figure-html/unnamed-chunk-8-1.png b/_book/r-session-03_files/figure-html/unnamed-chunk-8-1.png
index 9c7d9e0..31cc526 100644
Binary files a/_book/r-session-03_files/figure-html/unnamed-chunk-8-1.png and b/_book/r-session-03_files/figure-html/unnamed-chunk-8-1.png differ
diff --git a/_book/search.json b/_book/search.json
index c7bedea..78f8672 100644
--- a/_book/search.json
+++ b/_book/search.json
@@ -544,7 +544,7 @@
"href": "r-session-03.html#estimating-r_0-problem-background",
"title": "\n8 R Session 03\n",
"section": "\n8.2 Estimating \\(R_0\\) Problem Background",
- "text": "8.2 Estimating \\(R_0\\) Problem Background\nSo far in this class we have focused on the theory of infectious disease. Often, however, we will want to apply this theory to particular situations. One of the key applied problems in epidemic modeling is the estimation of \\(R_0\\) from outbreak data. In this session, we study two methods for estimating \\(R_0\\) from an epidemic curve. As a running example, we will use the data on influenza in a British boarding school.\n\nCodeggplot(flu, aes(x = day, y = flu)) +\n geom_line(color = \"slategray4\") +\n geom_point(shape = 21, size = 5, fill = \"slategray4\", alpha = 0.8) +\n labs(x = \"Day\", y = \"Active Influenza Cases\")",
+ "text": "8.2 Estimating \\(R_0\\) Problem Background\nSo far in this class we have focused on the theory of infectious disease. Often, however, we will want to apply this theory to particular situations. One of the key applied problems in epidemic modeling is the estimation of \\(R_0\\) from outbreak data. In this session, we study two methods for estimating \\(R_0\\) from an epidemic curve. As a running example, we will use the data on influenza in a British boarding school.\n\nCodeggplot(flu, aes(x = day, y = flu)) +\n geom_line(color = \"slategray4\") +\n geom_point(shape = 21, size = 5, fill = \"slategray4\", alpha = 0.8) +\n labs(x = \"Day\", y = \"Active Influenza Cases\")\n\n\n\n\n\n\n\nThe figure above, and the corresponding tabulated outbreak data, report the number of active influenza cases on each day. These are point-in-time counts, or prevalence: the number of students who were sick at that moment. By contrast, the final-size calculation below uses cumulative incidence: the total number of students who became cases at any point during the outbreak. For this boarding school outbreak, the total school size was \\(N = 764\\), and the total number of students who became cases over the full outbreak was \\(Z = 512\\). Because \\(512\\) is a cumulative incidence count, it is not obtained by reading a single value from the prevalence figure or table.",
"crumbs": [
"Day 2",
"8 R Session 03"
diff --git a/_freeze/r-session-03/execute-results/html.json b/_freeze/r-session-03/execute-results/html.json
index ae64fbd..daf358e 100644
--- a/_freeze/r-session-03/execute-results/html.json
+++ b/_freeze/r-session-03/execute-results/html.json
@@ -1,8 +1,8 @@
{
- "hash": "55bfc603c4d505ccf282c9d09ba5a3ca",
+ "hash": "6d03494b6000d56e4af9f6c1b125d774",
"result": {
"engine": "knitr",
- "markdown": "---\nsubtitle: \"Estimation\"\nabstract-title: \"\"\nabstract: |\n *Materials adapted from John Drake and Pej Rohani [@drakeEstimation2019]*\nexecute:\n warning: false\nmetadata-files:\n - metadata/matthewferrari.yml\n - metadata/mathjax-packages.yml\neditor:\n markdown:\n wrap: sentence\n---\n\n# R Session 03\n\n## Setup\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(here)\nlibrary(rio)\nlibrary(deSolve)\nlibrary(tidyverse)\nlibrary(ggtext)\nlibrary(gt)\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\ntheme_set(theme_minimal())\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Loads the datasets: flu, measles, niamey, plauge\nflu <- rio::import(\n \"https://raw.githubusercontent.com/arnold-c/SISMID-Module-02_2025/main/data/flu.csv\"\n)\nniamey <- rio::import(\n \"https://raw.githubusercontent.com/arnold-c/SISMID-Module-02_2025/main/data/niamey.csv\"\n)\n# flu <- rio::import(here::here(\"data\", \"flu.csv\"))\n# niamey <- rio::import(here::here(\"data\", \"niamey.csv\"))\n```\n:::\n\n\n## Estimating $R_0$ Problem Background\n\nSo far in this class we have focused on the **theory** of infectious disease.\nOften, however, we will want to apply this theory to particular situations.\nOne of the key applied problems in epidemic modeling is the estimation of $R_0$ from outbreak data.\nIn this session, we study two methods for estimating $R_0$ from an epidemic curve.\nAs a running example, we will use the data on influenza in a British boarding school.\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nggplot(flu, aes(x = day, y = flu)) +\n geom_line(color = \"slategray4\") +\n geom_point(shape = 21, size = 5, fill = \"slategray4\", alpha = 0.8) +\n labs(x = \"Day\", y = \"Active Influenza Cases\")\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\n## Estimating $R_0$ From The Final Outbreak Size\n\nOur first approach is to estimate $R_0$ from the final outbreak size.\nAlthough unhelpful at the early stages of an epidemic (before the final epidemic size is observed), this method is nonetheless a useful tool for *post hoc* analysis.\nThe method is general and can be motivated by the argument listed in [@keelingIntroductionSimpleEpidemic2008]:\n\nFirst, we assume that the epidemic is started by a single infectious individual in a completely susceptible population.\nOn average, this individual infects $R_0$ others.\nThe probability a particular individual escaped infection is therefore $e^{-R_0 / N}$.\n\nIf $Z$ individuals have been infected, the probability of an individual escaping infection from all potential sources is $e^{-Z R_0 / N}$.\nIt follows that at the end of the epidemic a proportion $R(\\infty) = Z / N$ have been infected and the fraction remaining susceptible is $S(\\infty) = e^{-R(\\infty) R_0}$, which is equal to $1 - R(\\infty)$.\n\n::: {.callout-note collapse=\"true\"}\n$S(\\infty) = e^{-R(\\infty) R_0}$ can be calculated by acknowledging that at equilibrium ($t = \\infty$), $S(\\infty) = 1 - R(\\infty) = Z / N$, so substituting $R(\\infty)$ into $1 - e^{-Z R_0 / N}$ gives the desired result.\n\nIt could also be calculated by dividing $\\frac{\\dd{S}}{\\dd{t}}$ by $\\frac{\\dd{R}}{\\dd{t}}$:\n\n```{=tex}\n\\begin{aligned}\n\\frac{\\dd{S}}{\\dd{R}} &= - \\frac{\\beta S}{\\gamma} \\\\\n&= - R_0 S\n\\end{aligned}\n```\nwhich is a [separable differential equation](https://tutorial.math.lamar.edu/classes/de/separable.aspx), so can be integrated as follows:\n\n```{=tex}\n\\begin{aligned}\n- \\int_{0}^{t} \\frac{1}{R_0 S} \\dd{S} &= \\int_{0}^{t} \\dd{R} \\\\\n- \\frac{1}{R_0} \\left(\\ln{S(t)} - \\ln{S(0)} \\right) &= R(t) - \\cancelto{0}{R(0)} \\\\\n\\ln{S(t)} &= \\ln{S(0)} - R_0 R(t) \\\\\nS(t) &= S(0) e^{-R_0 R(t)}\n\n\\end{aligned}\n```\n:::\n\nPutting this together, we get:\n\n$$\n1 - R(\\infty) - e^{-R(\\infty) R_0} = 0\n$$\n\nRearranging, we have the estimator\n\n$$\n \\hat{R_0} = \\frac{\\ln(1 - Z / N)}{-Z / N},\n$$\n\nwhich, in this case, evaluates to $\\frac{\\ln(1 - 512 / 764)}{-512 / 764} = 1.655$.\n\n::: exercise\n### Exercise 1\n:::\n\nThis equation shows the important one-to-one relationship between $R_0$ and the final epidemic size.\nPlot the relationship between the total epidemic size and $R_0$ for the complete range of values between 0 and 1.\n\n## Linear Approximation\n\nThe next method we introduce takes advantage of the fact that during the early stages of an outbreak, the number of infected individuals is given approximately as $I(t) \\approx I_0 e^{((R_0 - 1)(\\gamma + \\mu)t)}$.\nTaking logarithms of both sides, we have $\\ln(I(t)) \\approx \\ln(I_0) + (R_0 - 1)(\\gamma + \\mu)t$, showing that the log of the number of infected individuals is approximately linear in time with a slope that reflects both $R_0$ and the recovery rate.\n\nThis suggests that a simple linear regression fit to the first several data points on a log-scale, corrected to account for $\\gamma$ and $\\mu$, provides a rough and ready estimate of $R_0$.\nFor flu, we can assume $\\mu =0$ because the epidemic occurred over a time period during which natural mortality is negligible.\nFurther, assuming an infectious period of about 2.5 days, we use $\\gamma = (2.5)^{-1} = 0.4$ for the correction.\nFitting to the first four data points, we obtain the slope as follows.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Fit a linear model\nlinear_model <- lm(log(flu[1:4]) ~ day[1:4], data = flu)\n\n# Summary statistics for fit model\nsummary(linear_model)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n\nCall:\nlm(formula = log(flu[1:4]) ~ day[1:4], data = flu)\n\nResiduals:\n 1 2 3 4 \n 0.03073 -0.08335 0.07450 -0.02188 \n\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) -0.02703 0.10218 -0.265 0.81611 \nday[1:4] 1.09491 0.03731 29.346 0.00116 **\n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n\nResidual standard error: 0.08343 on 2 degrees of freedom\nMultiple R-squared: 0.9977,\tAdjusted R-squared: 0.9965 \nF-statistic: 861.2 on 1 and 2 DF, p-value: 0.001159\n```\n\n\n:::\n\n```{.r .cell-code}\n# Extract slope parameter\ncoef(linear_model)[2]\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\nday[1:4] \n1.094913 \n```\n\n\n:::\n:::\n\n\nRearranging the linear equation above and denoting the slope coefficient by $\\hat \\beta_1$ we have the estimator $\\hat R_0 = \\hat \\beta_1 / \\gamma + 1$ giving $\\hat R_0 = 1.094913 / 0.4 + 1 \\approx 3.7$.\n\n::: exercise\n### Exercise 2\n:::\n\nOur estimate assumes that boys remained infectious during the natural course of infection.\nThe original report on this epidemic indicates that boys found to have symptoms were immediately confined to bed in the infirmary.\nThe report also indicates that only 1 out of 130 adults at the school exhibited any symptoms.\nIt is reasonable, then, to suppose that transmission in each case ceased once he had been admitted to the infirmary.\nSupposing admission happened within 24 hours of the onset of symptoms.\nHow does this affect our estimate of $R_0$?\nTwelve hours?\n\n::: exercise\n### Exercise 3\n:::\n\nBiweekly data for outbreaks of measles in three communities in Niamey, Niger are provided in the dataframe `niamey`.\nUse this method to obtain estimates of $R_0$ for measles from the first community assuming that the infectious period is approximately two weeks or $\\frac{14}{365} \\approx 0.0384$ years.\n\n::: exercise\n### Exercise 4\n:::\n\nA defect with this method is that it uses only a small fraction of the information that might be available, i.e., the first few data points.\nIndeed, there is nothing in the method that tells one how many data points to use--this is a matter of judgment.\nFurther, there is a tradeoff in that as more and more data points are used the precision of the estimate increases, but this comes at a cost of additional bias.\nPlot the estimate of $R_0$ obtained from $n=3, 4, 5, ...$ data points against the standard error of the slope from the regression analysis to show this tradeoff.\n\n## Estimating dynamical parameters with least squares\n\nThe objective of the previous exercise was to estimate $R_0$.\nKnowing $R_0$ is critical to understanding the dynamics of any epidemic system.\nIt is, however, a composite quantity and is not sufficient to completely describe the epidemic trajectory.\nFor this, we require estimates for all parameters of the model.\nIn this exercise, we introduce a simple approach to model estimation called **least squares fitting**, sometimes called **trajectory matching**.\nThe basic idea is that we find the values of the model parameters that minimize the squared differences between model predictions and the observed data.\nTo demonstrate least squares fitting, we consider an outbreak of measles in Niamey, Niger, reported on by [@graisEstimatingTransmissionIntensity2006].\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Replace an \"NA\"\nniamey[5, 3] <- 0\n\nniamey_df <- niamey %>%\n # Rename columns to remove automatic \"V1\" etc columns names\n rename_with(., ~ paste0(\"Site_\", str_remove(.x, \"V\"))) %>%\n # Add a column for the biweekly time period\n mutate(biweek = 1:16) %>%\n # Convert to long format for plotting\n pivot_longer(\n cols = contains(\"Site\"),\n names_to = \"site\",\n values_to = \"cases\"\n )\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Create a vector of colors for each site in the Niamey dataset\nniamey_site_colors <- RColorBrewer::brewer.pal(3, \"Dark2\")\n# Assign names to the colors\nnames(niamey_site_colors) <- unique(niamey_df$site)\n\n# Create a vector of labels for each site for nicer plotting legends\nniamey_site_labels <- str_replace_all(names(niamey_site_colors), \"_\", \" \")\nnames(niamey_site_labels) <- names(niamey_site_colors)\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nggplot(\n niamey_df,\n aes(x = biweek, y = cases, color = site, fill = site, group = site)\n) +\n geom_line() +\n geom_point(shape = 21, size = 5, alpha = 0.8) +\n scale_color_manual(\n values = niamey_site_colors,\n aesthetics = c(\"color\", \"fill\"),\n labels = niamey_site_labels\n ) +\n guides(color = \"none\") +\n labs(x = \"Biweek\", y = \"Number of Cases\", fill = \"Site Number\") +\n theme(legend.position = \"bottom\")\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\n## Dynamical Model\n\nFirst, we write a specialized function for simulating the SIR model in a case where the removal rate is *\"hard-wired\"* and with no demography.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n#' Basic SIR model\n#'\n#' A basic SIR model with no demographic structure to be used in deSolve\n#'\n#' @param time deSolve passes the time parameter to the function.\n#' @param state A vector of states.\n#' @param params The beta parameter\n#' @param ... Other arguments passed by deSolve.\n#'\n#' @return A deSolve matrix of states at each time step.\n#' @examples\n#' sir_params <- 0.0005\n#' sir_init_states <- c(S = 5000, I = 1, R = 0)\n#' sim_times <- seq(0, 16 / 365, by = 0.1 / 365)\n#'\n#' sir_sol <- deSolve::ode(\n#' y = sir_init_states,\n#' times = sim_times,\n#' func = closed_sir_model,\n#' parms = sir_params\n#' ))\nclosed_sir_model <- function(time, state, params, ...) {\n # Unpack states\n S <- state[\"S\"]\n I <- state[\"I\"]\n\n # Unpack parameters\n beta <- params\n dur_inf <- 14 / 365\n gamma <- 1 / dur_inf\n\n new_inf <- beta * S * I\n\n # Calculate the ODEs\n dSdt <- -new_inf\n dIdt <- new_inf - (gamma * I)\n\n # Return the ODEs\n return(list(c(dSdt, dIdt, new_inf)))\n}\n```\n:::\n\n\n## Interactive Optimization\n\n\n::: {.cell}\n\n:::\n\n\n\n::: {.cell}\n\n:::\n\n\n```{ojs}\n//| echo: false\nfiltered_niamey_data = aq.table(niamey_data)\n .filter(aq.escape(d => d.site == site_select))\n```\n\n```{ojs}\n//| echo: false\nreset_S = 10000\nreset_I = 20\nreset_beta = 5.00\n```\n\n```{ojs}\n//| echo: false\nfunction set(input, value) {\n input.value = value;\n input.dispatchEvent(new Event(\"input\", {bubbles: true}));\n}\n```\n\n```{ojs}\n//| echo: false\nfunction sse(obs, preds) {\n\n if(obs.length == preds.length) {\n var squared_errs = obs.map((e, i) => (e - preds[i])**2 )\n return squared_errs.reduce((a, b) => a + b, 0)\n } else {\n return(\"lengths are not the same\")\n }\n}\n```\n\n```{ojs}\n//| panel: sidebar\n//| echo: false\nviewof reset = Inputs.button([\n [\"Reset all sliders\", () => {\n set(viewof S0, reset_S)\n set(viewof I0, reset_I)\n set(viewof beta_input, reset_beta)\n }]\n])\nviewof S0 = Inputs.range(\n [500, 15000],\n {value: reset_S, step: 1, label: md`${tex`S(0)`}`}\n)\n\nviewof I0 = Inputs.range(\n [0.001, 50],\n {value: reset_I, step: 0.001, label: md`${tex`I(0)`}`}\n)\n\nviewof beta_input = Inputs.range(\n [1, 100],\n {value: reset_beta, step: 0.001, label: md`${tex`\\beta (\\times 10^{-3})`}`}\n)\n\nviewof site_select = Inputs.select(\n [\"Site 1\", \"Site 2\", \"Site 3\"],\n {label: \"Select a site:\"}\n)\n\n// convert to daily time scale as easier to manipulate\nbeta = (beta_input / 365) * (10 ** (-3))\n\nmd`${tex`R_0 = ${R0_str}`}`\n```\n\n```{ojs}\n//| echo: false\ndur_inf = 14\ngamma = 1 / dur_inf\nR0 = beta * (S0 + I0)/ gamma\n\nR0_str = R0.toPrecision(2).toLocaleString()\n```\n\n```{ojs}\n//| echo: false\ndt = 0.01\ntmax = 16 * 14\n```\n\n```{ojs}\n//| echo: false\nimport {odeRK4} from '@rreusser/integration@3064'\nimport { aq, op } from '@uwdata/arquero'\n```\n\n```{ojs}\n//| echo: false\nfunction sir(dydt, y, t) {\n dydt[3] = beta * y[0] * y[1]\n\n dydt[0] = - dydt[3]\n dydt[1] = dydt[3] - gamma * y[1]\n dydt[2] = gamma * y[1]\n}\n```\n\n```{ojs}\n//| echo: false\nfunction simulate(f, t0, y0, dt, tmax) {\n var t = t0\n var y = y0\n var i = 0\n\n var ysim = [y0]\n\n for (t = t0 + dt; t <= tmax; t += dt) {\n ysim.push(odeRK4([], ysim[i], f, dt))\n i += 1\n }\n\n // return cumulative infections\n return ysim.map(d => d[3])\n}\n```\n\n```{ojs}\n//| echo: false\nsir_sol = simulate(sir, 0, [S0, I0, 0.0, 0.0], dt, tmax)\n```\n\n```{ojs}\n//| echo: false\nsiteColors = [\"#1b9e77\", \"#d95f02\", \"#7570b3\"]\n```\n\n```{ojs}\n//| echo: false\n// Create an array of time indices, dropping first as from returns index 0, which\n// doesn't exist in the real data\ntimes = Array.from({length: 17}, (_, i) => i * 14).slice(1)\ntindex = times.map((e, i) => e * (1 / dt))\n\ncum_inc = tindex.map((i) => sir_sol[i])\n\npreds = [cum_inc[0] + I0, ...cum_inc.map((e, i) => cum_inc[i] - cum_inc[i-1]).slice(1)]\n\nsir_tbl = aq.table({\n Biweek: times.map(t => t / 14),\n \"Cumulative Incidence\": cum_inc,\n \"Number of Individuals\": preds\n})\n\nsim_sse = [({\n sse: sse(\n filtered_niamey_data.array(\"Number of Individuals\"),\n preds\n ).toPrecision(4),\n Biweek: 3,\n \"Number of Individuals\": Math.max(\n ...sir_tbl.array(\"Number of Individuals\"),\n ...filtered_niamey_data.array(\"Number of Individuals\")\n ) * 0.9\n})]\n```\n\n```{ojs}\n//| echo: false\n//| panel: fill\nPlot.plot({\n color: {\n legend: true,\n domain: [\"Site 1\", \"Site 2\", \"Site 3\"],\n range: siteColors,\n },\n style: {fontSize: \"20px\"},\n marginLeft: 75,\n marginTop: 40,\n marginBottom: 55,\n grid: true,\n width: 800,\n height: 670,\n marks: [\n Plot.lineY(\n sir_tbl,\n {x: \"Biweek\", y: \"Number of Individuals\", stroke: \"#4d4d4dff\", strokeWidth: 6, strokeOpacity: 0.8}\n ),\n Plot.dot(\n filtered_niamey_data,\n {x: \"Biweek\", y: \"Number of Individuals\", stroke: \"site\", fill: \"site\", fillOpacity: 0.6, r: 12}\n ),\n Plot.lineY(\n filtered_niamey_data,\n {x: \"Biweek\", y: \"Number of Individuals\", stroke: \"site\"}\n ),\n Plot.text(\n sim_sse,\n {x: \"Biweek\", y: \"Number of Individuals\", text: (d) => `SSE = ${d.sse}`, dx: 0, dy: 0, fontWeight: \"bold\"}\n )\n ]\n})\n```\n\n## Objective Function\n\nNow we set up a function that will calculate the sum of the squared differences between the observations and the model at any parameterization (more commonly known as *\"sum of squared errors\"*).\nIn general, this is called the **objective function** because it is the quantity that optimization seeks to minimize.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n#' Calculate the Sum of Squared Errors\n#'\n#' A function to take in biweekly incidence data, and SIR parameters, and\n#' calculate the SSE\n#'\n#' @param params A vector of parameter values\n#' @param data A dataframe containing biweekly incidence data in the case column\n#'\n#' @return The SSE of type double\n#' @examples\nsse_sir <- function(params, data) {\n # Convert biweekly time series into annual time scale\n # Daily time scale has requires beta values to be too small - doesn't\n # optimize well\n dt <- 0.01\n max_biweek <- max(data$biweek)\n t <- seq(0, max_biweek * 14, dt) / 365\n\n # Extract the number of observed incidence\n obs_inc <- data$cases\n\n # Note the parameters are updated throughout the optimization process by\n # the optim() function\n # Unpack the transmission parameter and exponentiate to fit on ln scale\n beta <- exp(params[[\"beta\"]])\n\n # Unpack the initial states and exponentiate to fit on normal scale\n S_init <- exp(params[[\"S_init\"]])\n I_init <- exp(params[[\"I_init\"]])\n\n # Fit SIR model to the parameters\n sol <- deSolve::ode(\n y = c(S = S_init, I = I_init, new_inf = 0),\n times = t,\n func = closed_sir_model,\n parms = beta,\n # Use rk4 as fixed time steps, which is important for indexing\n method = \"rk4\"\n )\n\n # Extract the cumulative incidence\n cum_inc <- sol[, \"new_inf\"]\n\n # Find the indices of the cumulative incidence to extract\n biweek_index <- seq(1, max_biweek) * (14 / dt) + 1\n\n # Index cumulative incidence to get the values at the end of the biweeks\n biweek_cum_inc <- cum_inc[biweek_index]\n\n # Calculate the biweekly incidence by using the difference between\n # consecutive biweeks. Need to manually prepend first week's incidence\n # and add in the initial number of infectious individuals, as ODE model\n # only returns the cumulative differences, which is 0 at the start.\n biweek_inc <- c(biweek_cum_inc[1] + I_init, diff(biweek_cum_inc, lag = 1))\n\n # return SSE of predicted vs observed incidence\n return(sum((biweek_inc - obs_inc)^2))\n}\n```\n:::\n\n\nNotice that the code for `sse_sir()` makes use of the following modeling trick.\nWe know that $\\beta$, $S_0$, and $I_0$ must be positive, but our search to optimize these parameters will be over the entire number line.\nWe could constrain the search using a more sophisticated algorithm, but this might introduce other problems (i.e., stability at the boundaries).\nInstead, we parameterize our objective function (`sse_sir`) in terms of some alternative variables $\\ln(\\beta)$, $\\ln(S_0)$, and $\\ln(I_0)$.\nWhile these numbers range from $-\\infty$ to $\\infty$ (the range of our search) they map to our model parameters on a range from $0$ to $\\infty$ (the range that is biologically meaningful).\n\n## Optimization\n\nOur final step is to use the function `optim` to find the values of $\\beta$, $S_0$, and $I_0$ that minimize the sum of squared errors as calculated using our function.\n\nFinally, we plot these fits against the data.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Initial guess\nsse_optim_params <- c(beta = log(0.055), S_init = log(5000), I_init = log(1))\n\n# Create a dataframe of optimized parameters\nniamey_optims <- niamey_df %>%\n # Create a nested dataframe i.e. one row for each site, and the data column\n # now is a list column that contains a separate dataframe of times and\n # cases for each site\n nest(data = -site) %>%\n mutate(\n # Map the optim() function call to each of the separate dataframes\n # stored in the nested data column we just created\n fit = map(data, ~ optim(sse_optim_params, sse_sir, data = .x)),\n # Map the exp() function to each of the model fits just created, and\n # output to a dataframe instead of a list (like in map()), for easier\n # use in the plottinge predictions later\n map_dfr(fit, ~ exp(.x$par))\n )\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nniamey_optims %>%\n select(-c(data, fit)) %>%\n mutate(site = str_replace_all(site, \"_\", \" \")) %>%\n gt() %>%\n fmt_number(columns = -site, decimals = 3) %>%\n fmt_scientific(columns = beta, decimals = 3) %>%\n # Relabel the column headers\n cols_label(\n site = md(\"**Site**\"),\n beta = md(\"**Beta**\"),\n S_init = md(\"**Initial S**\"),\n I_init = md(\"**Initial I**\")\n ) %>%\n # Apply style to the table with gray alternating rows\n opt_stylize(style = 1, color = 'gray') %>%\n # Increate space between columns\n opt_horizontal_padding(scale = 3) %>%\n cols_align(\"center\")\n```\n\n::: {.cell-output-display}\n\n```{=html}\n
Code
times = Array.from({length: 17}, (_, i) => i * 14).slice(1)
-tindex = times.map((e, i) => e * (1 / dt))
-
-cum_inc = tindex.map((i) => sir_sol[i])
-
-preds = [cum_inc[0] + I0, ...cum_inc.map((e, i) => cum_inc[i] - cum_inc[i-1]).slice(1)]
-
-sir_tbl = aq.table({
- Biweek: times.map(t => t / 14),
- "Cumulative Incidence": cum_inc,
- "Number of Individuals": preds
-})
+Code
times = Array.from({length: 17}, (_, i) => i * 14).slice(1)
+tindex = times.map((e, i) => e * (1 / dt))
+
+cum_inc = tindex.map((i) => sir_sol[i])
+
+preds = [cum_inc[0] + I0, ...cum_inc.map((e, i) => cum_inc[i] - cum_inc[i-1]).slice(1)]
-sim_sse = [({
- sse: sse(
- filtered_niamey_data.array("Number of Individuals"),
- preds
- ).toPrecision(4),
- Biweek: 3,
- "Number of Individuals": Math.max(
- ...sir_tbl.array("Number of Individuals"),
- ...filtered_niamey_data.array("Number of Individuals")
- ) * 0.9
-})]
+sir_tbl = aq.table({
+ Biweek: times.map(t => t / 14),
+ "Cumulative Incidence": cum_inc,
+ "Number of Individuals": preds
+})
+
+sim_sse = [({
+ sse: sse(
+ filtered_niamey_data.array("Number of Individuals"),
+ preds
+ ).toPrecision(4),
+ Biweek: 3,
+ "Number of Individuals": Math.max(
+ ...sir_tbl.array("Number of Individuals"),
+ ...filtered_niamey_data.array("Number of Individuals")
+ ) * 0.9
+})]
@@ -1193,38 +1194,38 @@
-
+
Code
Plot.plot({
- color: {
- legend: true,
- domain: ["Site 1", "Site 2", "Site 3"],
- range: siteColors,
- },
- style: {fontSize: "20px"},
- marginLeft: 75,
- marginTop: 40,
- marginBottom: 55,
- grid: true,
- width: 800,
- height: 670,
- marks: [
- Plot.lineY(
- sir_tbl,
- {x: "Biweek", y: "Number of Individuals", stroke: "#4d4d4dff", strokeWidth: 6, strokeOpacity: 0.8}
- ),
- Plot.dot(
- filtered_niamey_data,
- {x: "Biweek", y: "Number of Individuals", stroke: "site", fill: "site", fillOpacity: 0.6, r: 12}
- ),
- Plot.lineY(
- filtered_niamey_data,
- {x: "Biweek", y: "Number of Individuals", stroke: "site"}
- ),
- Plot.text(
- sim_sse,
- {x: "Biweek", y: "Number of Individuals", text: (d) => `SSE = ${d.sse}`, dx: 0, dy: 0, fontWeight: "bold"}
- )
- ]
-})Code
Plot.plot({
+ color: {
+ legend: true,
+ domain: ["Site 1", "Site 2", "Site 3"],
+ range: siteColors,
+ },
+ style: {fontSize: "20px"},
+ marginLeft: 75,
+ marginTop: 40,
+ marginBottom: 55,
+ grid: true,
+ width: 800,
+ height: 670,
+ marks: [
+ Plot.lineY(
+ sir_tbl,
+ {x: "Biweek", y: "Number of Individuals", stroke: "#4d4d4dff", strokeWidth: 6, strokeOpacity: 0.8}
+ ),
+ Plot.dot(
+ filtered_niamey_data,
+ {x: "Biweek", y: "Number of Individuals", stroke: "site", fill: "site", fillOpacity: 0.6, r: 12}
+ ),
+ Plot.lineY(
+ filtered_niamey_data,
+ {x: "Biweek", y: "Number of Individuals", stroke: "site"}
+ ),
+ Plot.text(
+ sim_sse,
+ {x: "Biweek", y: "Number of Individuals", text: (d) => `SSE = ${d.sse}`, dx: 0, dy: 0, fontWeight: "bold"}
+ )
+ ]
+})
-
-
@@ -3452,1517 +3453,1523 @@
labs(x = "Day", y = "Active Influenza Cases")
```
-## Estimating $R_0$ From The Final Outbreak Size
-
-Our first approach is to estimate $R_0$ from the final outbreak size.
-Although unhelpful at the early stages of an epidemic (before the final epidemic size is observed), this method is nonetheless a useful tool for *post hoc* analysis.
-The method is general and can be motivated by the argument listed in [@keelingIntroductionSimpleEpidemic2008]:
+The figure above, and the corresponding tabulated outbreak data, report the number of **active influenza cases** on each day.
+These are point-in-time counts, or prevalence: the number of students who were sick at that moment.
+By contrast, the final-size calculation below uses cumulative **incidence**: the total number of students who became cases at any point during the outbreak.
+For this boarding school outbreak, the total school size was $N = 764$, and the total number of students who became cases over the full outbreak was $Z = 512$.
+Because $512$ is a cumulative incidence count, it is not obtained by reading a single value from the prevalence figure or table.
-First, we assume that the epidemic is started by a single infectious individual in a completely susceptible population.
-On average, this individual infects $R_0$ others.
-The probability a particular individual escaped infection is therefore $e^{-R_0 / N}$.
-
-If $Z$ individuals have been infected, the probability of an individual escaping infection from all potential sources is $e^{-Z R_0 / N}$.
-It follows that at the end of the epidemic a proportion $R(\infty) = Z / N$ have been infected and the fraction remaining susceptible is $S(\infty) = e^{-R(\infty) R_0}$, which is equal to $1 - R(\infty)$.
-
-::: {.callout-note collapse="true"}
-$S(\infty) = e^{-R(\infty) R_0}$ can be calculated by acknowledging that at equilibrium ($t = \infty$), $S(\infty) = 1 - R(\infty) = Z / N$, so substituting $R(\infty)$ into $1 - e^{-Z R_0 / N}$ gives the desired result.
+## Estimating $R_0$ From The Final Outbreak Size
+
+Our first approach is to estimate $R_0$ from the final outbreak size.
+Although unhelpful at the early stages of an epidemic (before the final epidemic size is observed), this method is nonetheless a useful tool for *post hoc* analysis.
+The method is general and can be motivated by the argument listed in [@keelingIntroductionSimpleEpidemic2008]:
+
+First, we assume that the epidemic is started by a single infectious individual in a completely susceptible population.
+On average, this individual infects $R_0$ others.
+The probability a particular individual escaped infection is therefore $e^{-R_0 / N}$.
-It could also be calculated by dividing $\frac{\dd{S}}{\dd{t}}$ by $\frac{\dd{R}}{\dd{t}}$:
-
-```{=tex}
-\begin{aligned}
-\frac{\dd{S}}{\dd{R}} &= - \frac{\beta S}{\gamma} \\
-&= - R_0 S
-\end{aligned}
-```
-which is a [separable differential equation](https://tutorial.math.lamar.edu/classes/de/separable.aspx), so can be integrated as follows:
-
-```{=tex}
-\begin{aligned}
-- \int_{0}^{t} \frac{1}{R_0 S} \dd{S} &= \int_{0}^{t} \dd{R} \\
-- \frac{1}{R_0} \left(\ln{S(t)} - \ln{S(0)} \right) &= R(t) - \cancelto{0}{R(0)} \\
-\ln{S(t)} &= \ln{S(0)} - R_0 R(t) \\
-S(t) &= S(0) e^{-R_0 R(t)}
-
-\end{aligned}
-```
-:::
-
-Putting this together, we get:
+If $Z$ individuals have been infected, the probability of an individual escaping infection from all potential sources is $e^{-Z R_0 / N}$.
+It follows that at the end of the epidemic a proportion $R(\infty) = Z / N$ have been infected and the fraction remaining susceptible is $S(\infty) = e^{-R(\infty) R_0}$, which is equal to $1 - R(\infty)$.
+
+::: {.callout-note collapse="true"}
+$S(\infty) = e^{-R(\infty) R_0}$ can be calculated by acknowledging that at equilibrium ($t = \infty$), $S(\infty) = 1 - R(\infty) = Z / N$, so substituting $R(\infty)$ into $1 - e^{-Z R_0 / N}$ gives the desired result.
+
+It could also be calculated by dividing $\frac{\dd{S}}{\dd{t}}$ by $\frac{\dd{R}}{\dd{t}}$:
+
+```{=tex}
+\begin{aligned}
+\frac{\dd{S}}{\dd{R}} &= - \frac{\beta S}{\gamma} \\
+&= - R_0 S
+\end{aligned}
+```
+which is a [separable differential equation](https://tutorial.math.lamar.edu/classes/de/separable.aspx), so can be integrated as follows:
+
+```{=tex}
+\begin{aligned}
+- \int_{0}^{t} \frac{1}{R_0 S} \dd{S} &= \int_{0}^{t} \dd{R} \\
+- \frac{1}{R_0} \left(\ln{S(t)} - \ln{S(0)} \right) &= R(t) - \cancelto{0}{R(0)} \\
+\ln{S(t)} &= \ln{S(0)} - R_0 R(t) \\
+S(t) &= S(0) e^{-R_0 R(t)}
-$$
-1 - R(\infty) - e^{-R(\infty) R_0} = 0
-$$
+\end{aligned}
+```
+:::
-Rearranging, we have the estimator
+Putting this together, we get:
$$
- \hat{R_0} = \frac{\ln(1 - Z / N)}{-Z / N},
+1 - R(\infty) - e^{-R(\infty) R_0} = 0
$$
-which, in this case, evaluates to $\frac{\ln(1 - 512 / 764)}{-512 / 764} = 1.655$.
+Rearranging, we have the estimator
-::: exercise
-### Exercise 1
-:::
+$$
+ \hat{R_0} = \frac{\ln(1 - Z / N)}{-Z / N},
+$$
-This equation shows the important one-to-one relationship between $R_0$ and the final epidemic size.
-Plot the relationship between the total epidemic size and $R_0$ for the complete range of values between 0 and 1.
-
-## Linear Approximation
-
-The next method we introduce takes advantage of the fact that during the early stages of an outbreak, the number of infected individuals is given approximately as $I(t) \approx I_0 e^{((R_0 - 1)(\gamma + \mu)t)}$.
-Taking logarithms of both sides, we have $\ln(I(t)) \approx \ln(I_0) + (R_0 - 1)(\gamma + \mu)t$, showing that the log of the number of infected individuals is approximately linear in time with a slope that reflects both $R_0$ and the recovery rate.
-
-This suggests that a simple linear regression fit to the first several data points on a log-scale, corrected to account for $\gamma$ and $\mu$, provides a rough and ready estimate of $R_0$.
-For flu, we can assume $\mu =0$ because the epidemic occurred over a time period during which natural mortality is negligible.
-Further, assuming an infectious period of about 2.5 days, we use $\gamma = (2.5)^{-1} = 0.4$ for the correction.
-Fitting to the first four data points, we obtain the slope as follows.
-
-```{r}
-# Fit a linear model
-linear_model <- lm(log(flu[1:4]) ~ day[1:4], data = flu)
+which, in this case, evaluates to $\frac{\ln(1 - 512 / 764)}{-512 / 764} = 1.655$.
+
+::: exercise
+### Exercise 1
+:::
+
+This equation shows the important one-to-one relationship between $R_0$ and the final epidemic size.
+Plot the relationship between the total epidemic size and $R_0$ for the complete range of values between 0 and 1.
+
+## Linear Approximation
+
+The next method we introduce takes advantage of the fact that during the early stages of an outbreak, the number of infected individuals is given approximately as $I(t) \approx I_0 e^{((R_0 - 1)(\gamma + \mu)t)}$.
+Taking logarithms of both sides, we have $\ln(I(t)) \approx \ln(I_0) + (R_0 - 1)(\gamma + \mu)t$, showing that the log of the number of infected individuals is approximately linear in time with a slope that reflects both $R_0$ and the recovery rate.
+
+This suggests that a simple linear regression fit to the first several data points on a log-scale, corrected to account for $\gamma$ and $\mu$, provides a rough and ready estimate of $R_0$.
+For flu, we can assume $\mu =0$ because the epidemic occurred over a time period during which natural mortality is negligible.
+Further, assuming an infectious period of about 2.5 days, we use $\gamma = (2.5)^{-1} = 0.4$ for the correction.
+Fitting to the first four data points, we obtain the slope as follows.
-# Summary statistics for fit model
-summary(linear_model)
-
-# Extract slope parameter
-coef(linear_model)[2]
-```
-
-Rearranging the linear equation above and denoting the slope coefficient by $\hat \beta_1$ we have the estimator $\hat R_0 = \hat \beta_1 / \gamma + 1$ giving $\hat R_0 = 1.094913 / 0.4 + 1 \approx 3.7$.
+```{r}
+# Fit a linear model
+linear_model <- lm(log(flu[1:4]) ~ day[1:4], data = flu)
+
+# Summary statistics for fit model
+summary(linear_model)
-::: exercise
-### Exercise 2
-:::
+# Extract slope parameter
+coef(linear_model)[2]
+```
-Our estimate assumes that boys remained infectious during the natural course of infection.
-The original report on this epidemic indicates that boys found to have symptoms were immediately confined to bed in the infirmary.
-The report also indicates that only 1 out of 130 adults at the school exhibited any symptoms.
-It is reasonable, then, to suppose that transmission in each case ceased once he had been admitted to the infirmary.
-Supposing admission happened within 24 hours of the onset of symptoms.
-How does this affect our estimate of $R_0$?
-Twelve hours?
-
-::: exercise
-### Exercise 3
-:::
-
-Biweekly data for outbreaks of measles in three communities in Niamey, Niger are provided in the dataframe `niamey`.
-Use this method to obtain estimates of $R_0$ for measles from the first community assuming that the infectious period is approximately two weeks or $\frac{14}{365} \approx 0.0384$ years.
-
-::: exercise
-### Exercise 4
-:::
-
-A defect with this method is that it uses only a small fraction of the information that might be available, i.e., the first few data points.
-Indeed, there is nothing in the method that tells one how many data points to use--this is a matter of judgment.
-Further, there is a tradeoff in that as more and more data points are used the precision of the estimate increases, but this comes at a cost of additional bias.
-Plot the estimate of $R_0$ obtained from $n=3, 4, 5, ...$ data points against the standard error of the slope from the regression analysis to show this tradeoff.
-
-## Estimating dynamical parameters with least squares
-
-The objective of the previous exercise was to estimate $R_0$.
-Knowing $R_0$ is critical to understanding the dynamics of any epidemic system.
-It is, however, a composite quantity and is not sufficient to completely describe the epidemic trajectory.
-For this, we require estimates for all parameters of the model.
-In this exercise, we introduce a simple approach to model estimation called **least squares fitting**, sometimes called **trajectory matching**.
-The basic idea is that we find the values of the model parameters that minimize the squared differences between model predictions and the observed data.
-To demonstrate least squares fitting, we consider an outbreak of measles in Niamey, Niger, reported on by [@graisEstimatingTransmissionIntensity2006].
-
-```{r}
-# Replace an "NA"
-niamey[5, 3] <- 0
+Rearranging the linear equation above and denoting the slope coefficient by $\hat \beta_1$ we have the estimator $\hat R_0 = \hat \beta_1 / \gamma + 1$ giving $\hat R_0 = 1.094913 / 0.4 + 1 \approx 3.7$.
+
+::: exercise
+### Exercise 2
+:::
+
+Our estimate assumes that boys remained infectious during the natural course of infection.
+The original report on this epidemic indicates that boys found to have symptoms were immediately confined to bed in the infirmary.
+The report also indicates that only 1 out of 130 adults at the school exhibited any symptoms.
+It is reasonable, then, to suppose that transmission in each case ceased once he had been admitted to the infirmary.
+Supposing admission happened within 24 hours of the onset of symptoms.
+How does this affect our estimate of $R_0$?
+Twelve hours?
+
+::: exercise
+### Exercise 3
+:::
+
+Biweekly data for outbreaks of measles in three communities in Niamey, Niger are provided in the dataframe `niamey`.
+Use this method to obtain estimates of $R_0$ for measles from the first community assuming that the infectious period is approximately two weeks or $\frac{14}{365} \approx 0.0384$ years.
+
+::: exercise
+### Exercise 4
+:::
+
+A defect with this method is that it uses only a small fraction of the information that might be available, i.e., the first few data points.
+Indeed, there is nothing in the method that tells one how many data points to use--this is a matter of judgment.
+Further, there is a tradeoff in that as more and more data points are used the precision of the estimate increases, but this comes at a cost of additional bias.
+Plot the estimate of $R_0$ obtained from $n=3, 4, 5, ...$ data points against the standard error of the slope from the regression analysis to show this tradeoff.
+
+## Estimating dynamical parameters with least squares
+
+The objective of the previous exercise was to estimate $R_0$.
+Knowing $R_0$ is critical to understanding the dynamics of any epidemic system.
+It is, however, a composite quantity and is not sufficient to completely describe the epidemic trajectory.
+For this, we require estimates for all parameters of the model.
+In this exercise, we introduce a simple approach to model estimation called **least squares fitting**, sometimes called **trajectory matching**.
+The basic idea is that we find the values of the model parameters that minimize the squared differences between model predictions and the observed data.
+To demonstrate least squares fitting, we consider an outbreak of measles in Niamey, Niger, reported on by [@graisEstimatingTransmissionIntensity2006].
-niamey_df <- niamey %>%
- # Rename columns to remove automatic "V1" etc columns names
- rename_with(., ~ paste0("Site_", str_remove(.x, "V"))) %>%
- # Add a column for the biweekly time period
- mutate(biweek = 1:16) %>%
- # Convert to long format for plotting
- pivot_longer(
- cols = contains("Site"),
- names_to = "site",
- values_to = "cases"
- )
-```
-
-```{r}
-# Create a vector of colors for each site in the Niamey dataset
-niamey_site_colors <- RColorBrewer::brewer.pal(3, "Dark2")
-# Assign names to the colors
-names(niamey_site_colors) <- unique(niamey_df$site)
-
-# Create a vector of labels for each site for nicer plotting legends
-niamey_site_labels <- str_replace_all(names(niamey_site_colors), "_", " ")
-names(niamey_site_labels) <- names(niamey_site_colors)
-```
-
-```{r}
-#| column: body
-#| out-width: 100%
-ggplot(
- niamey_df,
- aes(x = biweek, y = cases, color = site, fill = site, group = site)
-) +
- geom_line() +
- geom_point(shape = 21, size = 5, alpha = 0.8) +
- scale_color_manual(
- values = niamey_site_colors,
- aesthetics = c("color", "fill"),
- labels = niamey_site_labels
- ) +
- guides(color = "none") +
- labs(x = "Biweek", y = "Number of Cases", fill = "Site Number") +
- theme(legend.position = "bottom")
-```
-
-## Dynamical Model
-
-First, we write a specialized function for simulating the SIR model in a case where the removal rate is *"hard-wired"* and with no demography.
-
-```{r}
-#' Basic SIR model
-#'
-#' A basic SIR model with no demographic structure to be used in deSolve
-#'
-#' @param time deSolve passes the time parameter to the function.
-#' @param state A vector of states.
-#' @param params The beta parameter
-#' @param ... Other arguments passed by deSolve.
-#'
-#' @return A deSolve matrix of states at each time step.
-#' @examples
-#' sir_params <- 0.0005
-#' sir_init_states <- c(S = 5000, I = 1, R = 0)
-#' sim_times <- seq(0, 16 / 365, by = 0.1 / 365)
+```{r}
+# Replace an "NA"
+niamey[5, 3] <- 0
+
+niamey_df <- niamey %>%
+ # Rename columns to remove automatic "V1" etc columns names
+ rename_with(., ~ paste0("Site_", str_remove(.x, "V"))) %>%
+ # Add a column for the biweekly time period
+ mutate(biweek = 1:16) %>%
+ # Convert to long format for plotting
+ pivot_longer(
+ cols = contains("Site"),
+ names_to = "site",
+ values_to = "cases"
+ )
+```
+
+```{r}
+# Create a vector of colors for each site in the Niamey dataset
+niamey_site_colors <- RColorBrewer::brewer.pal(3, "Dark2")
+# Assign names to the colors
+names(niamey_site_colors) <- unique(niamey_df$site)
+
+# Create a vector of labels for each site for nicer plotting legends
+niamey_site_labels <- str_replace_all(names(niamey_site_colors), "_", " ")
+names(niamey_site_labels) <- names(niamey_site_colors)
+```
+
+```{r}
+#| column: body
+#| out-width: 100%
+ggplot(
+ niamey_df,
+ aes(x = biweek, y = cases, color = site, fill = site, group = site)
+) +
+ geom_line() +
+ geom_point(shape = 21, size = 5, alpha = 0.8) +
+ scale_color_manual(
+ values = niamey_site_colors,
+ aesthetics = c("color", "fill"),
+ labels = niamey_site_labels
+ ) +
+ guides(color = "none") +
+ labs(x = "Biweek", y = "Number of Cases", fill = "Site Number") +
+ theme(legend.position = "bottom")
+```
+
+## Dynamical Model
+
+First, we write a specialized function for simulating the SIR model in a case where the removal rate is *"hard-wired"* and with no demography.
+
+```{r}
+#' Basic SIR model
+#'
+#' A basic SIR model with no demographic structure to be used in deSolve
+#'
+#' @param time deSolve passes the time parameter to the function.
+#' @param state A vector of states.
+#' @param params The beta parameter
+#' @param ... Other arguments passed by deSolve.
#'
-#' sir_sol <- deSolve::ode(
-#' y = sir_init_states,
-#' times = sim_times,
-#' func = closed_sir_model,
-#' parms = sir_params
-#' ))
-closed_sir_model <- function(time, state, params, ...) {
- # Unpack states
- S <- state["S"]
- I <- state["I"]
-
- # Unpack parameters
- beta <- params
- dur_inf <- 14 / 365
- gamma <- 1 / dur_inf
-
- new_inf <- beta * S * I
-
- # Calculate the ODEs
- dSdt <- -new_inf
- dIdt <- new_inf - (gamma * I)
+#' @return A deSolve matrix of states at each time step.
+#' @examples
+#' sir_params <- 0.0005
+#' sir_init_states <- c(S = 5000, I = 1, R = 0)
+#' sim_times <- seq(0, 16 / 365, by = 0.1 / 365)
+#'
+#' sir_sol <- deSolve::ode(
+#' y = sir_init_states,
+#' times = sim_times,
+#' func = closed_sir_model,
+#' parms = sir_params
+#' ))
+closed_sir_model <- function(time, state, params, ...) {
+ # Unpack states
+ S <- state["S"]
+ I <- state["I"]
+
+ # Unpack parameters
+ beta <- params
+ dur_inf <- 14 / 365
+ gamma <- 1 / dur_inf
- # Return the ODEs
- return(list(c(dSdt, dIdt, new_inf)))
-}
-```
-
-## Interactive Optimization
-
-```{r}
-#| echo: false
-ojs_niamey <- niamey_df %>%
- rename("Biweek" = biweek, "Number of Individuals" = cases) %>%
- mutate(site = str_replace_all(site, "_", " "))
-```
-
-```{r}
-#| echo: false
-#| cache: false
-ojs_define(niamey_data = ojs_niamey)
-```
-
-```{ojs}
-//| echo: false
-filtered_niamey_data = aq.table(niamey_data)
- .filter(aq.escape(d => d.site == site_select))
-```
-
-```{ojs}
-//| echo: false
-reset_S = 10000
-reset_I = 20
-reset_beta = 5.00
-```
-
-```{ojs}
-//| echo: false
-function set(input, value) {
- input.value = value;
- input.dispatchEvent(new Event("input", {bubbles: true}));
-}
-```
-
-```{ojs}
-//| echo: false
-function sse(obs, preds) {
+ new_inf <- beta * S * I
+
+ # Calculate the ODEs
+ dSdt <- -new_inf
+ dIdt <- new_inf - (gamma * I)
+
+ # Return the ODEs
+ return(list(c(dSdt, dIdt, new_inf)))
+}
+```
+
+## Interactive Optimization
+
+```{r}
+#| echo: false
+ojs_niamey <- niamey_df %>%
+ rename("Biweek" = biweek, "Number of Individuals" = cases) %>%
+ mutate(site = str_replace_all(site, "_", " "))
+```
+
+```{r}
+#| echo: false
+#| cache: false
+ojs_define(niamey_data = ojs_niamey)
+```
+
+```{ojs}
+//| echo: false
+filtered_niamey_data = aq.table(niamey_data)
+ .filter(aq.escape(d => d.site == site_select))
+```
+
+```{ojs}
+//| echo: false
+reset_S = 10000
+reset_I = 20
+reset_beta = 5.00
+```
+
+```{ojs}
+//| echo: false
+function set(input, value) {
+ input.value = value;
+ input.dispatchEvent(new Event("input", {bubbles: true}));
+}
+```
- if(obs.length == preds.length) {
- var squared_errs = obs.map((e, i) => (e - preds[i])**2 )
- return squared_errs.reduce((a, b) => a + b, 0)
- } else {
- return("lengths are not the same")
- }
-}
-```
-
-```{ojs}
-//| panel: sidebar
-//| echo: false
-viewof reset = Inputs.button([
- ["Reset all sliders", () => {
- set(viewof S0, reset_S)
- set(viewof I0, reset_I)
- set(viewof beta_input, reset_beta)
- }]
-])
-viewof S0 = Inputs.range(
- [500, 15000],
- {value: reset_S, step: 1, label: md`${tex`S(0)`}`}
-)
-
-viewof I0 = Inputs.range(
- [0.001, 50],
- {value: reset_I, step: 0.001, label: md`${tex`I(0)`}`}
-)
-
-viewof beta_input = Inputs.range(
- [1, 100],
- {value: reset_beta, step: 0.001, label: md`${tex`\beta (\times 10^{-3})`}`}
-)
-
-viewof site_select = Inputs.select(
- ["Site 1", "Site 2", "Site 3"],
- {label: "Select a site:"}
-)
-
-// convert to daily time scale as easier to manipulate
-beta = (beta_input / 365) * (10 ** (-3))
-
-md`${tex`R_0 = ${R0_str}`}`
-```
+```{ojs}
+//| echo: false
+function sse(obs, preds) {
+
+ if(obs.length == preds.length) {
+ var squared_errs = obs.map((e, i) => (e - preds[i])**2 )
+ return squared_errs.reduce((a, b) => a + b, 0)
+ } else {
+ return("lengths are not the same")
+ }
+}
+```
+
+```{ojs}
+//| panel: sidebar
+//| echo: false
+viewof reset = Inputs.button([
+ ["Reset all sliders", () => {
+ set(viewof S0, reset_S)
+ set(viewof I0, reset_I)
+ set(viewof beta_input, reset_beta)
+ }]
+])
+viewof S0 = Inputs.range(
+ [500, 15000],
+ {value: reset_S, step: 1, label: md`${tex`S(0)`}`}
+)
+
+viewof I0 = Inputs.range(
+ [0.001, 50],
+ {value: reset_I, step: 0.001, label: md`${tex`I(0)`}`}
+)
+
+viewof beta_input = Inputs.range(
+ [1, 100],
+ {value: reset_beta, step: 0.001, label: md`${tex`\beta (\times 10^{-3})`}`}
+)
+
+viewof site_select = Inputs.select(
+ ["Site 1", "Site 2", "Site 3"],
+ {label: "Select a site:"}
+)
-```{ojs}
-//| echo: false
-dur_inf = 14
-gamma = 1 / dur_inf
-R0 = beta * (S0 + I0)/ gamma
-
-R0_str = R0.toPrecision(2).toLocaleString()
-```
-
-```{ojs}
-//| echo: false
-dt = 0.01
-tmax = 16 * 14
-```
-
-```{ojs}
-//| echo: false
-import {odeRK4} from '@rreusser/integration@3064'
-import { aq, op } from '@uwdata/arquero'
-```
-
-```{ojs}
-//| echo: false
-function sir(dydt, y, t) {
- dydt[3] = beta * y[0] * y[1]
-
- dydt[0] = - dydt[3]
- dydt[1] = dydt[3] - gamma * y[1]
- dydt[2] = gamma * y[1]
-}
-```
+// convert to daily time scale as easier to manipulate
+beta = (beta_input / 365) * (10 ** (-3))
+
+md`${tex`R_0 = ${R0_str}`}`
+```
+
+```{ojs}
+//| echo: false
+dur_inf = 14
+gamma = 1 / dur_inf
+R0 = beta * (S0 + I0)/ gamma
+
+R0_str = R0.toPrecision(2).toLocaleString()
+```
+
+```{ojs}
+//| echo: false
+dt = 0.01
+tmax = 16 * 14
+```
+
+```{ojs}
+//| echo: false
+import {odeRK4} from '@rreusser/integration@3064'
+import { aq, op } from '@uwdata/arquero'
+```
+
+```{ojs}
+//| echo: false
+function sir(dydt, y, t) {
+ dydt[3] = beta * y[0] * y[1]
-```{ojs}
-//| echo: false
-function simulate(f, t0, y0, dt, tmax) {
- var t = t0
- var y = y0
- var i = 0
-
- var ysim = [y0]
-
- for (t = t0 + dt; t <= tmax; t += dt) {
- ysim.push(odeRK4([], ysim[i], f, dt))
- i += 1
- }
-
- // return cumulative infections
- return ysim.map(d => d[3])
-}
-```
-
-```{ojs}
-//| echo: false
-sir_sol = simulate(sir, 0, [S0, I0, 0.0, 0.0], dt, tmax)
-```
-
-```{ojs}
-//| echo: false
-siteColors = ["#1b9e77", "#d95f02", "#7570b3"]
-```
-
-```{ojs}
-//| echo: false
-// Create an array of time indices, dropping first as from returns index 0, which
-// doesn't exist in the real data
-times = Array.from({length: 17}, (_, i) => i * 14).slice(1)
-tindex = times.map((e, i) => e * (1 / dt))
-
-cum_inc = tindex.map((i) => sir_sol[i])
-
-preds = [cum_inc[0] + I0, ...cum_inc.map((e, i) => cum_inc[i] - cum_inc[i-1]).slice(1)]
-
-sir_tbl = aq.table({
- Biweek: times.map(t => t / 14),
- "Cumulative Incidence": cum_inc,
- "Number of Individuals": preds
-})
+ dydt[0] = - dydt[3]
+ dydt[1] = dydt[3] - gamma * y[1]
+ dydt[2] = gamma * y[1]
+}
+```
+
+```{ojs}
+//| echo: false
+function simulate(f, t0, y0, dt, tmax) {
+ var t = t0
+ var y = y0
+ var i = 0
+
+ var ysim = [y0]
+
+ for (t = t0 + dt; t <= tmax; t += dt) {
+ ysim.push(odeRK4([], ysim[i], f, dt))
+ i += 1
+ }
+
+ // return cumulative infections
+ return ysim.map(d => d[3])
+}
+```
+
+```{ojs}
+//| echo: false
+sir_sol = simulate(sir, 0, [S0, I0, 0.0, 0.0], dt, tmax)
+```
+
+```{ojs}
+//| echo: false
+siteColors = ["#1b9e77", "#d95f02", "#7570b3"]
+```
+
+```{ojs}
+//| echo: false
+// Create an array of time indices, dropping first as from returns index 0, which
+// doesn't exist in the real data
+times = Array.from({length: 17}, (_, i) => i * 14).slice(1)
+tindex = times.map((e, i) => e * (1 / dt))
+
+cum_inc = tindex.map((i) => sir_sol[i])
+
+preds = [cum_inc[0] + I0, ...cum_inc.map((e, i) => cum_inc[i] - cum_inc[i-1]).slice(1)]
-sim_sse = [({
- sse: sse(
- filtered_niamey_data.array("Number of Individuals"),
- preds
- ).toPrecision(4),
- Biweek: 3,
- "Number of Individuals": Math.max(
- ...sir_tbl.array("Number of Individuals"),
- ...filtered_niamey_data.array("Number of Individuals")
- ) * 0.9
-})]
-```
-
-```{ojs}
-//| echo: false
-//| panel: fill
-Plot.plot({
- color: {
- legend: true,
- domain: ["Site 1", "Site 2", "Site 3"],
- range: siteColors,
- },
- style: {fontSize: "20px"},
- marginLeft: 75,
- marginTop: 40,
- marginBottom: 55,
- grid: true,
- width: 800,
- height: 670,
- marks: [
- Plot.lineY(
- sir_tbl,
- {x: "Biweek", y: "Number of Individuals", stroke: "#4d4d4dff", strokeWidth: 6, strokeOpacity: 0.8}
- ),
- Plot.dot(
- filtered_niamey_data,
- {x: "Biweek", y: "Number of Individuals", stroke: "site", fill: "site", fillOpacity: 0.6, r: 12}
- ),
- Plot.lineY(
- filtered_niamey_data,
- {x: "Biweek", y: "Number of Individuals", stroke: "site"}
- ),
- Plot.text(
- sim_sse,
- {x: "Biweek", y: "Number of Individuals", text: (d) => `SSE = ${d.sse}`, dx: 0, dy: 0, fontWeight: "bold"}
- )
- ]
-})
-```
-
-## Objective Function
-
-Now we set up a function that will calculate the sum of the squared differences between the observations and the model at any parameterization (more commonly known as *"sum of squared errors"*).
-In general, this is called the **objective function** because it is the quantity that optimization seeks to minimize.
-
-```{r}
-#' Calculate the Sum of Squared Errors
-#'
-#' A function to take in biweekly incidence data, and SIR parameters, and
-#' calculate the SSE
-#'
-#' @param params A vector of parameter values
-#' @param data A dataframe containing biweekly incidence data in the case column
+sir_tbl = aq.table({
+ Biweek: times.map(t => t / 14),
+ "Cumulative Incidence": cum_inc,
+ "Number of Individuals": preds
+})
+
+sim_sse = [({
+ sse: sse(
+ filtered_niamey_data.array("Number of Individuals"),
+ preds
+ ).toPrecision(4),
+ Biweek: 3,
+ "Number of Individuals": Math.max(
+ ...sir_tbl.array("Number of Individuals"),
+ ...filtered_niamey_data.array("Number of Individuals")
+ ) * 0.9
+})]
+```
+
+```{ojs}
+//| echo: false
+//| panel: fill
+Plot.plot({
+ color: {
+ legend: true,
+ domain: ["Site 1", "Site 2", "Site 3"],
+ range: siteColors,
+ },
+ style: {fontSize: "20px"},
+ marginLeft: 75,
+ marginTop: 40,
+ marginBottom: 55,
+ grid: true,
+ width: 800,
+ height: 670,
+ marks: [
+ Plot.lineY(
+ sir_tbl,
+ {x: "Biweek", y: "Number of Individuals", stroke: "#4d4d4dff", strokeWidth: 6, strokeOpacity: 0.8}
+ ),
+ Plot.dot(
+ filtered_niamey_data,
+ {x: "Biweek", y: "Number of Individuals", stroke: "site", fill: "site", fillOpacity: 0.6, r: 12}
+ ),
+ Plot.lineY(
+ filtered_niamey_data,
+ {x: "Biweek", y: "Number of Individuals", stroke: "site"}
+ ),
+ Plot.text(
+ sim_sse,
+ {x: "Biweek", y: "Number of Individuals", text: (d) => `SSE = ${d.sse}`, dx: 0, dy: 0, fontWeight: "bold"}
+ )
+ ]
+})
+```
+
+## Objective Function
+
+Now we set up a function that will calculate the sum of the squared differences between the observations and the model at any parameterization (more commonly known as *"sum of squared errors"*).
+In general, this is called the **objective function** because it is the quantity that optimization seeks to minimize.
+
+```{r}
+#' Calculate the Sum of Squared Errors
#'
-#' @return The SSE of type double
-#' @examples
-sse_sir <- function(params, data) {
- # Convert biweekly time series into annual time scale
- # Daily time scale has requires beta values to be too small - doesn't
- # optimize well
- dt <- 0.01
- max_biweek <- max(data$biweek)
- t <- seq(0, max_biweek * 14, dt) / 365
-
- # Extract the number of observed incidence
- obs_inc <- data$cases
-
- # Note the parameters are updated throughout the optimization process by
- # the optim() function
- # Unpack the transmission parameter and exponentiate to fit on ln scale
- beta <- exp(params[["beta"]])
-
- # Unpack the initial states and exponentiate to fit on normal scale
- S_init <- exp(params[["S_init"]])
- I_init <- exp(params[["I_init"]])
-
- # Fit SIR model to the parameters
- sol <- deSolve::ode(
- y = c(S = S_init, I = I_init, new_inf = 0),
- times = t,
- func = closed_sir_model,
- parms = beta,
- # Use rk4 as fixed time steps, which is important for indexing
- method = "rk4"
- )
-
- # Extract the cumulative incidence
- cum_inc <- sol[, "new_inf"]
-
- # Find the indices of the cumulative incidence to extract
- biweek_index <- seq(1, max_biweek) * (14 / dt) + 1
+#' A function to take in biweekly incidence data, and SIR parameters, and
+#' calculate the SSE
+#'
+#' @param params A vector of parameter values
+#' @param data A dataframe containing biweekly incidence data in the case column
+#'
+#' @return The SSE of type double
+#' @examples
+sse_sir <- function(params, data) {
+ # Convert biweekly time series into annual time scale
+ # Daily time scale has requires beta values to be too small - doesn't
+ # optimize well
+ dt <- 0.01
+ max_biweek <- max(data$biweek)
+ t <- seq(0, max_biweek * 14, dt) / 365
+
+ # Extract the number of observed incidence
+ obs_inc <- data$cases
+
+ # Note the parameters are updated throughout the optimization process by
+ # the optim() function
+ # Unpack the transmission parameter and exponentiate to fit on ln scale
+ beta <- exp(params[["beta"]])
+
+ # Unpack the initial states and exponentiate to fit on normal scale
+ S_init <- exp(params[["S_init"]])
+ I_init <- exp(params[["I_init"]])
+
+ # Fit SIR model to the parameters
+ sol <- deSolve::ode(
+ y = c(S = S_init, I = I_init, new_inf = 0),
+ times = t,
+ func = closed_sir_model,
+ parms = beta,
+ # Use rk4 as fixed time steps, which is important for indexing
+ method = "rk4"
+ )
- # Index cumulative incidence to get the values at the end of the biweeks
- biweek_cum_inc <- cum_inc[biweek_index]
+ # Extract the cumulative incidence
+ cum_inc <- sol[, "new_inf"]
- # Calculate the biweekly incidence by using the difference between
- # consecutive biweeks. Need to manually prepend first week's incidence
- # and add in the initial number of infectious individuals, as ODE model
- # only returns the cumulative differences, which is 0 at the start.
- biweek_inc <- c(biweek_cum_inc[1] + I_init, diff(biweek_cum_inc, lag = 1))
+ # Find the indices of the cumulative incidence to extract
+ biweek_index <- seq(1, max_biweek) * (14 / dt) + 1
+
+ # Index cumulative incidence to get the values at the end of the biweeks
+ biweek_cum_inc <- cum_inc[biweek_index]
- # return SSE of predicted vs observed incidence
- return(sum((biweek_inc - obs_inc)^2))
-}
-```
-
-Notice that the code for `sse_sir()` makes use of the following modeling trick.
-We know that $\beta$, $S_0$, and $I_0$ must be positive, but our search to optimize these parameters will be over the entire number line.
-We could constrain the search using a more sophisticated algorithm, but this might introduce other problems (i.e., stability at the boundaries).
-Instead, we parameterize our objective function (`sse_sir`) in terms of some alternative variables $\ln(\beta)$, $\ln(S_0)$, and $\ln(I_0)$.
-While these numbers range from $-\infty$ to $\infty$ (the range of our search) they map to our model parameters on a range from $0$ to $\infty$ (the range that is biologically meaningful).
+ # Calculate the biweekly incidence by using the difference between
+ # consecutive biweeks. Need to manually prepend first week's incidence
+ # and add in the initial number of infectious individuals, as ODE model
+ # only returns the cumulative differences, which is 0 at the start.
+ biweek_inc <- c(biweek_cum_inc[1] + I_init, diff(biweek_cum_inc, lag = 1))
+
+ # return SSE of predicted vs observed incidence
+ return(sum((biweek_inc - obs_inc)^2))
+}
+```
-## Optimization
-
-Our final step is to use the function `optim` to find the values of $\beta$, $S_0$, and $I_0$ that minimize the sum of squared errors as calculated using our function.
-
-Finally, we plot these fits against the data.
+Notice that the code for `sse_sir()` makes use of the following modeling trick.
+We know that $\beta$, $S_0$, and $I_0$ must be positive, but our search to optimize these parameters will be over the entire number line.
+We could constrain the search using a more sophisticated algorithm, but this might introduce other problems (i.e., stability at the boundaries).
+Instead, we parameterize our objective function (`sse_sir`) in terms of some alternative variables $\ln(\beta)$, $\ln(S_0)$, and $\ln(I_0)$.
+While these numbers range from $-\infty$ to $\infty$ (the range of our search) they map to our model parameters on a range from $0$ to $\infty$ (the range that is biologically meaningful).
-```{r}
-# Initial guess
-sse_optim_params <- c(beta = log(0.055), S_init = log(5000), I_init = log(1))
+## Optimization
+
+Our final step is to use the function `optim` to find the values of $\beta$, $S_0$, and $I_0$ that minimize the sum of squared errors as calculated using our function.
+
+Finally, we plot these fits against the data.
-# Create a dataframe of optimized parameters
-niamey_optims <- niamey_df %>%
- # Create a nested dataframe i.e. one row for each site, and the data column
- # now is a list column that contains a separate dataframe of times and
- # cases for each site
- nest(data = -site) %>%
- mutate(
- # Map the optim() function call to each of the separate dataframes
- # stored in the nested data column we just created
- fit = map(data, ~ optim(sse_optim_params, sse_sir, data = .x)),
- # Map the exp() function to each of the model fits just created, and
- # output to a dataframe instead of a list (like in map()), for easier
- # use in the plottinge predictions later
- map_dfr(fit, ~ exp(.x$par))
- )
-```
-
-```{r}
-niamey_optims %>%
- select(-c(data, fit)) %>%
- mutate(site = str_replace_all(site, "_", " ")) %>%
- gt() %>%
- fmt_number(columns = -site, decimals = 3) %>%
- fmt_scientific(columns = beta, decimals = 3) %>%
- # Relabel the column headers
- cols_label(
- site = md("**Site**"),
- beta = md("**Beta**"),
- S_init = md("**Initial S**"),
- I_init = md("**Initial I**")
- ) %>%
- # Apply style to the table with gray alternating rows
- opt_stylize(style = 1, color = 'gray') %>%
- # Increate space between columns
- opt_horizontal_padding(scale = 3) %>%
- cols_align("center")
-```
-
-::: callout-note
-You may have noticed that you can achieve slightly different results for the optimal parameter values using the interactive plot than are being presented here (though they are very similar).
-This is because while the optimization code is running in `R`, the interactive plot and the calculation of the SSE is implemented using JavaScript.
-Therefore, despite using the same underlying model structure, the answers will vary slightly, because the ODE solvers are different, resulting in different model simulations.
-The difference is not enough to be concerned with here, but it is a point that's worth being aware of when you build your own models - you may want to perform sensitivity to confirm that your model implementation is not driving the magnitude of the results you see, and the inferences you make.
-:::
-
-```{r}
-niamey_predictions <- niamey_optims %>%
- mutate(
- # For each of the different site's nested dataframes, fit the SIR model
- # with the optimal parameters to get best fit predictions
- predictions = pmap(
- .l = list(
- S_init = S_init,
- I_init = I_init,
- beta = beta,
- time_data = data
- ),
- .f = function(S_init, I_init, beta, time_data) {
- site_times <- time_data$biweek * 14 / 365
-
- # Return a dataframe of model solutions
- as_tibble(ode(
- y = c(S = S_init, I = I_init, new_inf = 0),
- times = site_times,
- func = closed_sir_model,
- parms = beta,
- hmax = 1 / 120
- )) %>%
- # Make sure all values are numeric for plotting purposes
- mutate(across(everything(), as.numeric)) %>%
- mutate(
- incidence = ifelse(
- row_number() == 1,
- new_inf[1],
- diff(new_inf, lag = 1)
- )
- )
- }
- )
- ) %>%
- unnest(c(data, predictions))
-```
-
-::: callout-important
-An important point to note is that our data is biweekly **incidence** (new cases in time period), whereas out SIR model produces **prevalence** (total cases at any time point).
-To account for this, our SIR model returns the cumulative incidence (line 32 of the model code), and our objective function extracts the biweekly incidence (lines 41-54), to ensure we are fitting the same data!
-This is a common source of error in interpretation when people fit models to data.
-:::
+```{r}
+# Initial guess
+sse_optim_params <- c(beta = log(0.055), S_init = log(5000), I_init = log(1))
+
+# Create a dataframe of optimized parameters
+niamey_optims <- niamey_df %>%
+ # Create a nested dataframe i.e. one row for each site, and the data column
+ # now is a list column that contains a separate dataframe of times and
+ # cases for each site
+ nest(data = -site) %>%
+ mutate(
+ # Map the optim() function call to each of the separate dataframes
+ # stored in the nested data column we just created
+ fit = map(data, ~ optim(sse_optim_params, sse_sir, data = .x)),
+ # Map the exp() function to each of the model fits just created, and
+ # output to a dataframe instead of a list (like in map()), for easier
+ # use in the plottinge predictions later
+ map_dfr(fit, ~ exp(.x$par))
+ )
+```
+
+```{r}
+niamey_optims %>%
+ select(-c(data, fit)) %>%
+ mutate(site = str_replace_all(site, "_", " ")) %>%
+ gt() %>%
+ fmt_number(columns = -site, decimals = 3) %>%
+ fmt_scientific(columns = beta, decimals = 3) %>%
+ # Relabel the column headers
+ cols_label(
+ site = md("**Site**"),
+ beta = md("**Beta**"),
+ S_init = md("**Initial S**"),
+ I_init = md("**Initial I**")
+ ) %>%
+ # Apply style to the table with gray alternating rows
+ opt_stylize(style = 1, color = 'gray') %>%
+ # Increate space between columns
+ opt_horizontal_padding(scale = 3) %>%
+ cols_align("center")
+```
+
+::: callout-note
+You may have noticed that you can achieve slightly different results for the optimal parameter values using the interactive plot than are being presented here (though they are very similar).
+This is because while the optimization code is running in `R`, the interactive plot and the calculation of the SSE is implemented using JavaScript.
+Therefore, despite using the same underlying model structure, the answers will vary slightly, because the ODE solvers are different, resulting in different model simulations.
+The difference is not enough to be concerned with here, but it is a point that's worth being aware of when you build your own models - you may want to perform sensitivity to confirm that your model implementation is not driving the magnitude of the results you see, and the inferences you make.
+:::
+
+```{r}
+niamey_predictions <- niamey_optims %>%
+ mutate(
+ # For each of the different site's nested dataframes, fit the SIR model
+ # with the optimal parameters to get best fit predictions
+ predictions = pmap(
+ .l = list(
+ S_init = S_init,
+ I_init = I_init,
+ beta = beta,
+ time_data = data
+ ),
+ .f = function(S_init, I_init, beta, time_data) {
+ site_times <- time_data$biweek * 14 / 365
+
+ # Return a dataframe of model solutions
+ as_tibble(ode(
+ y = c(S = S_init, I = I_init, new_inf = 0),
+ times = site_times,
+ func = closed_sir_model,
+ parms = beta,
+ hmax = 1 / 120
+ )) %>%
+ # Make sure all values are numeric for plotting purposes
+ mutate(across(everything(), as.numeric)) %>%
+ mutate(
+ incidence = ifelse(
+ row_number() == 1,
+ new_inf[1],
+ diff(new_inf, lag = 1)
+ )
+ )
+ }
+ )
+ ) %>%
+ unnest(c(data, predictions))
+```
-```{r}
-#| column: body
-#| out-width: 100%
-# Create a dataframe to store the positions of the text labels
-niamey_preds_labels <- tibble(
- site = c("Site_1", "Site_2"),
- x_label = c(6.5, 6.5),
- x_arrow_just = c(-0.5, -0.5),
- x_arrow_end = c(7, 7.75),
- y_label = c(900, 600),
- y_arrow_just = c(-80, -70),
- y_arrow_end = c(350, 290),
- commentary = c("**Predicted", "**Observed"),
- color = c("grey20", niamey_site_colors["Site_2"])
-)
-
-ggplot(niamey_predictions, aes(x = biweek, group = site)) +
- # Plot the actual data in color
- geom_line(aes(y = cases, color = site)) +
- geom_point(aes(y = cases, color = site), size = 4, alpha = 0.8) +
- # Plot the best-fit model predictions in black
- geom_line(aes(y = incidence), color = "black") +
- scale_color_manual(
- values = niamey_site_colors,
- aesthetics = c("color", "fill")
- ) +
- # Place each site on it's own subplot and change labels
- facet_wrap(
- ~site,
- ncol = 1,
- scales = "free_y",
- labeller = as_labeller(niamey_site_labels)
- ) +
- labs(x = "Biweek", y = "Number of Case") +
- theme(legend.position = "none") +
- ggtext::geom_textbox(
- data = niamey_preds_labels,
- aes(
- label = paste0(
- "<span style = \"color:",
- color,
- "\">",
- commentary,
- " Cases**",
- "</span>"
- ),
- x = x_label,
- y = y_label
- ),
- size = 4,
- fill = NA,
- box.colour = NA
- ) +
- geom_curve(
- data = niamey_preds_labels,
- aes(
- x = x_label + x_arrow_just,
- xend = x_arrow_end,
- y = y_label + y_arrow_just,
- yend = y_arrow_end
- ),
- linewidth = 0.75,
- arrow = arrow(length = unit(0.2, "cm")),
- curvature = list(0.25),
- color = "grey20"
- )
-```
-
-::: exercise
-### Exercise 5
-:::
-
-To make things easier, we have assumed the infectious period is known to be 14 days.
-In terms of years, $\text{D} = \frac{14}{365} \approx 0.0384$, and the recovery rate is the inverse i.e., $\gamma = \frac{14}{365}$.
-Now, modify the code above to estimate $\gamma$ and $\beta$ simultaneously.
-
-::: exercise
-### Exercise 6
-:::
-
-What happens if one or both of the other unknowns ($S_0$ and $I_0$) is fixed instead of $\gamma$?
+::: callout-important
+An important point to note is that our data is biweekly **incidence** (new cases in time period), whereas out SIR model produces **prevalence** (total cases at any time point).
+To account for this, our SIR model returns the cumulative incidence (line 32 of the model code), and our objective function extracts the biweekly incidence (lines 41-54), to ensure we are fitting the same data!
+This is a common source of error in interpretation when people fit models to data.
+:::
+
+```{r}
+#| column: body
+#| out-width: 100%
+# Create a dataframe to store the positions of the text labels
+niamey_preds_labels <- tibble(
+ site = c("Site_1", "Site_2"),
+ x_label = c(6.5, 6.5),
+ x_arrow_just = c(-0.5, -0.5),
+ x_arrow_end = c(7, 7.75),
+ y_label = c(900, 600),
+ y_arrow_just = c(-80, -70),
+ y_arrow_end = c(350, 290),
+ commentary = c("**Predicted", "**Observed"),
+ color = c("grey20", niamey_site_colors["Site_2"])
+)
+
+ggplot(niamey_predictions, aes(x = biweek, group = site)) +
+ # Plot the actual data in color
+ geom_line(aes(y = cases, color = site)) +
+ geom_point(aes(y = cases, color = site), size = 4, alpha = 0.8) +
+ # Plot the best-fit model predictions in black
+ geom_line(aes(y = incidence), color = "black") +
+ scale_color_manual(
+ values = niamey_site_colors,
+ aesthetics = c("color", "fill")
+ ) +
+ # Place each site on it's own subplot and change labels
+ facet_wrap(
+ ~site,
+ ncol = 1,
+ scales = "free_y",
+ labeller = as_labeller(niamey_site_labels)
+ ) +
+ labs(x = "Biweek", y = "Number of Case") +
+ theme(legend.position = "none") +
+ ggtext::geom_textbox(
+ data = niamey_preds_labels,
+ aes(
+ label = paste0(
+ "<span style = \"color:",
+ color,
+ "\">",
+ commentary,
+ " Cases**",
+ "</span>"
+ ),
+ x = x_label,
+ y = y_label
+ ),
+ size = 4,
+ fill = NA,
+ box.colour = NA
+ ) +
+ geom_curve(
+ data = niamey_preds_labels,
+ aes(
+ x = x_label + x_arrow_just,
+ xend = x_arrow_end,
+ y = y_label + y_arrow_just,
+ yend = y_arrow_end
+ ),
+ linewidth = 0.75,
+ arrow = arrow(length = unit(0.2, "cm")),
+ curvature = list(0.25),
+ color = "grey20"
+ )
+```
+
+::: exercise
+### Exercise 5
+:::
+
+To make things easier, we have assumed the infectious period is known to be 14 days.
+In terms of years, $\text{D} = \frac{14}{365} \approx 0.0384$, and the recovery rate is the inverse i.e., $\gamma = \frac{14}{365}$.
+Now, modify the code above to estimate $\gamma$ and $\beta$ simultaneously.
-## Solutions
-
-### Exercise 1
+::: exercise
+### Exercise 6
+:::
-```{r}
-#| code-fold: true
-p_infec <- (seq(0, 1, by = 0.001))
-r0_p <- (log(1 - p_infec)) / (-p_infec)
-plot(
- x = p_infec,
- y = r0_p,
- main = "Relationship between final proportion infected and R0",
- xlab = "Final proportion infected",
- ylab = "R0"
-)
-```
-
-### Exercise 2
-
-```{r}
-#| code-fold: true
-#A: If the cases were isolated after 24 hours, then gamma would be 1/1 = 1, and if the cases were isolated after 12 hours, gamma would be 1/0.5 = 2. R0 would be calculated as the beta coefficient over gamma, below:
+What happens if one or both of the other unknowns ($S_0$ and $I_0$) is fixed instead of $\gamma$?
+
+## Solutions
+
+### Exercise 1
+
+```{r}
+#| code-fold: true
+p_infec <- (seq(0, 1, by = 0.001))
+r0_p <- (log(1 - p_infec)) / (-p_infec)
+plot(
+ x = p_infec,
+ y = r0_p,
+ main = "Relationship between final proportion infected and R0",
+ xlab = "Final proportion infected",
+ ylab = "R0"
+)
+```
+
+### Exercise 2
-r0_g1 <- 1.094913 / 1 + 1
-r0_g1
-
-r0_g2 <- 1.094913 / 2 + 1
-r0_g2
-```
-
-### Exercise 3
+```{r}
+#| code-fold: true
+#A: If the cases were isolated after 24 hours, then gamma would be 1/1 = 1, and if the cases were isolated after 12 hours, gamma would be 1/0.5 = 2. R0 would be calculated as the beta coefficient over gamma, below:
+
+r0_g1 <- 1.094913 / 1 + 1
+r0_g1
-```{r}
-#| code-fold: true
-niamey[5, 3] <- 0 #replace a "NA"
-#the command below organizes the data so it can be plotted and analyzed
-niamey <- data.frame(
- biweek = rep(seq(1, 16), 3),
- site = c(rep(1, 16), rep(2, 16), rep(3, 16)),
- cases = c(niamey[, 1], niamey[, 2], niamey[, 3])
-) #define "biweeks"
-
-# As the data are reported every two weeks, this corresponds to the 10th observation. Let’s fit a linear model
-```
-
-```{r}
-#| code-fold: true
-
-# First let's see what the outbreak looks like for the first community on a linear scale
-plot(
- niamey$biweek[niamey$site == 1],
- niamey$cases[niamey$site == 1],
- type = 'p',
- col = niamey$site,
- xlab = 'Biweek',
- ylab = 'Cases'
-)
-lines(niamey$biweek[niamey$site == 1], niamey$cases[niamey$site == 1])
-```
-
-```{r}
-#| code-fold: true
-
-# Now let’s try it on a log scale to see until when the outbreak is roughly linear
-plot(
- niamey$biweek[niamey$site == 1],
- niamey$cases[niamey$site == 1],
- type = 'p',
- col = niamey$site,
- xlab = 'Biweek',
- ylab = 'Cases',
- log = 'y'
-)
-lines(niamey$biweek[niamey$site == 1], niamey$cases[niamey$site == 1])
-```
-
-```{r}
-#| code-fold: true
-# here we create a "week" variable to run the analysis on a weekly
-niamey$week <- niamey$biweek * 2
-# we use the `head` command to take the first N values of a vector. In this case, we're taking the first 10 values of our outcome (cases) and predictor (time) variables.
-model <- lm(
- log(head(niamey$cases[niamey$site == 1], 10)) ~
- (head(niamey$week[niamey$site == 1], 10))
-)
-summary(model) #summary statistics for fit model
-```
-
-```{r}
-#| code-fold: true
-# Now let's get the slope and display it
-slope <- coef(model)[2] #extract slope parameter
-slope #print to screen
-```
-
-As we ran the model by weeks, the $\gamma$ value is $2^{-1}$ and $\hat R_0 = \hat \beta_1 / \gamma +1$ giving $\hat R_0=0.2196041/0.5+1 \approx 1.44$.
-
-### Exercise 4
-
-Here, we can use a loop and repeat the regression procedure we used above for varying numbers of initial data points in our model.
+r0_g2 <- 1.094913 / 2 + 1
+r0_g2
+```
+
+### Exercise 3
+
+```{r}
+#| code-fold: true
+niamey[5, 3] <- 0 #replace a "NA"
+#the command below organizes the data so it can be plotted and analyzed
+niamey <- data.frame(
+ biweek = rep(seq(1, 16), 3),
+ site = c(rep(1, 16), rep(2, 16), rep(3, 16)),
+ cases = c(niamey[, 1], niamey[, 2], niamey[, 3])
+) #define "biweeks"
+
+# As the data are reported every two weeks, this corresponds to the 10th observation. Let’s fit a linear model
+```
+
+```{r}
+#| code-fold: true
+
+# First let's see what the outbreak looks like for the first community on a linear scale
+plot(
+ niamey$biweek[niamey$site == 1],
+ niamey$cases[niamey$site == 1],
+ type = 'p',
+ col = niamey$site,
+ xlab = 'Biweek',
+ ylab = 'Cases'
+)
+lines(niamey$biweek[niamey$site == 1], niamey$cases[niamey$site == 1])
+```
+
+```{r}
+#| code-fold: true
+
+# Now let’s try it on a log scale to see until when the outbreak is roughly linear
+plot(
+ niamey$biweek[niamey$site == 1],
+ niamey$cases[niamey$site == 1],
+ type = 'p',
+ col = niamey$site,
+ xlab = 'Biweek',
+ ylab = 'Cases',
+ log = 'y'
+)
+lines(niamey$biweek[niamey$site == 1], niamey$cases[niamey$site == 1])
+```
+
+```{r}
+#| code-fold: true
+# here we create a "week" variable to run the analysis on a weekly
+niamey$week <- niamey$biweek * 2
+# we use the `head` command to take the first N values of a vector. In this case, we're taking the first 10 values of our outcome (cases) and predictor (time) variables.
+model <- lm(
+ log(head(niamey$cases[niamey$site == 1], 10)) ~
+ (head(niamey$week[niamey$site == 1], 10))
+)
+summary(model) #summary statistics for fit model
+```
+
+```{r}
+#| code-fold: true
+# Now let's get the slope and display it
+slope <- coef(model)[2] #extract slope parameter
+slope #print to screen
+```
-```{r}
-#| code-fold: true
-slope <- NULL
-se <- NULL
-for (i in 3:18) {
- model <- lm(
- log(head(niamey$cases[niamey$site == 1], i)) ~
- (head(niamey$week[niamey$site == 1], i))
- )
- slope <- c(slope, as.numeric(coef(model)[2]))
- se <- c(se, summary(model)$coefficients[4])
-}
-R0 <- slope / 0.5 + 1
-plot(R0, se, ylab = 'Standard error')
-```
-
-### Exercise 5
-
-First, let's add in $\gamma$ estimation into the `sse.sir` function and create a new `sse.sir` function called `sse.sir.g`
-
-```{r}
-#' Basic SIR model
-#'
-#' A basic SIR model with no demographic structure to be used in deSolve
-#'
-#' @param time deSolve passes the time parameter to the function.
-#' @param state A vector of states.
-#' @param params The beta parameter
-#' @param ... Other arguments passed by deSolve.
-#'
-#' @return A deSolve matrix of states at each time step.
-#' @examples
-#' sir_params <- 0.0005
-#' sir_init_states <- c(S = 5000, I = 1, R = 0)
-#' sim_times <- seq(0, 16 / 365, by = 0.1 / 365)
+As we ran the model by weeks, the $\gamma$ value is $2^{-1}$ and $\hat R_0 = \hat \beta_1 / \gamma +1$ giving $\hat R_0=0.2196041/0.5+1 \approx 1.44$.
+
+### Exercise 4
+
+Here, we can use a loop and repeat the regression procedure we used above for varying numbers of initial data points in our model.
+
+```{r}
+#| code-fold: true
+slope <- NULL
+se <- NULL
+for (i in 3:18) {
+ model <- lm(
+ log(head(niamey$cases[niamey$site == 1], i)) ~
+ (head(niamey$week[niamey$site == 1], i))
+ )
+ slope <- c(slope, as.numeric(coef(model)[2]))
+ se <- c(se, summary(model)$coefficients[4])
+}
+R0 <- slope / 0.5 + 1
+plot(R0, se, ylab = 'Standard error')
+```
+
+### Exercise 5
+
+First, let's add in $\gamma$ estimation into the `sse.sir` function and create a new `sse.sir` function called `sse.sir.g`
+
+```{r}
+#' Basic SIR model
+#'
+#' A basic SIR model with no demographic structure to be used in deSolve
+#'
+#' @param time deSolve passes the time parameter to the function.
+#' @param state A vector of states.
+#' @param params The beta parameter
+#' @param ... Other arguments passed by deSolve.
#'
-#' sir_sol <- deSolve::ode(
-#' y = sir_init_states,
-#' times = sim_times,
-#' func = closed_sir_model,
-#' parms = sir_params
-#' )
-
-# Create a new SIR model to include gamma
-closed_sir_model_g <- function(time, state, params, ...) {
- # Unpack states
- S <- state["S"]
- I <- state["I"]
+#' @return A deSolve matrix of states at each time step.
+#' @examples
+#' sir_params <- 0.0005
+#' sir_init_states <- c(S = 5000, I = 1, R = 0)
+#' sim_times <- seq(0, 16 / 365, by = 0.1 / 365)
+#'
+#' sir_sol <- deSolve::ode(
+#' y = sir_init_states,
+#' times = sim_times,
+#' func = closed_sir_model,
+#' parms = sir_params
+#' )
- # Unpack parameters
- beta <- params[["beta"]]
- dur_inf <- 14 / 365
- gamma <- params[["gamma"]]
-
- new_inf <- beta * S * I
-
- # Calculate the ODEs
- dSdt <- -new_inf
- dIdt <- new_inf - (gamma * I)
+# Create a new SIR model to include gamma
+closed_sir_model_g <- function(time, state, params, ...) {
+ # Unpack states
+ S <- state["S"]
+ I <- state["I"]
+
+ # Unpack parameters
+ beta <- params[["beta"]]
+ dur_inf <- 14 / 365
+ gamma <- params[["gamma"]]
- # Return the ODEs
- return(list(c(dSdt, dIdt, new_inf)))
-}
-
-
-#' Calculate the Sum of Squared Errors
-#'
-#' A function to take in biweekly incidence data, and SIR parameters, and
-#' calculate the SSE
-#'
-#' @param params A vector of parameter values
-#' @param data A dataframe containing biweekly incidence data in the case column
+ new_inf <- beta * S * I
+
+ # Calculate the ODEs
+ dSdt <- -new_inf
+ dIdt <- new_inf - (gamma * I)
+
+ # Return the ODEs
+ return(list(c(dSdt, dIdt, new_inf)))
+}
+
+
+#' Calculate the Sum of Squared Errors
#'
-#' @return The SSE of type double
-#' @examples
-sse_sir_g <- function(params, data) {
- # Convert biweekly time series into annual time scale
- # Daily time scale has requires beta values to be too small - doesn't
- # optimize well
- dt <- 0.01
- max_biweek <- max(data$biweek)
- t <- seq(0, max_biweek * 14, dt) / 365
-
- # Extract the number of observed incidence
- obs_inc <- data$cases
-
- # Note the parameters are updated throughout the optimization process by
- # the optim() function
- # Unpack the transmission parameter and exponentiate to fit on ln scale
- beta <- exp(params[["beta"]])
- gamma <- exp(params[["gamma"]])
+#' A function to take in biweekly incidence data, and SIR parameters, and
+#' calculate the SSE
+#'
+#' @param params A vector of parameter values
+#' @param data A dataframe containing biweekly incidence data in the case column
+#'
+#' @return The SSE of type double
+#' @examples
+sse_sir_g <- function(params, data) {
+ # Convert biweekly time series into annual time scale
+ # Daily time scale has requires beta values to be too small - doesn't
+ # optimize well
+ dt <- 0.01
+ max_biweek <- max(data$biweek)
+ t <- seq(0, max_biweek * 14, dt) / 365
+
+ # Extract the number of observed incidence
+ obs_inc <- data$cases
- in_parms <- c("beta" = beta, "gamma" = gamma)
- # Unpack the initial states and exponentiate to fit on normal scale
- S_init <- exp(params[["S_init"]])
- I_init <- exp(params[["I_init"]])
-
- # Fit SIR model to the parameters
- sol <- deSolve::ode(
- y = c(S = S_init, I = I_init, new_inf = 0),
- times = t,
- func = closed_sir_model_g,
- parms = in_parms,
- # Use rk4 as fixed time steps, which is important for indexing
- method = "rk4"
- )
-
- # Extract the cumulative incidence
- cum_inc <- sol[, "new_inf"]
-
- # Find the indices of the cumulative incidence to extract
- biweek_index <- seq(1, max_biweek) * (14 / dt) + 1
+ # Note the parameters are updated throughout the optimization process by
+ # the optim() function
+ # Unpack the transmission parameter and exponentiate to fit on ln scale
+ beta <- exp(params[["beta"]])
+ gamma <- exp(params[["gamma"]])
+
+ in_parms <- c("beta" = beta, "gamma" = gamma)
+ # Unpack the initial states and exponentiate to fit on normal scale
+ S_init <- exp(params[["S_init"]])
+ I_init <- exp(params[["I_init"]])
+
+ # Fit SIR model to the parameters
+ sol <- deSolve::ode(
+ y = c(S = S_init, I = I_init, new_inf = 0),
+ times = t,
+ func = closed_sir_model_g,
+ parms = in_parms,
+ # Use rk4 as fixed time steps, which is important for indexing
+ method = "rk4"
+ )
- # Index cumulative incidence to get the values at the end of the biweeks
- biweek_cum_inc <- cum_inc[biweek_index]
+ # Extract the cumulative incidence
+ cum_inc <- sol[, "new_inf"]
- # Calculate the biweekly incidence by using the difference between
- # consecutive biweeks. Need to manually prepend first week's incidence
- # and add in the initial number of infectious individuals, as ODE model
- # only returns the cumulative differences, which is 0 at the start.
- biweek_inc <- c(biweek_cum_inc[1] + I_init, diff(biweek_cum_inc, lag = 1))
+ # Find the indices of the cumulative incidence to extract
+ biweek_index <- seq(1, max_biweek) * (14 / dt) + 1
+
+ # Index cumulative incidence to get the values at the end of the biweeks
+ biweek_cum_inc <- cum_inc[biweek_index]
- # return SSE of predicted vs observed incidence
- return(sum((biweek_inc - obs_inc)^2))
-}
-```
-
-Now we can run the code to optimize, beta, gamma, I0, and S0
-
-```{r}
-sse_optim_params_g <- c(
- beta = log(0.055),
- gamma = log(365 / 14),
- S_init = log(5000),
- I_init = log(1)
-)
-
-# Create a dataframe of optimized parameters
-niamey_optims_g <- niamey_df %>%
- # Create a nested dataframe i.e. one row for each site, and the data column
- # now is a list column that contains a separate dataframe of times and
- # cases for each site
- nest(data = -site) %>%
- mutate(
- # Map the optim() function call to each of the separate dataframes
- # stored in the nested data column we just created
- fit = map(data, ~ optim(sse_optim_params_g, sse_sir_g, data = .x)),
- # Map the exp() function to each of the model fits just created, and
- # output to a dataframe instead of a list (like in map()), for easier
- # use in the plottinge predictions later
- map_dfr(fit, ~ exp(.x$par))
- )
-```
-
-Now we can get our optimized parameters
-
-```{r}
-niamey_optims_g %>%
- select(-c(data, fit)) %>%
- mutate(site = str_replace_all(site, "_", " ")) %>%
- gt() %>%
- fmt_number(columns = -site, decimals = 3) %>%
- fmt_scientific(columns = beta, decimals = 3) %>%
- # Relabel the column headers
- cols_label(
- site = md("**Site**"),
- beta = md("**Beta**"),
- gamma = md("**Gamma**"),
- S_init = md("**Initial S**"),
- I_init = md("**Initial I**")
- ) %>%
- # Apply style to the table with gray alternating rows
- opt_stylize(style = 1, color = 'gray') %>%
- # Increate space between columns
- opt_horizontal_padding(scale = 3) %>%
- cols_align("center")
-```
-
-Getting new predictions from the optimized parameters and a new model run for each site.
-
-```{r}
-niamey_predictions_g <- niamey_optims_g %>%
- mutate(
- # For each of the different site's nested dataframes, fit the SIR model
- # with the optimal parameters to get best fit predictions
- predictions = pmap(
- .l = list(
- S_init = S_init,
- I_init = I_init,
- beta = beta,
- gamma = gamma,
- time_data = data
- ),
- .f = function(S_init, I_init, beta, gamma, time_data) {
- site_times <- time_data$biweek * 14 / 365
-
- in_parms <- c("beta" = beta, "gamma" = gamma)
- # Return a dataframe of model solutions
- as_tibble(ode(
- y = c(S = S_init, I = I_init, new_inf = 0),
- times = site_times,
- func = closed_sir_model_g,
- parms = in_parms,
- hmax = 1 / 120
- )) %>%
- # Make sure all values are numeric for plotting purposes
- mutate(across(everything(), as.numeric)) %>%
- mutate(
- incidence = ifelse(
- row_number() == 1,
- new_inf[1],
- diff(new_inf, lag = 1)
- )
- )
- }
- )
- ) %>%
- unnest(c(data, predictions))
-```
-
-Finally we can plot:
-
-```{r}
-# Create a dataframe to store the positions of the text labels
-niamey_preds_labels <- tibble(
- site = c("Site_1", "Site_2"),
- x_label = c(6.5, 6.5),
- x_arrow_just = c(-0.5, -0.5),
- x_arrow_end = c(7, 7.75),
- y_label = c(900, 600),
- y_arrow_just = c(-80, -70),
- y_arrow_end = c(350, 290),
- commentary = c("**Predicted", "**Observed"),
- color = c("grey20", niamey_site_colors["Site_2"])
-)
-
-ggplot(niamey_predictions_g, aes(x = biweek, group = site)) +
- # Plot the actual data in color
- geom_line(aes(y = cases, color = site)) +
- geom_point(aes(y = cases, color = site), size = 4, alpha = 0.8) +
- # Plot the best-fit model predictions in black
- geom_line(aes(y = incidence), color = "black") +
- scale_color_manual(
- values = niamey_site_colors,
- aesthetics = c("color", "fill")
- ) +
- # Place each site on it's own subplot and change labels
- facet_wrap(
- ~site,
- ncol = 1,
- scales = "free_y",
- labeller = as_labeller(niamey_site_labels)
- ) +
- labs(x = "Biweek", y = "Number of Case") +
- theme(legend.position = "none") +
- ggtext::geom_textbox(
- data = niamey_preds_labels,
- aes(
- label = paste0(
- "<span style = \"color:",
- color,
- "\">",
- commentary,
- " Cases**",
- "</span>"
- ),
- x = x_label,
- y = y_label
- ),
- size = 4,
- fill = NA,
- box.colour = NA
- ) +
- geom_curve(
- data = niamey_preds_labels,
- aes(
- x = x_label + x_arrow_just,
- xend = x_arrow_end,
- y = y_label + y_arrow_just,
- yend = y_arrow_end
- ),
- linewidth = 0.75,
- arrow = arrow(length = unit(0.2, "cm")),
- curvature = list(0.25),
- color = "grey20"
- )
-```
-
-```{r}
-#| code-fold: true
-#| echo: false
-#| eval: false
-closed.sir.model.g <- function(t, x, params) {
- #SIR model equations
- X <- x[1]
- Y <- x[2]
- beta <- params[1]
- gamma <- params[2]
- dX <- -beta * X * Y
- dY <- beta * X * Y - gamma * Y
- list(c(dX, dY))
-}
-
-sse.sir.g <- function(params0, data, site) {
- #function to calculate squared errors
- data <- data[data$site == site, ] #working dataset, based on site
- t <- data[, 1] * 14 / 365 #time in biweeks
- cases <- data[, 3] #number of cases
- beta <- exp(params0[1]) #parameter beta
- X0 <- exp(params0[2]) #initial susceptibles
- Y0 <- exp(params0[3]) #initial infected
- gamma <- exp(params0[4]) #initial gamma
- out <- as.data.frame(ode(
- c(X = X0, Y = Y0),
- times = t,
- closed.sir.model.g,
- c(beta, gamma),
- hmax = 1 / 120
- ))
- sum((out$Y - cases)^2) #sum of squared errors that is returned to the optim function
-}
-```
-
-```{r}
-#| code-fold: true
-#| echo: false
-#| eval: false
-params0 <- c(-3.2, 7.3, -2.6, log(14 / 365)) # initial guess
-print("Site 1")
-fit1.g <- optim(params0, sse.sir.g, data = niamey, site = 1, hessian = TRUE) #fit
-exp(fit1.g$par) #back-transform parameters
-print("Site 2")
-fit2.g <- optim(params0, sse.sir.g, data = niamey, site = 2, hessian = TRUE) #fit
-exp(fit2.g$par) #back-transform parameters
-print("Site 3")
-fit3.g <- optim(params0, sse.sir.g, data = niamey, site = 3, hessian = TRUE) #fit
-exp(fit3.g$par) #back-transform parameters
-```
-
-```{r}
-#| code-fold: true
-#| echo: false
-#| eval: false
-par(mfrow = c(2, 2)) #set up plotting area for multiple panels
-plot(cases ~ biweek, data = subset(niamey, site == 1), type = 'p', pch = 21) #plot site 1
-t <- subset(niamey, site == 1)[, 1] * 14 / 365
-mod.pred.g1 <- as.data.frame(ode(
- c(X = exp(fit1.g$par[2]), Y = exp(fit1.g$par[3])),
- times = t,
- closed.sir.model.g,
- c(exp(fit1.g$par[1]), exp(fit1.g$par[4])),
- hmax = 1 / 120
-))
-#obtain model predictions
-lines(mod.pred.g1$Y ~ subset(niamey, site == 1)[, 1]) #and plot as a line
-
-plot(cases ~ biweek, data = subset(niamey, site == 2), type = 'b', col = site) #site 2
-t <- subset(niamey, site == 2)[, 1] * 14 / 365
-mod.pred.g2 <- as.data.frame(ode(
- c(X = exp(fit2.g$par[2]), Y = exp(fit2.g$par[3])),
- times = t,
- closed.sir.model.g,
- c(exp(fit2.g$par[1]), exp(fit2.g$par[4])),
- hmax = 1 / 120
-))
-#obtain model predictions
-lines(mod.pred.g2$Y ~ subset(niamey, site == 2)[, 1]) #and plot as a line
-
-plot(cases ~ biweek, data = subset(niamey, site == 3), type = 'b', col = site) #site 3
-t <- subset(niamey, site == 3)[, 1] * 14 / 365
-mod.pred.g3 <- as.data.frame(ode(
- c(X = exp(fit3.g$par[2]), Y = exp(fit3.g$par[3])),
- times = t,
- closed.sir.model.g,
- c(exp(fit3.g$par[1]), exp(fit3.g$par[4])),
- hmax = 1 / 120
-))
-lines(mod.pred.g3$Y ~ subset(niamey, site == 3)[, 1])
-```
-
-### Exercise 6
-
-What happens if one or both of the other unknowns ($X_0$ and $Y_0$) is fixed instead of $\gamma$?
-
-#### Solutions coming soon!
+ # Calculate the biweekly incidence by using the difference between
+ # consecutive biweeks. Need to manually prepend first week's incidence
+ # and add in the initial number of infectious individuals, as ODE model
+ # only returns the cumulative differences, which is 0 at the start.
+ biweek_inc <- c(biweek_cum_inc[1] + I_init, diff(biweek_cum_inc, lag = 1))
+
+ # return SSE of predicted vs observed incidence
+ return(sum((biweek_inc - obs_inc)^2))
+}
+```
+
+Now we can run the code to optimize, beta, gamma, I0, and S0
+
+```{r}
+sse_optim_params_g <- c(
+ beta = log(0.055),
+ gamma = log(365 / 14),
+ S_init = log(5000),
+ I_init = log(1)
+)
+
+# Create a dataframe of optimized parameters
+niamey_optims_g <- niamey_df %>%
+ # Create a nested dataframe i.e. one row for each site, and the data column
+ # now is a list column that contains a separate dataframe of times and
+ # cases for each site
+ nest(data = -site) %>%
+ mutate(
+ # Map the optim() function call to each of the separate dataframes
+ # stored in the nested data column we just created
+ fit = map(data, ~ optim(sse_optim_params_g, sse_sir_g, data = .x)),
+ # Map the exp() function to each of the model fits just created, and
+ # output to a dataframe instead of a list (like in map()), for easier
+ # use in the plottinge predictions later
+ map_dfr(fit, ~ exp(.x$par))
+ )
+```
+
+Now we can get our optimized parameters
+
+```{r}
+niamey_optims_g %>%
+ select(-c(data, fit)) %>%
+ mutate(site = str_replace_all(site, "_", " ")) %>%
+ gt() %>%
+ fmt_number(columns = -site, decimals = 3) %>%
+ fmt_scientific(columns = beta, decimals = 3) %>%
+ # Relabel the column headers
+ cols_label(
+ site = md("**Site**"),
+ beta = md("**Beta**"),
+ gamma = md("**Gamma**"),
+ S_init = md("**Initial S**"),
+ I_init = md("**Initial I**")
+ ) %>%
+ # Apply style to the table with gray alternating rows
+ opt_stylize(style = 1, color = 'gray') %>%
+ # Increate space between columns
+ opt_horizontal_padding(scale = 3) %>%
+ cols_align("center")
+```
+
+Getting new predictions from the optimized parameters and a new model run for each site.
+
+```{r}
+niamey_predictions_g <- niamey_optims_g %>%
+ mutate(
+ # For each of the different site's nested dataframes, fit the SIR model
+ # with the optimal parameters to get best fit predictions
+ predictions = pmap(
+ .l = list(
+ S_init = S_init,
+ I_init = I_init,
+ beta = beta,
+ gamma = gamma,
+ time_data = data
+ ),
+ .f = function(S_init, I_init, beta, gamma, time_data) {
+ site_times <- time_data$biweek * 14 / 365
+
+ in_parms <- c("beta" = beta, "gamma" = gamma)
+ # Return a dataframe of model solutions
+ as_tibble(ode(
+ y = c(S = S_init, I = I_init, new_inf = 0),
+ times = site_times,
+ func = closed_sir_model_g,
+ parms = in_parms,
+ hmax = 1 / 120
+ )) %>%
+ # Make sure all values are numeric for plotting purposes
+ mutate(across(everything(), as.numeric)) %>%
+ mutate(
+ incidence = ifelse(
+ row_number() == 1,
+ new_inf[1],
+ diff(new_inf, lag = 1)
+ )
+ )
+ }
+ )
+ ) %>%
+ unnest(c(data, predictions))
+```
+
+Finally we can plot:
+
+```{r}
+# Create a dataframe to store the positions of the text labels
+niamey_preds_labels <- tibble(
+ site = c("Site_1", "Site_2"),
+ x_label = c(6.5, 6.5),
+ x_arrow_just = c(-0.5, -0.5),
+ x_arrow_end = c(7, 7.75),
+ y_label = c(900, 600),
+ y_arrow_just = c(-80, -70),
+ y_arrow_end = c(350, 290),
+ commentary = c("**Predicted", "**Observed"),
+ color = c("grey20", niamey_site_colors["Site_2"])
+)
+
+ggplot(niamey_predictions_g, aes(x = biweek, group = site)) +
+ # Plot the actual data in color
+ geom_line(aes(y = cases, color = site)) +
+ geom_point(aes(y = cases, color = site), size = 4, alpha = 0.8) +
+ # Plot the best-fit model predictions in black
+ geom_line(aes(y = incidence), color = "black") +
+ scale_color_manual(
+ values = niamey_site_colors,
+ aesthetics = c("color", "fill")
+ ) +
+ # Place each site on it's own subplot and change labels
+ facet_wrap(
+ ~site,
+ ncol = 1,
+ scales = "free_y",
+ labeller = as_labeller(niamey_site_labels)
+ ) +
+ labs(x = "Biweek", y = "Number of Case") +
+ theme(legend.position = "none") +
+ ggtext::geom_textbox(
+ data = niamey_preds_labels,
+ aes(
+ label = paste0(
+ "<span style = \"color:",
+ color,
+ "\">",
+ commentary,
+ " Cases**",
+ "</span>"
+ ),
+ x = x_label,
+ y = y_label
+ ),
+ size = 4,
+ fill = NA,
+ box.colour = NA
+ ) +
+ geom_curve(
+ data = niamey_preds_labels,
+ aes(
+ x = x_label + x_arrow_just,
+ xend = x_arrow_end,
+ y = y_label + y_arrow_just,
+ yend = y_arrow_end
+ ),
+ linewidth = 0.75,
+ arrow = arrow(length = unit(0.2, "cm")),
+ curvature = list(0.25),
+ color = "grey20"
+ )
+```
+
+```{r}
+#| code-fold: true
+#| echo: false
+#| eval: false
+closed.sir.model.g <- function(t, x, params) {
+ #SIR model equations
+ X <- x[1]
+ Y <- x[2]
+ beta <- params[1]
+ gamma <- params[2]
+ dX <- -beta * X * Y
+ dY <- beta * X * Y - gamma * Y
+ list(c(dX, dY))
+}
+
+sse.sir.g <- function(params0, data, site) {
+ #function to calculate squared errors
+ data <- data[data$site == site, ] #working dataset, based on site
+ t <- data[, 1] * 14 / 365 #time in biweeks
+ cases <- data[, 3] #number of cases
+ beta <- exp(params0[1]) #parameter beta
+ X0 <- exp(params0[2]) #initial susceptibles
+ Y0 <- exp(params0[3]) #initial infected
+ gamma <- exp(params0[4]) #initial gamma
+ out <- as.data.frame(ode(
+ c(X = X0, Y = Y0),
+ times = t,
+ closed.sir.model.g,
+ c(beta, gamma),
+ hmax = 1 / 120
+ ))
+ sum((out$Y - cases)^2) #sum of squared errors that is returned to the optim function
+}
+```
+
+```{r}
+#| code-fold: true
+#| echo: false
+#| eval: false
+params0 <- c(-3.2, 7.3, -2.6, log(14 / 365)) # initial guess
+print("Site 1")
+fit1.g <- optim(params0, sse.sir.g, data = niamey, site = 1, hessian = TRUE) #fit
+exp(fit1.g$par) #back-transform parameters
+print("Site 2")
+fit2.g <- optim(params0, sse.sir.g, data = niamey, site = 2, hessian = TRUE) #fit
+exp(fit2.g$par) #back-transform parameters
+print("Site 3")
+fit3.g <- optim(params0, sse.sir.g, data = niamey, site = 3, hessian = TRUE) #fit
+exp(fit3.g$par) #back-transform parameters
+```
+
+```{r}
+#| code-fold: true
+#| echo: false
+#| eval: false
+par(mfrow = c(2, 2)) #set up plotting area for multiple panels
+plot(cases ~ biweek, data = subset(niamey, site == 1), type = 'p', pch = 21) #plot site 1
+t <- subset(niamey, site == 1)[, 1] * 14 / 365
+mod.pred.g1 <- as.data.frame(ode(
+ c(X = exp(fit1.g$par[2]), Y = exp(fit1.g$par[3])),
+ times = t,
+ closed.sir.model.g,
+ c(exp(fit1.g$par[1]), exp(fit1.g$par[4])),
+ hmax = 1 / 120
+))
+#obtain model predictions
+lines(mod.pred.g1$Y ~ subset(niamey, site == 1)[, 1]) #and plot as a line
+
+plot(cases ~ biweek, data = subset(niamey, site == 2), type = 'b', col = site) #site 2
+t <- subset(niamey, site == 2)[, 1] * 14 / 365
+mod.pred.g2 <- as.data.frame(ode(
+ c(X = exp(fit2.g$par[2]), Y = exp(fit2.g$par[3])),
+ times = t,
+ closed.sir.model.g,
+ c(exp(fit2.g$par[1]), exp(fit2.g$par[4])),
+ hmax = 1 / 120
+))
+#obtain model predictions
+lines(mod.pred.g2$Y ~ subset(niamey, site == 2)[, 1]) #and plot as a line
+
+plot(cases ~ biweek, data = subset(niamey, site == 3), type = 'b', col = site) #site 3
+t <- subset(niamey, site == 3)[, 1] * 14 / 365
+mod.pred.g3 <- as.data.frame(ode(
+ c(X = exp(fit3.g$par[2]), Y = exp(fit3.g$par[3])),
+ times = t,
+ closed.sir.model.g,
+ c(exp(fit3.g$par[1]), exp(fit3.g$par[4])),
+ hmax = 1 / 120
+))
+lines(mod.pred.g3$Y ~ subset(niamey, site == 3)[, 1])
+```
-First we modify the `sse_sir_g` function again to fix $X_0$ and $Y_0$ (or both)
+### Exercise 6
-```{r}
-#| code-fold: true
-#| echo: false
-#| eval: false
-# fix S0
-sse_sir_S0 <- function(params, data) {
- # Convert biweekly time series into annual time scale
- # Daily time scale has requires beta values to be too small - doesn't
- # optimize well
- dt <- 0.01
- max_biweek <- max(data$biweek)
- t <- seq(0, max_biweek * 14, dt) / 365
-
- # Extract the number of observed incidence
- obs_inc <- data$cases
-
- # Note the parameters are updated throughout the optimization process by
- # the optim() function
- # Unpack the transmission parameter and exponentiate to fit on ln scale
- beta <- exp(params[["beta"]])
- gamma <- exp(params[["gamma"]])
+What happens if one or both of the other unknowns ($X_0$ and $Y_0$) is fixed instead of $\gamma$?
+
+#### Solutions coming soon!
+
+First we modify the `sse_sir_g` function again to fix $X_0$ and $Y_0$ (or both)
+
+```{r}
+#| code-fold: true
+#| echo: false
+#| eval: false
+# fix S0
+sse_sir_S0 <- function(params, data) {
+ # Convert biweekly time series into annual time scale
+ # Daily time scale has requires beta values to be too small - doesn't
+ # optimize well
+ dt <- 0.01
+ max_biweek <- max(data$biweek)
+ t <- seq(0, max_biweek * 14, dt) / 365
+
+ # Extract the number of observed incidence
+ obs_inc <- data$cases
- in_parms <- c("beta" = beta, "gamma" = gamma)
- # Unpack the initial states and exponentiate to fit on normal scale
- S_init <- exp(12500) #initial susceptibles (FIXED)
- I_init <- exp(params[["I_init"]])
-
- # Fit SIR model to the parameters
- sol <- deSolve::ode(
- y = c(I = I_init, new_inf = 0),
- times = t,
- func = closed_sir_model_g,
- parms = in_parms,
- # Use rk4 as fixed time steps, which is important for indexing
- method = "rk4"
- )
-
- # Extract the cumulative incidence
- cum_inc <- sol[, "new_inf"]
-
- # Find the indices of the cumulative incidence to extract
- biweek_index <- seq(1, max_biweek) * (14 / dt) + 1
+ # Note the parameters are updated throughout the optimization process by
+ # the optim() function
+ # Unpack the transmission parameter and exponentiate to fit on ln scale
+ beta <- exp(params[["beta"]])
+ gamma <- exp(params[["gamma"]])
+
+ in_parms <- c("beta" = beta, "gamma" = gamma)
+ # Unpack the initial states and exponentiate to fit on normal scale
+ S_init <- exp(12500) #initial susceptibles (FIXED)
+ I_init <- exp(params[["I_init"]])
+
+ # Fit SIR model to the parameters
+ sol <- deSolve::ode(
+ y = c(I = I_init, new_inf = 0),
+ times = t,
+ func = closed_sir_model_g,
+ parms = in_parms,
+ # Use rk4 as fixed time steps, which is important for indexing
+ method = "rk4"
+ )
- # Index cumulative incidence to get the values at the end of the biweeks
- biweek_cum_inc <- cum_inc[biweek_index]
+ # Extract the cumulative incidence
+ cum_inc <- sol[, "new_inf"]
- # Calculate the biweekly incidence by using the difference between
- # consecutive biweeks. Need to manually prepend first week's incidence
- # and add in the initial number of infectious individuals, as ODE model
- # only returns the cumulative differences, which is 0 at the start.
- biweek_inc <- c(biweek_cum_inc[1] + I_init, diff(biweek_cum_inc, lag = 1))
+ # Find the indices of the cumulative incidence to extract
+ biweek_index <- seq(1, max_biweek) * (14 / dt) + 1
+
+ # Index cumulative incidence to get the values at the end of the biweeks
+ biweek_cum_inc <- cum_inc[biweek_index]
- # return SSE of predicted vs observed incidence
- return(sum((biweek_inc - obs_inc)^2))
-}
-
-#fix I0
-sse_sir_I0 <- function(params, data) {
- # Convert biweekly time series into annual time scale
- # Daily time scale has requires beta values to be too small - doesn't
- # optimize well
- dt <- 0.01
- max_biweek <- max(data$biweek)
- t <- seq(0, max_biweek * 14, dt) / 365
-
- # Extract the number of observed incidence
- obs_inc <- data$cases
-
- # Note the parameters are updated throughout the optimization process by
- # the optim() function
- # Unpack the transmission parameter and exponentiate to fit on ln scale
- beta <- exp(params[["beta"]])
- gamma <- exp(params[["gamma"]])
+ # Calculate the biweekly incidence by using the difference between
+ # consecutive biweeks. Need to manually prepend first week's incidence
+ # and add in the initial number of infectious individuals, as ODE model
+ # only returns the cumulative differences, which is 0 at the start.
+ biweek_inc <- c(biweek_cum_inc[1] + I_init, diff(biweek_cum_inc, lag = 1))
+
+ # return SSE of predicted vs observed incidence
+ return(sum((biweek_inc - obs_inc)^2))
+}
+
+#fix I0
+sse_sir_I0 <- function(params, data) {
+ # Convert biweekly time series into annual time scale
+ # Daily time scale has requires beta values to be too small - doesn't
+ # optimize well
+ dt <- 0.01
+ max_biweek <- max(data$biweek)
+ t <- seq(0, max_biweek * 14, dt) / 365
+
+ # Extract the number of observed incidence
+ obs_inc <- data$cases
- in_parms <- c("beta" = beta, "gamma" = gamma)
- # Unpack the initial states and exponentiate to fit on normal scale
- S_init <- exp(params[["S_init"]])
- I_init <- exp(1000) # initial infecteds (FIXED)
-
- # Fit SIR model to the parameters
- sol <- deSolve::ode(
- y = c(S = S_init, I = I_init, new_inf = 0),
- times = t,
- func = closed_sir_model_g,
- parms = in_parms,
- # Use rk4 as fixed time steps, which is important for indexing
- method = "rk4"
- )
-
- # Extract the cumulative incidence
- cum_inc <- sol[, "new_inf"]
-
- # Find the indices of the cumulative incidence to extract
- biweek_index <- seq(1, max_biweek) * (14 / dt) + 1
+ # Note the parameters are updated throughout the optimization process by
+ # the optim() function
+ # Unpack the transmission parameter and exponentiate to fit on ln scale
+ beta <- exp(params[["beta"]])
+ gamma <- exp(params[["gamma"]])
+
+ in_parms <- c("beta" = beta, "gamma" = gamma)
+ # Unpack the initial states and exponentiate to fit on normal scale
+ S_init <- exp(params[["S_init"]])
+ I_init <- exp(1000) # initial infecteds (FIXED)
+
+ # Fit SIR model to the parameters
+ sol <- deSolve::ode(
+ y = c(S = S_init, I = I_init, new_inf = 0),
+ times = t,
+ func = closed_sir_model_g,
+ parms = in_parms,
+ # Use rk4 as fixed time steps, which is important for indexing
+ method = "rk4"
+ )
- # Index cumulative incidence to get the values at the end of the biweeks
- biweek_cum_inc <- cum_inc[biweek_index]
+ # Extract the cumulative incidence
+ cum_inc <- sol[, "new_inf"]
- # Calculate the biweekly incidence by using the difference between
- # consecutive biweeks. Need to manually prepend first week's incidence
- # and add in the initial number of infectious individuals, as ODE model
- # only returns the cumulative differences, which is 0 at the start.
- biweek_inc <- c(biweek_cum_inc[1] + I_init, diff(biweek_cum_inc, lag = 1))
+ # Find the indices of the cumulative incidence to extract
+ biweek_index <- seq(1, max_biweek) * (14 / dt) + 1
+
+ # Index cumulative incidence to get the values at the end of the biweeks
+ biweek_cum_inc <- cum_inc[biweek_index]
- # return SSE of predicted vs observed incidence
- return(sum((biweek_inc - obs_inc)^2))
-}
-
-# fix S0 and I0
-sse_sir_S0_I0 <- function(params, data) {
- # Convert biweekly time series into annual time scale
- # Daily time scale has requires beta values to be too small - doesn't
- # optimize well
- dt <- 0.01
- max_biweek <- max(data$biweek)
- t <- seq(0, max_biweek * 14, dt) / 365
-
- # Extract the number of observed incidence
- obs_inc <- data$cases
-
- # Note the parameters are updated throughout the optimization process by
- # the optim() function
- # Unpack the transmission parameter and exponentiate to fit on ln scale
- beta <- exp(params[["beta"]])
- gamma <- exp(params[["gamma"]])
+ # Calculate the biweekly incidence by using the difference between
+ # consecutive biweeks. Need to manually prepend first week's incidence
+ # and add in the initial number of infectious individuals, as ODE model
+ # only returns the cumulative differences, which is 0 at the start.
+ biweek_inc <- c(biweek_cum_inc[1] + I_init, diff(biweek_cum_inc, lag = 1))
+
+ # return SSE of predicted vs observed incidence
+ return(sum((biweek_inc - obs_inc)^2))
+}
+
+# fix S0 and I0
+sse_sir_S0_I0 <- function(params, data) {
+ # Convert biweekly time series into annual time scale
+ # Daily time scale has requires beta values to be too small - doesn't
+ # optimize well
+ dt <- 0.01
+ max_biweek <- max(data$biweek)
+ t <- seq(0, max_biweek * 14, dt) / 365
+
+ # Extract the number of observed incidence
+ obs_inc <- data$cases
- in_parms <- c("beta" = beta, "gamma" = gamma)
- # Unpack the initial states and exponentiate to fit on normal scale
- S_init <- exp(12500) # initial susceptibles (FIXED)
- I_init <- exp(1000) # initial infecteds (FIXED)
-
- # Fit SIR model to the parameters
- sol <- deSolve::ode(
- y = c(S = S_init, I = I_init, new_inf = 0),
- times = t,
- func = closed_sir_model_g,
- parms = in_parms,
- # Use rk4 as fixed time steps, which is important for indexing
- method = "rk4"
- )
-
- # Extract the cumulative incidence
- cum_inc <- sol[, "new_inf"]
-
- # Find the indices of the cumulative incidence to extract
- biweek_index <- seq(1, max_biweek) * (14 / dt) + 1
+ # Note the parameters are updated throughout the optimization process by
+ # the optim() function
+ # Unpack the transmission parameter and exponentiate to fit on ln scale
+ beta <- exp(params[["beta"]])
+ gamma <- exp(params[["gamma"]])
+
+ in_parms <- c("beta" = beta, "gamma" = gamma)
+ # Unpack the initial states and exponentiate to fit on normal scale
+ S_init <- exp(12500) # initial susceptibles (FIXED)
+ I_init <- exp(1000) # initial infecteds (FIXED)
+
+ # Fit SIR model to the parameters
+ sol <- deSolve::ode(
+ y = c(S = S_init, I = I_init, new_inf = 0),
+ times = t,
+ func = closed_sir_model_g,
+ parms = in_parms,
+ # Use rk4 as fixed time steps, which is important for indexing
+ method = "rk4"
+ )
- # Index cumulative incidence to get the values at the end of the biweeks
- biweek_cum_inc <- cum_inc[biweek_index]
+ # Extract the cumulative incidence
+ cum_inc <- sol[, "new_inf"]
- # Calculate the biweekly incidence by using the difference between
- # consecutive biweeks. Need to manually prepend first week's incidence
- # and add in the initial number of infectious individuals, as ODE model
- # only returns the cumulative differences, which is 0 at the start.
- biweek_inc <- c(biweek_cum_inc[1] + I_init, diff(biweek_cum_inc, lag = 1))
+ # Find the indices of the cumulative incidence to extract
+ biweek_index <- seq(1, max_biweek) * (14 / dt) + 1
+
+ # Index cumulative incidence to get the values at the end of the biweeks
+ biweek_cum_inc <- cum_inc[biweek_index]
- # return SSE of predicted vs observed incidence
- return(sum((biweek_inc - obs_inc)^2))
-}
-```
-
-Now we can run the optim algorithm and find the best parameters for the SIR model for each of the scenario.
-
-First, for S0 fixed:
-
-```{r}
-#| code-fold: true
-#| echo: false
-#| eval: false
-sse_optim_params_S0 <- c(
- beta = log(0.055),
- gamma = log(365 / 14),
- I_init = log(1)
-)
-
-# Create a dataframe of optimized parameters
-niamey_optims_S0 <- niamey_df %>%
- # Create a nested dataframe i.e. one row for each site, and the data column
- # now is a list column that contains a separate dataframe of times and
- # cases for each site
- nest(data = -site) %>%
- mutate(
- # Map the optim() function call to each of the separate dataframes
- # stored in the nested data column we just created
- fit = map(data, ~ optim(sse_optim_params_S0, sse_sir_S0, data = .x)),
- # Map the exp() function to each of the model fits just created, and
- # output to a dataframe instead of a list (like in map()), for easier
- # use in the plottinge predictions later
- map_dfr(fit, ~ exp(.x$par))
- )
-```
-
-Second, for I0 fixed:
-
-```{r}
-#| code-fold: true
-#| echo: false
-#| eval: false
-sse_optim_params_I0 <- c(
- beta = log(0.055),
- gamma = log(365 / 14),
- S_init = log(5000),
- I_init = log(1)
-)
-
-# Create a dataframe of optimized parameters
-niamey_optims_I0 <- niamey_df %>%
- # Create a nested dataframe i.e. one row for each site, and the data column
- # now is a list column that contains a separate dataframe of times and
- # cases for each site
- nest(data = -site) %>%
- mutate(
- # Map the optim() function call to each of the separate dataframes
- # stored in the nested data column we just created
- fit = map(data, ~ optim(sse_optim_params_I0, sse_sir_I0, data = .x)),
- # Map the exp() function to each of the model fits just created, and
- # output to a dataframe instead of a list (like in map()), for easier
- # use in the plottinge predictions later
- map_dfr(fit, ~ exp(.x$par))
- )
-```
-
-Third, for fixing both S0 and I0:
-
-```{r}
-#| code-fold: true
-#| echo: false
-#| eval: false
-sse_optim_params_S0_I0 <- c(
- beta = log(0.055),
- gamma = log(365 / 14),
- S_init = log(5000),
- I_init = log(1)
-)
-
-# Create a dataframe of optimized parameters
-niamey_optims_S0_I0 <- niamey_df %>%
- # Create a nested dataframe i.e. one row for each site, and the data column
- # now is a list column that contains a separate dataframe of times and
- # cases for each site
- nest(data = -site) %>%
- mutate(
- # Map the optim() function call to each of the separate dataframes
- # stored in the nested data column we just created
- fit = map(data, ~ optim(sse_optim_params_S0_I0, sse_sir_S0_I0, data = .x)),
- # Map the exp() function to each of the model fits just created, and
- # output to a dataframe instead of a list (like in map()), for easier
- # use in the plottinge predictions later
- map_dfr(fit, ~ exp(.x$par))
- )
-```
-
-```{r}
-#| code-fold: true
-#| echo: false
-#| eval: false
-sse.sir.x0 <- function(params0, data, site) {
- #function to calculate squared errors
- data <- data[data$site == site, ] #working dataset, based on site
- t <- data[, 1] * 14 / 365 #time in biweeks
- cases <- data[, 3] #number of cases
- beta <- exp(params0[1]) #parameter beta
- X0 <- 12500 #initial susceptibles (FIXED)
- Y0 <- exp(params0[2]) #initial infected
- gamma <- exp(params0[3])
- out <- as.data.frame(ode(
- c(X = X0, Y = Y0),
- times = t,
- closed.sir.model.g,
- c(beta, gamma),
- hmax = 1 / 120
- ))
- sum((out$Y - cases)^2) #sum of squared errors that is returned to the optim function
-}
-
-
-sse.sir.y0 <- function(params0, data, site) {
- #function to calculate squared errors
- data <- data[data$site == site, ] #working dataset, based on site
- t <- data[, 1] * 14 / 365 #time in biweeks
- cases <- data[, 3] #number of cases
- beta <- exp(params0[1]) #parameter beta
- X0 <- exp(params0[2]) #initial susceptibles (FIXED)
- Y0 <- 1000 #initial infected
- gamma <- exp(params0[3])
- out <- as.data.frame(ode(
- c(X = X0, Y = Y0),
- times = t,
- closed.sir.model.g,
- c(beta, gamma),
- hmax = 1 / 120
- ))
- sum((out$Y - cases)^2) #sum of squared errors that is returned to the optim function
-}
-
-
-sse.sir.x0y0 <- function(params0, data, site) {
- #function to calculate squared errors
- data <- data[data$site == site, ] #working dataset, based on site
- t <- data[, 1] * 14 / 365 #time in biweeks
- cases <- data[, 3] #number of cases
- beta <- exp(params0[1]) #parameter beta
- X0 <- 12500 #initial susceptibles (FIXED)
- Y0 <- 1000 #initial infected (FIXED)
- gamma <- exp(params0[2])
- out <- as.data.frame(ode(
- c(X = X0, Y = Y0),
- times = t,
- closed.sir.model.g,
- c(beta, gamma),
- hmax = 1 / 120
- ))
- sum((out$Y - cases)^2) #sum of squared errors that is returned to the optim function
-}
-```
-
-Now let's pass these models through `optim` and fit the models:
-
-```{r}
-#| code-fold: true
-#| echo: false
-#| eval: false
-# Fixing X0
-params0 <- c(-3.2, 1, log(14 / 365)) #initial guess of beta, Y0, and gamma
-fit1.x0 <- optim(params0, sse.sir.x0, data = niamey, site = 1, hessian = TRUE) #fit
-exp(fit1.x0$par) #back-transform parameters
-
-#Fixing Y0
-params0 <- c(-3.2, 3, log(14 / 365)) #initial guess of beta, X0, and gamma
-fit1.y0 <- optim(params0, sse.sir.y0, data = niamey, site = 1, hessian = TRUE) #fit
-exp(fit1.y0$par) #back-transform parameters
-
-#Fixing X0 and Y0
-params0 <- c(-3.2, log(14 / 365)) #initial guess of beta and gamma
-fit1.x0y0 <- optim(
- params0,
- sse.sir.x0y0,
- data = niamey,
- site = 1,
- hessian = TRUE
-) #fit
-exp(fit1.x0y0$par) #back-transform parameters
-```
+ # Calculate the biweekly incidence by using the difference between
+ # consecutive biweeks. Need to manually prepend first week's incidence
+ # and add in the initial number of infectious individuals, as ODE model
+ # only returns the cumulative differences, which is 0 at the start.
+ biweek_inc <- c(biweek_cum_inc[1] + I_init, diff(biweek_cum_inc, lag = 1))
+
+ # return SSE of predicted vs observed incidence
+ return(sum((biweek_inc - obs_inc)^2))
+}
+```
+
+Now we can run the optim algorithm and find the best parameters for the SIR model for each of the scenario.
+
+First, for S0 fixed:
+
+```{r}
+#| code-fold: true
+#| echo: false
+#| eval: false
+sse_optim_params_S0 <- c(
+ beta = log(0.055),
+ gamma = log(365 / 14),
+ I_init = log(1)
+)
+
+# Create a dataframe of optimized parameters
+niamey_optims_S0 <- niamey_df %>%
+ # Create a nested dataframe i.e. one row for each site, and the data column
+ # now is a list column that contains a separate dataframe of times and
+ # cases for each site
+ nest(data = -site) %>%
+ mutate(
+ # Map the optim() function call to each of the separate dataframes
+ # stored in the nested data column we just created
+ fit = map(data, ~ optim(sse_optim_params_S0, sse_sir_S0, data = .x)),
+ # Map the exp() function to each of the model fits just created, and
+ # output to a dataframe instead of a list (like in map()), for easier
+ # use in the plottinge predictions later
+ map_dfr(fit, ~ exp(.x$par))
+ )
+```
+
+Second, for I0 fixed:
+
+```{r}
+#| code-fold: true
+#| echo: false
+#| eval: false
+sse_optim_params_I0 <- c(
+ beta = log(0.055),
+ gamma = log(365 / 14),
+ S_init = log(5000),
+ I_init = log(1)
+)
+
+# Create a dataframe of optimized parameters
+niamey_optims_I0 <- niamey_df %>%
+ # Create a nested dataframe i.e. one row for each site, and the data column
+ # now is a list column that contains a separate dataframe of times and
+ # cases for each site
+ nest(data = -site) %>%
+ mutate(
+ # Map the optim() function call to each of the separate dataframes
+ # stored in the nested data column we just created
+ fit = map(data, ~ optim(sse_optim_params_I0, sse_sir_I0, data = .x)),
+ # Map the exp() function to each of the model fits just created, and
+ # output to a dataframe instead of a list (like in map()), for easier
+ # use in the plottinge predictions later
+ map_dfr(fit, ~ exp(.x$par))
+ )
+```
+
+Third, for fixing both S0 and I0:
+
+```{r}
+#| code-fold: true
+#| echo: false
+#| eval: false
+sse_optim_params_S0_I0 <- c(
+ beta = log(0.055),
+ gamma = log(365 / 14),
+ S_init = log(5000),
+ I_init = log(1)
+)
+
+# Create a dataframe of optimized parameters
+niamey_optims_S0_I0 <- niamey_df %>%
+ # Create a nested dataframe i.e. one row for each site, and the data column
+ # now is a list column that contains a separate dataframe of times and
+ # cases for each site
+ nest(data = -site) %>%
+ mutate(
+ # Map the optim() function call to each of the separate dataframes
+ # stored in the nested data column we just created
+ fit = map(data, ~ optim(sse_optim_params_S0_I0, sse_sir_S0_I0, data = .x)),
+ # Map the exp() function to each of the model fits just created, and
+ # output to a dataframe instead of a list (like in map()), for easier
+ # use in the plottinge predictions later
+ map_dfr(fit, ~ exp(.x$par))
+ )
+```
+
+```{r}
+#| code-fold: true
+#| echo: false
+#| eval: false
+sse.sir.x0 <- function(params0, data, site) {
+ #function to calculate squared errors
+ data <- data[data$site == site, ] #working dataset, based on site
+ t <- data[, 1] * 14 / 365 #time in biweeks
+ cases <- data[, 3] #number of cases
+ beta <- exp(params0[1]) #parameter beta
+ X0 <- 12500 #initial susceptibles (FIXED)
+ Y0 <- exp(params0[2]) #initial infected
+ gamma <- exp(params0[3])
+ out <- as.data.frame(ode(
+ c(X = X0, Y = Y0),
+ times = t,
+ closed.sir.model.g,
+ c(beta, gamma),
+ hmax = 1 / 120
+ ))
+ sum((out$Y - cases)^2) #sum of squared errors that is returned to the optim function
+}
+
+
+sse.sir.y0 <- function(params0, data, site) {
+ #function to calculate squared errors
+ data <- data[data$site == site, ] #working dataset, based on site
+ t <- data[, 1] * 14 / 365 #time in biweeks
+ cases <- data[, 3] #number of cases
+ beta <- exp(params0[1]) #parameter beta
+ X0 <- exp(params0[2]) #initial susceptibles (FIXED)
+ Y0 <- 1000 #initial infected
+ gamma <- exp(params0[3])
+ out <- as.data.frame(ode(
+ c(X = X0, Y = Y0),
+ times = t,
+ closed.sir.model.g,
+ c(beta, gamma),
+ hmax = 1 / 120
+ ))
+ sum((out$Y - cases)^2) #sum of squared errors that is returned to the optim function
+}
+
+
+sse.sir.x0y0 <- function(params0, data, site) {
+ #function to calculate squared errors
+ data <- data[data$site == site, ] #working dataset, based on site
+ t <- data[, 1] * 14 / 365 #time in biweeks
+ cases <- data[, 3] #number of cases
+ beta <- exp(params0[1]) #parameter beta
+ X0 <- 12500 #initial susceptibles (FIXED)
+ Y0 <- 1000 #initial infected (FIXED)
+ gamma <- exp(params0[2])
+ out <- as.data.frame(ode(
+ c(X = X0, Y = Y0),
+ times = t,
+ closed.sir.model.g,
+ c(beta, gamma),
+ hmax = 1 / 120
+ ))
+ sum((out$Y - cases)^2) #sum of squared errors that is returned to the optim function
+}
+```
+
+Now let's pass these models through `optim` and fit the models:
+
+```{r}
+#| code-fold: true
+#| echo: false
+#| eval: false
+# Fixing X0
+params0 <- c(-3.2, 1, log(14 / 365)) #initial guess of beta, Y0, and gamma
+fit1.x0 <- optim(params0, sse.sir.x0, data = niamey, site = 1, hessian = TRUE) #fit
+exp(fit1.x0$par) #back-transform parameters
+
+#Fixing Y0
+params0 <- c(-3.2, 3, log(14 / 365)) #initial guess of beta, X0, and gamma
+fit1.y0 <- optim(params0, sse.sir.y0, data = niamey, site = 1, hessian = TRUE) #fit
+exp(fit1.y0$par) #back-transform parameters
+
+#Fixing X0 and Y0
+params0 <- c(-3.2, log(14 / 365)) #initial guess of beta and gamma
+fit1.x0y0 <- optim(
+ params0,
+ sse.sir.x0y0,
+ data = niamey,
+ site = 1,
+ hessian = TRUE
+) #fit
+exp(fit1.x0y0$par) #back-transform parameters
+```
\n\n\n \n
\n
\n```\n\n:::\n:::\n\n\n::: callout-note\nYou may have noticed that you can achieve slightly different results for the optimal parameter values using the interactive plot than are being presented here (though they are very similar).\nThis is because while the optimization code is running in `R`, the interactive plot and the calculation of the SSE is implemented using JavaScript.\nTherefore, despite using the same underlying model structure, the answers will vary slightly, because the ODE solvers are different, resulting in different model simulations.\nThe difference is not enough to be concerned with here, but it is a point that's worth being aware of when you build your own models - you may want to perform sensitivity to confirm that your model implementation is not driving the magnitude of the results you see, and the inferences you make.\n:::\n\n\n::: {.cell}\n\n```{.r .cell-code}\nniamey_predictions <- niamey_optims %>%\n mutate(\n # For each of the different site's nested dataframes, fit the SIR model\n # with the optimal parameters to get best fit predictions\n predictions = pmap(\n .l = list(\n S_init = S_init,\n I_init = I_init,\n beta = beta,\n time_data = data\n ),\n .f = function(S_init, I_init, beta, time_data) {\n site_times <- time_data$biweek * 14 / 365\n\n # Return a dataframe of model solutions\n as_tibble(ode(\n y = c(S = S_init, I = I_init, new_inf = 0),\n times = site_times,\n func = closed_sir_model,\n parms = beta,\n hmax = 1 / 120\n )) %>%\n # Make sure all values are numeric for plotting purposes\n mutate(across(everything(), as.numeric)) %>%\n mutate(\n incidence = ifelse(\n row_number() == 1,\n new_inf[1],\n diff(new_inf, lag = 1)\n )\n )\n }\n )\n ) %>%\n unnest(c(data, predictions))\n```\n:::\n\n\n::: callout-important\nAn important point to note is that our data is biweekly **incidence** (new cases in time period), whereas out SIR model produces **prevalence** (total cases at any time point).\nTo account for this, our SIR model returns the cumulative incidence (line 32 of the model code), and our objective function extracts the biweekly incidence (lines 41-54), to ensure we are fitting the same data!\nThis is a common source of error in interpretation when people fit models to data.\n:::\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\n# Create a dataframe to store the positions of the text labels\nniamey_preds_labels <- tibble(\n site = c(\"Site_1\", \"Site_2\"),\n x_label = c(6.5, 6.5),\n x_arrow_just = c(-0.5, -0.5),\n x_arrow_end = c(7, 7.75),\n y_label = c(900, 600),\n y_arrow_just = c(-80, -70),\n y_arrow_end = c(350, 290),\n commentary = c(\"**Predicted\", \"**Observed\"),\n color = c(\"grey20\", niamey_site_colors[\"Site_2\"])\n)\n\nggplot(niamey_predictions, aes(x = biweek, group = site)) +\n # Plot the actual data in color\n geom_line(aes(y = cases, color = site)) +\n geom_point(aes(y = cases, color = site), size = 4, alpha = 0.8) +\n # Plot the best-fit model predictions in black\n geom_line(aes(y = incidence), color = \"black\") +\n scale_color_manual(\n values = niamey_site_colors,\n aesthetics = c(\"color\", \"fill\")\n ) +\n # Place each site on it's own subplot and change labels\n facet_wrap(\n ~site,\n ncol = 1,\n scales = \"free_y\",\n labeller = as_labeller(niamey_site_labels)\n ) +\n labs(x = \"Biweek\", y = \"Number of Case\") +\n theme(legend.position = \"none\") +\n ggtext::geom_textbox(\n data = niamey_preds_labels,\n aes(\n label = paste0(\n \"\",\n commentary,\n \" Cases**\",\n \"\"\n ),\n x = x_label,\n y = y_label\n ),\n size = 4,\n fill = NA,\n box.colour = NA\n ) +\n geom_curve(\n data = niamey_preds_labels,\n aes(\n x = x_label + x_arrow_just,\n xend = x_arrow_end,\n y = y_label + y_arrow_just,\n yend = y_arrow_end\n ),\n linewidth = 0.75,\n arrow = arrow(length = unit(0.2, \"cm\")),\n curvature = list(0.25),\n color = \"grey20\"\n )\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\n::: exercise\n### Exercise 5\n:::\n\nTo make things easier, we have assumed the infectious period is known to be 14 days.\nIn terms of years, $\\text{D} = \\frac{14}{365} \\approx 0.0384$, and the recovery rate is the inverse i.e., $\\gamma = \\frac{14}{365}$.\nNow, modify the code above to estimate $\\gamma$ and $\\beta$ simultaneously.\n\n::: exercise\n### Exercise 6\n:::\n\nWhat happens if one or both of the other unknowns ($S_0$ and $I_0$) is fixed instead of $\\gamma$?\n\n## Solutions\n\n### Exercise 1\n\n\n::: {.cell}\n\n```{.r .cell-code code-fold=\"true\"}\np_infec <- (seq(0, 1, by = 0.001))\nr0_p <- (log(1 - p_infec)) / (-p_infec)\nplot(\n x = p_infec,\n y = r0_p,\n main = \"Relationship between final proportion infected and R0\",\n xlab = \"Final proportion infected\",\n ylab = \"R0\"\n)\n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n:::\n\n\n### Exercise 2\n\n\n::: {.cell}\n\n```{.r .cell-code code-fold=\"true\"}\n#A: If the cases were isolated after 24 hours, then gamma would be 1/1 = 1, and if the cases were isolated after 12 hours, gamma would be 1/0.5 = 2. R0 would be calculated as the beta coefficient over gamma, below:\n\nr0_g1 <- 1.094913 / 1 + 1\nr0_g1\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] 2.094913\n```\n\n\n:::\n\n```{.r .cell-code code-fold=\"true\"}\nr0_g2 <- 1.094913 / 2 + 1\nr0_g2\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] 1.547457\n```\n\n\n:::\n:::\n\n\n### Exercise 3\n\n\n::: {.cell}\n\n```{.r .cell-code code-fold=\"true\"}\nniamey[5, 3] <- 0 #replace a \"NA\"\n#the command below organizes the data so it can be plotted and analyzed\nniamey <- data.frame(\n biweek = rep(seq(1, 16), 3),\n site = c(rep(1, 16), rep(2, 16), rep(3, 16)),\n cases = c(niamey[, 1], niamey[, 2], niamey[, 3])\n) #define \"biweeks\"\n\n# As the data are reported every two weeks, this corresponds to the 10th observation. Let’s fit a linear model\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code code-fold=\"true\"}\n# First let's see what the outbreak looks like for the first community on a linear scale\nplot(\n niamey$biweek[niamey$site == 1],\n niamey$cases[niamey$site == 1],\n type = 'p',\n col = niamey$site,\n xlab = 'Biweek',\n ylab = 'Cases'\n)\nlines(niamey$biweek[niamey$site == 1], niamey$cases[niamey$site == 1])\n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code code-fold=\"true\"}\n# Now let’s try it on a log scale to see until when the outbreak is roughly linear\nplot(\n niamey$biweek[niamey$site == 1],\n niamey$cases[niamey$site == 1],\n type = 'p',\n col = niamey$site,\n xlab = 'Biweek',\n ylab = 'Cases',\n log = 'y'\n)\nlines(niamey$biweek[niamey$site == 1], niamey$cases[niamey$site == 1])\n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code code-fold=\"true\"}\n# here we create a \"week\" variable to run the analysis on a weekly\nniamey$week <- niamey$biweek * 2\n# we use the `head` command to take the first N values of a vector. In this case, we're taking the first 10 values of our outcome (cases) and predictor (time) variables.\nmodel <- lm(\n log(head(niamey$cases[niamey$site == 1], 10)) ~\n (head(niamey$week[niamey$site == 1], 10))\n)\nsummary(model) #summary statistics for fit model\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n\nCall:\nlm(formula = log(head(niamey$cases[niamey$site == 1], 10)) ~ \n (head(niamey$week[niamey$site == 1], 10)))\n\nResiduals:\n Min 1Q Median 3Q Max \n-0.51796 -0.07364 0.00790 0.10727 0.36805 \n\nCoefficients:\n Estimate Std. Error t value Pr(>|t|)\n(Intercept) 2.59679 0.17563 14.79 4.31e-07\nhead(niamey$week[niamey$site == 1], 10) 0.21960 0.01415 15.52 2.96e-07\n \n(Intercept) ***\nhead(niamey$week[niamey$site == 1], 10) ***\n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n\nResidual standard error: 0.2571 on 8 degrees of freedom\nMultiple R-squared: 0.9678,\tAdjusted R-squared: 0.9638 \nF-statistic: 240.8 on 1 and 8 DF, p-value: 2.963e-07\n```\n\n\n:::\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code code-fold=\"true\"}\n# Now let's get the slope and display it\nslope <- coef(model)[2] #extract slope parameter\nslope #print to screen\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\nhead(niamey$week[niamey$site == 1], 10) \n 0.2196041 \n```\n\n\n:::\n:::\n\n\nAs we ran the model by weeks, the $\\gamma$ value is $2^{-1}$ and $\\hat R_0 = \\hat \\beta_1 / \\gamma +1$ giving $\\hat R_0=0.2196041/0.5+1 \\approx 1.44$.\n\n### Exercise 4\n\nHere, we can use a loop and repeat the regression procedure we used above for varying numbers of initial data points in our model.\n\n\n::: {.cell}\n\n```{.r .cell-code code-fold=\"true\"}\nslope <- NULL\nse <- NULL\nfor (i in 3:18) {\n model <- lm(\n log(head(niamey$cases[niamey$site == 1], i)) ~\n (head(niamey$week[niamey$site == 1], i))\n )\n slope <- c(slope, as.numeric(coef(model)[2]))\n se <- c(se, summary(model)$coefficients[4])\n}\nR0 <- slope / 0.5 + 1\nplot(R0, se, ylab = 'Standard error')\n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n:::\n\n\n### Exercise 5\n\nFirst, let's add in $\\gamma$ estimation into the `sse.sir` function and create a new `sse.sir` function called `sse.sir.g`\n\n\n::: {.cell}\n\n```{.r .cell-code}\n#' Basic SIR model\n#'\n#' A basic SIR model with no demographic structure to be used in deSolve\n#'\n#' @param time deSolve passes the time parameter to the function.\n#' @param state A vector of states.\n#' @param params The beta parameter\n#' @param ... Other arguments passed by deSolve.\n#'\n#' @return A deSolve matrix of states at each time step.\n#' @examples\n#' sir_params <- 0.0005\n#' sir_init_states <- c(S = 5000, I = 1, R = 0)\n#' sim_times <- seq(0, 16 / 365, by = 0.1 / 365)\n#'\n#' sir_sol <- deSolve::ode(\n#' y = sir_init_states,\n#' times = sim_times,\n#' func = closed_sir_model,\n#' parms = sir_params\n#' )\n\n# Create a new SIR model to include gamma\nclosed_sir_model_g <- function(time, state, params, ...) {\n # Unpack states\n S <- state[\"S\"]\n I <- state[\"I\"]\n\n # Unpack parameters\n beta <- params[[\"beta\"]]\n dur_inf <- 14 / 365\n gamma <- params[[\"gamma\"]]\n\n new_inf <- beta * S * I\n\n # Calculate the ODEs\n dSdt <- -new_inf\n dIdt <- new_inf - (gamma * I)\n\n # Return the ODEs\n return(list(c(dSdt, dIdt, new_inf)))\n}\n\n\n#' Calculate the Sum of Squared Errors\n#'\n#' A function to take in biweekly incidence data, and SIR parameters, and\n#' calculate the SSE\n#'\n#' @param params A vector of parameter values\n#' @param data A dataframe containing biweekly incidence data in the case column\n#'\n#' @return The SSE of type double\n#' @examples\nsse_sir_g <- function(params, data) {\n # Convert biweekly time series into annual time scale\n # Daily time scale has requires beta values to be too small - doesn't\n # optimize well\n dt <- 0.01\n max_biweek <- max(data$biweek)\n t <- seq(0, max_biweek * 14, dt) / 365\n\n # Extract the number of observed incidence\n obs_inc <- data$cases\n\n # Note the parameters are updated throughout the optimization process by\n # the optim() function\n # Unpack the transmission parameter and exponentiate to fit on ln scale\n beta <- exp(params[[\"beta\"]])\n gamma <- exp(params[[\"gamma\"]])\n\n in_parms <- c(\"beta\" = beta, \"gamma\" = gamma)\n # Unpack the initial states and exponentiate to fit on normal scale\n S_init <- exp(params[[\"S_init\"]])\n I_init <- exp(params[[\"I_init\"]])\n\n # Fit SIR model to the parameters\n sol <- deSolve::ode(\n y = c(S = S_init, I = I_init, new_inf = 0),\n times = t,\n func = closed_sir_model_g,\n parms = in_parms,\n # Use rk4 as fixed time steps, which is important for indexing\n method = \"rk4\"\n )\n\n # Extract the cumulative incidence\n cum_inc <- sol[, \"new_inf\"]\n\n # Find the indices of the cumulative incidence to extract\n biweek_index <- seq(1, max_biweek) * (14 / dt) + 1\n\n # Index cumulative incidence to get the values at the end of the biweeks\n biweek_cum_inc <- cum_inc[biweek_index]\n\n # Calculate the biweekly incidence by using the difference between\n # consecutive biweeks. Need to manually prepend first week's incidence\n # and add in the initial number of infectious individuals, as ODE model\n # only returns the cumulative differences, which is 0 at the start.\n biweek_inc <- c(biweek_cum_inc[1] + I_init, diff(biweek_cum_inc, lag = 1))\n\n # return SSE of predicted vs observed incidence\n return(sum((biweek_inc - obs_inc)^2))\n}\n```\n:::\n\n\nNow we can run the code to optimize, beta, gamma, I0, and S0\n\n\n::: {.cell}\n\n```{.r .cell-code}\nsse_optim_params_g <- c(\n beta = log(0.055),\n gamma = log(365 / 14),\n S_init = log(5000),\n I_init = log(1)\n)\n\n# Create a dataframe of optimized parameters\nniamey_optims_g <- niamey_df %>%\n # Create a nested dataframe i.e. one row for each site, and the data column\n # now is a list column that contains a separate dataframe of times and\n # cases for each site\n nest(data = -site) %>%\n mutate(\n # Map the optim() function call to each of the separate dataframes\n # stored in the nested data column we just created\n fit = map(data, ~ optim(sse_optim_params_g, sse_sir_g, data = .x)),\n # Map the exp() function to each of the model fits just created, and\n # output to a dataframe instead of a list (like in map()), for easier\n # use in the plottinge predictions later\n map_dfr(fit, ~ exp(.x$par))\n )\n```\n:::\n\n\nNow we can get our optimized parameters\n\n\n::: {.cell}\n\n```{.r .cell-code}\nniamey_optims_g %>%\n select(-c(data, fit)) %>%\n mutate(site = str_replace_all(site, \"_\", \" \")) %>%\n gt() %>%\n fmt_number(columns = -site, decimals = 3) %>%\n fmt_scientific(columns = beta, decimals = 3) %>%\n # Relabel the column headers\n cols_label(\n site = md(\"**Site**\"),\n beta = md(\"**Beta**\"),\n gamma = md(\"**Gamma**\"),\n S_init = md(\"**Initial S**\"),\n I_init = md(\"**Initial I**\")\n ) %>%\n # Apply style to the table with gray alternating rows\n opt_stylize(style = 1, color = 'gray') %>%\n # Increate space between columns\n opt_horizontal_padding(scale = 3) %>%\n cols_align(\"center\")\n```\n\n::: {.cell-output-display}\n\n```{=html}\n| Site | \nBeta | \nInitial S | \nInitial I | \n
|---|---|---|---|
| Site 1 | \n5.370 × 10−3 | \n8,566.648 | \n1.349 |
| Site 2 | \n8.317 × 10−3 | \n5,961.388 | \n0.201 |
| Site 3 | \n7.134 × 10−2 | \n792.990 | \n0.001 |
\n\n\n \n
\n
\n```\n\n:::\n:::\n\n\nGetting new predictions from the optimized parameters and a new model run for each site.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nniamey_predictions_g <- niamey_optims_g %>%\n mutate(\n # For each of the different site's nested dataframes, fit the SIR model\n # with the optimal parameters to get best fit predictions\n predictions = pmap(\n .l = list(\n S_init = S_init,\n I_init = I_init,\n beta = beta,\n gamma = gamma,\n time_data = data\n ),\n .f = function(S_init, I_init, beta, gamma, time_data) {\n site_times <- time_data$biweek * 14 / 365\n\n in_parms <- c(\"beta\" = beta, \"gamma\" = gamma)\n # Return a dataframe of model solutions\n as_tibble(ode(\n y = c(S = S_init, I = I_init, new_inf = 0),\n times = site_times,\n func = closed_sir_model_g,\n parms = in_parms,\n hmax = 1 / 120\n )) %>%\n # Make sure all values are numeric for plotting purposes\n mutate(across(everything(), as.numeric)) %>%\n mutate(\n incidence = ifelse(\n row_number() == 1,\n new_inf[1],\n diff(new_inf, lag = 1)\n )\n )\n }\n )\n ) %>%\n unnest(c(data, predictions))\n```\n:::\n\n\nFinally we can plot:\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Create a dataframe to store the positions of the text labels\nniamey_preds_labels <- tibble(\n site = c(\"Site_1\", \"Site_2\"),\n x_label = c(6.5, 6.5),\n x_arrow_just = c(-0.5, -0.5),\n x_arrow_end = c(7, 7.75),\n y_label = c(900, 600),\n y_arrow_just = c(-80, -70),\n y_arrow_end = c(350, 290),\n commentary = c(\"**Predicted\", \"**Observed\"),\n color = c(\"grey20\", niamey_site_colors[\"Site_2\"])\n)\n\nggplot(niamey_predictions_g, aes(x = biweek, group = site)) +\n # Plot the actual data in color\n geom_line(aes(y = cases, color = site)) +\n geom_point(aes(y = cases, color = site), size = 4, alpha = 0.8) +\n # Plot the best-fit model predictions in black\n geom_line(aes(y = incidence), color = \"black\") +\n scale_color_manual(\n values = niamey_site_colors,\n aesthetics = c(\"color\", \"fill\")\n ) +\n # Place each site on it's own subplot and change labels\n facet_wrap(\n ~site,\n ncol = 1,\n scales = \"free_y\",\n labeller = as_labeller(niamey_site_labels)\n ) +\n labs(x = \"Biweek\", y = \"Number of Case\") +\n theme(legend.position = \"none\") +\n ggtext::geom_textbox(\n data = niamey_preds_labels,\n aes(\n label = paste0(\n \"\",\n commentary,\n \" Cases**\",\n \"\"\n ),\n x = x_label,\n y = y_label\n ),\n size = 4,\n fill = NA,\n box.colour = NA\n ) +\n geom_curve(\n data = niamey_preds_labels,\n aes(\n x = x_label + x_arrow_just,\n xend = x_arrow_end,\n y = y_label + y_arrow_just,\n yend = y_arrow_end\n ),\n linewidth = 0.75,\n arrow = arrow(length = unit(0.2, \"cm\")),\n curvature = list(0.25),\n color = \"grey20\"\n )\n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n:::\n\n\n\n::: {.cell}\n\n:::\n\n\n\n::: {.cell}\n\n:::\n\n\n\n::: {.cell}\n\n:::\n\n\n### Exercise 6\n\nWhat happens if one or both of the other unknowns ($X_0$ and $Y_0$) is fixed instead of $\\gamma$?\n\n#### Solutions coming soon!\n\nFirst we modify the `sse_sir_g` function again to fix $X_0$ and $Y_0$ (or both)\n\n\n::: {.cell}\n\n:::\n\n\nNow we can run the optim algorithm and find the best parameters for the SIR model for each of the scenario.\n\nFirst, for S0 fixed:\n\n\n::: {.cell}\n\n:::\n\n\nSecond, for I0 fixed:\n\n\n::: {.cell}\n\n:::\n\n\nThird, for fixing both S0 and I0:\n\n\n::: {.cell}\n\n:::\n\n\n\n::: {.cell}\n\n:::\n\n\nNow let's pass these models through `optim` and fit the models:\n\n\n::: {.cell}\n\n:::\n\n",
+ "markdown": "---\nsubtitle: \"Estimation\"\nabstract-title: \"\"\nabstract: |\n *Materials adapted from John Drake and Pej Rohani [@drakeEstimation2019]*\nexecute:\n warning: false\nmetadata-files:\n - metadata/matthewferrari.yml\n - metadata/mathjax-packages.yml\neditor:\n markdown:\n wrap: sentence\n---\n\n# R Session 03\n\n## Setup\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(here)\nlibrary(rio)\nlibrary(deSolve)\nlibrary(tidyverse)\nlibrary(ggtext)\nlibrary(gt)\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\ntheme_set(theme_minimal())\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Loads the datasets: flu, measles, niamey, plauge\nflu <- rio::import(\n \"https://raw.githubusercontent.com/arnold-c/SISMID-Module-02_2025/main/data/flu.csv\"\n)\nniamey <- rio::import(\n \"https://raw.githubusercontent.com/arnold-c/SISMID-Module-02_2025/main/data/niamey.csv\"\n)\n# flu <- rio::import(here::here(\"data\", \"flu.csv\"))\n# niamey <- rio::import(here::here(\"data\", \"niamey.csv\"))\n```\n:::\n\n\n## Estimating $R_0$ Problem Background\n\nSo far in this class we have focused on the **theory** of infectious disease.\nOften, however, we will want to apply this theory to particular situations.\nOne of the key applied problems in epidemic modeling is the estimation of $R_0$ from outbreak data.\nIn this session, we study two methods for estimating $R_0$ from an epidemic curve.\nAs a running example, we will use the data on influenza in a British boarding school.\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nggplot(flu, aes(x = day, y = flu)) +\n geom_line(color = \"slategray4\") +\n geom_point(shape = 21, size = 5, fill = \"slategray4\", alpha = 0.8) +\n labs(x = \"Day\", y = \"Active Influenza Cases\")\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\nThe figure above, and the corresponding tabulated outbreak data, report the number of **active influenza cases** on each day.\nThese are point-in-time counts, or prevalence: the number of students who were sick at that moment.\nBy contrast, the final-size calculation below uses cumulative **incidence**: the total number of students who became cases at any point during the outbreak.\nFor this boarding school outbreak, the total school size was $N = 764$, and the total number of students who became cases over the full outbreak was $Z = 512$.\nBecause $512$ is a cumulative incidence count, it is not obtained by reading a single value from the prevalence figure or table.\n\n## Estimating $R_0$ From The Final Outbreak Size\n\nOur first approach is to estimate $R_0$ from the final outbreak size.\nAlthough unhelpful at the early stages of an epidemic (before the final epidemic size is observed), this method is nonetheless a useful tool for *post hoc* analysis.\nThe method is general and can be motivated by the argument listed in [@keelingIntroductionSimpleEpidemic2008]:\n\nFirst, we assume that the epidemic is started by a single infectious individual in a completely susceptible population.\nOn average, this individual infects $R_0$ others.\nThe probability a particular individual escaped infection is therefore $e^{-R_0 / N}$.\n\nIf $Z$ individuals have been infected, the probability of an individual escaping infection from all potential sources is $e^{-Z R_0 / N}$.\nIt follows that at the end of the epidemic a proportion $R(\\infty) = Z / N$ have been infected and the fraction remaining susceptible is $S(\\infty) = e^{-R(\\infty) R_0}$, which is equal to $1 - R(\\infty)$.\n\n::: {.callout-note collapse=\"true\"}\n$S(\\infty) = e^{-R(\\infty) R_0}$ can be calculated by acknowledging that at equilibrium ($t = \\infty$), $S(\\infty) = 1 - R(\\infty) = Z / N$, so substituting $R(\\infty)$ into $1 - e^{-Z R_0 / N}$ gives the desired result.\n\nIt could also be calculated by dividing $\\frac{\\dd{S}}{\\dd{t}}$ by $\\frac{\\dd{R}}{\\dd{t}}$:\n\n```{=tex}\n\\begin{aligned}\n\\frac{\\dd{S}}{\\dd{R}} &= - \\frac{\\beta S}{\\gamma} \\\\\n&= - R_0 S\n\\end{aligned}\n```\nwhich is a [separable differential equation](https://tutorial.math.lamar.edu/classes/de/separable.aspx), so can be integrated as follows:\n\n```{=tex}\n\\begin{aligned}\n- \\int_{0}^{t} \\frac{1}{R_0 S} \\dd{S} &= \\int_{0}^{t} \\dd{R} \\\\\n- \\frac{1}{R_0} \\left(\\ln{S(t)} - \\ln{S(0)} \\right) &= R(t) - \\cancelto{0}{R(0)} \\\\\n\\ln{S(t)} &= \\ln{S(0)} - R_0 R(t) \\\\\nS(t) &= S(0) e^{-R_0 R(t)}\n\n\\end{aligned}\n```\n:::\n\nPutting this together, we get:\n\n$$\n1 - R(\\infty) - e^{-R(\\infty) R_0} = 0\n$$\n\nRearranging, we have the estimator\n\n$$\n \\hat{R_0} = \\frac{\\ln(1 - Z / N)}{-Z / N},\n$$\n\nwhich, in this case, evaluates to $\\frac{\\ln(1 - 512 / 764)}{-512 / 764} = 1.655$.\n\n::: exercise\n### Exercise 1\n:::\n\nThis equation shows the important one-to-one relationship between $R_0$ and the final epidemic size.\nPlot the relationship between the total epidemic size and $R_0$ for the complete range of values between 0 and 1.\n\n## Linear Approximation\n\nThe next method we introduce takes advantage of the fact that during the early stages of an outbreak, the number of infected individuals is given approximately as $I(t) \\approx I_0 e^{((R_0 - 1)(\\gamma + \\mu)t)}$.\nTaking logarithms of both sides, we have $\\ln(I(t)) \\approx \\ln(I_0) + (R_0 - 1)(\\gamma + \\mu)t$, showing that the log of the number of infected individuals is approximately linear in time with a slope that reflects both $R_0$ and the recovery rate.\n\nThis suggests that a simple linear regression fit to the first several data points on a log-scale, corrected to account for $\\gamma$ and $\\mu$, provides a rough and ready estimate of $R_0$.\nFor flu, we can assume $\\mu =0$ because the epidemic occurred over a time period during which natural mortality is negligible.\nFurther, assuming an infectious period of about 2.5 days, we use $\\gamma = (2.5)^{-1} = 0.4$ for the correction.\nFitting to the first four data points, we obtain the slope as follows.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Fit a linear model\nlinear_model <- lm(log(flu[1:4]) ~ day[1:4], data = flu)\n\n# Summary statistics for fit model\nsummary(linear_model)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n\nCall:\nlm(formula = log(flu[1:4]) ~ day[1:4], data = flu)\n\nResiduals:\n 1 2 3 4 \n 0.03073 -0.08335 0.07450 -0.02188 \n\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) -0.02703 0.10218 -0.265 0.81611 \nday[1:4] 1.09491 0.03731 29.346 0.00116 **\n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n\nResidual standard error: 0.08343 on 2 degrees of freedom\nMultiple R-squared: 0.9977,\tAdjusted R-squared: 0.9965 \nF-statistic: 861.2 on 1 and 2 DF, p-value: 0.001159\n```\n\n\n:::\n\n```{.r .cell-code}\n# Extract slope parameter\ncoef(linear_model)[2]\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\nday[1:4] \n1.094913 \n```\n\n\n:::\n:::\n\n\nRearranging the linear equation above and denoting the slope coefficient by $\\hat \\beta_1$ we have the estimator $\\hat R_0 = \\hat \\beta_1 / \\gamma + 1$ giving $\\hat R_0 = 1.094913 / 0.4 + 1 \\approx 3.7$.\n\n::: exercise\n### Exercise 2\n:::\n\nOur estimate assumes that boys remained infectious during the natural course of infection.\nThe original report on this epidemic indicates that boys found to have symptoms were immediately confined to bed in the infirmary.\nThe report also indicates that only 1 out of 130 adults at the school exhibited any symptoms.\nIt is reasonable, then, to suppose that transmission in each case ceased once he had been admitted to the infirmary.\nSupposing admission happened within 24 hours of the onset of symptoms.\nHow does this affect our estimate of $R_0$?\nTwelve hours?\n\n::: exercise\n### Exercise 3\n:::\n\nBiweekly data for outbreaks of measles in three communities in Niamey, Niger are provided in the dataframe `niamey`.\nUse this method to obtain estimates of $R_0$ for measles from the first community assuming that the infectious period is approximately two weeks or $\\frac{14}{365} \\approx 0.0384$ years.\n\n::: exercise\n### Exercise 4\n:::\n\nA defect with this method is that it uses only a small fraction of the information that might be available, i.e., the first few data points.\nIndeed, there is nothing in the method that tells one how many data points to use--this is a matter of judgment.\nFurther, there is a tradeoff in that as more and more data points are used the precision of the estimate increases, but this comes at a cost of additional bias.\nPlot the estimate of $R_0$ obtained from $n=3, 4, 5, ...$ data points against the standard error of the slope from the regression analysis to show this tradeoff.\n\n## Estimating dynamical parameters with least squares\n\nThe objective of the previous exercise was to estimate $R_0$.\nKnowing $R_0$ is critical to understanding the dynamics of any epidemic system.\nIt is, however, a composite quantity and is not sufficient to completely describe the epidemic trajectory.\nFor this, we require estimates for all parameters of the model.\nIn this exercise, we introduce a simple approach to model estimation called **least squares fitting**, sometimes called **trajectory matching**.\nThe basic idea is that we find the values of the model parameters that minimize the squared differences between model predictions and the observed data.\nTo demonstrate least squares fitting, we consider an outbreak of measles in Niamey, Niger, reported on by [@graisEstimatingTransmissionIntensity2006].\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Replace an \"NA\"\nniamey[5, 3] <- 0\n\nniamey_df <- niamey %>%\n # Rename columns to remove automatic \"V1\" etc columns names\n rename_with(., ~ paste0(\"Site_\", str_remove(.x, \"V\"))) %>%\n # Add a column for the biweekly time period\n mutate(biweek = 1:16) %>%\n # Convert to long format for plotting\n pivot_longer(\n cols = contains(\"Site\"),\n names_to = \"site\",\n values_to = \"cases\"\n )\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Create a vector of colors for each site in the Niamey dataset\nniamey_site_colors <- RColorBrewer::brewer.pal(3, \"Dark2\")\n# Assign names to the colors\nnames(niamey_site_colors) <- unique(niamey_df$site)\n\n# Create a vector of labels for each site for nicer plotting legends\nniamey_site_labels <- str_replace_all(names(niamey_site_colors), \"_\", \" \")\nnames(niamey_site_labels) <- names(niamey_site_colors)\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nggplot(\n niamey_df,\n aes(x = biweek, y = cases, color = site, fill = site, group = site)\n) +\n geom_line() +\n geom_point(shape = 21, size = 5, alpha = 0.8) +\n scale_color_manual(\n values = niamey_site_colors,\n aesthetics = c(\"color\", \"fill\"),\n labels = niamey_site_labels\n ) +\n guides(color = \"none\") +\n labs(x = \"Biweek\", y = \"Number of Cases\", fill = \"Site Number\") +\n theme(legend.position = \"bottom\")\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\n## Dynamical Model\n\nFirst, we write a specialized function for simulating the SIR model in a case where the removal rate is *\"hard-wired\"* and with no demography.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n#' Basic SIR model\n#'\n#' A basic SIR model with no demographic structure to be used in deSolve\n#'\n#' @param time deSolve passes the time parameter to the function.\n#' @param state A vector of states.\n#' @param params The beta parameter\n#' @param ... Other arguments passed by deSolve.\n#'\n#' @return A deSolve matrix of states at each time step.\n#' @examples\n#' sir_params <- 0.0005\n#' sir_init_states <- c(S = 5000, I = 1, R = 0)\n#' sim_times <- seq(0, 16 / 365, by = 0.1 / 365)\n#'\n#' sir_sol <- deSolve::ode(\n#' y = sir_init_states,\n#' times = sim_times,\n#' func = closed_sir_model,\n#' parms = sir_params\n#' ))\nclosed_sir_model <- function(time, state, params, ...) {\n # Unpack states\n S <- state[\"S\"]\n I <- state[\"I\"]\n\n # Unpack parameters\n beta <- params\n dur_inf <- 14 / 365\n gamma <- 1 / dur_inf\n\n new_inf <- beta * S * I\n\n # Calculate the ODEs\n dSdt <- -new_inf\n dIdt <- new_inf - (gamma * I)\n\n # Return the ODEs\n return(list(c(dSdt, dIdt, new_inf)))\n}\n```\n:::\n\n\n## Interactive Optimization\n\n\n::: {.cell}\n\n:::\n\n\n\n::: {.cell}\n\n:::\n\n\n```{ojs}\n//| echo: false\nfiltered_niamey_data = aq.table(niamey_data)\n .filter(aq.escape(d => d.site == site_select))\n```\n\n```{ojs}\n//| echo: false\nreset_S = 10000\nreset_I = 20\nreset_beta = 5.00\n```\n\n```{ojs}\n//| echo: false\nfunction set(input, value) {\n input.value = value;\n input.dispatchEvent(new Event(\"input\", {bubbles: true}));\n}\n```\n\n```{ojs}\n//| echo: false\nfunction sse(obs, preds) {\n\n if(obs.length == preds.length) {\n var squared_errs = obs.map((e, i) => (e - preds[i])**2 )\n return squared_errs.reduce((a, b) => a + b, 0)\n } else {\n return(\"lengths are not the same\")\n }\n}\n```\n\n```{ojs}\n//| panel: sidebar\n//| echo: false\nviewof reset = Inputs.button([\n [\"Reset all sliders\", () => {\n set(viewof S0, reset_S)\n set(viewof I0, reset_I)\n set(viewof beta_input, reset_beta)\n }]\n])\nviewof S0 = Inputs.range(\n [500, 15000],\n {value: reset_S, step: 1, label: md`${tex`S(0)`}`}\n)\n\nviewof I0 = Inputs.range(\n [0.001, 50],\n {value: reset_I, step: 0.001, label: md`${tex`I(0)`}`}\n)\n\nviewof beta_input = Inputs.range(\n [1, 100],\n {value: reset_beta, step: 0.001, label: md`${tex`\\beta (\\times 10^{-3})`}`}\n)\n\nviewof site_select = Inputs.select(\n [\"Site 1\", \"Site 2\", \"Site 3\"],\n {label: \"Select a site:\"}\n)\n\n// convert to daily time scale as easier to manipulate\nbeta = (beta_input / 365) * (10 ** (-3))\n\nmd`${tex`R_0 = ${R0_str}`}`\n```\n\n```{ojs}\n//| echo: false\ndur_inf = 14\ngamma = 1 / dur_inf\nR0 = beta * (S0 + I0)/ gamma\n\nR0_str = R0.toPrecision(2).toLocaleString()\n```\n\n```{ojs}\n//| echo: false\ndt = 0.01\ntmax = 16 * 14\n```\n\n```{ojs}\n//| echo: false\nimport {odeRK4} from '@rreusser/integration@3064'\nimport { aq, op } from '@uwdata/arquero'\n```\n\n```{ojs}\n//| echo: false\nfunction sir(dydt, y, t) {\n dydt[3] = beta * y[0] * y[1]\n\n dydt[0] = - dydt[3]\n dydt[1] = dydt[3] - gamma * y[1]\n dydt[2] = gamma * y[1]\n}\n```\n\n```{ojs}\n//| echo: false\nfunction simulate(f, t0, y0, dt, tmax) {\n var t = t0\n var y = y0\n var i = 0\n\n var ysim = [y0]\n\n for (t = t0 + dt; t <= tmax; t += dt) {\n ysim.push(odeRK4([], ysim[i], f, dt))\n i += 1\n }\n\n // return cumulative infections\n return ysim.map(d => d[3])\n}\n```\n\n```{ojs}\n//| echo: false\nsir_sol = simulate(sir, 0, [S0, I0, 0.0, 0.0], dt, tmax)\n```\n\n```{ojs}\n//| echo: false\nsiteColors = [\"#1b9e77\", \"#d95f02\", \"#7570b3\"]\n```\n\n```{ojs}\n//| echo: false\n// Create an array of time indices, dropping first as from returns index 0, which\n// doesn't exist in the real data\ntimes = Array.from({length: 17}, (_, i) => i * 14).slice(1)\ntindex = times.map((e, i) => e * (1 / dt))\n\ncum_inc = tindex.map((i) => sir_sol[i])\n\npreds = [cum_inc[0] + I0, ...cum_inc.map((e, i) => cum_inc[i] - cum_inc[i-1]).slice(1)]\n\nsir_tbl = aq.table({\n Biweek: times.map(t => t / 14),\n \"Cumulative Incidence\": cum_inc,\n \"Number of Individuals\": preds\n})\n\nsim_sse = [({\n sse: sse(\n filtered_niamey_data.array(\"Number of Individuals\"),\n preds\n ).toPrecision(4),\n Biweek: 3,\n \"Number of Individuals\": Math.max(\n ...sir_tbl.array(\"Number of Individuals\"),\n ...filtered_niamey_data.array(\"Number of Individuals\")\n ) * 0.9\n})]\n```\n\n```{ojs}\n//| echo: false\n//| panel: fill\nPlot.plot({\n color: {\n legend: true,\n domain: [\"Site 1\", \"Site 2\", \"Site 3\"],\n range: siteColors,\n },\n style: {fontSize: \"20px\"},\n marginLeft: 75,\n marginTop: 40,\n marginBottom: 55,\n grid: true,\n width: 800,\n height: 670,\n marks: [\n Plot.lineY(\n sir_tbl,\n {x: \"Biweek\", y: \"Number of Individuals\", stroke: \"#4d4d4dff\", strokeWidth: 6, strokeOpacity: 0.8}\n ),\n Plot.dot(\n filtered_niamey_data,\n {x: \"Biweek\", y: \"Number of Individuals\", stroke: \"site\", fill: \"site\", fillOpacity: 0.6, r: 12}\n ),\n Plot.lineY(\n filtered_niamey_data,\n {x: \"Biweek\", y: \"Number of Individuals\", stroke: \"site\"}\n ),\n Plot.text(\n sim_sse,\n {x: \"Biweek\", y: \"Number of Individuals\", text: (d) => `SSE = ${d.sse}`, dx: 0, dy: 0, fontWeight: \"bold\"}\n )\n ]\n})\n```\n\n## Objective Function\n\nNow we set up a function that will calculate the sum of the squared differences between the observations and the model at any parameterization (more commonly known as *\"sum of squared errors\"*).\nIn general, this is called the **objective function** because it is the quantity that optimization seeks to minimize.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n#' Calculate the Sum of Squared Errors\n#'\n#' A function to take in biweekly incidence data, and SIR parameters, and\n#' calculate the SSE\n#'\n#' @param params A vector of parameter values\n#' @param data A dataframe containing biweekly incidence data in the case column\n#'\n#' @return The SSE of type double\n#' @examples\nsse_sir <- function(params, data) {\n # Convert biweekly time series into annual time scale\n # Daily time scale has requires beta values to be too small - doesn't\n # optimize well\n dt <- 0.01\n max_biweek <- max(data$biweek)\n t <- seq(0, max_biweek * 14, dt) / 365\n\n # Extract the number of observed incidence\n obs_inc <- data$cases\n\n # Note the parameters are updated throughout the optimization process by\n # the optim() function\n # Unpack the transmission parameter and exponentiate to fit on ln scale\n beta <- exp(params[[\"beta\"]])\n\n # Unpack the initial states and exponentiate to fit on normal scale\n S_init <- exp(params[[\"S_init\"]])\n I_init <- exp(params[[\"I_init\"]])\n\n # Fit SIR model to the parameters\n sol <- deSolve::ode(\n y = c(S = S_init, I = I_init, new_inf = 0),\n times = t,\n func = closed_sir_model,\n parms = beta,\n # Use rk4 as fixed time steps, which is important for indexing\n method = \"rk4\"\n )\n\n # Extract the cumulative incidence\n cum_inc <- sol[, \"new_inf\"]\n\n # Find the indices of the cumulative incidence to extract\n biweek_index <- seq(1, max_biweek) * (14 / dt) + 1\n\n # Index cumulative incidence to get the values at the end of the biweeks\n biweek_cum_inc <- cum_inc[biweek_index]\n\n # Calculate the biweekly incidence by using the difference between\n # consecutive biweeks. Need to manually prepend first week's incidence\n # and add in the initial number of infectious individuals, as ODE model\n # only returns the cumulative differences, which is 0 at the start.\n biweek_inc <- c(biweek_cum_inc[1] + I_init, diff(biweek_cum_inc, lag = 1))\n\n # return SSE of predicted vs observed incidence\n return(sum((biweek_inc - obs_inc)^2))\n}\n```\n:::\n\n\nNotice that the code for `sse_sir()` makes use of the following modeling trick.\nWe know that $\\beta$, $S_0$, and $I_0$ must be positive, but our search to optimize these parameters will be over the entire number line.\nWe could constrain the search using a more sophisticated algorithm, but this might introduce other problems (i.e., stability at the boundaries).\nInstead, we parameterize our objective function (`sse_sir`) in terms of some alternative variables $\\ln(\\beta)$, $\\ln(S_0)$, and $\\ln(I_0)$.\nWhile these numbers range from $-\\infty$ to $\\infty$ (the range of our search) they map to our model parameters on a range from $0$ to $\\infty$ (the range that is biologically meaningful).\n\n## Optimization\n\nOur final step is to use the function `optim` to find the values of $\\beta$, $S_0$, and $I_0$ that minimize the sum of squared errors as calculated using our function.\n\nFinally, we plot these fits against the data.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Initial guess\nsse_optim_params <- c(beta = log(0.055), S_init = log(5000), I_init = log(1))\n\n# Create a dataframe of optimized parameters\nniamey_optims <- niamey_df %>%\n # Create a nested dataframe i.e. one row for each site, and the data column\n # now is a list column that contains a separate dataframe of times and\n # cases for each site\n nest(data = -site) %>%\n mutate(\n # Map the optim() function call to each of the separate dataframes\n # stored in the nested data column we just created\n fit = map(data, ~ optim(sse_optim_params, sse_sir, data = .x)),\n # Map the exp() function to each of the model fits just created, and\n # output to a dataframe instead of a list (like in map()), for easier\n # use in the plottinge predictions later\n map_dfr(fit, ~ exp(.x$par))\n )\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nniamey_optims %>%\n select(-c(data, fit)) %>%\n mutate(site = str_replace_all(site, \"_\", \" \")) %>%\n gt() %>%\n fmt_number(columns = -site, decimals = 3) %>%\n fmt_scientific(columns = beta, decimals = 3) %>%\n # Relabel the column headers\n cols_label(\n site = md(\"**Site**\"),\n beta = md(\"**Beta**\"),\n S_init = md(\"**Initial S**\"),\n I_init = md(\"**Initial I**\")\n ) %>%\n # Apply style to the table with gray alternating rows\n opt_stylize(style = 1, color = 'gray') %>%\n # Increate space between columns\n opt_horizontal_padding(scale = 3) %>%\n cols_align(\"center\")\n```\n\n::: {.cell-output-display}\n\n```{=html}\n| Site | \nBeta | \nGamma | \nInitial S | \nInitial I | \n
|---|---|---|---|---|
| Site 1 | \n5.854 × 10−3 | \n81.150 | \n17,113.450 | \n0.786 |
| Site 2 | \n9.066 × 10−3 | \n77.613 | \n11,030.375 | \n0.142 |
| Site 3 | \n8.517 × 10−2 | \n351.697 | \n4,462.286 | \n0.000 |
\n\n\n \n
\n
\n```\n\n:::\n:::\n\n\n::: callout-note\nYou may have noticed that you can achieve slightly different results for the optimal parameter values using the interactive plot than are being presented here (though they are very similar).\nThis is because while the optimization code is running in `R`, the interactive plot and the calculation of the SSE is implemented using JavaScript.\nTherefore, despite using the same underlying model structure, the answers will vary slightly, because the ODE solvers are different, resulting in different model simulations.\nThe difference is not enough to be concerned with here, but it is a point that's worth being aware of when you build your own models - you may want to perform sensitivity to confirm that your model implementation is not driving the magnitude of the results you see, and the inferences you make.\n:::\n\n\n::: {.cell}\n\n```{.r .cell-code}\nniamey_predictions <- niamey_optims %>%\n mutate(\n # For each of the different site's nested dataframes, fit the SIR model\n # with the optimal parameters to get best fit predictions\n predictions = pmap(\n .l = list(\n S_init = S_init,\n I_init = I_init,\n beta = beta,\n time_data = data\n ),\n .f = function(S_init, I_init, beta, time_data) {\n site_times <- time_data$biweek * 14 / 365\n\n # Return a dataframe of model solutions\n as_tibble(ode(\n y = c(S = S_init, I = I_init, new_inf = 0),\n times = site_times,\n func = closed_sir_model,\n parms = beta,\n hmax = 1 / 120\n )) %>%\n # Make sure all values are numeric for plotting purposes\n mutate(across(everything(), as.numeric)) %>%\n mutate(\n incidence = ifelse(\n row_number() == 1,\n new_inf[1],\n diff(new_inf, lag = 1)\n )\n )\n }\n )\n ) %>%\n unnest(c(data, predictions))\n```\n:::\n\n\n::: callout-important\nAn important point to note is that our data is biweekly **incidence** (new cases in time period), whereas out SIR model produces **prevalence** (total cases at any time point).\nTo account for this, our SIR model returns the cumulative incidence (line 32 of the model code), and our objective function extracts the biweekly incidence (lines 41-54), to ensure we are fitting the same data!\nThis is a common source of error in interpretation when people fit models to data.\n:::\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\n# Create a dataframe to store the positions of the text labels\nniamey_preds_labels <- tibble(\n site = c(\"Site_1\", \"Site_2\"),\n x_label = c(6.5, 6.5),\n x_arrow_just = c(-0.5, -0.5),\n x_arrow_end = c(7, 7.75),\n y_label = c(900, 600),\n y_arrow_just = c(-80, -70),\n y_arrow_end = c(350, 290),\n commentary = c(\"**Predicted\", \"**Observed\"),\n color = c(\"grey20\", niamey_site_colors[\"Site_2\"])\n)\n\nggplot(niamey_predictions, aes(x = biweek, group = site)) +\n # Plot the actual data in color\n geom_line(aes(y = cases, color = site)) +\n geom_point(aes(y = cases, color = site), size = 4, alpha = 0.8) +\n # Plot the best-fit model predictions in black\n geom_line(aes(y = incidence), color = \"black\") +\n scale_color_manual(\n values = niamey_site_colors,\n aesthetics = c(\"color\", \"fill\")\n ) +\n # Place each site on it's own subplot and change labels\n facet_wrap(\n ~site,\n ncol = 1,\n scales = \"free_y\",\n labeller = as_labeller(niamey_site_labels)\n ) +\n labs(x = \"Biweek\", y = \"Number of Case\") +\n theme(legend.position = \"none\") +\n ggtext::geom_textbox(\n data = niamey_preds_labels,\n aes(\n label = paste0(\n \"\",\n commentary,\n \" Cases**\",\n \"\"\n ),\n x = x_label,\n y = y_label\n ),\n size = 4,\n fill = NA,\n box.colour = NA\n ) +\n geom_curve(\n data = niamey_preds_labels,\n aes(\n x = x_label + x_arrow_just,\n xend = x_arrow_end,\n y = y_label + y_arrow_just,\n yend = y_arrow_end\n ),\n linewidth = 0.75,\n arrow = arrow(length = unit(0.2, \"cm\")),\n curvature = list(0.25),\n color = \"grey20\"\n )\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\n::: exercise\n### Exercise 5\n:::\n\nTo make things easier, we have assumed the infectious period is known to be 14 days.\nIn terms of years, $\\text{D} = \\frac{14}{365} \\approx 0.0384$, and the recovery rate is the inverse i.e., $\\gamma = \\frac{14}{365}$.\nNow, modify the code above to estimate $\\gamma$ and $\\beta$ simultaneously.\n\n::: exercise\n### Exercise 6\n:::\n\nWhat happens if one or both of the other unknowns ($S_0$ and $I_0$) is fixed instead of $\\gamma$?\n\n## Solutions\n\n### Exercise 1\n\n\n::: {.cell}\n\n```{.r .cell-code code-fold=\"true\"}\np_infec <- (seq(0, 1, by = 0.001))\nr0_p <- (log(1 - p_infec)) / (-p_infec)\nplot(\n x = p_infec,\n y = r0_p,\n main = \"Relationship between final proportion infected and R0\",\n xlab = \"Final proportion infected\",\n ylab = \"R0\"\n)\n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n:::\n\n\n### Exercise 2\n\n\n::: {.cell}\n\n```{.r .cell-code code-fold=\"true\"}\n#A: If the cases were isolated after 24 hours, then gamma would be 1/1 = 1, and if the cases were isolated after 12 hours, gamma would be 1/0.5 = 2. R0 would be calculated as the beta coefficient over gamma, below:\n\nr0_g1 <- 1.094913 / 1 + 1\nr0_g1\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] 2.094913\n```\n\n\n:::\n\n```{.r .cell-code code-fold=\"true\"}\nr0_g2 <- 1.094913 / 2 + 1\nr0_g2\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] 1.547457\n```\n\n\n:::\n:::\n\n\n### Exercise 3\n\n\n::: {.cell}\n\n```{.r .cell-code code-fold=\"true\"}\nniamey[5, 3] <- 0 #replace a \"NA\"\n#the command below organizes the data so it can be plotted and analyzed\nniamey <- data.frame(\n biweek = rep(seq(1, 16), 3),\n site = c(rep(1, 16), rep(2, 16), rep(3, 16)),\n cases = c(niamey[, 1], niamey[, 2], niamey[, 3])\n) #define \"biweeks\"\n\n# As the data are reported every two weeks, this corresponds to the 10th observation. Let’s fit a linear model\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code code-fold=\"true\"}\n# First let's see what the outbreak looks like for the first community on a linear scale\nplot(\n niamey$biweek[niamey$site == 1],\n niamey$cases[niamey$site == 1],\n type = 'p',\n col = niamey$site,\n xlab = 'Biweek',\n ylab = 'Cases'\n)\nlines(niamey$biweek[niamey$site == 1], niamey$cases[niamey$site == 1])\n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code code-fold=\"true\"}\n# Now let’s try it on a log scale to see until when the outbreak is roughly linear\nplot(\n niamey$biweek[niamey$site == 1],\n niamey$cases[niamey$site == 1],\n type = 'p',\n col = niamey$site,\n xlab = 'Biweek',\n ylab = 'Cases',\n log = 'y'\n)\nlines(niamey$biweek[niamey$site == 1], niamey$cases[niamey$site == 1])\n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code code-fold=\"true\"}\n# here we create a \"week\" variable to run the analysis on a weekly\nniamey$week <- niamey$biweek * 2\n# we use the `head` command to take the first N values of a vector. In this case, we're taking the first 10 values of our outcome (cases) and predictor (time) variables.\nmodel <- lm(\n log(head(niamey$cases[niamey$site == 1], 10)) ~\n (head(niamey$week[niamey$site == 1], 10))\n)\nsummary(model) #summary statistics for fit model\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n\nCall:\nlm(formula = log(head(niamey$cases[niamey$site == 1], 10)) ~ \n (head(niamey$week[niamey$site == 1], 10)))\n\nResiduals:\n Min 1Q Median 3Q Max \n-0.51796 -0.07364 0.00790 0.10727 0.36805 \n\nCoefficients:\n Estimate Std. Error t value Pr(>|t|)\n(Intercept) 2.59679 0.17563 14.79 4.31e-07\nhead(niamey$week[niamey$site == 1], 10) 0.21960 0.01415 15.52 2.96e-07\n \n(Intercept) ***\nhead(niamey$week[niamey$site == 1], 10) ***\n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n\nResidual standard error: 0.2571 on 8 degrees of freedom\nMultiple R-squared: 0.9678,\tAdjusted R-squared: 0.9638 \nF-statistic: 240.8 on 1 and 8 DF, p-value: 2.963e-07\n```\n\n\n:::\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code code-fold=\"true\"}\n# Now let's get the slope and display it\nslope <- coef(model)[2] #extract slope parameter\nslope #print to screen\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\nhead(niamey$week[niamey$site == 1], 10) \n 0.2196041 \n```\n\n\n:::\n:::\n\n\nAs we ran the model by weeks, the $\\gamma$ value is $2^{-1}$ and $\\hat R_0 = \\hat \\beta_1 / \\gamma +1$ giving $\\hat R_0=0.2196041/0.5+1 \\approx 1.44$.\n\n### Exercise 4\n\nHere, we can use a loop and repeat the regression procedure we used above for varying numbers of initial data points in our model.\n\n\n::: {.cell}\n\n```{.r .cell-code code-fold=\"true\"}\nslope <- NULL\nse <- NULL\nfor (i in 3:18) {\n model <- lm(\n log(head(niamey$cases[niamey$site == 1], i)) ~\n (head(niamey$week[niamey$site == 1], i))\n )\n slope <- c(slope, as.numeric(coef(model)[2]))\n se <- c(se, summary(model)$coefficients[4])\n}\nR0 <- slope / 0.5 + 1\nplot(R0, se, ylab = 'Standard error')\n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n:::\n\n\n### Exercise 5\n\nFirst, let's add in $\\gamma$ estimation into the `sse.sir` function and create a new `sse.sir` function called `sse.sir.g`\n\n\n::: {.cell}\n\n```{.r .cell-code}\n#' Basic SIR model\n#'\n#' A basic SIR model with no demographic structure to be used in deSolve\n#'\n#' @param time deSolve passes the time parameter to the function.\n#' @param state A vector of states.\n#' @param params The beta parameter\n#' @param ... Other arguments passed by deSolve.\n#'\n#' @return A deSolve matrix of states at each time step.\n#' @examples\n#' sir_params <- 0.0005\n#' sir_init_states <- c(S = 5000, I = 1, R = 0)\n#' sim_times <- seq(0, 16 / 365, by = 0.1 / 365)\n#'\n#' sir_sol <- deSolve::ode(\n#' y = sir_init_states,\n#' times = sim_times,\n#' func = closed_sir_model,\n#' parms = sir_params\n#' )\n\n# Create a new SIR model to include gamma\nclosed_sir_model_g <- function(time, state, params, ...) {\n # Unpack states\n S <- state[\"S\"]\n I <- state[\"I\"]\n\n # Unpack parameters\n beta <- params[[\"beta\"]]\n dur_inf <- 14 / 365\n gamma <- params[[\"gamma\"]]\n\n new_inf <- beta * S * I\n\n # Calculate the ODEs\n dSdt <- -new_inf\n dIdt <- new_inf - (gamma * I)\n\n # Return the ODEs\n return(list(c(dSdt, dIdt, new_inf)))\n}\n\n\n#' Calculate the Sum of Squared Errors\n#'\n#' A function to take in biweekly incidence data, and SIR parameters, and\n#' calculate the SSE\n#'\n#' @param params A vector of parameter values\n#' @param data A dataframe containing biweekly incidence data in the case column\n#'\n#' @return The SSE of type double\n#' @examples\nsse_sir_g <- function(params, data) {\n # Convert biweekly time series into annual time scale\n # Daily time scale has requires beta values to be too small - doesn't\n # optimize well\n dt <- 0.01\n max_biweek <- max(data$biweek)\n t <- seq(0, max_biweek * 14, dt) / 365\n\n # Extract the number of observed incidence\n obs_inc <- data$cases\n\n # Note the parameters are updated throughout the optimization process by\n # the optim() function\n # Unpack the transmission parameter and exponentiate to fit on ln scale\n beta <- exp(params[[\"beta\"]])\n gamma <- exp(params[[\"gamma\"]])\n\n in_parms <- c(\"beta\" = beta, \"gamma\" = gamma)\n # Unpack the initial states and exponentiate to fit on normal scale\n S_init <- exp(params[[\"S_init\"]])\n I_init <- exp(params[[\"I_init\"]])\n\n # Fit SIR model to the parameters\n sol <- deSolve::ode(\n y = c(S = S_init, I = I_init, new_inf = 0),\n times = t,\n func = closed_sir_model_g,\n parms = in_parms,\n # Use rk4 as fixed time steps, which is important for indexing\n method = \"rk4\"\n )\n\n # Extract the cumulative incidence\n cum_inc <- sol[, \"new_inf\"]\n\n # Find the indices of the cumulative incidence to extract\n biweek_index <- seq(1, max_biweek) * (14 / dt) + 1\n\n # Index cumulative incidence to get the values at the end of the biweeks\n biweek_cum_inc <- cum_inc[biweek_index]\n\n # Calculate the biweekly incidence by using the difference between\n # consecutive biweeks. Need to manually prepend first week's incidence\n # and add in the initial number of infectious individuals, as ODE model\n # only returns the cumulative differences, which is 0 at the start.\n biweek_inc <- c(biweek_cum_inc[1] + I_init, diff(biweek_cum_inc, lag = 1))\n\n # return SSE of predicted vs observed incidence\n return(sum((biweek_inc - obs_inc)^2))\n}\n```\n:::\n\n\nNow we can run the code to optimize, beta, gamma, I0, and S0\n\n\n::: {.cell}\n\n```{.r .cell-code}\nsse_optim_params_g <- c(\n beta = log(0.055),\n gamma = log(365 / 14),\n S_init = log(5000),\n I_init = log(1)\n)\n\n# Create a dataframe of optimized parameters\nniamey_optims_g <- niamey_df %>%\n # Create a nested dataframe i.e. one row for each site, and the data column\n # now is a list column that contains a separate dataframe of times and\n # cases for each site\n nest(data = -site) %>%\n mutate(\n # Map the optim() function call to each of the separate dataframes\n # stored in the nested data column we just created\n fit = map(data, ~ optim(sse_optim_params_g, sse_sir_g, data = .x)),\n # Map the exp() function to each of the model fits just created, and\n # output to a dataframe instead of a list (like in map()), for easier\n # use in the plottinge predictions later\n map_dfr(fit, ~ exp(.x$par))\n )\n```\n:::\n\n\nNow we can get our optimized parameters\n\n\n::: {.cell}\n\n```{.r .cell-code}\nniamey_optims_g %>%\n select(-c(data, fit)) %>%\n mutate(site = str_replace_all(site, \"_\", \" \")) %>%\n gt() %>%\n fmt_number(columns = -site, decimals = 3) %>%\n fmt_scientific(columns = beta, decimals = 3) %>%\n # Relabel the column headers\n cols_label(\n site = md(\"**Site**\"),\n beta = md(\"**Beta**\"),\n gamma = md(\"**Gamma**\"),\n S_init = md(\"**Initial S**\"),\n I_init = md(\"**Initial I**\")\n ) %>%\n # Apply style to the table with gray alternating rows\n opt_stylize(style = 1, color = 'gray') %>%\n # Increate space between columns\n opt_horizontal_padding(scale = 3) %>%\n cols_align(\"center\")\n```\n\n::: {.cell-output-display}\n\n```{=html}\n| Site | \nBeta | \nInitial S | \nInitial I | \n
|---|---|---|---|
| Site 1 | \n5.370 × 10−3 | \n8,566.648 | \n1.349 |
| Site 2 | \n8.317 × 10−3 | \n5,961.388 | \n0.201 |
| Site 3 | \n7.134 × 10−2 | \n792.990 | \n0.001 |
\n\n\n \n
\n
\n```\n\n:::\n:::\n\n\nGetting new predictions from the optimized parameters and a new model run for each site.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nniamey_predictions_g <- niamey_optims_g %>%\n mutate(\n # For each of the different site's nested dataframes, fit the SIR model\n # with the optimal parameters to get best fit predictions\n predictions = pmap(\n .l = list(\n S_init = S_init,\n I_init = I_init,\n beta = beta,\n gamma = gamma,\n time_data = data\n ),\n .f = function(S_init, I_init, beta, gamma, time_data) {\n site_times <- time_data$biweek * 14 / 365\n\n in_parms <- c(\"beta\" = beta, \"gamma\" = gamma)\n # Return a dataframe of model solutions\n as_tibble(ode(\n y = c(S = S_init, I = I_init, new_inf = 0),\n times = site_times,\n func = closed_sir_model_g,\n parms = in_parms,\n hmax = 1 / 120\n )) %>%\n # Make sure all values are numeric for plotting purposes\n mutate(across(everything(), as.numeric)) %>%\n mutate(\n incidence = ifelse(\n row_number() == 1,\n new_inf[1],\n diff(new_inf, lag = 1)\n )\n )\n }\n )\n ) %>%\n unnest(c(data, predictions))\n```\n:::\n\n\nFinally we can plot:\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Create a dataframe to store the positions of the text labels\nniamey_preds_labels <- tibble(\n site = c(\"Site_1\", \"Site_2\"),\n x_label = c(6.5, 6.5),\n x_arrow_just = c(-0.5, -0.5),\n x_arrow_end = c(7, 7.75),\n y_label = c(900, 600),\n y_arrow_just = c(-80, -70),\n y_arrow_end = c(350, 290),\n commentary = c(\"**Predicted\", \"**Observed\"),\n color = c(\"grey20\", niamey_site_colors[\"Site_2\"])\n)\n\nggplot(niamey_predictions_g, aes(x = biweek, group = site)) +\n # Plot the actual data in color\n geom_line(aes(y = cases, color = site)) +\n geom_point(aes(y = cases, color = site), size = 4, alpha = 0.8) +\n # Plot the best-fit model predictions in black\n geom_line(aes(y = incidence), color = \"black\") +\n scale_color_manual(\n values = niamey_site_colors,\n aesthetics = c(\"color\", \"fill\")\n ) +\n # Place each site on it's own subplot and change labels\n facet_wrap(\n ~site,\n ncol = 1,\n scales = \"free_y\",\n labeller = as_labeller(niamey_site_labels)\n ) +\n labs(x = \"Biweek\", y = \"Number of Case\") +\n theme(legend.position = \"none\") +\n ggtext::geom_textbox(\n data = niamey_preds_labels,\n aes(\n label = paste0(\n \"\",\n commentary,\n \" Cases**\",\n \"\"\n ),\n x = x_label,\n y = y_label\n ),\n size = 4,\n fill = NA,\n box.colour = NA\n ) +\n geom_curve(\n data = niamey_preds_labels,\n aes(\n x = x_label + x_arrow_just,\n xend = x_arrow_end,\n y = y_label + y_arrow_just,\n yend = y_arrow_end\n ),\n linewidth = 0.75,\n arrow = arrow(length = unit(0.2, \"cm\")),\n curvature = list(0.25),\n color = \"grey20\"\n )\n```\n\n::: {.cell-output-display}\n{width=672}\n:::\n:::\n\n\n\n::: {.cell}\n\n:::\n\n\n\n::: {.cell}\n\n:::\n\n\n\n::: {.cell}\n\n:::\n\n\n### Exercise 6\n\nWhat happens if one or both of the other unknowns ($X_0$ and $Y_0$) is fixed instead of $\\gamma$?\n\n#### Solutions coming soon!\n\nFirst we modify the `sse_sir_g` function again to fix $X_0$ and $Y_0$ (or both)\n\n\n::: {.cell}\n\n:::\n\n\nNow we can run the optim algorithm and find the best parameters for the SIR model for each of the scenario.\n\nFirst, for S0 fixed:\n\n\n::: {.cell}\n\n:::\n\n\nSecond, for I0 fixed:\n\n\n::: {.cell}\n\n:::\n\n\nThird, for fixing both S0 and I0:\n\n\n::: {.cell}\n\n:::\n\n\n\n::: {.cell}\n\n:::\n\n\nNow let's pass these models through `optim` and fit the models:\n\n\n::: {.cell}\n\n:::\n\n",
"supporting": [],
"filters": [
"rmarkdown/pagebreak.lua"
diff --git a/_freeze/r-session-03/figure-html/unnamed-chunk-30-1.png b/_freeze/r-session-03/figure-html/unnamed-chunk-30-1.png
index dbcc19f..3848a80 100644
Binary files a/_freeze/r-session-03/figure-html/unnamed-chunk-30-1.png and b/_freeze/r-session-03/figure-html/unnamed-chunk-30-1.png differ
diff --git a/_freeze/r-session-03/figure-html/unnamed-chunk-31-1.png b/_freeze/r-session-03/figure-html/unnamed-chunk-31-1.png
index 2f71477..0f0976b 100644
Binary files a/_freeze/r-session-03/figure-html/unnamed-chunk-31-1.png and b/_freeze/r-session-03/figure-html/unnamed-chunk-31-1.png differ
diff --git a/_freeze/r-session-03/figure-html/unnamed-chunk-34-1.png b/_freeze/r-session-03/figure-html/unnamed-chunk-34-1.png
index df4ff47..2e230d4 100644
Binary files a/_freeze/r-session-03/figure-html/unnamed-chunk-34-1.png and b/_freeze/r-session-03/figure-html/unnamed-chunk-34-1.png differ
diff --git a/_freeze/r-session-03/figure-html/unnamed-chunk-35-1.png b/_freeze/r-session-03/figure-html/unnamed-chunk-35-1.png
index 899cbab..93133ae 100644
Binary files a/_freeze/r-session-03/figure-html/unnamed-chunk-35-1.png and b/_freeze/r-session-03/figure-html/unnamed-chunk-35-1.png differ
diff --git a/_freeze/r-session-03/figure-html/unnamed-chunk-38-1.png b/_freeze/r-session-03/figure-html/unnamed-chunk-38-1.png
index 3c276b2..b98e979 100644
Binary files a/_freeze/r-session-03/figure-html/unnamed-chunk-38-1.png and b/_freeze/r-session-03/figure-html/unnamed-chunk-38-1.png differ
diff --git a/_freeze/r-session-03/figure-html/unnamed-chunk-4-1.png b/_freeze/r-session-03/figure-html/unnamed-chunk-4-1.png
index 7f19d74..d6856ea 100644
Binary files a/_freeze/r-session-03/figure-html/unnamed-chunk-4-1.png and b/_freeze/r-session-03/figure-html/unnamed-chunk-4-1.png differ
diff --git a/_freeze/r-session-03/figure-html/unnamed-chunk-43-1.png b/_freeze/r-session-03/figure-html/unnamed-chunk-43-1.png
index cec9d20..45a0d9c 100644
Binary files a/_freeze/r-session-03/figure-html/unnamed-chunk-43-1.png and b/_freeze/r-session-03/figure-html/unnamed-chunk-43-1.png differ
diff --git a/_freeze/r-session-03/figure-html/unnamed-chunk-8-1.png b/_freeze/r-session-03/figure-html/unnamed-chunk-8-1.png
index 9c7d9e0..31cc526 100644
Binary files a/_freeze/r-session-03/figure-html/unnamed-chunk-8-1.png and b/_freeze/r-session-03/figure-html/unnamed-chunk-8-1.png differ
diff --git a/r-session-03.qmd b/r-session-03.qmd
index 8949c4b..91bda71 100644
--- a/r-session-03.qmd
+++ b/r-session-03.qmd
@@ -59,6 +59,12 @@ ggplot(flu, aes(x = day, y = flu)) +
labs(x = "Day", y = "Active Influenza Cases")
```
+The figure above, and the corresponding tabulated outbreak data, report the number of **active influenza cases** on each day.
+These are point-in-time counts, or prevalence: the number of students who were sick at that moment.
+By contrast, the final-size calculation below uses cumulative **incidence**: the total number of students who became cases at any point during the outbreak.
+For this boarding school outbreak, the total school size was $N = 764$, and the total number of students who became cases over the full outbreak was $Z = 512$.
+Because $512$ is a cumulative incidence count, it is not obtained by reading a single value from the prevalence figure or table.
+
## Estimating $R_0$ From The Final Outbreak Size
Our first approach is to estimate $R_0$ from the final outbreak size.
diff --git a/r-session-03_files/figure-html/unnamed-chunk-30-1.png b/r-session-03_files/figure-html/unnamed-chunk-30-1.png
index dbcc19f..3848a80 100644
Binary files a/r-session-03_files/figure-html/unnamed-chunk-30-1.png and b/r-session-03_files/figure-html/unnamed-chunk-30-1.png differ
diff --git a/r-session-03_files/figure-html/unnamed-chunk-31-1.png b/r-session-03_files/figure-html/unnamed-chunk-31-1.png
index 2f71477..0f0976b 100644
Binary files a/r-session-03_files/figure-html/unnamed-chunk-31-1.png and b/r-session-03_files/figure-html/unnamed-chunk-31-1.png differ
diff --git a/r-session-03_files/figure-html/unnamed-chunk-34-1.png b/r-session-03_files/figure-html/unnamed-chunk-34-1.png
index df4ff47..2e230d4 100644
Binary files a/r-session-03_files/figure-html/unnamed-chunk-34-1.png and b/r-session-03_files/figure-html/unnamed-chunk-34-1.png differ
diff --git a/r-session-03_files/figure-html/unnamed-chunk-35-1.png b/r-session-03_files/figure-html/unnamed-chunk-35-1.png
index 899cbab..93133ae 100644
Binary files a/r-session-03_files/figure-html/unnamed-chunk-35-1.png and b/r-session-03_files/figure-html/unnamed-chunk-35-1.png differ
diff --git a/r-session-03_files/figure-html/unnamed-chunk-38-1.png b/r-session-03_files/figure-html/unnamed-chunk-38-1.png
index 3c276b2..b98e979 100644
Binary files a/r-session-03_files/figure-html/unnamed-chunk-38-1.png and b/r-session-03_files/figure-html/unnamed-chunk-38-1.png differ
diff --git a/r-session-03_files/figure-html/unnamed-chunk-4-1.png b/r-session-03_files/figure-html/unnamed-chunk-4-1.png
index 7f19d74..d6856ea 100644
Binary files a/r-session-03_files/figure-html/unnamed-chunk-4-1.png and b/r-session-03_files/figure-html/unnamed-chunk-4-1.png differ
diff --git a/r-session-03_files/figure-html/unnamed-chunk-43-1.png b/r-session-03_files/figure-html/unnamed-chunk-43-1.png
index cec9d20..45a0d9c 100644
Binary files a/r-session-03_files/figure-html/unnamed-chunk-43-1.png and b/r-session-03_files/figure-html/unnamed-chunk-43-1.png differ
diff --git a/r-session-03_files/figure-html/unnamed-chunk-8-1.png b/r-session-03_files/figure-html/unnamed-chunk-8-1.png
index 9c7d9e0..31cc526 100644
Binary files a/r-session-03_files/figure-html/unnamed-chunk-8-1.png and b/r-session-03_files/figure-html/unnamed-chunk-8-1.png differ
diff --git a/wrangler.toml b/wrangler.toml
new file mode 100644
index 0000000..efabbcd
--- /dev/null
+++ b/wrangler.toml
@@ -0,0 +1,5 @@
+name = "sismid-module-01-2026"
+compatibility_date = "2026-06-22"
+
+[assets]
+directory = "./_book"
| Site | \nBeta | \nGamma | \nInitial S | \nInitial I | \n
|---|---|---|---|---|
| Site 1 | \n5.854 × 10−3 | \n81.150 | \n17,113.450 | \n0.786 |
| Site 2 | \n9.066 × 10−3 | \n77.613 | \n11,030.375 | \n0.142 |
| Site 3 | \n8.517 × 10−2 | \n351.697 | \n4,462.286 | \n0.000 |