6 Applications of SIR Models

This chapter brings the full arc of the module together. We have equations, a solver, and a method for estimating parameters. Now we use all three on real data, working through a complete example from raw observations to a fitted model to biological interpretation.

6.1 Epidemic SIR Model: A Complete Worked Example

Learning Objectives

By the end of Section 6.1, you should be able to:

Synthesize the information from Chapters 1 through 6 to construct, solve, and analyze an SIR model using real-world data.

Let’s apply all the information and skills we just learned about to a complete worked example of an epidemic SIR model using real-world data from the 2017-2018 US influenza season.

Step 1: Introducing and Contextualizing Our Data

We will use a modified version of the data on the CDC’s FluView tool, which tracks the weekly number of positive influenza tests reported to the CDC during the 2017–18 flu season in the US (CDC 2026a). The original data reports the number of tests sent into public labs (which report to the CDC) that came back positive for each strain of influenza. We aren’t really interested in what specific strains of influenza are present, so we combined this information to just tell us how many new positive cases of influenza the labs reported each week. We also modified the time reporting of the data set. The original data set reported tests weekly starting from the 40th week of 2017 to the 39th week of 2018 (inclusive). To make the calculations and code easier, we modified the data to only report weeks since the 39th week of 2017; for simplicity, instead of saying “x weeks since the 39th week of 2017,” we will just refer to it as week x. You can view the modified data set below.

You might have noticed that this is actually incidence data, but so far, we have only gone over how to use prevalence data. We are actually going to use this data as an approximation to the true influenza prevalence. We can do this because prevalence is approximately equal to the incidence times the duration of infection (TAMU 2025). According to the CDC, people are infectious with the flu for up to one day before symptoms appear and 5-7 days after (CDC 2026b). Using this information, we can estimate the average duration of infection for influenza to be about 7 days (1 week). Our incidence data is reported weekly, so using our approximation for prevalence, we can say that the prevalence is approximately equal to our incidence data.

We are going to make another approximation about the size of the population which this data comes from ($N$). It might seem obvious to use the population size of the United States for $N$, but that doesn’t account for the fact that this data only represents a very small portion of the total influenza cases seen in the 2017-18 flu season. Very few flu cases are sent to labs for testing, and even fewer of those cases are sent to the public labs which report their data to the CDC. We will assume a population size of $N=150000$ for this example simply because it produces a nice model fit!

Important Note

In a true research context, these approximations for the prevalence and population size would not hold and more exact measures would have to be used. We are just using this data as a learning example to build experience creating SIR models and apply the knowledge gained from this module, so these assumptions are okay to make.

Step 2: Visualizing Our Data and Building Model Equations

Before we construct our SIR model equations, let’s get an idea of what our data looks like plotted.

Our data follows a very nice epidemic-like curve: the disease is growing slowly from weeks 0 through 7 then grows rapidly, reaching a peak around weeks 14 to 16 before decreasing to pre-outbreak levels at a similar rate. Now that we have confirmed our data follows a typical epidemic trend, let’s build our model equations. Since we are using count instead of proportions, our model equations will look like so:

\[ \begin{aligned} \frac{dS}{dt} &= -\beta \left(\frac{I}{N}\right) S \\[10pt] \frac{dI}{dt} &= \beta \left(\frac{I}{N}\right) S - \gamma I \\[10pt] \frac{dR}{dt} &= \gamma I \end{aligned} \]

$\beta$: Transmission rate
$\gamma$: Recovery rate
$N$: Population size

Read through Section 2.1 and Section 2.3 if you need a refresher on how we derive these model equations.

Step 3: Estimating Our Model Parameters

Our model only has two parameters: the recovery rate ($\gamma$) and the transmission rate ($\beta$). We should estimate the recovery rate first because we will need it later when we estimate the transmission rate using MLE.

Recalling from Section 2.2, we can estimate the recovery rate as $\gamma=\frac{1}{D}$ where $D$ is the average duration of infection (i.e. how long individuals are contagious with the disease). Like we mentioned in Step 1, influenza is contagious up to 1 day before symptoms appear and up to 5-7 days after. This gives us a total range of 6-8 days that a person with the flu is infectious. For the sake of simplicity, we will use the median value of this range and suppose that $D=7 \text{ days}$. This gives us an estimated recovery rate of $\gamma=\frac{1}{7} \text{ day}^{-1}$. Our data is measured in units of weeks, not days, so using a simple conversion our final estimated recovery rate is $\gamma=1 \text{ week}^{-1}$.

