model_default() is merging rows when 4 populations or more

[avallecam]

Bugs may include behaviour that does not impact package functionality, but which seems as though it might not be intended by the package developers.

Please place an "x" in all the boxes that apply

• I have the most recent version of epidemics and R
• I have found a bug
• I have a reproducible example
• I want to request a new feature

---

Please include a brief description of the problem with a code example:

library(epidemics)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(ggplot2)

# --- Contact data (UK, as in vignette) ---
polymod <- socialmixr::polymod

contact_data1 <- socialmixr::contact_matrix(
  polymod,
  countries = "United Kingdom",
  age.limits = c(0, 20, 40),
  symmetric = TRUE
)
#> Removing participants that have contacts without age information. To change this behaviour, set the 'missing.contact.age' option
contact_matrix1 <- t(contact_data1$matrix)
demography_vector1 <- contact_data1$demography$population
names(demography_vector1) <- rownames(contact_matrix1)

# --- Build n populations (first n-1 with no infected) ---
n <- 4

all_population <- list()
for (i in seq_len(n - 1)) {
  all_population[i] <- list(population(
    contact_matrix = contact_matrix1,
    demography_vector = demography_vector1 / 100,
    initial_conditions = rbind(
      c(S = 1, E = 0, I = 0, R = 0, V = 0),
      c(S = 1, E = 0, I = 0, R = 0, V = 0),
      c(S = 1, E = 0, I = 0, R = 0, V = 0)
    ),
    name = paste0("Region ", i)
  ))
}

# --- Add last population with a few infectious individuals ---
all_population[n] <- list(population(
  contact_matrix = contact_matrix1,
  demography_vector = demography_vector1 / 100,
  initial_conditions = rbind(
    c(S = 1 - 1e-6, E = 0, I = 1e-6, R = 0, V = 0),
    c(S = 1 - 1e-6, E = 0, I = 1e-6, R = 0, V = 0),
    c(S = 1 - 1e-6, E = 0, I = 1e-6, R = 0, V = 0)
  ),
  name = paste0("Region ", n)
))

# --- Distance matrix: consecutive populations 150km apart ---
distance_matrix_n <-
  matrix(150 * abs(rep(seq_len(n), n) - rep(seq_len(n), each = n)), nrow = n)

# --- Combine and run ---
combined_n_populations <-
  combine_populations(
    populations = all_population,
    connectivity_matrix = distance_matrix_n,
    method = "gravity", name = "combine_gravity_n"
  )

output_gravity_n <- model_default(
  population = combined_n_populations,
  time_end = 700, increment = 1.0
)

# --- Test ---
output_gravity_n %>% dplyr::count(demography_group)
#>    demography_group     n
#>              <char> <int>
#> 1:     Region 1_40+  3505
#> 2:  Region 1_[0,20) 14020
#> 3: Region 1_[20,40)  3505
#> 4:     Region 2_40+  3505
#> 5:  Region 2_[0,20)  3505
#> 6: Region 2_[20,40)  3505
#> 7:     Region 3_40+  3505
#> 8:  Region 3_[0,20)  3505
#> 9: Region 3_[20,40)  3505

^{Created on 2026-06-02 with reprex v2.1.1}

---

Expected solution: the output from output_gravity_n %>% dplyr::count(demography_group) should have the same number of rows for each group. Currently region 4 is merged in region 1.

Where this was found?

Directly diagnosed from output in vignette, last figure: https://epiverse-trace.github.io/epidemics/dev/articles/modelling_populations.html#combining-n-populations-using-a-gravity-model

[avallecam]

Root cause

The bug originates at line 447:

temp_demography = substring(temp, 3L, 3L), # e.g. S[1] -> 1

substring(temp, 3L, 3L) always extracts exactly one character at position 3. This is the index digit inside the odin variable name (e.g. S[1], E[2]). It works correctly for single-digit indices (1–9) but silently truncates multi-digit ones:

odin variable	extracted index	correct index
S[1]	"1" ✓	"1"
S[9]	"9" ✓	"9"
S[10]	"1" ✗	"10"
S[11]	"1" ✗	"11"
S[12]	"1" ✗	"12"

Trigger condition

The total number of demographic groups = n populations × age groups per population. With the standard 3-age-group UK/Germany contact matrices:

• n = 2 → 6 groups (indices 1–6) — works fine
• n = 3 → 9 groups (indices 1–9) — works fine
• n = 4 → 12 groups (indices 1–12) — breaks: groups 10, 11, 12 all map to index "1", silently merging their rows with Region 1

This is why the issue is not visible in the two-population vignette examples but appears in the n-populations example.

Proposed fix

temp_demography = substring(temp, 3L, nchar(as.character(temp)) - 1L), # e.g. S[10] -> 10

This strips the leading X[ (positions 1–2) and the trailing ] (last character), correctly extracting the full index regardless of digit count.

[avallecam]

Follow-up: there is a second related bug in the same code path. After looking at the output more closely, values for groups 10+ appeared as 0Region 4_[0,20), 1Region 4_[20,40), etc.

The cause: the sorting prefix in age_group_mappings uses seq_len(n_age), which becomes two digits for n_age ≥ 10. The strip at line 464 was hardcoded to substring(demography_group, 2L), removing only 1 character, leaving the second digit behind.

Fix: use a zero-padded prefix of width nchar(n_age) via formatC, and strip n_demography_prefix + 1L characters instead of a hardcoded 2L. This is included as a separate commit on the fix branch.