Previously

• Markov models: the current state of the process fixes the distribution of the future; the past is irrelevant \[ \mathrm{Pr}\left( X(t+1), X(t+2), \ldots | X(t), X(t-1), \ldots X(0) \right) = \mathrm{Pr}\left( X(t+1), X(t+2), \ldots| X(t) \right) \]
• Everything boils down to the transition probabilities, \(\mathrm{Pr}\left( X(t+1) = y|X(t) = x \right)\)
• Inference for Markov processes:
  • Based on the likelihood, which we can write in terms of the transition probabilities
  • Or based on aggregated counts

“SIR” Epidemic Model

• Models the spread of a disease through a population
• Compartments: Every person can be in one of three states:
  • Susceptible (\(S\)): not sick, could get sick
  • Infectious (\(I\)): sick, can make others sick
  • Removed (\(R\)): either recovered or dead; not sick, can’t get sick, can’t make others sick
• Contagion: When an \(S\) meets an \(I\), some probability of the \(S\) also becoming \(I\)
• Removal: \(I\)s spontaneously change into \(R\)s
  • Maybe after some random time
  • Maybe at some steady rate per unit time
• Mixing or mass action: probability of an \(S\) meeting an \(I\) just depends on total numbers of \(S\)s, \(I\)s and \(R\)s

Stochastic form of the SIR model (1)

• \(S(t) =\) number of \(S\)s at time t, similarly for \(I(t)\) and \(R(t)\)
• Total (initial) population is \(n\)
• Small time increment \(h\), during which each \(S\) person encounters \(k h\) people
• The probability of contagion from an encounter between an \(S\) and an \(I\) is \(c\)
• Probability that a random person is infected is \(I(t)/n\)
• So the probability of an \(S\) becoming \(I\) is \(kch I(t)/n\)
  • Strictly it’s \(1-(1-c I(t)/n)^{hk} \approx kch I(t)/n\) for small \(h\)
• Number of new infections in time \(h\) is \(\mathrm{Binom}(S(t), kch I(t)/n)\)
  • If \(h\) is very small but \(S(t)\) is very large, this is approximately \(\mathrm{Poisson}(\frac{kch}{n} S(t) I(t))\)
• Every \(I\) gets removed with probability \(\gamma h\) so the number of removals in time \(h\) is \(\mathrm{Binom}(I(t), \gamma h)\)
  • If \(h\) is very small but \(I(t)\) is very large, this is approximately \(\mathrm{Poisson}(\gamma h I(t))\)

Stochastic form of the SIR model (2)

• Over-all equations: \[\begin{eqnarray} C(t+h) & = & \mathrm{Binom}(S(t), \frac{\beta}{n} h I(t))\\ D(t+h) & = & \mathrm{Binom}(I(t), \gamma h)\\ S(t+h) & = & S(t) - C(t+h)\\ I(t+h) & = & I(t) + C(t+h) - D(t+h)\\ R(t+h) & = & R(t) + D(t+h) \end{eqnarray}\]

• This bundles up multiple parameters we don’t get to measure separately into some “effective” parameters \(\beta\) and \(\gamma\)
  • I left in the factors of \(h\) and \(n\) for reasons which will become clear later

The SIR model is a Markov model

• The state variable \(X(t) = (S(t), I(t), R(t))\)
  • Because \(S(t) + I(t) + R(t) = n\), could eliminate one variable from the state (usually pick \(R(t)\))
• The equations on the previous slide implicitly give the transition probabilities
  • Explicit expressions are possible but tedious
• Markov property holds by assumption

SIR is easy to simulate

sim.SIR <- function(n, beta, gamma, s.initial = n - 1, i.initial = 1, r.initial = 0, 
    T) {
    stopifnot(s.initial + i.initial + r.initial == n)
    states <- matrix(NA, nrow = 3, ncol = T)
    rownames(states) <- c("S", "I", "R")
    states[, 1] <- c(s.initial, i.initial, r.initial)
    for (t in 2:T) {
        contagions <- rbinom(n = 1, size = states["S", t - 1], prob = beta * states["I", 
            t - 1]/n)
        removals <- rbinom(n = 1, size = states["I", t - 1], prob = gamma)
        states["S", t] <- states["S", t - 1] - contagions
        states["I", t] <- states["I", t - 1] + contagions - removals
        states["R", t] <- states["R", t - 1] + removals
    }
    return(states)
}

• The code doesn’t need \(h\) here because it’s (implicitly) taking \(h=1\)

• We’ll start by simulating the model to see what our assumptions imply (as one should always do)

What One Simulation Looks Like

Multiple Simulations with the Same Settings

• Most of these show similar behavior, but not quite the same
  • Why do some runs have later outbreaks than others?
  • What’s up with the horizontal black line at the top?

What If We Make Contagion Harder?

• The one initial \(I\) gets removed before it can start an epidemic

