DiseasyRegions

library(diseasy)
#> Loading required package: diseasystore
#> 
#> Attaching package: 'diseasy'
#> The following object is masked from 'package:diseasystore':
#> 
#>     diseasyoption

Overview

The DiseasyRegions module defines the regional scope for the models. In addition, it provides a common place to define the demography available for those areas, and the adjacency between regions within the areas.

Data structure

DiseasyRegions uses three inputs:

• area: a character vector defining the selected regional scope.
• demography: a data.frame with one “region” column, optional stratification columns, and one population column.
• adjacency: a data.frame describing connectedness between pairs of regions (see the Regional adjacency section for further details).

The demography data must include the columns region and population. Any additional columns are treated as stratification variables.

demography_nordic
#> # A tibble: 505 × 3
#>   region   age population
#>   <chr>  <dbl>      <dbl>
#> 1 DK         0      60965
#> 2 DK         1      61357
#> 3 DK         2      61663
#> 4 DK         3      62128
#> 5 DK         4      61148
#> # ℹ 500 more rows

The adjacency data must include exactly the columns from, to, and adjacency. The from and to columns should contain the same set of region identifiers.

adjacency_meta_nordic
#> # A tibble: 25 × 3
#>   from  to    adjacency
#>   <chr> <chr>     <dbl>
#> 1 NO    NO       6413. 
#> 2 NO    FI         13.7
#> 3 FI    NO         13.7
#> 4 SE    NO         52.1
#> 5 NO    SE         52.1
#> # ℹ 20 more rows

Regional adjacency

DiseasyRegions builds on the concept of a degree of “connectedness” between the different regions. Here, this “connectedness” is referred to as “adjacency” in alignment with other models.

Different measures for adjacency can be extracted from real-world data which have very different interpretations. Adjacency is often described from information about “movement” such as mobility data or surveys on how people move between regions. Alternatively adjacency can be derived from “infection” data such as genetic phylogeny.

While these seem very similar, there are substantial differences in how these adjacencies are expressed numerically, and how they should be implemented in disease spread models. Our framework assumes a reality where disease is transported between regions by people returning to their own region (travel or commuting), that is: there is no immigration or emigration. If the supplied adjacency describes a movement matrix it is furthermore assumed that people travelling to a place transmits disease equally with people travelling from a place. Which in turns means that the resulting infection matrix used in the models will be symmetric.

In diseasy, the adjacency of DiseasyRegions are implemented in SEIR-like models and it helps to be clear on how this adjacency is used to illustrate these differences.

To begin, let’s define the “movement matrix” $\Phi$ between $R$ connected regions:

$$\overline{\overline{\Phi}} = \begin{bmatrix} \phi_{1,1} & \cdots & \phi_{1,R} \\ \vdots & \ddots & \vdots \\ \phi_{R,1} & \cdots & \phi_{R,R} \end{bmatrix}$$

If we interpret these adjacencies as “the average fraction of contacts between regions”, we can diagram the SEIR compartments of a two-region model as:

SEIR model overview for a two region model

With the following equations for the infectious compartments in region $1$$$\dot{I}_1 = \beta\left( \underbrace{S_1 \phi_{1,1}}_\textrm{stay} \underbrace{\phi_{1,1} I_1}_\textrm{stay} + \underbrace{S_1 \phi_{1,1}}_\textrm{stay} \underbrace{\phi_{1,2} I_2}_\textrm{travel} + \underbrace{S_1 \phi_{2,1}}_\textrm{travel}\underbrace{\phi_{2,2} I_2}_\textrm{stay} + \underbrace{S_1 \phi_{2,1}}_\textrm{travel}\underbrace{\phi_{2,1} I_1}_\textrm{travel}\right) - r_I I_1$$ where we divide the types of contacts into whether individuals stay in their own region or travel¹ to other regions.

In general, we can express the dynamics in matrix form assuming $R$ different regions $$\dot{\overline{I}} = \beta\overline{S}\circ\overline{\overline{\Theta}}\cdot\overline{I} - r_I \overline{I}$$

Where $\circ$ is the Hadamard product (element-wise multiplication).

The matrix elements of $\overline{\overline{\Theta}}$ contains the probability that individual interact across regions, with elements $\Theta_{x,y} = \sum_z \phi_{z,x} \phi_{z,y}$.

These elements compute the probability that individuals in region $x$ travel to region $z$ AND individuals in region $y$ travel to region $z$.

For the 2-region system, the mixing can be expressed in matrix form as: $$\overline{\overline{\Theta}} = \begin{bmatrix} \phi_{1,1}^2 + \phi_{2,1}^2 & \phi_{1,1}\phi_{1,2} + \phi_{2,1}\phi_{2,2} \\ \phi_{2,2}\phi_{2,1} + \phi_{1,2}\phi_{1,1} & \phi_{2,2}^2 + \phi_{1,2}^2 \end{bmatrix}$$

