finalsize provides methods to calculate the final proportion of individuals in different age groups infected during an epidemic, as predicted by an age-stratified SIR model. In this vignette we lay out the mathematical concepts behind the functionality available in the package.
The dynamics of the proportion of susceptible and infectious individuals over time in the standard susceptible-infectious-recovered (SIR) epidemic model are given by the following ordinary differential equations: $$ \frac{d S(t)}{d t} = -\beta S(t) I(t) \\ \frac{d I(t)}{d t} = \beta S(t) I(t) -\gamma I(t) $$ where β is the transmission rate and γ is the rate of recovery.
To calculate the proportion infected, we need to derive an expression for the proportion of susceptibles who got infected during the epidemic, i.e. ϕ = 1 − S(∞)/S(0). Rather than simulating the model over a very long time period (i.e. t → ∞), we can instead calculate dI/dS and integrate:
$$
\frac{d I}{d S} = - 1 + \frac{\gamma}{\beta S} \\
I(t) = - S(t) + \frac{\gamma}{\beta} \log S(t) + c
$$ where c depends on
initial conditions. Under the assumption the population is fully
susceptible initially (i.e. S(0) = 1), and I(t) → 0 as t → ∞, we obtain the following:
$$
\log S(\infty) = - R_0 [1-S(\infty) ] \\
\mathrm{log}~\phi -R_0 (1-\phi) = 0 = F(\phi)
$$ This equation doesn’t have an analytical solution, but there
are various methods available to allow it to be solved numerically to
find the unique solution in (0,1]. For example, we can use Newton’s
method (solver = "newton" in the final_size()
function) to find the final epidemic size:
We can extend the above formulation to account for variable susceptibility among individuals, following the methods outlined in Miller (2012). First, we define p(x) to be the probability density for a randomly chosen individual to be of susceptibility type x. If we define the probability an individual of type x to be ultimately infected as ϕ(x), then the overall proportion of the population to be infected is defined as follows, depending on whether p(x) is a continuous or discrete distribution:
$$ \hat\phi = \int_0^\infty \phi(x) p(x) dx \\ \text{or} \\ \hat\phi = \sum_x \phi(x) p(x) dx \\ $$
If we define xϕ̂ to be the expected number of transmissions to a given individual of type x, we can therefore calculate the probability that this individual remains susceptible. This gives the following relationship between ϕ and ϕ̂:
$$ 1-\phi(x) = S(x,0) e^{-x \hat\phi } \\ $$ And hence we can combine the above to give: $$ \hat\phi = \int_0^\infty (1- S(x,0) e^{-x \hat\phi }) p(x) dx \\ \hat\phi = 1- \int_0^\infty S(x,0) e^{-x \hat\phi } p(x) dx $$
Under the assumption that S(x, 0) is near 1 initially (i.e. f), we therefore obtain: ϕ̂ = 1 − ∫0∞p(x)e−xϕ̂dx
We can extend the same derivation for the above homogenous model for populations with multiple age groups, following Andreasen (2011). The age-dependent ODEs are defined as follows:
$$ \frac{d S(a,t)}{d t} = -S(a,t) \sum_b \beta_{ab} I(b,t) \\ \frac{d I(a,t)}{d t} = S(a,t) \sum_b \beta_{ab} I(b,t) -\gamma I(a,t) $$ Where βa, b is the transmission rate from group b to group a. If we define Ra, b = βa, b/γ to be group specific reproduction number, we can rearrange the above and integrate to derive an expression for the final size in each age group:
$$ - \gamma \sum_b R_{a,b} I(b,t) = \frac{d S(a,t)}{d t} \frac{1}{S(a,t)} \\ - \gamma \sum_b R_{a,b} \int_0^\infty I(b,t) ~dt = \int_0^\infty \frac{d S(a,t)}{d t} \frac{1}{S(a,t)} ~dt = \mathrm{log}~S(a,\infty) -\mathrm{log}~S(a,0) = \mathrm{log}~\frac{S(a,\infty)}{S(a,0)} \\ -\gamma \int_0^\infty I(a,t) ~dt = \int_0^\infty \frac{d S(a,t)}{d t} +\frac{dI(a,t)}{dt} ~dt = S(a,\infty) -S(a,0) ~. $$ and hence:
$$ \mathrm{log}~\frac{S(a,\infty)}{S(a,0)} = - \sum_b R_{a,b} [ S(b,\infty) -S(b,0)] $$
If we define $\underline\phi = (\phi_1, ... \phi_n)$ to be a vector of final sizes in each of n age group, we can rewrite the above as
$$ \mathrm{log}~\underline\phi -\mathbf{R} (1-\underline\phi) = 0 = F(\underline\phi) $$ where R is a matrix with entries Ri, j. Once again, we can use Newton’s method, this time with $\underline\phi$ as a vector, and hence we need to solve a set of linear equations in step 2.