Maybe Not Quite That Much Harder?

• Note the change on the horizontal axis
• Sometimes the epidemic never gets started, sometimes it takes off but doesn’t get everyone before dying out…
  • Notice that when it does break out, it stops when about 80% of the population has been infected, not nearly-100% like in the first set of simulations

Go Back to the Original Ease of Contagion, But Make Removal Faster

• Same long-run number of infected as the last simulations, though we get there faster

Some Suggestions from Those Simulations

• Epidemic outbreak when \(\beta\) is sufficiently large, compared to \(\gamma\)
  • At least, outbreaks with high probability
• Epidemics don’t have to affect everyone before they die out
• It may be the ratio \(\beta/\gamma\) is what matters for the probability of an outbreak and for the ultimate fraction ever infected
  • Rather than, say, the difference \(\beta-\gamma\) or even \(\beta^2-\gamma^2\)

Deterministic limit

• Imagine making the population size really big (\(n \rightarrow \infty\)), and the time step really small (\(h \rightarrow 0\))
• The fluctuations get (relatively) negligible, and we’re left with deterministic equations (see backup)

• When people talk about “the SIR model”, it’s often ambiguous between the stochastic model we started with and these ordinary differential equations (ODEs)
• See backup for manipulations to these equations which show how to re-write them so the only parameter is \(\beta/\gamma\)

Simulating the Deterministic Limit

Go back to the original parameter values:

Is There Going to be an Epidemic?

• Epidemic transition: Typically, the contagion either dies out rapidly or it takes over a large fraction of the population
  • We saw evidence for this in the simulation
  • More mathematically: As \(n\) grows, does the total number who are ever infected grow with \(n\) or not?
  • I.e., is \(R(\infty) = O(n)\) (an epidemic) or \(R(\infty) = o(n)\) (no epidemic)?
• One way to approach this is through the basic reproductive number \(\mathcal{R}_0\): \[ \mathcal{R}_0\equiv \mathbf{E}\left[ \text{number of new infections directly produced by adding one infected into an all-susceptible population} \right] \]
  • Not the same as \(R(0) =\) number removed at time \(0\)!
  • (Sometimes written \(R_0\), \(r_0\), etc., but I’ve used \(\mathcal{R}_0\) to lessen confusion)
• Each of those \(\mathcal{R}_0\) new infections also faces an (almost) entirely susceptible population, so they’ll cause \(\mathcal{R}_0\) new infections on average as well
  • Obviously can’t be true for ever, but is true initially
• After \(d\) rounds of contagion, we expect \(\mathcal{R}_0^d\) infected
  • Again, initially
• The early stage of an epidemic has exponential growth
• There is an epidemic threshold at \(\mathcal{R}_0=1\)
  • \(\mathcal{R}_0< 1\): contagion dies out before getting very far (sub-critical)
  • \(\mathcal{R}_0> 1\): contagion will grow exponentially, epidemic outbreak (super-criticial)
  • \(\mathcal{R}_0= 1\): Will die out, but can live a long time before doing so (critical)

The Epidemic Threshold, Illustrated

(Details: tracing out to a maximum of 6 rounds of growth; assuming a geometric distribution for the number of new infections, with the mean given by \(\mathcal{R}_0\); color-coded by generation)

Relating \(\mathcal{R}_0\) to SIR Parameters (1)

• Start by guestimating
• \(I(t)=1\), \(S(t) \approx n\), so expected number of infections per time step (of length \(h\)) is \(\approx \beta\)
• Expected number of time steps before the initial \(I\) gets removed is \(1/\gamma\)
• So (roughly!) \[ \mathcal{R}_0= \frac{\beta}{\gamma} \]
  • Our initial simulations had \(\beta=10^{-1}\), \(\gamma=10^{-2}\) so \(\mathcal{R}_0= 10\)
  • For the second set of simulation, \(\mathcal{R}_0= 1\)
  • For the third and fourth sets, \(\mathcal{R}_0= 2\)

Relating \(\mathcal{R}_0\) to SIR Parameters (2)

• There are more rigorous ways to come to the same conclusion…
  • For the stochastic model: approximate it as a branching process
  • For the deterministic model: look at the stability of the all-susceptible state
• Upshot: In the SIR model, \[ \mathcal{R}_0= \frac{\beta}{\gamma} \]
  • Holds in both the stochastic and deterministic versions

Qualitative Results for the Deterministic Model

• \(dS/dt \leq 0\) so the number of susceptibles never goes up
• \(dR/dt \geq 0\) so the number removed never goes down
• If we begin with \(I(0) \ll n\), \(R(0)=0\), \(S(0) \approx n\), then at small times \(t\), \[ I(t) \approx I(0) e^{(\beta-\gamma) t} \]
  • See backup for derivation
