In [1]:

# load finalsize
library(finalsize)
library(socialmixr)

Attaching package: 'socialmixr'

The following object is masked from 'package:utils':

    cite

library(ggplot2)

In [2]:

# get UK polymod data from socialmixr
polymod <- socialmixr::polymod
contact_data <- socialmixr::contact_matrix(
  polymod,
  countries = "United Kingdom",
  age.limits = c(0, 5, 18, 40, 65),
  symmetric = TRUE
)

Using POLYMOD social contact data. To cite this in a publication, use the 'get_citation()' function

Removing participants that have contacts without age information. To change this behaviour, set the 'missing.contact.age' option

Warning in pop_age(survey.pop, part.age.group.present, ...): Not all age groups represented in population data (5-year age band).
  Linearly estimating age group sizes from the 5-year bands.

# get the contact matrix and demography data
contact_matrix <- t(contact_data$matrix)
demography_vector <- contact_data$demography$population

# scale the contact matrix so the largest eigenvalue is 1.0
contact_matrix <- contact_matrix / max(Re(eigen(contact_matrix)$values))

# divide each row of the contact matrix by the corresponding demography
contact_matrix <- contact_matrix / demography_vector

n_demo_grps <- length(demography_vector)

In [3]:

# mean R0 is 1.5
r0_mean <- 1.5

In [4]:

# susceptibility is uniform
susc_uniform <- matrix(
  data = 1,
  nrow = n_demo_grps,
  ncol = 1L
)

# p_susceptibility is uniform
p_susc_uniform <- susc_uniform

In [5]:

# create an R0 samples vector
r0_samples <- rnorm(n = 1000, mean = r0_mean, sd = 0.1)

In [6]:

# run final size on each sample with the same data
final_size_data <- Map(
  r0_samples, seq_along(r0_samples),
  f = function(r0, i) {
    # the i-th final size estimate
    fs <- final_size(
      r0 = r0,
      contact_matrix = contact_matrix,
      demography_vector = demography_vector,
      susceptibility = susc_uniform,
      p_susceptibility = p_susc_uniform
    )

    fs$replicate <- i
    fs$r0_estimate <- r0
    fs
  }
)

# combine data
final_size_data <- Reduce(x = final_size_data, f = rbind)

# order age groups
final_size_data$demo_grp <- factor(
  final_size_data$demo_grp,
  levels = contact_data$demography$age.group
)

# examine some replicates in the data
head(final_size_data, 15)

   demo_grp   susc_grp susceptibility p_infected replicate r0_estimate
1     [0,5) susc_grp_1              1  0.4378500         1    1.536741
2    [5,18) susc_grp_1              1  0.7111826         1    1.536741
3   [18,40) susc_grp_1              1  0.5625151         1    1.536741
4   [40,65) susc_grp_1              1  0.4875332         1    1.536741
5       65+ susc_grp_1              1  0.3227313         1    1.536741
6     [0,5) susc_grp_1              1  0.4374186         2    1.535997
7    [5,18) susc_grp_1              1  0.7107502         2    1.535997
8   [18,40) susc_grp_1              1  0.5620324         2    1.535997
9   [40,65) susc_grp_1              1  0.4870674         2    1.535997
10      65+ susc_grp_1              1  0.3223626         2    1.535997
11    [0,5) susc_grp_1              1  0.4528238         3    1.563006
12   [5,18) susc_grp_1              1  0.7259570         3    1.563006
13  [18,40) susc_grp_1              1  0.5791672         3    1.563006
14  [40,65) susc_grp_1              1  0.5036707         3    1.563006
15      65+ susc_grp_1              1  0.3355975         3    1.563006