Note

Although people are technically infectious with the flu for 6-8 days, the contagiousness of the flu is at its highest during the first 3-4 days after symptoms start and dramatically reduces thereafter (Clinic 2022). More complex SIR models intended for research purposes would take this into account, but our simple estimation for $\gamma$ is good enough for this application.

Now that we have an estimate for our recovery rate, let’s shift our focus to estimate our transmission rate $\beta$ using MLE. The code we will use to do this is nearly identical to the code in Section 4.3, so we won’t go over it in detail here. Some variable names have been changed to reflect the new data set, but besides that, comments explain wherever the code below differs significantly.

If you still are confused about the process of coding an MLE algorithm, go back and read Section 4.2 and Section 4.3 before continuing!

Step 1: Extract data and code a model function.

Step 2: Code a function to return model predictions.

Step 3: Code a function to return negative log-likelihoods.

Step 4: Call mle2 and output estimate for $\beta$.

Keep this estimate for $\beta$ in mind because we will use it in the next section to solve and visualize our model!

Step 4: Visualizing Our Model

Now that we have estimates for our model parameters, it’s time to write our full SIR equations and solve them:

\[ \begin{aligned} \frac{dS}{dt} &= -(1.31) \left(\frac{I}{150,000}\right) S \\[10pt] \frac{dI}{dt} &= (1.31) \left(\frac{I}{150,000}\right) S - I \\[10pt] \frac{dR}{dt} &= I \end{aligned} \]

How did we get these equations?

We got these equations by plugging in our estimates we found before ($\beta=1.31$, $N=150000$, and $\gamma=1$) to the SIR equations from Step 1.

Now that we have our full equations, let’s go ahead and solve them using rk4sys just like we did in Section 3.3.

Step 1: Use rk4sys to produce a data frame of solutions.

Step 2: Plotting the curves.

Plotting all three SIR curves on a single plot will not be very helpful to us since the number of susceptibles is in the hundreds of thousands while the number of infected is only in the thousands. A base R plot will be too zoomed out to show any useful details. To produce a useful plot with all three curves, a special plotting package must be used. Step 5 has a plot of all three curves that can be panned and zoomed into. In this section, we will just plot each curve separately using base R, starting with the susceptible curve:

Now let’s plot the infectious curve. We will also overlay our data to see how good of a fit our model is:

Finally, let’s plot our recovered curve:

Looking at our second plot, it appears that our model fits quite well! It’s not exactly the same as the original data, but that’s not the point: the point is to get an idea of the general trends present in the data, not to capture every small detail of it.

Step 5: Analyzing and Evaluating Our Model

Like we mentioned in Step 1, our model uses data modified from its original form, assumes an approximate equivalence between prevalence and incidence (an assumption that does not always hold in real life), and a placeholder value for the population size. We did these things for the sake of producing a nice, clean example that applies the knowledge from this module, not to produce a model with real-life applications and implications. Any conclusions we draw from this model are purely hypothetical and just for practice; it doesn’t necessarily reflect how influenza acts in real life!

Models as approximations

Every model we build is a simplification. The question is not whether the simplification is perfect, but whether it is useful. A model that captures the general shape of an epidemic, even with rough parameter estimates, can still tell us something real about transmission dynamics and vaccination thresholds.

Now that we have that disclaimer out of the way, let’s analyze our model!

The basic reproduction number of our model is $R_0=\frac{\beta}{\gamma}=1.31$. This means that (according to our model) people with the flu, on average, spread it to a little more than one person throughout the entire time they are sick. In line with the threshold phenomenon, our epidemic was able to spread because the initial proportion of susceptible people $\tilde{S}(0)=\frac{150000-161}{150000}\approx99.9\%$ is much greater than the reciprocal of our basic reproduction number $\left(\frac{1}{R_0}\approx76.3\%\right)$.

The plot below shows all three SIR curves. Hovering over any of the three curves displays the effective reproduction number at that time. Like we learned in Section 5.2, the epidemic grows when $R_e>1$, meets its peak at $R_e=1$, and shrinks when $R_e<1$.