• That is, early in an outbreak we expect to see exponential growth at the rate \(\beta-\gamma\)
  • Notice: positive growth rate iff \(\mathcal{R}_0> 1\)
  • Notice: at this early stage, \(I(t)\) is doesn’t care about total population size \(n\), so don’t compare outbreaks by their per-capita infection rates
• Once \(S\) gets small enough, \(dI/dt < 0\)
  • Specifically, \(dI/dt < 0\) once \(\frac{S}{n} \leq \frac{\gamma}{\beta} = \frac{1}{\mathcal{R}_0}\)
• So \(I(t)\) will grow exponentially and then eventually decline
  • Implication: stop an outbreak by removing enough susceptibles
• \(I(t)\) is going to be a one-humped curve, usually skewed to the right
• \(R(t) = \gamma \int_{0}^{t}{I(s)ds}\), it’s \(I(t)\) integrated over time, so it’s going to look like the logistic function \(\frac{e^t}{1+e^t}\)
  • Or some other similar curve — not necessarily that functional form
• There are ways of calculating \(R(\infty)\) but they involve solving a transcendental equation numerically (Miller 2012)
  • The upshot is that \(R(\infty)/n\) is an increasing function of \(\mathcal{R}_0\)

Back-of-the-envelope SIR Modeling for Covid-19 in the USA

• Did this back in early March (file last updated 11 March)
• Initial state:
  • \(S(0) = 3\times {10}^8\) (approx. US population)
  • \(I(0) = 3\times {10}^3\) (number of known US cases at the time, \(\times\) a fudge factor of about 30)
  • \(R(0) = 1\) just because I wanted to plot everything on a log scale
• Parameters:
  • \(1/\gamma =\) average length of time someone is infectious, say 12 days
    • An early estimate; we could adjust that now
  • \(\mathcal{R}_0=\beta/\gamma\), and one plausible estimate back then was \(\mathcal{R}_0\approx 2.6\) (Imai et al. 2020)
  • I backed out \(\beta\) from \(\mathcal{R}_0\) and \(\gamma\)

What Did This Predict?

• About 85 days to reach 50% of the population, or early-mid June
  • More exactly, 85.06 days
• About 150 days to reach 90% of the population, or early-mid September
  • More exactly, 148.87 days
• Level off at a bit over 90% of the population
  • More exactly, 0.9048903
• Assuming a death rate of 1% among those who contract it, \(\approx 2.7\times 10^{6}\) deaths

How Well Did This Do?

• Use Kieran Healy’s covdata package to get data on reported cases and deaths for the US

• The 100th confirmed case in the US was reported on 2020-03-03, about when I did the simple model

• Remember I assumed 30 times more true infections than known infections
  • 30 was pretty arbitrary, so we should only expect to get the shape right
  • I.e., parallel curves on the log-scaled plot
  • I also show deaths because those are (arguably) less under-counted than cases; they’d lag infections but should show parallel behavior
• The model is not doing too badly up to about day 20, or late March, by which point everything was locked down
• N.B.: This is the simplest possible model and the parameter estimation was very crude

Extending the Basic SIR Model (1): More Compartments

• Reinfection: Change the \(I\rightarrow R\) transition to be \(I \rightarrow S\) (SIS model)
  • or \(I \rightarrow R \rightarrow S\) (SIRS model) if there’s a period of immunity after infection
• Can add an “exposed” compartment \(E\) for “have the disease and will become \(I\), but not yet”
  • SEIR models are very standard when the disease takes some time to develop
• Can add an “infectious but not showing symptoms” compartment that acts like \(I\)
  • Might be less (or more!) effective at spreading than the classic \(I\)s

Extending the Basic SIR Model (2): Time and Place

• Can make the basic coefficients \(\beta\) and \(\gamma\) functions of time, to account for seasonal effects
• Can divide space into multiple regions, give each region its own count of each compartment, and include movement between regions
  • E.g., we might have two regions \(A\) and \(B\), with \(S_A(t)\), \(S_B(t)\), \(I_A(t)\), \(I_B(t)\), etc., and equations like \[ \frac{dI_A}{dt} = \frac{\beta}{n}S_A(t) I_A(t) - \gamma I_A(t) - \mu_{AB}I_A(t) + \mu_{BA}I_B(t) \] where the \(\mu\) coefficients represent movement from region to region
  • Or use \(\frac{d I_A}{dt} = \frac{\beta_{AA}}{n}S_A(t) I_A(t) + \frac{\beta_{AB}}{n}S_A(t) I_B(t)\) if there’s cross-contagion without migration, etc.
• The spatial versions go back to Soviet work on flu in the 1970s (Sattenspiel 2009)

Social Network Structure

• When we derived the SIR model, we said each \(S\) encounters \(k\) other individuals in a time-step of length \(h\)
• Some people have a lot more contacts than others!
• Use a social network to represent this
• \(a_{ij} = 1\) if person \(i\) knows person \(j\) (this \(i \neq I\)!), otherwise \(a_{ij}=0\)
• \(k_i = \sum_{j=1}^{n}{a_{ij}} =\) number of people \(i\) knows, degree of \(i\)
• \(p(k) =\) fraction of people with degree \(k\), the degree distribution

