Move contact matrix transpose inside model functions

DESCRIPTION

@@ -8,11 +8,11 @@ Authors@R: c(
            comment = c(ORCID = "0000-0002-0362-6717")),
     person("Edwin", "Van Leeuwen", , "edwin.vanleeuwen@ukhsa.gov.uk", role = c("aut", "cph"),
            comment = c(ORCID = "0000-0002-2383-5305")),
+    person("Bankole", "Ahadzie", , "bankole.ahadzie@lshtm.ac.uk", role = c("aut")),
     person("Adam", "Kucharski", , "adam.kucharski@lshtm.ac.uk", role = c("ctb", "rev"),
            comment = c(ORCID = "0000-0001-8814-9421")),
     person("Tim", "Taylor", , "tim.taylor@hiddenelephants.co.uk", role = c("ctb", "rev"),
            comment = c(ORCID = "0000-0002-8587-7113")),
-    person("Banky", "Ahadzie", , "bahadzie@gmail.com", role = "ctb"),
     person("Alexis", "Robert", , "alexis.robert@lshtm.ac.uk", role = "ctb",
            comment = c(ORCID = "0000-0002-4516-2965")),
     person("Hugo", "Gruson", , "hugo.gruson@data.org", role = "rev",

NAMESPACE

@@ -26,16 +26,14 @@ export(new_infections)
 export(outcomes_averted)
 export(population)
 export(vaccination)
-importFrom(data.table,
-  ":=",
-  .BY,
-  .EACHI,
-  .GRP,
-  .I,
-  .N,
-  .NGRP,
-  .SD,
-  data.table
-)
+importFrom(data.table,":=")
+importFrom(data.table,.BY)
+importFrom(data.table,.EACHI)
+importFrom(data.table,.GRP)
+importFrom(data.table,.I)
+importFrom(data.table,.N)
+importFrom(data.table,.NGRP)
+importFrom(data.table,.SD)
+importFrom(data.table,data.table)
 importFrom(odin,odin)
 useDynLib(epidemics, .registration = TRUE)

NEWS.md

@@ -1,5 +1,9 @@
 # epidemics (development version)
 
+## Breaking changes
+
+1. Social contact matrices from _socialmixr_ must no longer be transposed before use. `model_default()` and `model_vacamole()` now transpose the contact matrix internally, so matrices should be passed to `population()` exactly as `socialmixr::contact_matrix()` returns them (#259, @bahadzie). Existing code that calls `t()` on the contact matrix will continue to run without error, but will silently give incorrect results because the matrix is then transposed twice; remove the `t()` call when upgrading. Vignettes, examples, tests and the README have been updated accordingly.
+
 ## Helper functions
 
 1. Added the `combine_populations()` function to combine age-structured populations into a new population object.

R/model_default.R

@@ -272,7 +272,7 @@ model_default <- function(population,
         types = c("vaccination", "null"), null.ok = TRUE
       )
   )
-
+  population[["contact_matrix"]] <- t(population[["contact_matrix"]])
   # make lists if not lists
   population <- list(population) # NOTE: currently not list, but see issue #181
   if (is_lofints) {

R/model_vacamole.R

@@ -366,6 +366,7 @@ model_vacamole <- function(population, # nolint: cyclocomp_linter.
         types = c("vaccination", "null"), null.ok = TRUE
       )
   )
+  population[["contact_matrix"]] <- t(population[["contact_matrix"]])
 
   # make lists if not lists
   population <- list(population)

R/population.R

@@ -3,6 +3,10 @@
 #' @param name Optional string for the population name.
 #' @param contact_matrix A matrix giving the contacts between the demographic
 #' groups in the population. Must be a square matrix.
+#' Contact matrices returned by `socialmixr::contact_matrix()` should be passed
+#' as they are, without transposing them: model functions transpose the matrix
+#' internally to the convention used by the model ODEs, in which
+#' \eqn{M_{ij}} represents contacts to group \eqn{i} from group \eqn{j}.
 #' @param demography_vector A vector of the sizes of each demographic group.
 #' Must have the same length as the dimensions of the contact matrix.
 #' @param initial_conditions Matrix representing the initial proportions of each
@@ -36,6 +40,10 @@ new_population <- function(name = NA_character_,
 #' @param name Optional string for the population name.
 #' @param contact_matrix A matrix giving the contacts between the demographic
 #' groups in the population. Must be a square matrix.
+#' Contact matrices returned by `socialmixr::contact_matrix()` should be passed
+#' as they are, without transposing them: model functions transpose the matrix
+#' internally to the convention used by the model ODEs, in which
+#' \eqn{M_{ij}} represents contacts to group \eqn{i} from group \eqn{j}.
 #' @param demography_vector A vector of the sizes of each demographic group.
 #' Must have the same length as the dimensions of the contact matrix.
 #' @param initial_conditions Matrix representing the initial proportions of each

R/scenario_comparison.R

@@ -64,11 +64,11 @@
 #' )
 #'
 #' # prepare contact matrix
-#' contact_matrix <- t(contact_data$matrix)
+#' contact_matrix <- contact_data[["matrix"]]
 #'
 #' # prepare the demography vector
 #' demography_vector <- contact_data$demography$population
-#' names(demography_vector) <- rownames(contact_matrix)
+#' names(demography_vector) <- colnames(contact_matrix)
 #'
 #' # initial conditions
 #' initial_i <- 1e-6

README.Rmd

@@ -89,11 +89,11 @@ contact_data <- socialmixr::contact_matrix(
 )
 
 # prepare contact matrix
-contact_matrix <- t(contact_data[["matrix"]])
+contact_matrix <- contact_data[["matrix"]]
 
 # prepare the demography vector
 demography_vector <- contact_data[["demography"]][["population"]]
-names(demography_vector) <- rownames(contact_matrix)
+names(demography_vector) <- colnames(contact_matrix)
 ```
 
 Prepare the initial conditions for the population by age group --- here, one in every million individuals is infected at the start of the epidemic (for a total of about 60 infections).
@@ -111,7 +111,7 @@ initial_conditions <- rbind(
   initial_conditions,
   initial_conditions
 )
-rownames(initial_conditions) <- rownames(contact_matrix)
+rownames(initial_conditions) <- colnames(contact_matrix)
 ```
 
 Prepare an object of the class `<population>`, using the function `population()`.
@@ -181,7 +181,7 @@ filter(output, compartment == "infectious") |>
   ) +
   scale_colour_brewer(
     palette = "Dark2",
-    labels = rownames(contact_matrix),
+    labels = colnames(contact_matrix),
     name = "Age group"
   ) +
   scale_y_continuous(

README.md

@@ -115,11 +115,11 @@ contact_data <- socialmixr::contact_matrix(
 )
 
 # prepare contact matrix
-contact_matrix <- t(contact_data[["matrix"]])
+contact_matrix <- contact_data[["matrix"]]
 
 # prepare the demography vector
 demography_vector <- contact_data[["demography"]][["population"]]
-names(demography_vector) <- rownames(contact_matrix)
+names(demography_vector) <- colnames(contact_matrix)
 ```
 
 Prepare the initial conditions for the population by age group — here,
@@ -139,7 +139,7 @@ initial_conditions <- rbind(
   initial_conditions,
   initial_conditions
 )
-rownames(initial_conditions) <- rownames(contact_matrix)
+rownames(initial_conditions) <- colnames(contact_matrix)
 ```
 
 Prepare an object of the class `<population>`, using the function

tests/testthat/_snaps/combine_populations.md

@@ -21,26 +21,26 @@
 
        Contact matrix 
                             Population 1_[0,40) Population 1_[40,65)
-      Population 1_[0,40)            9.69631902           4.78432690
-      Population 1_[40,65)           2.95709036           4.17821782
-      Population 1_[65,Inf)          0.68637984           1.22853242
-      Population 2_[0,40)            0.48481595           0.23921635
-      Population 2_[40,65)           0.14785452           0.20891089
-      Population 2_[65,Inf)          0.03431899           0.06142662
+      Population 1_[0,40)             9.6963190            2.9570904
+      Population 1_[40,65)            4.7843269            4.1782178
+      Population 1_[65,Inf)           2.2403417            2.4784498
+      Population 2_[0,40)             0.4848160            0.1478545
+      Population 2_[40,65)            0.2392163            0.2089109
+      Population 2_[65,Inf)           0.1120171            0.1239225
                             Population 1_[65,Inf) Population 2_[0,40)
-      Population 1_[0,40)              2.24034171          0.48481595
-      Population 1_[40,65)             2.47844976          0.14785452
-      Population 1_[65,Inf)            1.71428571          0.03431899
-      Population 2_[0,40)              0.11201709          9.69631902
-      Population 2_[40,65)             0.12392249          2.95709036
-      Population 2_[65,Inf)            0.08571429          0.68637984
+      Population 1_[0,40)              0.68637984           0.4848160
+      Population 1_[40,65)             1.22853242           0.2392163
+      Population 1_[65,Inf)            1.71428571           0.1120171
+      Population 2_[0,40)              0.03431899           9.6963190
+      Population 2_[40,65)             0.06142662           4.7843269
+      Population 2_[65,Inf)            0.08571429           2.2403417
                             Population 2_[40,65) Population 2_[65,Inf)
-      Population 1_[0,40)             0.23921635            0.11201709
-      Population 1_[40,65)            0.20891089            0.12392249
-      Population 1_[65,Inf)           0.06142662            0.08571429
-      Population 2_[0,40)             4.78432690            2.24034171
-      Population 2_[40,65)            4.17821782            2.47844976
-      Population 2_[65,Inf)           1.22853242            1.71428571
+      Population 1_[0,40)              0.1478545            0.03431899
+      Population 1_[40,65)             0.2089109            0.06142662
+      Population 1_[65,Inf)            0.1239225            0.08571429
+      Population 2_[0,40)              2.9570904            0.68637984
+      Population 2_[40,65)             4.1782178            1.22853242
+      Population 2_[65,Inf)            2.4784498            1.71428571
 
        Initial Conditions 
                               [,1]  [,2]  [,3] [,4] [,5]

tests/testthat/_snaps/population.md

@@ -17,11 +17,11 @@
       [65,Inf): 9,756,203 (20%)
 
        Contact matrix 
-                       age.group
-      contact.age.group    [0,40)  [40,65) [65,Inf)
-               [0,40)   9.6963190 4.784327 2.240342
-               [40,65)  2.9570904 4.178218 2.478450
-               [65,Inf) 0.6863798 1.228532 1.714286
+                contact.age.group
+      age.group    [0,40)  [40,65)  [65,Inf)
+        [0,40)   9.696319 2.957090 0.6863798
+        [40,65)  4.784327 4.178218 1.2285324
+        [65,Inf) 2.240342 2.478450 1.7142857
 
        Initial Conditions 
              [,1]  [,2]  [,3] [,4] [,5]
@@ -34,29 +34,29 @@
     Code
       uk_population
     Message
-      <population>object
+      <population> object
     Output
 
-      Populationname:
+       Population name: 
     Message
-      "UKpopulation"
+      "UK population"
     Output
 
-      Demography
-      Dem.grp.1:31,844,217(50%)
-      Dem.grp.2:19,682,231(30%)
-      Dem.grp.3:9,756,203(20%)
-
-      Contactmatrix
-      age.group
-      contact.age.group[0,40)[40,65)[65,Inf)
-      Dem.grp.1:9.69631904.7843272.240342
-      Dem.grp.2:2.95709044.1782182.478450
-      Dem.grp.3:0.68637981.2285321.714286
-
-      InitialConditions
-      [,1][,2][,3][,4][,5]
-      [1,]0.99995e-055e-0500
-      [2,]0.99995e-055e-0500
-      [3,]0.99995e-055e-0500
+       Demography 
+      Dem. grp. 1: 31,844,217 (50%)
+      Dem. grp. 2: 19,682,231 (30%)
+      Dem. grp. 3: 9,756,203 (20%)
+
+       Contact matrix 
+                    contact.age.group
+      age.group        [0,40)  [40,65)  [65,Inf)
+        Dem. grp. 1: 9.696319 2.957090 0.6863798
+        Dem. grp. 2: 4.784327 4.178218 1.2285324
+        Dem. grp. 3: 2.240342 2.478450 1.7142857
+
+       Initial Conditions 
+             [,1]  [,2]  [,3] [,4] [,5]
+      [1,] 0.9999 5e-05 5e-05    0    0
+      [2,] 0.9999 5e-05 5e-05    0    0
+      [3,] 0.9999 5e-05 5e-05    0    0
 

tests/testthat/test-combine_populations.R

@@ -17,7 +17,7 @@ contact_data <- socialmixr::contact_matrix(
   symmetric = TRUE,
   return_demography = TRUE
 )
-contact_matrix <- t(contact_data$matrix)
+contact_matrix <- contact_data[["matrix"]]
 demography_vector <- contact_data$demography$population
 
 # Create the first population

tests/testthat/test-epidemic_peak.R

@@ -20,11 +20,11 @@ test_data_default_model <- function() {
   )
 
   # prepare contact matrix
-  contact_matrix <- t(contact_data[["matrix"]])
+  contact_matrix <- contact_data[["matrix"]]
 
   # prepare the demography vector
   demography_vector <- contact_data[["demography"]][["population"]]
-  names(demography_vector) <- rownames(contact_matrix)
+  names(demography_vector) <- colnames(contact_matrix)
 
 
   # initial conditions: one in every 1 million is infected
@@ -39,7 +39,7 @@ test_data_default_model <- function() {
     initial_conditions,
     initial_conditions
   )
-  rownames(initial_conditions) <- rownames(contact_matrix)
+  rownames(initial_conditions) <- colnames(contact_matrix)
 
   # prepare the population to model as affected by the epidemic
   uk_population <- population(

tests/testthat/test-interventions.R

@@ -17,7 +17,7 @@ contact_data <- socialmixr::contact_matrix(
   symmetric = TRUE,
   return_demography = TRUE
 )
-contact_matrix <- t(contact_data$matrix)
+contact_matrix <- contact_data[["matrix"]]
 demography_vector <- contact_data$demography$population
 
 # Prepare some initial objects

tests/testthat/test-model_default.R

@@ -17,7 +17,7 @@ contact_data <- socialmixr::contact_matrix(
   symmetric = TRUE,
   return_demography = TRUE
 )
-contact_matrix <- t(contact_data$matrix)
+contact_matrix <- contact_data[["matrix"]]
 demography_vector <- contact_data$demography$population
 
 # make initial conditions - order is important
@@ -76,7 +76,7 @@ test_that("Default model: basic expectations, scalar arguments", {
     )
   )
   expect_identical(
-    unique(data$demography_group), rownames(contact_matrix)
+    unique(data$demography_group), colnames(contact_matrix)
   )
 
   # expect no individuals are vaccinated as vaccination is optional
@@ -650,9 +650,9 @@ test_that("Default model: demography groups with index >= 10 are labelled correc
     age_limits = c(0, 20, 40),
     symmetric = TRUE
   )
-  cm <- t(contact_data_3ag$matrix)
+  cm <- contact_data_3ag$matrix
   dv <- contact_data_3ag$demography$population
-  names(dv) <- rownames(cm)
+  names(dv) <- colnames(cm)
 
   n <- 4L
   n_age <- nrow(cm) # 3

tests/testthat/test-model_vacamole.R

@@ -17,7 +17,7 @@ contact_data <- socialmixr::contact_matrix(
   symmetric = TRUE,
   return_demography = TRUE
 )
-contact_matrix <- t(contact_data$matrix)
+contact_matrix <- contact_data[["matrix"]]
 demography_vector <- contact_data$demography$population
 
 # // 0| 1| 2|3| 4|5| 6|7| 8|9|10
@@ -100,7 +100,7 @@ test_that("Vacamole model: basic expectations, scalar arguments", {
     )
   )
   expect_identical(
-    unique(data$demography_group), rownames(contact_matrix)
+    unique(data$demography_group), colnames(contact_matrix)
   )
 
   # expect no individuals are vaccinated as vaccination is optional