The final size equation is derived analytically from a compartmental epidemiological model, the susceptible-infectious-recovered (SIR) model. Yet it can only be solved numerically — thus with perfect solvers the final epidemic size from an SIR model and a final size calculation should be the same.

This vignette compares the final epidemic sizes from a simple SIR model with those estimated using `finalsize::final_size()`

.

**New to finalsize?** It may help to read the “Get started” vignette first!

First, define a simple SIR model and a helper function to replicate the contact matrix by rows and columns. Here, the model is defined as an iteration of transitions between the epidemiological compartments. These transitions form a system of ordinary differential equations (ODEs), whose solution through time could also have been obtained by using a solver such as `deSolve::lsoda()`

.

**Interested in compartmental modelling to obtain timeseries of epidemic outcomes?**

The *epidemics* package might be useful!

```
# Replicate the contact matrix
repmat <- function(X, m, n) {
mx <- dim(X)[1]
nx <- dim(X)[2]
matrix(t(matrix(X, mx, nx * n)), mx * m, nx * n, byrow = TRUE)
}
# Calculate the proportion of individuals infected
final_size_r <- function(r0, contact_matrix,
demography_vector, susceptibility, p_susceptibility) {
# check p_susceptibility
stopifnot(all(rowSums(p_susceptibility) == 1))
p_susc <- as.vector(p_susceptibility)
susc <- as.vector(susceptibility)
demo_v <- rep(demography_vector, ncol(p_susceptibility)) * p_susc
nsteps <- 10000
# epidemiological parameters
h <- 0.1
beta <- r0
gamma <- 1
# initial conditions
I <- 1e-8 * demo_v
S <- demo_v - I
R <- 0
# prepare contact matrix
m <- repmat(contact_matrix, ncol(p_susceptibility), ncol(p_susceptibility))
# multiply contact matrix by susceptibility
m <- m * susc
# iterate over compartmental transitions
for (i in seq_len(nsteps)) {
force <- (m %*% I) * S * h * beta
S <- S - force
g <- I * gamma * h
I <- I + force - g
R <- R + g
}
# return proportion of individuals recovered
as.vector(R / demo_v)
}
```

Prepare \(R_0\).

Here, the assumption is of a uniformly mixing population where social contacts do not differ among demographic groups (here, age groups).

```
# estimated population size for the UK is 67 million
population_size <- 67e6
# estimate using finalsize::final_size
final_size_finalsize <- final_size(
r0 = r0,
contact_matrix = matrix(1) / population_size,
demography_vector = population_size,
susceptibility = matrix(1),
p_susceptibility = matrix(1),
solver = "newton"
)
# estimate from SIR model
final_size_sir <- final_size_r(
r0 = r0,
contact_matrix = 0.99 * matrix(1) / population_size,
demography_vector = population_size,
susceptibility = matrix(1),
p_susceptibility = matrix(1)
)
# View the estimates
final_size_finalsize
#> demo_grp susc_grp susceptibility p_infected
#> 1 demo_grp_1 susc_grp_1 1 0.7968121
final_size_sir
#> [1] 0.7968781
```

There is a small discrepancy between the final epidemics sizes returned by `final_size()`

and the SIR model. One quick check to identify the source of the error would be to calculate the difference between the two sides of the finalsize equation, i.e. \(x\), the estimated final size, against \(1 - e^{-R_0x}\) (for a difference of \(x - (1 - e^{-R_0x})\)) to see which of the final size estimates is incorrect — a small error, below the error tolerance value (\(1^{-6}\)), would indicate that the method being checked is more or less correct, while a larger error would suggest the method is not correct.

This can be written as a temporary function.

This suggests that it is the SIR model which makes a small error in the final size estimate, while `finalsize::final_size()`

provides an estimate that is within the specified solver tolerance (\(1^{-6}\)).

Here, the comparison considers scenarios in which social contacts vary by demographic groups, here, age groups. See the vignette on “Modelling heterogeneous social contacts” for more guidance on how to implement such scenarios.

First, prepare the population data using the dataset provided with *finalsize*. This includes preparing the population’s age-specific contact matrix and demography distribution vector, as well as the susceptibility and demography-in-susceptibility distribution vectors. See the “Guide to constructing susceptibility matrices” for help understanding this latter step, and see the the vignette on “Modelling heterogeneous susceptibility” for worked out examples.

Examine the contact matrix and demography vector.

In this scenario, susceptibility to infection only varies between groups, with older age groups (20 and above) only half as susceptible as individuals aged (0 – 19).

Examine the final size estimates.

```
# from {finalsize} using finalsize::final_size
final_size(
r0 = r0,
contact_matrix = contact_matrix,
demography_vector = demography_vector,
susceptibility = susceptibility,
p_susceptibility = p_susceptibility,
solver = "newton"
)
#> demo_grp susc_grp susceptibility p_infected
#> 1 [0,20) susc_grp_1 1.0 0.7119330
#> 2 [20,40) susc_grp_1 0.5 0.3252457
#> 3 40+ susc_grp_1 0.5 0.2334420
# using SIR model
final_size_r(
r0 = r0,
contact_matrix = contact_matrix * 0.99,
demography_vector = demography_vector,
susceptibility = susceptibility,
p_susceptibility = p_susceptibility
)
#> [1] 0.7057450 0.3187776 0.2281473
```