\(\mathcal{R}_0\) on Networks

• If you get infectious, you’re only infectious for so long, and you have some limited probability of passing on the disease to your contacts, call this \(\tau\) (for “transmission”)
  • For the SIR model, you can calculate this from \(\beta\) and \(\gamma\) (exercise!)
• We start with one transmission event, where Irene \(i\) gives the disease to Joey \(j\)
• \(\mathcal{R}_0=\) How many new infections does Joey cause, on average? \[ \mathcal{R}_0= \tau \mathbf{E}\left[ \text{degree of Joey}-1 \right] \]
  • Why \(-1\)?

The Degree Distribution for a Random Friend

• Might think \(\mathbf{E}\left[ \text{degree of Joey} \right] = \mathbf{E}\left[ \text{degree of random node} \right] = \mathbf{E}\left[ K \right] = \sum_{k=0}^{\infty}{k p(k)}\)
• Joey is not a random node
• Joey is a node reached by a random edge
• High degree nodes are extra likely to be reached by following edges
  • If Joey has twice the degree of James, Joey is twice as likely to be reached as James
  • But there are only so many nodes of high degree
• So \(\mathrm{Pr}\left( K_{\text{Joey}} = k \right) \propto k p(k)\), and, to normalize, \(\mathrm{Pr}\left( K_{\text{Joey}} = k \right) = k p(k)/\sum_{k=1}^{\infty}{k p(k)}\)
• The expected value of this distribution is going to be \(> \mathbf{E}\left[ K \right]\)
• The friendship paradox: “On average, your friends have more friends than you do”

Back to \(\mathcal{R}_0\) on Networks

• Plug in to the equation from two slides ago: \[\begin{eqnarray} \mathcal{R}_0& = & \tau \mathbf{E}\left[ \text{degree of Joey}-1 \right]\\ & = & \tau \frac{\sum_{k=1}^{\infty}{(k-1)k p(k)}}{\sum_{k=1}^{\infty}{k p(k)}}\\ & = & \tau \frac{\mathbf{E}\left[ K^2 - K \right]}{\mathbf{E}\left[ K \right]}\\ & = & \tau \frac{\mathbf{E}\left[ K^2 \right]-\mathbf{E}\left[ K \right]}{\mathbf{E}\left[ K \right]}\\ & = & \tau \frac{\mathrm{Var}\left[ K \right] + \mathbf{E}\left[ K \right]^2 - \mathbf{E}\left[ K \right]}{\mathbf{E}\left[ K \right]} \end{eqnarray}\]

(because \(\mathrm{Var}\left[ K \right] = \mathbf{E}\left[ K^2 \right] - \mathbf{E}\left[ K \right]^2\) for any variable)

Epidemic Threshold on Networks

• We’re above the epidemic threshold when \(\mathcal{R}_0> 1\), or \[ \tau\left(\frac{\mathrm{Var}\left[ K \right]}{\mathbf{E}\left[ K \right]} + \mathbf{E}\left[ K \right] - 1\right) > 1 \]
• So 3 things help the epidemic take off:
  • High transmissibility (large \(\tau\))
  • High average degree (\(\mathbf{E}\left[ K \right]\) big)
  • High variance in degree (\(\mathrm{Var}\left[ K \right]\) big)
• Lots of social networks have very high degree variance

Implications for Disease Control

• Removing the highest-degree nodes is especially effective for blocking epidemics
  • Drives down \(\mathbf{E}\left[ K \right]\)
  • Even more effective at driving down \(\mathrm{Var}\left[ K \right]\)
  • Idea goes back to Pastor-Satorras and Vespignani (2002), building off Cohen et al. (2000);Cohen et al. (2001)
• “Remove” could mean vaccinating
• “Remove” could mean isolating

Complications to the Basic Network Analysis

• Community structure: Different sorts of people have different probabilities of knowing each other
  • Usually high probability of ties within communities compared to across communities
  • Tight communities tend to favor localized outbreaks (Smith et al. 2013)
    • lots of those in the community gets sick
    • other communities not as affected
    • still some transmission across communities so still pan-demic outbreaks
• Higher-order interactions: Everyone in the same cafeteria, airplane, classroom, office… is exposed to each other, even if they don’t directly know each other
  • One approach: treat these locations as special nodes, with ties to the people who go there (Saint-Onge et al. 2020)
  • The locations will have very high degree
• More detailed dynamics instead of just “what’s \(\mathcal{R}_0\)”
  • “Agent-based” or “individual-based” models: Simulate!
  • “Mean-field” approximations: see backup

Some implications for the present situation

