epiverse-trace / simulist

An R package for simulating line lists
https://epiverse-trace.github.io/simulist/
Other
4 stars 0 forks source link

Issue when using infectious period which can sample negative values #141

Closed joshwlambert closed 2 weeks ago

joshwlambert commented 2 weeks ago

There an issue when simulating a line list with sim_linelist() when the distribution of the infectious period can generate negative values, reported by @CarmenTamayo.

Negative infectious periods do not make epidemiological sense (in the model the time an infector makes contact with their sampled contacts is uniformally distributed within the infectious period).

The sim_linelist() function runs and creates a line list <data.frame> as well as several warnings generated by stats::runif() from within .sim_network_bp() when the infect_period produces negative values. The columns of the line list get populated with NaNs. Here is an example:

library(simulist)
library(epiparameter)

# load data required to simulate line list
contact_distribution <- epiparameter::epidist(
  disease = "Ebola",
  epi_dist = "contact distribution",
  prob_distribution = "pois",
  prob_distribution_params = c(mean = 2)
)
#> Citation cannot be created as author, year, journal or title is missing

infect_period <- epiparameter::epidist(
  disease = "Ebola",
  epi_dist = "infectious period",
  prob_distribution = "norm",
  prob_distribution_params = c(mean = 9.4, sd = 5.5)
)
#> Citation cannot be created as author, year, journal or title is missing

onset_to_hosp <- epidist(
  disease = "ebola",
  epi_dist = "onset to hospitalisation",
  prob_distribution = "gamma",
  prob_distribution_params = c(shape = 1.6, scale = 3.0)
)
#> Citation cannot be created as author, year, journal or title is missing

onset_to_death <- epidist_db(
  disease = "ebola",
  epi_dist = "onset to death",
  single_epidist = TRUE
)
#> Using WHO Ebola Response Team, Agua-Agum J, Ariyarajah A, Aylward B, Blake I,
#> Brennan R, Cori A, Donnelly C, Dorigatti I, Dye C, Eckmanns T, Ferguson
#> N, Formenty P, Fraser C, Garcia E, Garske T, Hinsley W, Holmes D,
#> Hugonnet S, Iyengar S, Jombart T, Krishnan R, Meijers S, Mills H,
#> Mohamed Y, Nedjati-Gilani G, Newton E, Nouvellet P, Pelletier L,
#> Perkins D, Riley S, Sagrado M, Schnitzler J, Schumacher D, Shah A, Van
#> Kerkhove M, Varsaneux O, Kannangarage N (2015). "West African Ebola
#> Epidemic after One Year — Slowing but Not Yet under Control." _The New
#> England Journal of Medicine_. doi:10.1056/NEJMc1414992
#> <https://doi.org/10.1056/NEJMc1414992>.. 
#> To retrieve the citation use the 'get_citation' function

onset_to_recovery <- epidist(
  disease = "ebola",
  epi_dist = "onset to recovery",
  prob_distribution = "gamma",
  prob_distribution_params = c(shape = 19.8, scale = 0.9)
)
#> Citation cannot be created as author, year, journal or title is missing

