Open IzzyRou opened 2 years ago
Code to reproduce bug (using develop branch)
` library(tidyverse) library(devtools)
setwd("~/GitHub/SIMPLEGEN") devtools::load_all() check_SIMPLEGEN_loaded()
set.seed(1) n_demes <- 1 s <- simplegen_project() %>% define_epi_model_parameters( H = 1e3, M = 1e4, seed_infections = 10, mig_mat = diag(1, 1), prob_infection = 0.3, prob_acute = seq(1, 0.0000000000, l = 30), prob_AC = 1.0, duration_acute = dgeom(1:250, 1 / 100), duration_chronic = dgeom(1:250, 1 / 400), detectability_microscopy_acute = 0.9, detectability_microscopy_chronic = 0.45, treatment_seeking_mean = 0, infectivity_acute = 0.9, infectivity_chronic = 0.45 )
daily_dataframe <- rbind(data.frame(name = "acute_prev", deme = 1, state = "A", measure = "prevalence", diagnostic = c("microscopy", "PCR"), age_min = 2, age_max = 10), data.frame(name = "chronic_prev", deme = 1, state = "C", measure = "prevalence", diagnostic = c("microscopy", "PCR"), age_min = 2, age_max = 10), data.frame(name = "acute_inc", deme = 1, state = "A", measure = "incidence", diagnostic = "microscopy", age_min = 0, age_max = 5), data.frame(name = "EIR", deme = 1, state = "H", measure = "EIR", diagnostic = NA, age_min = 0, age_max = 110), make.row.names = FALSE)
sweep_dataframe <- rbind(data.frame(name = "prev_acute", time = 365, deme = 1, measure = "prevalence", state = "A", diagnostic = "PCR", age_min = seq(0, 100, 5), age_max = seq(0, 100, 5) + 4), data.frame(name = "inc_acute", time = 365, deme = 1, measure = "incidence", state = "A", diagnostic = "microscopy", age_min = seq(0, 100, 5), age_max = seq(0, 100, 5) + 4), data.frame(name = "EIR",time = 365, deme = 1, state = "H", measure = "EIR", diagnostic = NA, age_min = 0, age_max = 110), make.row.names = FALSE)
df_survey<-data.frame( t_start = 365, t_end= 365+29, reporting_interval = 30, deme = 1 , measure = c("prevalence", "prevalence", "incidence","prevalence", "prevalence", "incidence"), diagnostic = c("PCR","microscopy","microscopy","PCR","microscopy","microscopy") , sampling = c(NA,NA,"PCD",NA,NA,"PCD"), age_min = c(0,0,0,2,2,0), age_max = c(110,110,110,10,10,5), sample_size = 1000
)
s <- define_epi_sampling_parameters(project = s, surveys = df_survey, daily = daily_dataframe, sweep = sweep_dataframe)
s <- sim_epi( s, max_time = 2 * 365, save_transmission_record = FALSE, transmission_record_location = "ignore/trans_record.txt", overwrite_transmission_record = TRUE )
s$epi_output$sweeps %>% dplyr::filter(measure == "EIR")
s$epi_output$surveys %>% dplyr::filter(measure == "incidence")
s$epi_output$daily %>% dplyr::filter((measure == "incidence" | measure == "EIR") & time == unique(s$epi_output$sweeps$time))
`
Code to reproduce bug (using develop branch)
` library(tidyverse) library(devtools)
Set working directory to location of SIMPLEGEN and load
setwd("~/GitHub/SIMPLEGEN") devtools::load_all() check_SIMPLEGEN_loaded()
Define model parameters
set.seed(1) n_demes <- 1 s <- simplegen_project() %>% define_epi_model_parameters( H = 1e3, M = 1e4, seed_infections = 10, mig_mat = diag(1, 1), prob_infection = 0.3, prob_acute = seq(1, 0.0000000000, l = 30), prob_AC = 1.0, duration_acute = dgeom(1:250, 1 / 100), duration_chronic = dgeom(1:250, 1 / 400), detectability_microscopy_acute = 0.9, detectability_microscopy_chronic = 0.45, treatment_seeking_mean = 0, infectivity_acute = 0.9, infectivity_chronic = 0.45 )
Define sampling strategy
daily_dataframe <- rbind(data.frame(name = "acute_prev", deme = 1, state = "A", measure = "prevalence", diagnostic = c("microscopy", "PCR"), age_min = 2, age_max = 10), data.frame(name = "chronic_prev", deme = 1, state = "C", measure = "prevalence", diagnostic = c("microscopy", "PCR"), age_min = 2, age_max = 10), data.frame(name = "acute_inc", deme = 1, state = "A", measure = "incidence", diagnostic = "microscopy", age_min = 0, age_max = 5), data.frame(name = "EIR", deme = 1, state = "H", measure = "EIR", diagnostic = NA, age_min = 0, age_max = 110), make.row.names = FALSE)
sweep_dataframe <- rbind(data.frame(name = "prev_acute", time = 365, deme = 1, measure = "prevalence", state = "A", diagnostic = "PCR", age_min = seq(0, 100, 5), age_max = seq(0, 100, 5) + 4), data.frame(name = "inc_acute", time = 365, deme = 1, measure = "incidence", state = "A", diagnostic = "microscopy", age_min = seq(0, 100, 5), age_max = seq(0, 100, 5) + 4), data.frame(name = "EIR",time = 365, deme = 1, state = "H", measure = "EIR", diagnostic = NA, age_min = 0, age_max = 110), make.row.names = FALSE)
df_survey<-data.frame( t_start = 365, t_end= 365+29, reporting_interval = 30, deme = 1 , measure = c("prevalence", "prevalence", "incidence","prevalence", "prevalence", "incidence"), diagnostic = c("PCR","microscopy","microscopy","PCR","microscopy","microscopy") , sampling = c(NA,NA,"PCD",NA,NA,"PCD"), age_min = c(0,0,0,2,2,0), age_max = c(110,110,110,10,10,5), sample_size = 1000
)
s <- define_epi_sampling_parameters(project = s, surveys = df_survey, daily = daily_dataframe, sweep = sweep_dataframe)
Simulate from transmission model
s <- sim_epi( s, max_time = 2 * 365, save_transmission_record = FALSE, transmission_record_location = "ignore/trans_record.txt", overwrite_transmission_record = TRUE )
Note that EIR is listed as 0 when sampled at time = 365 days
s$epi_output$sweeps %>% dplyr::filter(measure == "EIR")
Note that incidence in 0-5 and 0-110 years over the 30 day period beginning on time = 365 is listed as zero
s$epi_output$surveys %>% dplyr::filter(measure == "incidence")
Yet at time = 365 the daily output lists EIR as 110.6 and incidence in 0 to 5 as 0.33
s$epi_output$daily %>% dplyr::filter((measure == "incidence" | measure == "EIR") & time == unique(s$epi_output$sweeps$time))
`