In this formalism, the $\Theta$ which satisfies: $$\dot{I} = \dot{I}_1 + \dot{I}_2 = \beta S I - r_I I$$

requires $$\overline{S}\cdot\overline{\overline{\Theta}}\cdot\overline{I} = S I$$

That leads to the following expression: $$S_1 (\phi_{1,1}^2 + \phi_{2,1}^2) I_1 + S_1 (\phi_{1,1}\phi_{1,2} + \phi_{2,1}\phi_{2,2}) I_2 + S_2 (\phi_{2,2}^2 + \phi_{1,2}^2) I_2 + S_2 (\phi_{2,2}\phi_{2,1} + \phi_{1,2}\phi_{1,1}) I_1$$

If we assume there is no spatial structure and therefore all matrix elements are identical: $\phi_{x,y} = \phi$

Then we can reduce the expression $$\Rightarrow S_1 (2\phi^2) I_1 + S_1 (2\phi^2) I_2 + S_2 (2\phi^2) I_2 + S_2 (2\phi^2) I_1$$

$$\Rightarrow S_1 (2\phi^2) (I_1 + I_2) + S_2 (2\phi^2) (I_1 + I_2)$$

$$\Rightarrow (2\phi^2) (S_1 + S_2) (I_1 + I_2)$$

$$\Rightarrow (2\phi^2) S I$$

Which means that the neutral matrix elements for the 2-region system are: $$\phi = \frac{1}{\sqrt{2}}$$

and in general for the n-region system²: $$\phi = \frac{1}{\sqrt{n}}$$

This “movement matrix” $\overline{\overline{\Phi}}$ can be estimated from information on how individuals travel, such as from survey or mobility data.

Alternatively, we can consider the $\overline{\overline{\Theta}}$ matrix which we can think of as a “infection-flow” matrix. This matrix identifies the flow of infections from individuals residing in regions $x$ to individuals residing in regions $y$. Such matrix may be estimated from more direct measurements of disease flow such as from genetic data or contact tracing.

The $\Phi$ matrix will be normalised such that its row sums are 1.

DiseasyRegions can take either type of adjacency matrix.

The stored adjacency input is long form, but the module also exposes an $infection_flow_matrix field. This field returns a matrix representation of adjacency for the currently selected regions.

Regional differences in infection risk

In DiseasyRegions we implement a regional modifier for the infection risk as a named vector of relative risks for each region in the model.

This modifier can conceptually be thought of as either a modifier on the human behaviour of individual residing in the region (e.g. socio-economic differences in behaviour) or a modifier relating to the physical location itself (e.g. related to infrastructure of the region).

If the modifier $\Gamma_x$ is related to the location, the elements of the infection-flow matrix is modified to $$\Theta_{x,y} = \sum_z \Gamma_z \phi_{z,x} \phi_{z,y}$$

Whereas if the modifier is related to behavioural differences, the elements are modified to: $$\Theta_{x,y} = \sqrt{\Gamma_x\Gamma_y} \sum_z \phi_{z,x} \phi_{z,y}$$

Notice that if the infection-matrix is directly specified, these modifies may already (implicitly) included.

The region-specific modifiers can be set via the set_regional_risks() method, and the interpretation of the modifier can be toggled via the risk_type argument.

Creating a regional module

A DiseasyRegions object can be created and then configured with the selected area, demography, and adjacency data. Each setter validates the candidate input before storing it.

regions <- DiseasyRegions$new()

regions
#> # DiseasyRegions #############################################
#> Area: No area has been specified
#> Total population: No population data loaded
#> Theta matrix: No adjacency data loaded

regions$set_area(area = c("DK", "SE", "NO", "FI"))
regions$set_demography(demography = demography_nordic)
regions$set_adjacency(adjacency = adjacency_meta_nordic)

regions
#> # DiseasyRegions #############################################
#> Area: DK, FI, NO, SE
#> Total population: 27,124,322
#> Theta matrix: Max eigenvalue 9754.49

The $demography field returns demography filtered to the currently selected regions.

regions %.% demography
#> # A tibble: 404 × 3
#>   region   age population
#>   <chr>  <dbl>      <dbl>
#> 1 DK         0      60965
#> 2 DK         1      61357
#> 3 DK         2      61663
#> 4 DK         3      62128
#> 5 DK         4      61148
#> # ℹ 399 more rows

Population totals can be computed directly from the filtered demography.

aggregate(
  population ~ region,
  data = regions %.% demography,
  FUN = sum
)
#>   region population
#> 1     DK    5837230
#> 2     FI    5522887
#> 3     NO    5381527
#> 4     SE   10382678

The active bindings $adjacency and $infection_flow_matrix expose derived views of adjacency input after filtering to the selected regions.