• Do things to reduce \(\tau\)
  • Reduce the amount of time people spend near each other
  • Increase physical distancing
  • Wear masks (which make it harder for what comes out of you when you cough or sneeze to spread to others)
• Do things to reduce \(\mathbf{E}\left[ K \right]\) and \(\mathrm{Var}\left[ K \right]\)
  • Close down events and public places
  • Stay home
• “Contact tracing” (= isolating those who are tied to those known to have been infected) could be effective, if it could be implemented well

Data issues

• Ideally we’d see \(S(t), I(t), R(t)\) at all times
  • Or any two of them, really
• Easier to observe new infections, \(I(t+h) - I(t)\), than current total infections, \(I(t)\)
• Removals by death are easier to observe than removals by recovery, so we mostly see \((R(t+h)-R(t))\times (\mathrm{death\ rate})\)
  • OK if we know the death rate and it doesn’t change
  • The death rate often changes
• The interval between measurements, say \(\Delta\), is often \(\gg h\)
• Measuring \(I(t)\) and \(R(t)\) (or their rates of change) is more complicated if testing is sporadic and/or error prone
  • Need to model test error (false positives, false negatives) and who gets tested
  • A good early example for COVID-19: Javan, Fox, and Meyers (2020)
• Parameters (especially but not just \(\beta\)) can change during the epidemic
  • Changing behavior, changing policy, environmental factors…
  • Sub-divide the data? Model behavior changes?

Connecting to Data

• Likelihood calculations are straightforward if we can measure \(I(t), R(t)\) at all times \(0, h, 2h \ldots T\)
  • Or \(I(0), R(0)\) and all the increments \(I(t+h)-I(t), R(t+h) - R(t)\)
  • May still have to optimize numerically
• Likelihood calculations already become difficult if the time between observations \(\Delta \gg h\)
  • In principle, this just defines another Markov process, with a longer interval \(\Delta\) between steps, but to get the likelihood of a \(\Delta\) step we have to sum over all possible paths of \(h\) steps adding up to it
• Other complications if we don’t observe all the compartments, and/or have a lot of noise in our observations
• Often more tractable to avoid likelihood:
  1. “Conditional least squares” (see backup)
  2. Simulation-based inference: match summary statistics, use indirect inference, etc.
• Intrinsic issues:
  • Initially, everything just looks like exponential growth, so it’s hard to discriminate between distinct models
  • Initially, everything just looks like exponential growth, so even assuming an SIR model, it’s easier to estimate \(\beta - \gamma\) than \((\beta, \gamma)\)
• Can sometimes calibrate the parameters from other sources
  • E.g., \(1/\gamma =\) average time someone is infectious, which could be determined from clinical studies / observations

Summary

• The SIR model divides the population among three “compartments”, Susceptible, Infectious, Removed
• Transitions between compartments happen at rates which vary with the current state of the population
• Even this simplest epidemic models tell us a lot about how contagious diseases work
• Refinements (more compartments, space, network structure) can yield additional insights
• Connecting to data is hard when the data is bad, but we already have the estimation tools we need
• Next time: the wider world of compartment models for demography, social phenomena, chemistry, etc.

Backup: Non-basic reproductive numbers

• If \(\mathcal{R}_0\) is the basic reproductive number, what are the other reproductive numbers?
• Remember that \(\mathcal{R}_0\) is the number of new infections directly caused by introducing one infected individual at time 0 of the epidemic
• \(\mathcal{R}_t =\) the number of extra infections directly caused by increasing \(I(t)\) by 1
• Generally, \(\mathcal{R}_t\) declines with \(t\) as there are fewer susceptibles to infect

Backup: From the stochastic to the deterministic SIR model

• Rate of change in susceptibles is \[ \frac{S(t+h) - S(t)}{h} = \frac{-C(t+h)}{h} \sim -\frac{1}{h}\mathrm{Binom}(S(t), \beta h I(t)/n) \]
• Similarly \[ \frac{I(t+h) - I(t)}{h} = \frac{C(t+h) - D(t+h)}{h} \sim \frac{1}{h} \mathrm{Binom}(S(t), \beta h I(t)/n) - \frac{1}{h}\mathrm{Binom}(I(t), \gamma h) \]
• So the expected rates of change are \[\begin{eqnarray} \mathbf{E}\left[ \frac{S(t+h) - S(t)}{h} | S(t), I(t) \right] & = & - \frac{\beta}{n} S(t) I(t)\\ \mathbf{E}\left[ \frac{I(t+h) - I(t)}{h} | S(t), I(t) \right] & = & \frac{\beta}{n} S(t) I(t) - \gamma I(t) \end{eqnarray}\]
• The variance of a binomial, \(\mathrm{Binom}(n,p)\) is \(n p(1-p)\), so as \(n\rightarrow\infty\), the fluctuations around the expected values become negligible in comparison
  • This is easy to say, but the real work of a proof goes here