Warning: package 'deSolve' was built under R version 4.5.2

Warning: package 'dplyr' was built under R version 4.5.2

Warning: package 'tidyr' was built under R version 4.5.2

Warning: package 'plotly' was built under R version 4.5.2

Warning: package 'ggplot2' was built under R version 4.5.2

Tip

If the curves are hard to see, try isolating a curve by double clicking on its entry in the legend then click on the “Autoscale” tool in the toolbar (tool to the left of the home icon). You can also use the “Pan” tool (icon that looks like a four-way arrow) to move the view window around.

Recalling the information from Section 5.3, we should have aimed for at least $1-\frac{1}{R_0}\approx24\%$ of the population to be successfully vaccinated against influenza to prevent this epidemic from occurring. The influenza vaccine for the 2017-18 season had an estimated $38\%$ efficacy rate (CDC 2026c), so the actual percentage of the population we should aim to vaccinate is closer to $\frac{0.24}{0.38}=63\%$.

It might seem counterintuitive to estimate these values using information after the outbreak already occurred, but it’s important to note that influenza outbreaks like the one we are analyzing occur every year during a time known as “flu season” (CDC 2025). These influenza outbreaks act similarly every year. We highly encourage you to visit the CDC’s Interactive FluView tool and look at the data they have on positive influenza cases and hospitalization rates from year to year: a pattern quickly emerges showing that flu numbers peak around the end of December and early January then drop drastically from February onward.

This means that the ideal percentage of people to vaccinate for previous years’ influenza outbreaks does a decent job at informing public health measures for upcoming outbreaks. In fact, studies have estimated the basic reproduction number for seasonal influenza outbreaks to be around $1.28$, which is not too far off from the basic reproduction number of our model ($R_0=1.31$) (Biggerstaff et al. 2014). Even though our model isn’t exactly reflective of how influenza acts in reality, it’s a simple yet decent approximation for how more complex SIR models of influenza would behave.

Chapter References

Biggerstaff, Matthew, Simon Cauchemez, Carrie Reed, Manoj Gambhir, and Lyn Finelli. 2014. “Estimates of the Reproduction Number for Seasonal, Pandemic, and Zoonotic Influenza: A Systematic Review of the Literature.” BMC Infectious Diseases 14 (1). https://doi.org/https://doi.org/10.1186/1471-2334-14-480.

CDC. 2025. “CDC Seasonal Flu Vaccine Effectiveness Studies.” Flu Vaccines Work. https://www.cdc.gov/flu-vaccines-work/php/effectiveness-studies/index.html.

———. 2026a. “National, Regional, and State Level Outpatient Illness and Viral Surveillance.” FluView Interactive. CDC. https://gis.cdc.gov/grasp/fluview/fluportaldashboard.html.

———. 2026b. “Influenza (Flu).” Centers for Disease Control and Prevention. https://www.cdc.gov/flu/about/index.html.

———. 2026c. “CDC Seasonal Flu Vaccine Effectiveness Studies.” Flu Vaccines Work. https://www.cdc.gov/flu-vaccines-work/php/effectiveness-studies/index.html.

Clinic, Cleveland. 2022. “Influenza (Flu): Treatment, Prevention, Symptoms Vs Cold.” Cleveland Clinic. https://my.clevelandclinic.org/health/diseases/4335-influenza-flu.

TAMU. 2025. “Epidemiology Incidence Vs. Prevalence: Exploring Two of the Most Impactful Concepts in Public Health.” Texas A&M University School of Public Health. Tamu.edu. https://public-health.tamu.edu/degrees/mph/blog/epidemiology-incidence-vs-prevalence-explained.html.

