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

Misspecification of time indexation in log-scale time-varying Stan model #208

Closed ntorresd closed 2 months ago

ntorresd commented 2 months ago

Running some tests I realized the time indexation of the log-scale time-varying Stan model is incorrect. See the following example:

image

I'm simulating both surveys using the following FOI:

foi <- data.frame(
  year = seq(1990, 2009, 1),
  foi = c(rep(0, 10), rep(0.1, 5), rep(0,5))
)

image

The difference between the two is that the one on the right is simulated using seroreversion_rate = 0.1 (for the other is just 0). I implemented the time-varying models with and without seroreversion specifying the correct FOI indexation (to reproduce the input of the simulation):

seromodel_no_serorev = fit_seromodel(
  serosurvey,
  model_type = "time",
  is_log_foi = TRUE,
  foi_prior = sf_normal(0, 1e-4),
  foi_index = c(rep(1, 10), rep(2, 5), rep(3, 5)),
  seed = "123"
)

seromodel_serorev <- fit_seromodel(
  serosurvey,
  model_type = "time",
  is_log_foi = TRUE,
  foi_prior = sf_normal(0, 1e-4),
  foi_index = c(rep(1, 10), rep(2, 5), rep(3, 5)),
  is_seroreversion = TRUE,
  seed = "123"
)

Note that the FOI on the right plot is indexed in the reversed order it should be.

ntorresd commented 2 months ago

Finding the error was kind of tricky. Comparing inst/stan/time_log_seroreversion.stan and inst/stan/time_log_no_seroreversion.stan I realized the indexation in the first case was being done like:

for(i in 1:age_max) {
  int time_index = age_max - i + 1;
  foi_expanded[time_index] = foi_vector[foi_index[i]];
}

time_index also appears on other stan files, but it isn't being used:

  for(i in 1:age_max) {
    int time_index = age_max - i + 1;
    foi_expanded[i] = foi_vector[foi_index[i]];
  }

I corrected the indexation in the time-log-foi model with seroreversion and removed the unnecessary lines from the others. Now the example is working as expected.