• Take the limit as \(n\rightarrow \infty\), \(h\rightarrow 0\): \[\begin{eqnarray} \frac{dS}{dt} & = & - \frac{\beta}{n} S(t) I(t)\\ \frac{dI}{dt} & = & \frac{\beta}{n} S(t) I(t) - \gamma I(t)\\ \frac{dR}{dt} & = & \gamma I(t) \end{eqnarray}\]

Backup: Exponential growth at the start of the epidemic in the deterministic SIR model

• Suppose that \(I(0) \ll n\), \(R(0)=0\) and so \(S(0)=n-1\), and \(n\) is big
• Then initially \[\begin{eqnarray} \frac{dI}{dt} & = & \frac{\beta}{n}S(t)I(t) - \gamma I(t)\\ & = & \left(\beta\frac{S(t)}{n} - \gamma\right) I(t)\\ & \approx & (\beta - \gamma) I(t) \end{eqnarray}\] and we know how to solve this differential equation: \[ I(t) \approx I(0) e^{(\beta-\gamma) t} \]
• So we’ll get exponential growth just when \(\beta > \gamma\), i.e., when \(\mathcal{R}_0> 1\)
  • If \(\mathcal{R}_0< 1\), we’ll get exponential decay
• Notice that \(n\) doesn’t matter to \(I(t)\)!
  • Early in epidemic, there’s no shortage of new susceptibles for the contagion to spread to
• If we have data on \(I(t)\) over time, we can use that to estimate \(\beta-\gamma\)
  • … but not \(\beta\) or \(\gamma\) separately
  • Works best if \(I(0)\) is big enough that the deterministic approximation is holding well, so you often need to wait until there are (say) 100 cases, rather than 1

Backup: Eliminating parameters from the deterministic SIR model

• We can change our units for time to whatever we like, i.e., define \(u=t\gamma\), so \(t=u/\gamma\)
• Then (chain rule of derivatives) \[\begin{eqnarray} \frac{dR}{du} & = & \left(\frac{dR}{dt}\right) \left(\frac{dt}{du} \right)\\ & = & \left( \gamma I \right) \left(\frac{1}{\gamma} \right)\\ & = & I \end{eqnarray}\]
• Doing the same thing for all the variables, \[\begin{eqnarray} \frac{dS}{du} & = & -\frac{\beta}{n\gamma} SI\\ \frac{dI}{du} & = & \frac{\beta}{n\gamma} SI - I\\ \frac{dR}{du} & = & I \end{eqnarray}\]
• Now re-scale our counts so they’re relative to the total population, so \(s=S/n\), \(i=I/n\), \(r=R/n\): \[\begin{eqnarray} \frac{ds}{du} & = & \left(\frac{ds}{dS}\right)\left(\frac{dS}{du}\right)\\ & = & \left(\frac{1}{n}\right)\left(-\frac{\beta}{\gamma}\frac{S}{n}I\right)\\ & = & - \frac{\beta}{\gamma} si \end{eqnarray}\]
• Doing the same thing for all the variables, \[\begin{eqnarray} \frac{ds}{du} & = & -\frac{\beta}{\gamma} si\\ \frac{di}{du} & = & \frac{\beta}{\gamma} si - i\\ \frac{dr}{du} & = & i \end{eqnarray}\]
• So the only parameter left is \(\beta/\gamma=\mathcal{R}_0\)
  • Units: \(dR/dt\) has units of \([\text{people}]/[\text{time}]\), so \(\gamma\) must have units of \(1/[\text{time}]\)
  • \(dS/dt\) also has units of \([\text{people}]/[\text{time}]\), but \(SI/n\) has units of \([\text{people}]^2/[\text{people}] = [\text{people}]\), so \(\beta\) must also have units of \(1/[\text{time}]\)
  • \(\beta/\gamma\) is a unitless ratio, which makes sense for interpreting it as the number of new infections directly caused by an initial infection
• The trick of manipulating variables so that parameters are combined into unitless ratios is a common way of simplifying dynamical systems
  • Often called “dimensional analysis”
  • Boccara (2004) has many good examples

Backup: Why diseases do not always evolve to be less deadly

• Many biologists used to think diseases would tend to evolve to be more benign on their hosts
  • Killing the hosts seems bad for the pathogen’s own reproduction!
  • It’s explicit in Zinsser (1935) (a wonderful, if eccentric, old book about epidemics, along with many other things)
    • Zinsser had a lot of unusual-for-his-time ideas but this was not one of them
  • You still find non-biologists (often economists or legal academics influenced by economists…) making this argument
• More mathematically:
  • if we have two variant strains of the disease, the one with the higher \(\mathcal{R}_0\) is going to dominate
    • Because: exponential growth
  • The longer the disease leaves the host alive, the smaller \(\gamma\) is
  • Smaller \(\gamma\) means larger \(\mathcal{R}_0\) (all else being equal)
  • So, shouldn’t evolution favor less-lethal variants of the disease?