set.seed(1234)
linelist <- sim_linelist(
  contact_distribution = contact_distribution,
  infect_period = infect_period,
  onset_to_hosp = onset_to_hosp,
  onset_to_death = onset_to_death,
  onset_to_recovery = onset_to_recovery,
  prob_infect = 0.6,
  outbreak_size = c(500, 2000),
  non_hosp_death_risk = 0.4
)
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning in stats::runif(n = length(vec_idx), min = 0, max =
#> contact_infect_period): NAs produced
#> Warning: Number of cases exceeds maximum outbreak size. 
#> Returning data early with 2107 cases and 3514 total contacts (including cases).
head(linelist$date_onset, n = 100)
#>   [1] "2023-01-01" "2023-01-09" "2023-01-11" "2023-01-15" "2023-01-18"
#>   [6] "2023-01-13" "2023-01-13" "2023-01-18" "2023-01-18" "2023-01-22"
#>  [11] "2023-01-20" "2023-01-22" "2023-01-28" "2023-01-21" "2023-01-25"
#>  [16] "2023-01-23" "2023-01-22" "2023-02-04" "2023-02-06" "2023-02-05"
#>  [21] "2023-01-25" "2023-01-28" "2023-01-31" "2023-02-15" "2023-02-08"
#>  [26] "2023-01-29" "2023-02-01" "2023-02-14" "2023-02-17" "2023-02-17"
#>  [31] "2023-02-17" "2023-02-16" "2023-02-14" "2023-01-31" "2023-02-16"
#>  [36] "2023-02-19" "2023-03-03" "2023-03-02" "2023-02-18" "2023-02-03"
#>  [41] "2023-02-21" "2023-03-03" "2023-03-05" "2023-02-07" "2023-02-23"
#>  [46] "2023-02-21" "2023-02-24" "2023-02-07" "2023-02-11" "2023-02-12"
#>  [51] "2023-02-10" "2023-02-21" "2023-02-28" "2023-02-09" "2023-02-17"
#>  [56] "2023-02-15" "2023-02-17" "2023-02-19" "2023-02-14" "2023-02-17"
#>  [61] "2023-02-22" "2023-03-01" "2023-03-02" "2023-02-15" "2023-02-09"
#>  [66] "2023-02-24" "2023-02-17" "2023-02-28" "2023-02-18" "2023-02-20"
#>  [71] "2023-02-20" "2023-02-21" "2023-02-20" "2023-02-24" "2023-02-17"
#>  [76] "2023-02-20" "2023-03-10" "2023-02-09" "2023-02-10" "2023-03-13"
#>  [81] "2023-03-10" "2023-02-27" "2023-02-18" "2023-03-01" "2023-02-28"
#>  [86] "2023-02-23" "2023-02-23" "NaN"        "2023-02-21" "2023-02-21"
#>  [91] "2023-03-09" "2023-03-04" "2023-03-04" "2023-03-01" "2023-02-25"
#>  [96] "NaN"        "NaN"        "2023-02-22" "2023-02-23" "2023-03-16"

Created on 2024-06-17 with reprex v2.1.0

joshwlambert commented 2 weeks ago

Reproducing the reprex above using the updates made in PR #142.

library(simulist)
library(epiparameter)

# load data required to simulate line list
contact_distribution <- epiparameter::epidist(
  disease = "Ebola",
  epi_dist = "contact distribution",
  prob_distribution = "pois",
  prob_distribution_params = c(mean = 2)
)
#> Citation cannot be created as author, year, journal or title is missing

infect_period <- epiparameter::epidist(
  disease = "Ebola",
  epi_dist = "infectious period",
  prob_distribution = "norm",
  prob_distribution_params = c(mean = 9.4, sd = 5.5)
)
#> Citation cannot be created as author, year, journal or title is missing

onset_to_hosp <- epidist(
  disease = "ebola",
  epi_dist = "onset to hospitalisation",
  prob_distribution = "gamma",
  prob_distribution_params = c(shape = 1.6, scale = 3.0)
)
#> Citation cannot be created as author, year, journal or title is missing

onset_to_death <- epidist_db(
  disease = "ebola",
  epi_dist = "onset to death",
  single_epidist = TRUE
)
#> Using WHO Ebola Response Team, Agua-Agum J, Ariyarajah A, Aylward B, Blake I,
#> Brennan R, Cori A, Donnelly C, Dorigatti I, Dye C, Eckmanns T, Ferguson
#> N, Formenty P, Fraser C, Garcia E, Garske T, Hinsley W, Holmes D,
#> Hugonnet S, Iyengar S, Jombart T, Krishnan R, Meijers S, Mills H,
#> Mohamed Y, Nedjati-Gilani G, Newton E, Nouvellet P, Pelletier L,
#> Perkins D, Riley S, Sagrado M, Schnitzler J, Schumacher D, Shah A, Van
#> Kerkhove M, Varsaneux O, Kannangarage N (2015). "West African Ebola
#> Epidemic after One Year — Slowing but Not Yet under Control." _The New
#> England Journal of Medicine_. doi:10.1056/NEJMc1414992
#> <https://doi.org/10.1056/NEJMc1414992>.. 
#> To retrieve the citation use the 'get_citation' function

onset_to_recovery <- epidist(
  disease = "ebola",
  epi_dist = "onset to recovery",
  prob_distribution = "gamma",
  prob_distribution_params = c(shape = 19.8, scale = 0.9)
)
#> Citation cannot be created as author, year, journal or title is missing

set.seed(1234)
linelist <- sim_linelist(
  contact_distribution = contact_distribution,
  infect_period = infect_period,
  onset_to_hosp = onset_to_hosp,
  onset_to_death = onset_to_death,
  onset_to_recovery = onset_to_recovery,
  prob_infect = 0.6,
  outbreak_size = c(500, 2000),
  non_hosp_death_risk = 0.4
)
#> Error: Negative infectious period sampled. The infectious period must be strictly positive.

Created on 2024-06-17 with reprex v2.1.0