Compute a reproduction number from a daily contact matrix with group-specific infectiousness and susceptibility

Compute the dominant eigenvalue of a next-generation matrix built from a daily contact matrix, a baseline per-contact transmission probability, an average infectious period, and optional group-specific infectiousness and susceptibility multipliers.

Usage

compute_reproduction_number(
  contact_matrix,
  transmission_prob,
  infectious_period_days = 1,
  infectiousness = NULL,
  susceptibility = NULL,
  check_reciprocity = FALSE
)

Arguments

contact_matrix: A square numeric matrix. Entry [i, j] must be the average number of daily contacts made by an infectious individual in group i with individuals in group j.
transmission_prob: Numeric scalar in [0, 1]. Baseline infection probability per contact.
infectious_period_days: Positive numeric scalar. Mean infectious period, in days. Since the contact matrix is daily, this scales daily contact opportunities to expected secondary infections over the infectious period.
infectiousness: Optional numeric vector of length nrow(contact_matrix). Entry i is the relative infectiousness multiplier for source group i. If NULL, a vector of ones is used.
susceptibility: Optional numeric vector of length nrow(contact_matrix). Entry j is the relative susceptibility multiplier for recipient group j. If NULL, a vector of ones is used.
check_reciprocity: Logical. If TRUE, the function issues a warning when the contact matrix is not symmetric. This is only a diagnostic and does not alter the calculation. Note that reciprocity for empirical contact matrices is often assessed after adjusting for group sizes, so lack of symmetry is not automatically an error.

Value

A list with the following elements:

R: Numeric scalar. Dominant eigenvalue (spectral radius) of the next-generation matrix.
type: Character string, either "R0" for the homogeneous case or "Reff" when group-specific modifiers are used.
next_generation_matrix: Numeric matrix used in the calculation.
infectiousness: Numeric vector of source-group infectiousness multipliers.
susceptibility: Numeric vector of recipient-group susceptibility multipliers.
eigenvalues: Complex vector of eigenvalues of the next-generation matrix.

Details

Let $C$ be a daily contact matrix where $c_{ij}$ is the average number of daily contacts made by one infectious individual in group $i$ with individuals in group $j$. If $p$ is a baseline infection probability per contact, $D$ is the mean infectious period in days, $b_i$ is the infectiousness multiplier for source group $i$, and $a_j$ is the susceptibility multiplier for recipient group $j$, then the expected next-generation matrix is

$$K = p \cdot D \cdot \mathrm{diag}(b_1, \ldots, b_n) \cdot C \cdot \mathrm{diag}(a_1, \ldots, a_n)$$

or element-wise

$$k_{ij} = p \cdot D \cdot b_i \cdot c_{ij} \cdot a_j$$

The reproduction number is the spectral radius of $K$, that is, its dominant eigenvalue:

$$R = \rho(K)$$

When infectiousness and susceptibility are both vectors of ones, the function reduces to the homogeneous case

$$R = \rho(p \cdot D \cdot C)$$

Vaccination can be represented indirectly through these inputs. For example, if vaccination reduces susceptibility only, a common choice is $a_j = 1 - e_j v_j$, where $v_j$ is coverage and $e_j$ is vaccine efficacy against susceptibility in group $j$.

Because the contact matrix is assumed to be daily, infectious_period_days should represent the mean number of days an infected individual remains infectious enough to transmit. In simple memoryless models this is often approximated by $D = 1 / \gamma$, where $\gamma$ is a per-day recovery rate in continuous time or, approximately, a daily recovery probability in a discrete-time model.

This function implements a next-generation-matrix calculation. It is most appropriate near the start of an outbreak, when the susceptible composition and contact structure are approximately stable over the time window of interest.

The calculation assumes:

the contact matrix is already on a per-day basis,
transmission_prob is a baseline per-contact infection probability,
infectiousness modifies the source side of transmission,
susceptibility modifies the recipient side of transmission,
a mean infectious period is an adequate summary of infectious duration.

Interpretation note: if infectiousness and susceptibility are both equal to one for all groups and the population is fully susceptible, the returned value corresponds to $R_0$. Otherwise, it is better interpreted as an initial or effective reproduction number under the supplied group-specific transmission modifiers.

For the per-contact interpretation to remain meaningful, the implied group-pair-specific transmission probabilities

$$p_{ij} = p \cdot b_i \cdot a_j$$

should ideally remain in $[0, 1]$.

References

Diekmann O, Heesterbeek JAP, Roberts MG (2010). "The construction of next-generation matrices for compartmental epidemic models." Journal of the Royal Society Interface, 7(47), 873-885. doi:10.1098/rsif.2009.0386

Wallinga J, Teunis P, Kretzschmar M (2006). "Using data on social contacts to estimate age-specific transmission parameters for respiratory-spread infectious agents." American Journal of Epidemiology, 164(10), 936-944. doi:10.1093/aje/kwj317

van den Driessche P (2017). "Reproduction numbers of infectious disease models." Infectious Disease Modelling, 2(3), 288-303. doi:10.1016/j.idm.2017.06.002

Examples

C <- matrix(c(
  8, 2, 1,
  3, 7, 2,
  1, 2, 5
), nrow = 3, byrow = TRUE)

# Homogeneous case
compute_reproduction_number(
  contact_matrix = C,
  transmission_prob = 0.04,
  infectious_period_days = 5
)
#> $R
#> [1] 2.149744
#> 
#> $type
#> [1] "R0"
#> 
#> $next_generation_matrix
#>      [,1] [,2] [,3]
#> [1,]  1.6  0.4  0.2
#> [2,]  0.6  1.4  0.4
#> [3,]  0.2  0.4  1.0
#> 
#> $infectiousness
#> [1] 1 1 1
#> 
#> $susceptibility
#> [1] 1 1 1
#> 
#> $eigenvalues
#> [1] 2.1497444 1.1079887 0.7422669
#> 

# Group-specific infectiousness and susceptibility
compute_reproduction_number(
  contact_matrix = C,
  transmission_prob = 0.04,
  infectious_period_days = 5,
  infectiousness = c(1.0, 0.8, 1.2),
  susceptibility = c(0.5, 0.9, 0.7)
)
#> $R
#> [1] 1.409109
#> 
#> $type
#> [1] "Reff"
#> 
#> $next_generation_matrix
#>      [,1]  [,2]  [,3]
#> [1,] 0.80 0.360 0.140
#> [2,] 0.24 1.008 0.224
#> [3,] 0.12 0.432 0.840
#> 
#> $infectiousness
#> [1] 1.0 0.8 1.2
#> 
#> $susceptibility
#> [1] 0.5 0.9 0.7
#> 
#> $eigenvalues
#> [1] 1.4091085 0.6849710 0.5539205
#>