• Biologists don’t believe this any more for a very simple reason: \[ \mathcal{R}_0= \frac{\beta}{\gamma} \]
  • A mutation which kills off hosts faster (increases \(\gamma\)) but also makes the disease easier to spread (increases \(\beta\)) can raise \(\mathcal{R}_0\)!
  • Lots of what diseases do to their hosts look like adaptations which increase spread at the cost of killing off the host faster (Ewald 1996)
  • I.e., they are adaptations which increase both \(\beta\) and \(\gamma\)
• Moral: pathogen evolution is not our friend

Backup: Slightly more rigorous calculation of the epidemic threshold on networks (1)

(Adapted from Newman (2002), with fewer generating functions)

• Start by infecting a single random node, Irene
• Irene is connected to a susceptible node, Joey
• Define \(\alpha_d\) as the probability that this connection starts a chain of infections of length \(\leq d\)
• So \(1-\alpha_d\) is the probability that the chain goes for more than \(d\) steps
• The probability that there’s an epidemic = the probability that the chain goes forever = \(\lim_{d\rightarrow\infty}{(1-\alpha_d})\)
• We’ll figure this out by recursively writing \(\alpha_d\) in terms of \(\alpha_{d-1}\) and taking the limit

Backup: Slightly more rigorous calculation of the epidemic threshold on networks (2)

• So let’s try to come up with the recursion for \(\alpha_d\)
• Irene’s connection to Joey starts a chain of length \(\leq d\) in one of two ways:
  1. Irene doesn’t infect Joey
  2. Irene does infect Joey, but all chains from Joey have length \(\leq d-1\)
• The average probability of an \(S\) being infected by a neighboring \(I\) is \(\tau\), so \[ \alpha_d = (1-\tau) + \mathrm{Pr}\left( \text{Joey gets infected but only starts chains of length } \leq d-1 \right) \]
• 2nd term: Joey has \(K_{\text{Joey}}-1\) friends other than Irene (Kalpana, Karl, … Kryztof), and none of those potential transmissions can start a chain more than \(d-1\) steps long
• Those events each have probability \(\alpha_{d-1}\)
• We worked out Joey’s degree distribution earlier as \(\mathrm{Pr}\left( \text{Joey has degree} k \right) = \frac{k p(k)}{\mathbf{E}\left[ K \right]}\)
• Plugging in, \[ \alpha_d = (1-\tau) + \tau\sum_{k=1}^{\infty}{\frac{k p(k)}{\mathbf{E}\left[ K \right]}\alpha^{k-1}_{d-1}} \]
• Now define the right-hand side as \(f(\alpha_{d-1})\), i.e,, \[ f(x) \equiv (1-\tau) + \tau\sum_{k=1}^{\infty}{\frac{k p(k)}{\mathbf{E}\left[ K \right]}x^{k-1}} \] so we can say that \[ \alpha_d = f(\alpha_{d-1}) \]

Backup: Slightly more rigorous calculation of the epidemic threshold on networks (3)

How does any of this help???

• “Starts a chain of length \(\leq d-1\)” implies “Starts a chain of length \(\leq d\)”, so \(\alpha_{d-1} \leq \alpha_d\), and \(\alpha_d\) is non-decreasing as \(d\rightarrow\infty\)
  • Alternately: use algebra to check that \(f(x) \geq x\)
• But \(\alpha_d\) is a probability so \(\alpha_d \leq 1\) for all \(d\)
• Since \(\alpha_d\) is non-decreasing and bounded above, \(\lim_{d\rightarrow\infty}{\alpha_d} \equiv \alpha_{\infty}\) must exist and be \(\leq 1\)
• The limit must be a fixed point of \(f\): \[ \alpha_{\infty}=f(\alpha_{\infty}) \]
• It’s easy to check that \(f(1)=1\) always
  • If \(\alpha_{\infty}=1\), then any chain of infections must stop in a finite number of steps, and the probability of an infinite epidemic is 0
• If \(f^{\prime}(1) < 1\), then the \(f(1)=1\) solution is stable, and this is what \(\alpha_{\infty}\) will converge to
• If \(f^{\prime}(1) > 1\), then the fixed point \(f(1)=1\) is unstable, and \(\alpha_d\) will converge to the other solution of \(f(x)=x\)
  • Proof that there is another solution: what is this, a math class?
• The derivative is \[ \frac{df}{dx} = \tau\sum_{k=1}^{\infty}{(k-1)\frac{k p(k)}{\mathbf{E}\left[ K \right]} x^{k-2}} \]
• So the epidemic criterion is \[ f^{\prime}(1) = \tau \sum_{k=1}^{\infty}{(k-1)\frac{k p(k)}{\mathbf{E}\left[ K \right]}} > 1 \] or \[ \tau \left(\frac{\mathbf{E}\left[ K^2 \right] - \mathbf{E}\left[ K \right]}{\mathbf{E}\left[ K \right]}\right) > 1 \] exactly as we got before