regions %.% adjacency
#> # A tibble: 16 × 3
#>    from  to    adjacency
#>    <chr> <chr>     <dbl>
#>  1 DK    DK       9753. 
#>  2 DK    FI         11.5
#>  3 DK    NO         75.9
#>  4 DK    SE         39.1
#>  5 FI    DK         11.5
#>  6 FI    FI       3512. 
#>  7 FI    NO         13.7
#>  8 FI    SE         24.2
#>  9 NO    DK         75.9
#> 10 NO    FI         13.7
#> 11 NO    NO       6413. 
#> 12 NO    SE         52.1
#> 13 SE    DK         39.1
#> 14 SE    FI         24.2
#> 15 SE    NO         52.1
#> 16 SE    SE       1360.

regions %.% infection_flow_matrix
#>            DK         FI         NO         SE
#> DK 9752.55170   11.48846   75.93018   39.14193
#> FI   11.48846 3511.90746   13.73980   24.17292
#> NO   75.93018   13.73980 6413.46763   52.11436
#> SE   39.14193   24.17292   52.11436 1360.30436

Regional modifiers of infection risk can be set via the $set_regional_risks() method

… tied to location:

regions$set_regional_risks(
  regional_risks = c("DK" = 1.1, "SE" = 0.9, "NO" = 0.8, "FI" = 0.5),
  regional_risks_type = "location"
)

# Our example produces an error since the provided `adjacency` is not movement data
try(regions %.% infection_flow_matrix)
#> Error in base::tryCatch(base::withCallingHandlers({ : 
#>   `regional_risks_type` can only be "location" if `adjacency_type = "movement"`

… or tied to population:

regions$set_regional_risks(
  regional_risks = c("DK" = 1.1, "SE" = 0.9, "NO" = 0.8, "FI" = 0.5),
  regional_risks_type = "behaviour"
)

regions %.% infection_flow_matrix
#>              DK          FI          NO         SE
#> DK 10727.806867    8.520073   71.228818   38.94573
#> FI     8.520073 1755.953729    8.689812   16.21569
#> NO    71.228818    8.689812 5130.774105   44.22050
#> SE    38.945732   16.215690   44.220503 1224.27392

The module configuration can also be visualised via the $plot() method with or plot().

plot(regions)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.4.0; sf_use_s2() is TRUE

Regional module with hierarchical structure (NUTS)

If you have demography and adjacency data for structured regions (such as NUTS), then DiseasyRegionsNuts can handle the nested data.

regions_nuts <- DiseasyRegionsNuts$new()

regions_nuts
#> # DiseasyRegions #############################################
#> Area: No area has been specified
#> Total population: No population data loaded
#> Theta matrix: No adjacency data loaded

regions_nuts$set_area(area = "DK")
regions_nuts$set_demography(demography = demography_nordic_nuts3)
regions_nuts$set_adjacency(adjacency = adjacency_meta_nordic_nuts3)

regions_nuts
#> # DiseasyRegions #############################################
#> Area: DK
#> Total population: 5,992,734
#> Theta matrix: Max eigenvalue 90130217

And again the module configuration can also be visualised with plot().

plot(regions_nuts)

Changing the selected regions

The selected regional scope can be changed with $set_area(). The new regions must be available in the stored adjacency and demography data.

regions$set_area(area = c("DK", "SE"))

regions %.% area
#> [1] "DK" "SE"

The active fields now reflect the updated scope.

regions %.% demography
#> # A tibble: 202 × 3
#>   region   age population
#>   <chr>  <dbl>      <dbl>
#> 1 DK         0      60965
#> 2 DK         1      61357
#> 3 DK         2      61663
#> 4 DK         3      62128
#> 5 DK         4      61148
#> # ℹ 197 more rows

regions %.% adjacency
#> # A tibble: 4 × 3
#>   from  to    adjacency
#>   <chr> <chr>     <dbl>
#> 1 DK    DK       9753. 
#> 2 DK    SE         39.1
#> 3 SE    DK         39.1
#> 4 SE    SE       1360.

regions %.% infection_flow_matrix
#>             DK         SE
#> DK 10727.80687   38.94573
#> SE    38.94573 1224.27392

Validation

The module checks that selected regions are available in the stored demography and adjacency. For example, trying to select a region that is not present in the data will fail.

try(DiseasyRegionsNuts$set_area(regions = "DE"))
#> Error in try(DiseasyRegionsNuts$set_area(regions = "DE")) : 
#>   attempt to apply non-function

The adjacency data must also contain the same region identifiers in the from and to columns.

bad_adjacency <- data.frame(
  from = c("DE", "BE"),
  to = c("DE", "BE"),
  adjacency = c(1, 1)
)

try(regions$set_adjacency(adjacency = bad_adjacency))
#> Error in base::tryCatch(base::withCallingHandlers({ : 
#>   `area` and `adjacency` must contain at least one common region.