Open avallecam opened 6 months ago
# load packages ----------------------------------------------------------- library(epidemics) library(socialmixr) #> #> Attaching package: 'socialmixr' #> The following object is masked from 'package:utils': #> #> cite library(tidyverse) # parameters -------------------------------------------------------------- preinfectious_period <- 3.0 infectious_period <- 7.0 basic_reproduction <- 1.46 infectiousness_rate <- 1.0/preinfectious_period recovery_rate <- 1.0/infectious_period transmission_rate <- basic_reproduction/infectious_period # paste from tutorial ----------------------------------------------------- polymod <- socialmixr::polymod contact_data <- socialmixr::contact_matrix( polymod, countries = "Italy", age.limits = c(0, 20, 40), symmetric = TRUE ) #> Using POLYMOD social contact data. To cite this in a publication, use the 'get_citation()' function #> Removing participants without age information. To change this behaviour, set the 'missing.participant.age' option #> Removing participants that have contacts without age information. To change this behaviour, set the 'missing.contact.age' option # prepare contact matrix contact_matrix <- t(contact_data$matrix) # prepare the demography vector demography_vector <- contact_data$demography$population names(demography_vector) <- rownames(contact_matrix) # initial conditions: one in every 1 million is infected initial_i <- 1e-6 initial_conditions <- c( S = 1 - initial_i, E = 0, I = initial_i, R = 0, V = 0 ) # build for all age groups initial_conditions <- matrix( rep(initial_conditions, dim(contact_matrix)[1]), ncol = 5, byrow = TRUE ) rownames(initial_conditions) <- rownames(contact_matrix) # prepare the population to model as affected by the epidemic uk_population <- population( name = "UK", contact_matrix = contact_matrix, demography_vector = demography_vector, initial_conditions = initial_conditions ) uk_population #> <population> object #> #> Population name: #> "UK" #> #> Demography #> [0,20): 11,204,261 (20%) #> [20,40): 16,305,622 (30%) #> 40+: 31,298,598 (50%) #> #> Contact matrix #> #> contact.age.group [0,20) [20,40) 40+ #> [0,20) 14.306502 2.255442 2.897234 #> [20,40) 3.282356 9.126904 4.182137 #> 40+ 8.093293 8.027601 8.900621 # intervention: close school ------------------------------------------------------------ rownames(contact_matrix) #> [1] "[0,20)" "[20,40)" "40+" close_schools <- intervention( name = "School closure", type = "contacts", time_begin = 200, time_end = 200 + 100, reduction = matrix(c(0.5, 0.01, 0.01)) # group specific effect ) close_schools #> <contacts_intervention> object #> #> Intervention name: #> "School closure" #> #> Begins at: #> [1] 200 #> #> Ends at: #> [1] 300 #> #> Reduction: #> Interv. 1 #> Demo. grp. 1 0.50 #> Demo. grp. 2 0.01 #> Demo. grp. 3 0.01 # intervention: mask mandate ---------------------------------------------- mask_mandate <- intervention( name = "Mask mandate", type = "rate", time_begin = 200, time_end = 200 + 200, reduction = 0.163 # rate impact all infectious (reduce infectiousness) ) mask_mandate #> <rate_intervention> object #> #> Intervention name: #> "Mask mandate" #> #> Begins at: #> [,1] #> [1,] 200 #> #> Ends at: #> [,1] #> [1,] 400 #> #> Reduction: #> Interv. 1 #> 0.163 # intervention: vaccination ----------------------------------------------- vaccinate <- vaccination( name = "vaccinate all", time_begin = matrix(200, nrow(contact_matrix)), time_end = matrix(200 + 150, nrow(contact_matrix)), nu = matrix(c(0.01, 0.0, 0.01)) ) # run baseline ------------------------------------------------------------ baseline <- model_default( # population population = uk_population, # rates transmission_rate = transmission_rate, infectiousness_rate = infectiousness_rate, recovery_rate = recovery_rate, # time time_end = 600,increment = 1 ) # run intervention -------------------------------------------------------- intervention_school <- model_default( # population population = uk_population, # rates transmission_rate = transmission_rate, infectiousness_rate = infectiousness_rate, recovery_rate = recovery_rate, # intervention intervention = list(contacts = close_schools), # time time_end = 600,increment = 1 ) intervention_mask <- model_default( # population population = uk_population, # rates transmission_rate = transmission_rate, infectiousness_rate = infectiousness_rate, recovery_rate = recovery_rate, # intervention intervention = list(transmission_rate = mask_mandate), # time time_end = 600,increment = 1 ) intervention_vaccine <- model_default( # population population = uk_population, # rates transmission_rate = transmission_rate, infectiousness_rate = infectiousness_rate, recovery_rate = recovery_rate, # intervention vaccination = vaccinate, # time time_end = 600,increment = 1 ) # intervention: combine --------------------------------------------------- intervention_mask_vaccine <- model_default( # population population = uk_population, # rates transmission_rate = transmission_rate, infectiousness_rate = infectiousness_rate, recovery_rate = recovery_rate, # intervention intervention = list(contacts = close_schools), # intervention = list(transmission_rate = mask_mandate), vaccination = vaccinate, # time time_end = 600,increment = 1 ) # # get new infections ------------------------------------------------------ # # baseline_infections <- new_infections(baseline,by_group = FALSE) # intervention_infections <- new_infections(intervention,by_group = FALSE) # # # # get finalsize ----------------------------------------------------------- # # epidemic_size(baseline) # epidemic_size(intervention) # visualize --------------------------------------------------------------- baseline$intervention_type <- "baseline" intervention_school$intervention_type <- "close_schools" intervention_mask$intervention_type <- "mask_mandate" intervention_vaccine$intervention_type <- "vaccinate" intervention_mask_vaccine$intervention_type <- "mask + vaccine" output <- bind_rows( baseline, intervention_school, intervention_mask, intervention_vaccine, intervention_mask_vaccine ) output %>% filter(compartment == "infectious") %>% ggplot( aes(x = time, y = value, linetype = intervention_type, colour = intervention_type ) ) + stat_summary( fun = "sum", geom = "line", linewidth = 1 ) + scale_y_continuous(labels = scales::comma) + geom_vline( xintercept = c( close_schools$time_begin, close_schools$time_end ), linetype = 2 ) + geom_vline( xintercept = c( mask_mandate$time_begin, mask_mandate$time_end ), linetype = 3 ) + geom_vline( xintercept = c( vaccinate$time_begin, vaccinate$time_end ), linetype = 4 ) + labs( x = "Simulation time (days)", colour = "Age group", y = "Individuals" ) + theme_bw()
Created on 2024-04-26 with reprex v2.1.0
Created on 2024-04-26 with reprex v2.1.0