diff --git a/scripts/create_simple_homog_pop.R b/scripts/create_simple_homog_pop.R index 00097da..796b139 100644 --- a/scripts/create_simple_homog_pop.R +++ b/scripts/create_simple_homog_pop.R @@ -1,9 +1,13 @@ # generate simple, unstructured synthetic population of arbitrary size n <- 1e3 -simple_people <- data.frame("age" = rep(1, times = n), - "homeId" = "10000000000") +simple_people <- data.frame( + "age" = rep(1, times = n), + "homeId" = "10000000000" +) -write.csv(x = simple_people, - file = "./input/simple_people.csv", - row.names = FALSE) \ No newline at end of file +write.csv( + x = simple_people, + file = "./input/simple_people.csv", + row.names = FALSE +) diff --git a/scripts/make_ode_and_compare.R b/scripts/make_ode_and_compare.R index 334cddf..b2d6ec2 100644 --- a/scripts/make_ode_and_compare.R +++ b/scripts/make_ode_and_compare.R @@ -12,23 +12,24 @@ sir_model <- function(time, state, parameters) { S <- state[1] I <- state[2] R <- state[3] - + # Unpack parameters beta <- parameters["beta"] gamma <- parameters["gamma"] - + # Calculate derivatives dS <- -beta * S * I / (S + I + R) dI <- beta * S * I / (S + I + R) - gamma * I dR <- gamma * I - + return(list(c(dS, dI, dR))) } # Get input params input_params <- fromJSON("./input/input.json") pop_size <- nrow(read.csv( - file = input_params$epi_isolation.Parameters$synth_population_file)) + file = input_params$epi_isolation.Parameters$synth_population_file +)) # Obtain via simulation an estimate of the mean duration of infectiousness @@ -36,14 +37,15 @@ R0 <- input_params$epi_isolation.Parameters$r_0 gi <- input_params$epi_isolation.Parameters$generation_interval reps <- 1e6 -ixa_style_duration_I_fx <- function(R0, gi){ +ixa_style_duration_I_fx <- function(R0, gi) { infect_attempts <- rpois(n = 1, lambda = R0) ifelse(test = infect_attempts == 0, - yes = 0, - no = max(rexp(n = infect_attempts, rate = 1/gi))) + yes = 0, + no = max(rexp(n = infect_attempts, rate = 1 / gi)) + ) } -sim_I_dur_ixa <- +sim_I_dur_ixa <- replicate(n = reps, expr = ixa_style_duration_I_fx(R0 = R0, gi = gi)) mean_sim_I_dur_ixa <- mean(sim_I_dur_ixa) @@ -52,25 +54,31 @@ mean_sim_I_dur_ixa <- mean(sim_I_dur_ixa) initial_state <- c(S = pop_size - 1, I = 1, R = 0) # Define parameters -gamma = 1 / mean_sim_I_dur_ixa -beta = input_params$epi_isolation.Parameters$r_0 * gamma +gamma <- 1 / mean_sim_I_dur_ixa +beta <- input_params$epi_isolation.Parameters$r_0 * gamma parameters <- c(beta = beta, gamma = gamma) # Time sequence for simulation (e.g., from day 0 to day 100) time_sequence <- seq(0, input_params$epi_isolation.Parameters$max_time, - by = input_params$epi_isolation.Parameters$report_period) + by = input_params$epi_isolation.Parameters$report_period +) # Run the ODE solver -ode_results <- ode(y = initial_state, - times = time_sequence, - func = sir_model, - parms = parameters) - -ode_results_df <- ode_results |> as.data.frame() |> +ode_results <- ode( + y = initial_state, + times = time_sequence, + func = sir_model, + parms = parameters +) + +ode_results_df <- ode_results |> + as.data.frame() |> rename(t = time, Susceptible = S, Infectious = I, Recovered = R) |> - pivot_longer(cols = c(Susceptible, Infectious, Recovered), - names_to = "InfectiousStatus", - values_to = "count") |> + pivot_longer( + cols = c(Susceptible, Infectious, Recovered), + names_to = "InfectiousStatus", + values_to = "count" + ) |> mutate(model = "ode_sir") # make a model where there are multiple infectious compartments @@ -84,18 +92,18 @@ s3ir_model <- function(time, state, parameters) { I2 <- state[3] I3 <- state[4] R <- state[5] - + # Unpack parameters beta <- parameters["beta"] gamma <- parameters["gamma"] * 3 - + # Calculate derivatives dS <- -beta * S * (I1 + I2 + I3) / (S + I1 + I2 + I3 + R) dI1 <- beta * S * (I1 + I2 + I3) / (S + I1 + I2 + I3 + R) - gamma * I1 dI2 <- gamma * I1 - gamma * I2 dI3 <- gamma * I2 - gamma * I3 - dR <- gamma * I3 # Change in recovered population - + dR <- gamma * I3 # Change in recovered population + return(list(c(dS, dI1, dI2, dI3, dR))) } @@ -103,44 +111,54 @@ s3ir_model <- function(time, state, parameters) { initial_state_s3ir <- c(S = pop_size - 1, I1 = 1, I2 = 0, I3 = 0, R = 0) # Run the ODE solver -ode_s3iR_results <- ode(y = initial_state_s3ir, - times = time_sequence, - func = s3ir_model, - parms = parameters) - -ode_s3ir_results_df <- ode_s3iR_results |> as.data.frame() |> - mutate(I = I1 + I2 + I3) |> select(!c(I1, I2, I3)) |> +ode_s3iR_results <- ode( + y = initial_state_s3ir, + times = time_sequence, + func = s3ir_model, + parms = parameters +) + +ode_s3ir_results_df <- ode_s3iR_results |> + as.data.frame() |> + mutate(I = I1 + I2 + I3) |> + select(!c(I1, I2, I3)) |> rename(t = time, Susceptible = S, Infectious = I, Recovered = R) |> - pivot_longer(cols = c(Susceptible, Infectious, Recovered), - names_to = "InfectiousStatus", - values_to = "count") |> + pivot_longer( + cols = c(Susceptible, Infectious, Recovered), + names_to = "InfectiousStatus", + values_to = "count" + ) |> mutate(model = "ode_s3ir") # Load simulation results from Ixa model ixa_results_df <- read.csv( - file = "./output/person_property_count_multi_rep.csv") + file = "./output/person_property_count_multi_rep.csv" +) # combine ODE and stochastic results and make plots -ode_results_df <- ode_results_df |> +ode_results_df <- ode_results_df |> mutate(ixa_rep = max(ixa_results_df$ixa_rep) + 1) -ode_s3ir_results_df <- ode_s3ir_results_df |> +ode_s3ir_results_df <- ode_s3ir_results_df |> mutate(ixa_rep = max(ixa_results_df$ixa_rep) + 2) results_df <- rbind(ixa_results_df, ode_results_df, ode_s3ir_results_df) -results_df |> - ggplot(aes(x = t, y = count, group = paste(InfectiousStatus,ixa_rep))) + +results_df |> + ggplot(aes(x = t, y = count, group = paste(InfectiousStatus, ixa_rep))) + geom_line( - aes(color = InfectiousStatus, linewidth = model, linetype = model)) + + aes(color = InfectiousStatus, linewidth = model, linetype = model) + ) + xlab("Day") + ylab("Number of people") + scale_linetype_manual(values = c(1, 2, 3)) + scale_linewidth_manual(values = c(0.1, 1, 1)) + theme_minimal() -results_df |> filter(InfectiousStatus == "Infectious") |> - ggplot(aes(x = t, y = count, group = paste(InfectiousStatus,ixa_rep))) + +results_df |> + filter(InfectiousStatus == "Infectious") |> + ggplot(aes(x = t, y = count, group = paste(InfectiousStatus, ixa_rep))) + geom_line( - aes(color = InfectiousStatus, linewidth = model, linetype = model)) + + aes(color = InfectiousStatus, linewidth = model, linetype = model) + ) + xlab("Day") + ylab("Number of people") + scale_linetype_manual(values = c(1, 2, 3)) + @@ -149,16 +167,18 @@ results_df |> filter(InfectiousStatus == "Infectious") |> # simulate durations of infectiousness and plot cdfs for different assumptions -s3ir_duration_I_fx <- function(n_I_comp, mean_ixa){ - sum(rexp(n = n_I_comp, rate = 1/mean_ixa * n_I_comp)) +s3ir_duration_I_fx <- function(n_I_comp, mean_ixa) { + sum(rexp(n = n_I_comp, rate = 1 / mean_ixa * n_I_comp)) } -sim_I_dur_ode_sir <- - rexp(n = reps, rate = 1/mean_sim_I_dur_ixa) +sim_I_dur_ode_sir <- + rexp(n = reps, rate = 1 / mean_sim_I_dur_ixa) -sim_I_dur_ode_s3ir <- - replicate(n = reps, expr = s3ir_duration_I_fx(n_I_comp = 3, - mean_ixa = mean_sim_I_dur_ixa)) +sim_I_dur_ode_s3ir <- + replicate(n = reps, expr = s3ir_duration_I_fx( + n_I_comp = 3, + mean_ixa = mean_sim_I_dur_ixa + )) plot(ecdf(sim_I_dur_ixa), main = "CDF of duration of infectiousness") plot(ecdf(sim_I_dur_ode_sir), col = "orange", add = TRUE) diff --git a/scripts/run_ixa_repeatedly.R b/scripts/run_ixa_repeatedly.R index 40aa6ad..95affb2 100644 --- a/scripts/run_ixa_repeatedly.R +++ b/scripts/run_ixa_repeatedly.R @@ -3,7 +3,7 @@ # so run the script directly from a terminal # function that changes the seed, runs the model, and gets the output -run_ixa_rep_fx <- function(ixa_rep){ +run_ixa_rep_fx <- function(ixa_rep) { input_params$epi_isolation.Parameters$seed <- ixa_rep jsonlite::toJSON(x = input_params, pretty = TRUE, auto_unbox = TRUE) |> writeLines("./input/input.json") @@ -11,7 +11,8 @@ run_ixa_rep_fx <- function(ixa_rep){ infectious_report <- readr::read_csv(file.path( "output", "person_property_count.csv" - )) |> dplyr::mutate(model = "ixa", ixa_rep = ixa_rep) |> + )) |> + dplyr::mutate(model = "ixa", ixa_rep = ixa_rep) |> dplyr::select(!c(Age, CensusTract)) return(infectious_report) } @@ -29,6 +30,8 @@ output <- lapply(X = seq_len(ixa_reps), FUN = run_ixa_rep_fx) output_df <- do.call(rbind, output) # export results -write.csv(x = output_df, - file = "./output/person_property_count_multi_rep.csv", - row.names = FALSE) \ No newline at end of file +write.csv( + x = output_df, + file = "./output/person_property_count_multi_rep.csv", + row.names = FALSE +)