--- title: "Applications of SIR Models" engine: knitr webr: packages: - pracma - bbmle --- This chapter brings the full arc of the module together. We have equations, a solver, and a method for estimating parameters. Now we use all three on real data, working through a complete example from raw observations to a fitted model to biological interpretation. ## Epidemic SIR Model: A Complete Worked Example {#sec-worked-example} ::: {.callout-tip title="Learning Objectives"} By the end of @sec-worked-example, you should be able to: - Synthesize the information from Chapters 1 through 6 to construct, solve, and analyze an SIR model using real-world data. ::: Let's apply all the information and skills we just learned about to a complete worked example of an epidemic SIR model using real-world data from the 2017-2018 US influenza season. ### Step 1: Introducing and Contextualizing Our Data {#step-1 .unnumbered} We will use a modified version of the data on the CDC's FluView tool, which tracks the weekly number of positive influenza tests reported to the CDC during the 2017–18 flu season in the US [@cdc_flu_view_data]. The original data reports the number of tests sent into public labs (which report to the CDC) that came back positive for each strain of influenza. We aren't really interested in what specific strains of influenza are present, so we combined this information to just tell us how many new positive cases of influenza the labs reported each week. We also modified the time reporting of the data set. The original data set reported tests weekly starting from the 40th week of 2017 to the 39th week of 2018 (inclusive). To make the calculations and code easier, we modified the data to only report weeks since the 39th week of 2017; for simplicity, instead of saying "x weeks since the 39th week of 2017," we will just refer to it as week x. You can view the modified data set below. ```{r, echo=FALSE} library(readxl) library(DT) data <- read_excel("../data/disease_flu_2017_MODIFIED.xlsx") DT::datatable(data) ``` You might have noticed that this is actually incidence data, but so far, we have only gone over how to use prevalence data. We are actually going to use this data as an approximation to the true influenza prevalence. We can do this because prevalence is approximately equal to the incidence times the duration of infection [@tamu_inc]. According to the CDC, people are infectious with the flu for up to one day before symptoms appear and 5-7 days after [@flu_cdc]. Using this information, we can estimate the average duration of infection for influenza to be about 7 days (1 week). Our incidence data is reported weekly, so using our approximation for prevalence, we can say that the prevalence is approximately equal to our incidence data. We are going to make another approximation about the size of the population which this data comes from ($N$). It might seem obvious to use the population size of the United States for $N$, but that doesn't account for the fact that this data only represents a very small portion of the total influenza cases seen in the 2017-18 flu season. Very few flu cases are sent to labs for testing, and even fewer of those cases are sent to the public labs which report their data to the CDC. We will assume a population size of $N=150000$ for this example simply because it produces a nice model fit! :::{.callout-important title="Important Note"} In a true research context, these approximations for the prevalence and population size would not hold and more exact measures would have to be used. We are just using this data as a learning example to build experience creating SIR models and apply the knowledge gained from this module, so these assumptions are okay to make. ::: ### Step 2: Visualizing Our Data and Building Model Equations {#step-2 .unnumbered} Before we construct our SIR model equations, let's get an idea of what our data looks like plotted. ```{webr-r} # Importing our data set from GitHub data <- read.csv("https://raw.githubusercontent.com/javirudolph/ecoSIR-modules/refs/heads/main/modules/data/disease_flu_2017_MODIFIED.csv") # Plotting the data plot(data$week, data$count, xlab = "week", ylab = "count") ``` Our data follows a very nice epidemic-like curve: the disease is growing slowly from weeks 0 through 7 then grows rapidly, reaching a peak around weeks 14 to 16 before decreasing to pre-outbreak levels at a similar rate. Now that we have confirmed our data follows a typical epidemic trend, let's build our model equations. Since we are using count instead of proportions, our model equations will look like so: <br> $$ \begin{aligned} \frac{dS}{dt} &= -\beta \left(\frac{I}{N}\right) S \\[10pt] \frac{dI}{dt} &= \beta \left(\frac{I}{N}\right) S - \gamma I \\[10pt] \frac{dR}{dt} &= \gamma I \end{aligned} $$ - $\beta$: Transmission rate - $\gamma$: Recovery rate - $N$: Population size Read through [@sec-the-classic-sir-epidemic-model] and [@sec-formulating-foi] if you need a refresher on how we derive these model equations. ### Step 3: Estimating Our Model Parameters {#step-3 .unnumbered} Our model only has two parameters: the recovery rate ($\gamma$) and the transmission rate ($\beta$). We should estimate the recovery rate first because we will need it later when we estimate the transmission rate using MLE. Recalling from [Section @sec-formulating-the-recovery-rate], we can estimate the recovery rate as $\gamma=\frac{1}{D}$ where $D$ is the average duration of infection (i.e. how long individuals are contagious with the disease). Like we mentioned in [Step 1](#step-1), influenza is contagious up to 1 day before symptoms appear and up to 5-7 days after. This gives us a total range of 6-8 days that a person with the flu is infectious. For the sake of simplicity, we will use the median value of this range and suppose that $D=7 \text{ days}$. This gives us an estimated recovery rate of $\gamma=\frac{1}{7} \text{ day}^{-1}$. Our data is measured in units of weeks, not days, so using a simple conversion our final estimated recovery rate is $\gamma=1 \text{ week}^{-1}$. :::{.callout-note} Although people are technically infectious with the flu for 6-8 days, the contagiousness of the flu is at its highest during the first 3-4 days after symptoms start and dramatically reduces thereafter [@cleveland_clinic_flu]. More complex SIR models intended for research purposes would take this into account, but our simple estimation for $\gamma$ is good enough for this application. ::: Now that we have an estimate for our recovery rate, let's shift our focus to estimate our transmission rate $\beta$ using MLE. The code we will use to do this is nearly identical to the code in [Section @sec-bbmle], so we won't go over it in detail here. Some variable names have been changed to reflect the new data set, but besides that, comments explain wherever the code below differs significantly. If you still are confused about the process of coding an MLE algorithm, go back and read [@sec-coding-scratch] and [@sec-bbmle] before continuing! **Step 1: Extract data and code a model function.** ```{webr-r} library(pracma) data <- read.csv("https://raw.githubusercontent.com/javirudolph/ecoSIR-modules/refs/heads/main/modules/data/disease_flu_2017_MODIFIED.csv") ModelFunc <- function(time, state, parameters) { S <- state[1] I <- state[2] R <- state[3] beta <- parameters["beta"] gamma <- parameters["gamma"] N <- parameters["N"] dS <- -(beta * I * S / N) dI <- (beta * I * S / N) - (gamma * I) dR <- (gamma * I) return(c(dS, dI, dR)) } print("Data has been imported and model function has been defined.") ``` --- **Step 2: Code a function to return model predictions.** ```{webr-r} library(pracma) pred <- function(parameters, t) { # CHANGE: Extracting the initial infected count from our data init_inf <- data[1, 2] # CHANGE: Different initial values initial_vals <- c(S = parameters["N"] - init_inf, I = init_inf, R = 0) t_start <- 0 t_end <- max(t) step_size <- t_end - t_start output <- rk4sys( f = function(t, y) ModelFunc(t, y, parameters), a = t_start, b = t_end, y0 = initial_vals, n = step_size ) output_df <- as.data.frame(output) prevalences <- output_df$y.I[t] prevalences[is.nan(prevalences) | prevalences < 0] <- 0 return(prevalences) } print("Model prediction function has been defined.") ``` --- **Step 3: Code a function to return negative log-likelihoods.** ```{webr-r} library(bbmle) binomial_ll <- function(b) { times <- data$week # CHANGE: Different parameters parameters <- c(beta = b, gamma = 1, N = 150000) predicted <- pred(parameters, times) -sum(dbinom( x = data$count, size = parameters["N"], prob = predicted / parameters["N"], log = TRUE )) } print("Negative log-likelihood function has been defined.") ``` --- **Step 4: Call `mle2` and output estimate for $\beta$.** ```{webr-r} library(bbmle) start <- list(b = 2) result <- mle2( minuslogl = binomial_ll, start = start, method = "L-BFGS-B", lower = c(b = 0.01) ) coef(result) ``` --- Keep this estimate for $\beta$ in mind because we will use it in the next section to solve and visualize our model! ### Step 4: Visualizing Our Model {#step-4 .unnumbered} Now that we have estimates for our model parameters, it's time to write our full SIR equations and solve them: $$ \begin{aligned} \frac{dS}{dt} &= -(1.31) \left(\frac{I}{150,000}\right) S \\[10pt] \frac{dI}{dt} &= (1.31) \left(\frac{I}{150,000}\right) S - I \\[10pt] \frac{dR}{dt} &= I \end{aligned} $$ :::{.callout-note collapse="true"} ## How did we get these equations? We got these equations by plugging in our estimates we found before ($\beta=1.31$, $N=150000$, and $\gamma=1$) to the SIR equations from [Step 1](#step-1). ::: Now that we have our full equations, let's go ahead and solve them using `rk4sys` just like we did in [@sec-pracma]. **Step 1: Use `rk4sys` to produce a data frame of solutions.** ```{webr-r} library(pracma) # Importing our data data <- read.csv("https://raw.githubusercontent.com/javirudolph/ecoSIR-modules/refs/heads/main/modules/data/disease_flu_2017_MODIFIED.csv") # SIR model function ModelFunc <- function(time, state) { S <- state[1] I <- state[2] R <- state[3] # The estimates found in Step 3 beta <- 1.31 gamma <- 1 dS <- -(beta * I * S / 150000) dI <- (beta * I * S / 150000) - (gamma * I) dR <- (gamma * I) return(c(dS, dI, dR)) } # Other values to be passed into rk4sys start <- 0 end <- 51 initial_vals <- c(S = 150000 - data[1, 2], I = data[1, 2], R = 0) step_size <- 10000 # Calling rk4sys output <- rk4sys( f = ModelFunc, a = start, b = end, y0 = initial_vals, n = step_size ) # Converting output to data frame for future plotting output <- as.data.frame(output) head(output, 5) ``` --- **Step 2: Plotting the curves.** Plotting all three SIR curves on a single plot will not be very helpful to us since the number of susceptibles is in the hundreds of thousands while the number of infected is only in the thousands. A base R plot will be too zoomed out to show any useful details. To produce a useful plot with all three curves, a special plotting package must be used. [Step 5](#step-5) has a plot of all three curves that can be panned and zoomed into. In this section, we will just plot each curve separately using base R, starting with the susceptible curve: ```{webr-r} times <- output$x S <- output$y.S beta <- 1.31 gamma <- 1 plot(times, S, type = "l", col = "steelblue", lwd = 2, xlab = "Time", ylab = "Prevalence (counts)", main = paste0( "SIR Susceptible Curve (Influenza, 2017-18) | β=", beta, " γ=", gamma, " R₀=", round(beta / gamma, 2) ) ) ``` Now let's plot the infectious curve. We will also overlay our data to see how good of a fit our model is: ```{webr-r} I <- output$y.I plot(times, I, type = "l", col = "firebrick", lwd = 2, xlab = "Time", ylab = "Prevalence (counts)", main = paste0( "SIR Infectious Curve (Influenza, 2017-18) | β=", beta, " γ=", gamma, " R₀=", round(beta / gamma, 2) ) ) points(data$week, data$count) ``` Finally, let's plot our recovered curve: ```{webr-r} R <- output$y.R plot(times, R, type = "l", col = "seagreen", lwd = 2, xlab = "Time", ylab = "Prevalence (counts)", main = paste0( "SIR Recovered Curve (Influenza, 2017-18) | β=", beta, " γ=", gamma, " R₀=", round(beta / gamma, 2) ) ) ``` Looking at our second plot, it appears that our model fits quite well! It's not exactly the same as the original data, but that's not the point: the point is to get an idea of the general trends present in the data, not to capture every small detail of it. ### Step 5: Analyzing and Evaluating Our Model {#step-5 .unnumbered} Like we mentioned in [Step 1](#step-1), our model uses data modified from its original form, assumes an approximate equivalence between prevalence and incidence (an assumption that does not always hold in real life), and a placeholder value for the population size. We did these things for the sake of producing a nice, clean example that applies the knowledge from this module, not to produce a model with real-life applications and implications. **Any conclusions we draw from this model are purely hypothetical and just for practice; it doesn't necessarily reflect how influenza acts in real life!** ::: {.callout-note title="Models as approximations"} Every model we build is a simplification. The question is not whether the simplification is perfect, but whether it is useful. A model that captures the general shape of an epidemic, even with rough parameter estimates, can still tell us something real about transmission dynamics and vaccination thresholds. ::: Now that we have that disclaimer out of the way, let's analyze our model! The basic reproduction number of our model is $R_0=\frac{\beta}{\gamma}=1.31$. This means that (according to our model) people with the flu, on average, spread it to a little more than one person throughout the entire time they are sick. In line with the threshold phenomenon, our epidemic was able to spread because the initial proportion of susceptible people $\tilde{S}(0)=\frac{150000-161}{150000}\approx99.9\%$ is much greater than the reciprocal of our basic reproduction number $\left(\frac{1}{R_0}\approx76.3\%\right)$. The plot below shows all three SIR curves. Hovering over any of the three curves displays the effective reproduction number at that time. Like we learned in [@sec-the-threshold-phenomenon], the epidemic grows when $R_e>1$, meets its peak at $R_e=1$, and shrinks when $R_e<1$. ```{r, echo=FALSE, message=FALSE} library(deSolve) library(dplyr) library(tidyr) library(plotly) data <- read.csv("https://raw.githubusercontent.com/javirudolph/ecoSIR-modules/refs/heads/main/modules/data/disease_flu_2017_MODIFIED.csv") solve_sir <- function(S0, I0, R0, b, g, n, len, dt) { sir <- function(time, state, parameters) { with(as.list(c(state, parameters)), { dS <- -beta * S * I / pop_size dI <- beta * S * I / pop_size - gamma * I dR <- gamma * I return(list(c(dS, dI, dR))) }) } init <- c(S=S0, I=I0, R=R0) parameters <- c(beta=b, gamma=g, pop_size=n) times <- seq(0, len, by = dt) out <- as.data.frame(ode(y = init, times = times, func = sir, parms = parameters)) out$Re <- (b/g)*(out$S/n) return(out) } # Change these if needed parameters <- c(b=1.31, g=1, N=150000) df <- solve_sir(S0=150000-data[1,2], I0=data[1,2], R0=0, b=as.numeric(parameters["b"]), g=as.numeric(parameters["g"]), n=as.numeric(parameters["N"]), len=51, dt=0.1) df_long <- pivot_longer(df, cols = c(S, I, R), names_to = "Compartment", values_to = "Value") df_long$Compartment <- factor(df_long$Compartment, levels = c("S", "I", "R"), labels = c("Susceptible", "Infected", "Recovered")) p <- ggplot(df_long, aes(x = time, y = Value, color = Compartment, group = Compartment, text = paste0("Time (weeks): ", round(time, 2), "<br>", Compartment, ": ", round(Value, 6), "<br>Effective Reproduction Number: ", round(Re, 3)))) + geom_line(linewidth = 0.6) + scale_color_manual(values=c("Susceptible"="steelblue", "Infected"="firebrick", "Recovered"="seagreen")) + labs(x = "Time (weeks)", y = "Prevalence (count)", color = "Compartment") + ggtitle(paste0("Epidemic SIR Curves (Influenza, 2017-18) | β=", parameters["b"], " γ=", parameters["g"], " R₀=", round(parameters["b"] / parameters["g"], 2))) ggplotly(p, tooltip = "text") ``` :::{.callout-tip} If the curves are hard to see, try isolating a curve by double clicking on its entry in the legend then click on the "Autoscale" tool in the toolbar (tool to the left of the home icon). You can also use the "Pan" tool (icon that looks like a four-way arrow) to move the view window around. ::: Recalling the information from [@sec-choosing-how-many-people-to-vaccinate], we should have aimed for at least $1-\frac{1}{R_0}\approx24\%$ of the population to be successfully vaccinated against influenza to prevent this epidemic from occurring. The influenza vaccine for the 2017-18 season had an estimated $38\%$ efficacy rate [@cdc_flu_vaccine_eff], so the actual percentage of the population we should aim to vaccinate is closer to $\frac{0.24}{0.38}=63\%$. It might seem counterintuitive to estimate these values using information after the outbreak already occurred, but it's important to note that influenza outbreaks like the one we are analyzing occur every year during a time known as "flu season" [@cdc_flu_season]. These influenza outbreaks act similarly every year. We highly encourage you to visit the [CDC's Interactive FluView tool](https://www.cdc.gov/fluview/overview/fluview-interactive.html) and look at the data they have on positive influenza cases and hospitalization rates from year to year: a pattern quickly emerges showing that flu numbers peak around the end of December and early January then drop drastically from February onward. This means that the ideal percentage of people to vaccinate for previous years' influenza outbreaks does a decent job at informing public health measures for upcoming outbreaks. In fact, studies have estimated the basic reproduction number for seasonal influenza outbreaks to be around $1.28$, which is not too far off from the basic reproduction number of our model ($R_0=1.31$) [@flu_reprod_num]. Even though our model isn't exactly reflective of how influenza acts in reality, it's a simple yet decent approximation for how more complex SIR models of influenza would behave. ## Chapter References {.unnumbered .unlisted}