Backup: “Mean-field” approximations to epidemic models on networks

• Keep track of number in each compartment with each degree, so \(S_2(t)\) or \(I_1(t)\)
  • Idea goes back to May and Anderson (1988) (if not before)
• Similarly say \(S_k I_l(t)\) is the number of \(S\)s who have degree \(k\) and are connected to \(I\)s of degree \(l\) at time \(t\)
• Then \(\frac{dS_k}{dt} = -\sum_{l}{\beta S_kI_l(t)}\) and similarly for the other first-order variables
• But now we need to know \(S_k I_l(t)\)
• Option 1, the “(pure) mean-field” approximation: \(S_k I_l \approx \frac{kl}{n\mathbf{E}\left[ K \right]} S_k I_l\)
  • The susceptibles of degree \(k\) “send out” \(k S_k\) edges at random
  • There are \(n \mathbf{E}\left[ K \right]\) incoming edges those could match
  • \(l I_l\) of those incoming edges belong to infectious of degree \(l\)
• Option 2, the “pairwise mean-field” approximation
  • Introduce \(S_k S_l I_m\) (and so forth), so \[ \frac{d S_k I_l}{dt} = -\gamma S_k I_l + \tau\sum_{m}{S_k S_l I_m - I_m S_k I_l - S_k I_l} \]
  • Now approximate triples like \(S_k S_l I_m\) in terms of pairs and singles
    • similar, if more elaborate, counting arguments as in the simple mean-field approach above
• The pairwise mean-field approximation often gives pretty good matches to individual-level simulations, with a lot less computational work
• Details: Kiss, Miller, and Simon (2017)

Backup: Conditional least squares

• Say we observe \(X(0), X(\Delta), X(2\Delta), \ldots\)
• Everything observed before time \(t\) is \(\mathcal{H}(t)\) (“history”)
• Work out \(\mathbf{E}\left[ X(t)|\mathcal{H}(t); \theta \right] \equiv \mu(t;\theta)\)
  • Usually a lot more tractable than finding the full distribution for \(X(t)\) given history
  • For instance, use the limiting deterministic system, started from whatever state \(\mathcal{H}(t)\) left the system in
• \(MSE = T^{-1} \sum_{t}{(X(t) - \mu(t;\theta))^2}\)
• Adjust \(\theta\) to minimize MSE
• Even better: also work out \(V(t;\theta) \equiv \mathrm{Var}\left[ X(t)|\mathcal{H}(t) \right]\), form \(M = T^{-1}\sum_{t}{(X(t) - \mu(t;\theta))^2/V(t;\theta)}\) and minimize \(M\)
  • Sometimes called “quasi-likelihood” (Heyde 1997)

Further Reading

• SIR or “compartment” type models of epidemics go back to the 1920s, especially to work on malaria
  • See Bacaër (2011), chapters 12 and 16, for some of the history
• For a modern textbook treatment of epidemic models, see Daley and Gani (1999)
  • There’s also a good chapter in Ellner and Guckenheimer (2006), which covers models of many other biological processes
  • Britton and Pardoux (2019) is a state-of-the-art review of the rigorous theory of epidemic models and their statistical inference, but presumes comfort with graduate-level probability and stochastic processes (like Fristedt and Gray (1997) or Kallenberg (2002), or at least Øksendal (1995))
• I have only touched on the spatial aspects, which are obviously important; see Sattenspiel (2009) for a really clear introduction and literature review, which does a nice job of balancing the theory against the applications
• Branching processes (which I used to derive \(\mathcal{R}_0\) from the SIR model) are an important part of applied probability theory; there’s a good treatment in Grimmett and Stirzaker (1992), among many other places
  • See also Bacaër (2011), chapter 9, on their origin in the 19th century
  • Guttorp (1995), chapter 2, goes over some statistical inference for branching processes
• The way I derived \(\mathcal{R}_0\) for the network model follows Newman (2002), but without (explicitly) using generating functions
  • It’s an example of using percolation theory, another important topic in applied probability (Grimmett 1999), with many surprising applications to epidemiology (Davis et al. 2008)
  • See Bacaër (2011), chapter 22 for some of the 20th century history
• Kiss, Miller, and Simon (2017) is an excellent textbook on epidemics on networks, but assumes more knowledge of differential equations than is usual for statistics students
  • You could always read, say, Arnol’d (1973)
• Simulating epidemic models at the individual level is a really good way to build insight and/or to handle complicated situations which defy exact solutions
  • Bartlett (1960) was one of the first to appreciate the power of simulations in this area
  • This is now a huge area of active research, e.g., the recent paper Gallagher et al. (2018) from our department, just on the (important) topic of initializing the simulation