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!

`library(finalsize)`

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

Prepare \(R_0\).

`<- 2.0 r0 `

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
<- 67e6
population_size
# estimate using finalsize::final_size
<- final_size(
final_size_finalsize 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_r(
final_size_sir 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.

```
# function to check error in final size estimates
<- function(x, r0, tolerance = 1e-6) {
final_size_error <- abs(x - (1.0 - exp(-r0 * x)))
error if (any(error > tolerance)) {
print("Final size estimate error is greater than tolerance!")
}# return error
error }
```

```
# error for estimate using finalsize::final_size()
final_size_error(final_size_finalsize$p_infected, r0)
#> [1] 0
# error for estimate from SIR model
final_size_error(final_size_sir, r0)
#> [1] "Final size estimate error is greater than tolerance!"
#> [1] 3.915923e-05
```

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.

```
# Load example POLYMOD data included with the package
data(polymod_uk)
# Define contact matrix (entry {ij} is contacts in group i
# reported by group j)
<- polymod_uk$contact_matrix
contact_matrix
# Define population in each age group
<- polymod_uk$demography_vector demography_vector
```

Examine the contact matrix and demography vector.

```
contact_matrix#>
#> contact.age.group [,1] [,2] [,3]
#> [0,20) 4.514575e-08 1.600071e-08 8.965770e-09
#> [20,40) 1.600071e-08 2.489595e-08 1.346051e-08
#> 40+ 8.965770e-09 1.346051e-08 1.464763e-08
demography_vector#> [1] 14799290 16526302 28961159
```

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).

```
# Define susceptibility of each group
<- matrix(
susceptibility data = c(1.0, 0.5, 0.5),
nrow = length(demography_vector),
ncol = 1
)
# Assume uniform susceptibility within age groups
<- matrix(
p_susceptibility data = 1.0,
nrow = length(demography_vector),
ncol = 1
)
```

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
<- matrix(
susceptibility 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
<- 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]
p_susceptibility[,
# 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
<- matrix(1, nrow = length(demography_vector), ncol = 2)
susceptibility 2] <- 0.0
susceptibility[,
# 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
<- 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]
p_susceptibility[,
# 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.