The final size estimates are very close, and the difference does not fundamentally alter our understanding of the potential outcomes of an epidemic in this population.

In this scenario, susceptibility to infection varies within age groups, with a fraction of each age group much less susceptible. Additionally, the proportion of individuals with low susceptibility varies between groups, with individuals aged (0 – 19) much more likely to have low susceptibility than older age groups.

```
# Define susceptibility of each group
# Here
susceptibility <- matrix(
c(1.0, 0.1),
nrow = length(demography_vector), ncol = 2,
byrow = TRUE
)
# view the susceptibility matrix
susceptibility
#> [,1] [,2]
#> [1,] 1 0.1
#> [2,] 1 0.1
#> [3,] 1 0.1
# Assume variation in susceptibility within age groups
# A higher proportion of 0 -- 20s are in the lower susceptibility group
p_susceptibility <- matrix(1, nrow = length(demography_vector), ncol = 2)
p_susceptibility[, 2] <- c(0.6, 0.5, 0.3)
p_susceptibility[, 1] <- 1 - p_susceptibility[, 2]
# view p_susceptibility
p_susceptibility
#> [,1] [,2]
#> [1,] 0.4 0.6
#> [2,] 0.5 0.5
#> [3,] 0.7 0.3
```

Examine the final size estimates for this scenario.

```
# estimate group-specific final sizes using finalsize::final_size
final_size(
r0 = r0,
contact_matrix = contact_matrix,
demography_vector = demography_vector,
susceptibility = susceptibility,
p_susceptibility = p_susceptibility
)
#> demo_grp susc_grp susceptibility p_infected
#> 1 [0,20) susc_grp_1 1.0 0.27211694
#> 2 [20,40) susc_grp_1 1.0 0.24819087
#> 3 40+ susc_grp_1 1.0 0.19449495
#> 4 [0,20) susc_grp_2 0.1 0.03126239
#> 5 [20,40) susc_grp_2 0.1 0.02812422
#> 6 40+ susc_grp_2 0.1 0.02139636
# estimate group-specific final sizes from the SIR model
final_size_r(
r0 = r0,
contact_matrix = contact_matrix,
demography_vector = demography_vector,
susceptibility = susceptibility,
p_susceptibility = p_susceptibility
)
#> [1] 0.27261723 0.24864863 0.19485488 0.03132008 0.02817620 0.02143593
```

Here too, the final size estimates are very similar, and provide essentially the same understanding of the potential outcomes of an epidemic in this population.

In this scenario, a fraction of each age group is completely immune to infection. Additionally, individuals aged (0 – 19) are much more likely to have complete immunity than older age groups.

```
# Define susceptibility of each group
susceptibility <- matrix(1, nrow = length(demography_vector), ncol = 2)
susceptibility[, 2] <- 0.0
# view susceptibility
susceptibility
#> [,1] [,2]
#> [1,] 1 0
#> [2,] 1 0
#> [3,] 1 0
# Assume that some proportion are completely immune
# complete immunity is more common among younger individuals
p_susceptibility <- matrix(1, nrow = length(demography_vector), ncol = 2)
p_susceptibility[, 2] <- c(0.6, 0.5, 0.3)
p_susceptibility[, 1] <- 1 - p_susceptibility[, 2]
# view p_susceptibility
p_susceptibility
#> [,1] [,2]
#> [1,] 0.4 0.6
#> [2,] 0.5 0.5
#> [3,] 0.7 0.3
```

```
final_size(
r0 = r0,
contact_matrix = contact_matrix,
demography_vector = demography_vector,
susceptibility = susceptibility,
p_susceptibility = p_susceptibility
)
#> demo_grp susc_grp susceptibility p_infected
#> 1 [0,20) susc_grp_1 1 0.07787362
#> 2 [20,40) susc_grp_1 1 0.07227733
#> 3 40+ susc_grp_1 1 0.05594389
#> 4 [0,20) susc_grp_2 0 0.00000000
#> 5 [20,40) susc_grp_2 0 0.00000000
#> 6 40+ susc_grp_2 0 0.00000000
final_size_r(
r0 = r0,
contact_matrix = contact_matrix,
demography_vector = demography_vector,
susceptibility = susceptibility,
p_susceptibility = p_susceptibility
)
#> [1] 0.07791032 0.07231121 0.05597004 0.00000001 0.00000001 0.00000001
```

Here, the method from `finalsize::final_size()`

correctly provides estimates of 0.0 for the susceptibility group (in each age group) that is completely immune, while the SIR model method estimates that a very small fraction of the completely immune groups are still infected. This discrepancy is due to the assumption of initial conditions for the SIR model, wherein the same fraction (\(1^{-8}\)) are assumed to be initially infected.