epiverse-trace / serofoi

Estimates the Force-of-Infection of a given pathogen from population based sero-prevalence studies
https://epiverse-trace.github.io/serofoi/
Other
17 stars 4 forks source link

Bug in matrix calculation `get_prev_expanded` #145

Closed ntorresd closed 7 months ago

ntorresd commented 7 months ago

Currently the package is not giving the right results for the seroprevalence due to a bug that went unnoticed (my bad...) in #135 . The reason is a change in how the expanded prevalence is computed (see commit ca885da). The original code:

age_class <- seq_len(NCOL(foi_expanded))
ly <- NCOL(foi_expanded)
exposure <- matrix(0, nrow = length(age_class), ncol = ly)
for (k in seq_along(age_class)) {
  exposure[k, (ly - age_class[k] + 1):ly] <- 1
}
exposure_expanded <- exposure

Outputs a lower triangular matrix with respect to the antidiagonal:

exposure_expanded
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17]
 [1,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0
 [2,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0
 [3,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0
 [4,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0
 [5,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0
 [6,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0
 [7,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0
 [8,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0
 [9,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0
[10,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0
[11,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0
[12,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0
...      
 [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77]
 [1,]     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
 [2,]     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
 [3,]     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
 [4,]     0     0     0     0     0     0     0     0     0     0     0     0     0     0     1
 [5,]     0     0     0     0     0     0     0     0     0     0     0     0     0     1     1
 [6,]     0     0     0     0     0     0     0     0     0     0     0     0     1     1     1
 [7,]     0     0     0     0     0     0     0     0     0     0     0     1     1     1     1
 [8,]     0     0     0     0     0     0     0     0     0     0     1     1     1     1     1
 [9,]     0     0     0     0     0     0     0     0     0     1     1     1     1     1     1
[10,]     0     0     0     0     0     0     0     0     1     1     1     1     1     1     1
[11,]     0     0     0     0     0     0     0     1     1     1     1     1     1     1     1
[12,]     0     0     0     0     0     0     1     1     1     1     1     1     1     1     1
      [,78] [,79] [,80]
 [1,]     0     0     1
 [2,]     0     1     1
 [3,]     1     1     1
 [4,]     1     1     1
 [5,]     1     1     1
 [6,]     1     1     1
 [7,]     1     1     1
 [8,]     1     1     1
 [9,]     1     1     1
[10,]     1     1     1
[11,]     1     1     1
[12,]     1     1     1

Meanwhile this simplification avoiding the loop:

ly <- NCOL(foi_expanded)
exposure_expanded <- matrix(0, nrow = ly, ncol = ly)
exposure_expanded[lower.tri(exposure_expanded, diag = TRUE)] <- 1

yields the regular lower triangular matrix.

exposure_expanded
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
 [1,]    1    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0
 [2,]    1    1    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0
 [3,]    1    1    1    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0
 [4,]    1    1    1    1    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0
 [5,]    1    1    1    1    1    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0
 [6,]    1    1    1    1    1    1    0    0    0     0     0     0     0     0     0     0     0     0     0     0
 [7,]    1    1    1    1    1    1    1    0    0     0     0     0     0     0     0     0     0     0     0     0
 [8,]    1    1    1    1    1    1    1    1    0     0     0     0     0     0     0     0     0     0     0     0
 [9,]    1    1    1    1    1    1    1    1    1     0     0     0     0     0     0     0     0     0     0     0
[10,]    1    1    1    1    1    1    1    1    1     1     0     0     0     0     0     0     0     0     0     0
[11,]    1    1    1    1    1    1    1    1    1     1     1     0     0     0     0     0     0     0     0     0
[12,]    1    1    1    1    1    1    1    1    1     1     1     1     0     0     0     0     0     0     0     0
      [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39]
 [1,]     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
 [2,]     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
 [3,]     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
 [4,]     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
 [5,]     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
 [6,]     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
 [7,]     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
 [8,]     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
 [9,]     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
[10,]     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
[11,]     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
[12,]     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0