Evaluation of management alternatives across systems
Analysis referenced in Cook et al. In Prep
Source:vignettes/Management_Alternatives_Systems.Rmd
Management_Alternatives_Systems.Rmd
Introduction:
In this vignette, we run outbreak simulations to estimate the effects of various management alternatives on the dynamics of a SARS-CoV-2 outbreak in wild and captive white-tailed deer. We considered captive:wild systems that are separated by a fence. These management alternatives may be implemented in both wild and captive populations, or in one scenario. We focus on 11 alternatives, ranging across agriculture, public health, and wildlife sectors. These alternatives are detailed in Cook et al. In Prep.. These alternatives (and parameter changes) include:
Local reductions in wild deer populations through agency culling or hunter harvest (Reduce wild deer densities by 10% and 25%)
Prohibit deer feeding and baiting to reduce congregations of wild deer (Remove inflation factor of proximity rates due to bait and supplemental feeding)
Pause/reduce research permits for work on wild deer or other susceptible wildlife/animal species (Reduce human-wild deer interaction)
Encourage limiting human-deer interactions including participating in conduct that attracts wildlife (Proportionally reduce human-deer proximity rates reflecting rural:suburban rates)
Encourage (and train on the proper) use of enhanced personal protective equipment during interactions with captive and wild deer (Integrate effects of masks in aerosol transmission from humans to deer)
Animal vaccination requirement (All captive deer start simulations as recovered; All captive deer start simulations as recovered and are boosted when immunity wanes)
Upgrade fences to double fences that prevent transmission between captive and wild deer (Transmission between captive and wild deer is set to zero)
Increase ventilation in indoor, intensive facilities (Increase air exchange rate)
We detail the alternatives and the systems where they can be implemented:
ActionTable <- data.frame(Alternative = c("Reduce wild density by 10%", "Reduce wild density by 25%", "Eliminate baiting", "Pause Research", "Restrict human interactions with deer","PPE, Wild" , "PPE, Captive", "PPE, Both", "Vaccinate captive deer", "Vaccinate and boost captive deer", "Double-fence captive facilities", "Increase air exchange rate in intensive facilities"), Sector = c("Wildlife", "Wildlife", "Wildlife", "Wildlife", "Public Health","Public Health","Public Health","Public Health", "Agriculture", "Agriculture","Agriculture", "Agriculture"), Ranch_Rural_System = c(TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,FALSE), Ranch_Suburban_System = c(TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,FALSE), Intensive_Facility_Rural_System = c(TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE), Intensive_Facility_Suburban_System = c(TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE))
ActionTable %>%
kbl() %>%
kable_paper("hover", full_width = F)
Alternative | Sector | Ranch_Rural_System | Ranch_Suburban_System | Intensive_Facility_Rural_System | Intensive_Facility_Suburban_System |
---|---|---|---|---|---|
Reduce wild density by 10% | Wildlife | TRUE | TRUE | TRUE | TRUE |
Reduce wild density by 25% | Wildlife | TRUE | TRUE | TRUE | TRUE |
Eliminate baiting | Wildlife | TRUE | TRUE | TRUE | TRUE |
Pause Research | Wildlife | TRUE | TRUE | TRUE | TRUE |
Restrict human interactions with deer | Public Health | TRUE | TRUE | TRUE | TRUE |
PPE, Wild | Public Health | TRUE | TRUE | TRUE | TRUE |
PPE, Captive | Public Health | TRUE | TRUE | TRUE | TRUE |
PPE, Both | Public Health | TRUE | TRUE | TRUE | TRUE |
Vaccinate captive deer | Agriculture | TRUE | TRUE | TRUE | TRUE |
Vaccinate and boost captive deer | Agriculture | TRUE | TRUE | TRUE | TRUE |
Double-fence captive facilities | Agriculture | TRUE | TRUE | TRUE | TRUE |
Increase air exchange rate in intensive facilities | Agriculture | FALSE | FALSE | TRUE | TRUE |
To evaluate the effects of these interventions, we follow the same analytical framework as in other vignettes. Because of we are looking at multiple systems, our work flow will differ a little in this vignette. Steps include:
Set number of simulations to run
Sample parameters for simulations
Establish systems and test alternatives
Combine and summarize results across alternatives and systems
Impact of cross-sector collaborations (multiple alternatives implemented in concert).
1. Set number of simulations to run
Like in other vignettes, we set a seed and sample size. We use a small sample size for ease of vignette loading. We also define the length of the simulation here.
2. Sample parameters for simulations
Common Parameters
We bring in elicitation data for various parameters included in these simulations.
We can then define a number of parameters that will be used in every system in which we want to test alternatives.
Scenario parameters
Here, we define parameter values used for the baseline projection for each system:
Captive deer in outdoor ranch facilities
#Ranch Parameters
#Deer-to-deer proximity rate in ranch context, mimicing wild proximity rates with the influence of baiting (events per day).
omega_cc_ranch <- calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild)*(get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17)
#Human-to-deer proximity rate in ranch context, mimicing suburban proximitty rates (events per day).
omega_hc_ranch <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Suburban (per 120 days)") /120
#Estimate duration of human-deer proximity event
t_contact_deer_human_ranch <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Duration, Captive (minutes)")
#Estimate probability of infection from deer-deer proximity given that duration of proximity event.
sigma_aero_deer_deer_ranch <- calc_sigma_aero(C_nu = C_nu_deer,
t_contact = t_contact_deer_deer_null / 60,
r = r_deer, nsamples = nsamples)
#Calculate infection probability given the duration of a human-deer proximity event.
sigma_aero_deer_human_ranch <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human,
t_contact = t_contact_deer_human_ranch / 60,
r = r_deer, nsamples = nsamples)
Captive deer in intensive facilities
#Deer-to-deer proximity rate in intensive context (events per day).
omega_cc_intensive <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Deer Proximity Rate, Captive (per day)")
#Human-to-deer proximity rate in intensive context (events per day).
omega_hc_intensive <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Captive (per 120 days)") /120
#Estimate duration of human-deer proximity event
t_contact_deer_human_intensive <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Duration, Captive (minutes)")
#Estimate probability of infection from deer-deer proximity given that duration of proximity event.
sigma_aero_deer_deer_intensive <- calc_sigma_aero(C_nu = C_nu_deer,
t_contact = t_contact_deer_deer_null / 60,
r = r_deer, nsamples = nsamples, AER = rep(1, nsamples))
#Calculate infection probability given the duration of a human-deer proximity event.
sigma_aero_deer_human_intensive <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human,
t_contact = t_contact_deer_human_intensive / 60,
r = r_deer, nsamples = nsamples, AER = rep(1, nsamples))
Wild deer in a suburban setting
#Deer-to-deer proximity rate (events per day).
omega_ww_suburban <- calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild)*(get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17)
#Human-to-deer proximity rate (events per day).
omega_hw_suburban <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Suburban (per 120 days)") /120
#Estimate duration of human-deer proximity event
t_contact_deer_human_suburban <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Duration, Suburban (minutes)")
#Estimate probability of infection from deer-deer proximity given that duration of proximity event.
sigma_aero_deer_deer_suburban <- calc_sigma_aero(C_nu = C_nu_deer,
t_contact = t_contact_deer_deer_null / 60,
r = r_deer, nsamples = nsamples, AER = rep(4, nsamples))
#Calculate infection probability given the duration of a human-deer proximity event.
sigma_aero_deer_human_suburban <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human,
t_contact = t_contact_deer_human_suburban / 60,
r = r_deer, nsamples = nsamples, AER = rep(4, nsamples))
Wild deer in a rural setting
#Deer-to-deer proximity rate (events per day).
omega_ww_rural <- calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild)*(get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17)
#Human-to-deer proximity rate (events per day).
omega_hw_rural <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Rural (per 120 days)") /120
#Estimate duration of human-deer proximity event
t_contact_deer_human_rural <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Duration, Rural (minutes)")
#Estimate probability of infection from deer-deer proximity given that duration of proximity event.
sigma_aero_deer_deer_rural <- calc_sigma_aero(C_nu = C_nu_deer,
t_contact = t_contact_deer_deer_null / 60,
r = r_deer, nsamples = nsamples, AER = rep(4, nsamples))
#Calculate infection probability given the duration of a human-deer proximity event.
sigma_aero_deer_human_rural <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human,
t_contact = t_contact_deer_human_rural / 60,
r = r_deer, nsamples = nsamples, AER = rep(4, nsamples))
3.Establish Systems and Test Alternatives
Below we define a captive-wild system, simulate SARS-CoV-2 outbreaks with no interventions, and then apply the 12 candidate interventions by altering a relevant parameter. All 13 simulations (1 no action + 12 interventions) are combined into one data frame at the end of each code chunk.
Combined system of outdoor ranch facilities and wild deer in a suburban setting
ranch_suburban.inits.fall <- initial_compartments(draws = nsamples)
ranch_suburban.inits.steady <- initial_compartments(draws = nsamples, steady = TRUE)
ranch_suburban.params <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_suburban, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_suburban, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_ranch, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_ranch_suburban_null <- run(iter = nsamples, initial_compartments = ranch_suburban.inits.fall, initial_compartments_steady = ranch_suburban.inits.steady, params = ranch_suburban.params, times = times, name = "Baseline")
#Alternatives
#1. 10% density reduction
omega_ww_alt1_density.red.10 <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.9)*(get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17)
ranch_suburban.params_alt1 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_alt1_density.red.10, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_suburban, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_ranch, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_ranch_suburban_alt1 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits.fall, initial_compartments_steady = ranch_suburban.inits.steady, params = ranch_suburban.params_alt1, times = times, name = "Reduce wild density by 10%")
#2. 25% density reduction
omega_ww_alt2_density.red.25 <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.75)*(get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17)
ranch_suburban.params_alt2 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_alt2_density.red.25, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_suburban, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_ranch, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_ranch_suburban_alt2 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits.fall, initial_compartments_steady = ranch_suburban.inits.steady, params = ranch_suburban.params_alt2, times = times, name = "Reduce wild density by 25%")
#3. Eliminate baiting
omega_ww_alt3_no.baiting <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.75)
ranch_suburban.params_alt3 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_alt3_no.baiting, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_suburban, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_ranch, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_ranch_suburban_alt3 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits.fall, initial_compartments_steady = ranch_suburban.inits.steady, params = ranch_suburban.params_alt3, times = times, name = "Eliminate baiting")
#4. Pause research (assuming random 1-50% of human-wild deer interaction)
omega_hw_alt4_pause.research <- omega_hw_suburban*runif(nsamples, min = 0.5, max = 0.99) #reduce human deer proximity rate in the wild up to 50%
ranch_suburban.params_alt4 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_suburban, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_alt4_pause.research, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_ranch, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_ranch_suburban_alt4 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits.fall, initial_compartments_steady = ranch_suburban.inits.steady, params = ranch_suburban.params_alt4, times = times, name = "Pause Research")
#5. Require PPE for humans in wild settings
sigma_aero_deer_human_wild_alt5_PPE <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_suburban/60, r = r_deer, V_d = 0.009*runif(nsamples, min = (1-0.994), max = (1-0.509)), nsamples = nsamples)
ranch_suburban.params_alt5 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_suburban, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_suburban, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, sigma_aero_deer_human_capt = sigma_aero_deer_human_ranch, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_ranch_suburban_alt5 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits.fall, initial_compartments_steady = ranch_suburban.inits.steady, params = ranch_suburban.params_alt5, times = times, name = "PPE, Wild")
#6. Require PPE for humans in captive settings
sigma_aero_deer_human_capt_alt6_PPE <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_ranch/60, r = r_deer, V_d = 0.009*runif(nsamples, min = (1-0.994), max = (1-0.509)), nsamples = nsamples)
ranch_suburban.params_alt6 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_suburban, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_suburban, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_ranch_suburban_alt6 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits.fall, initial_compartments_steady = ranch_suburban.inits.steady, params = ranch_suburban.params_alt6, times = times, name = "PPE, Captive")
#7. Require PPE for humans in wild and captive settings (use derived parameters from 5 and 6)
ranch_suburban.params_alt7 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_suburban, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_suburban, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_ranch_suburban_alt7 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits.fall, initial_compartments_steady = ranch_suburban.inits.steady, params = ranch_suburban.params_alt7, times = times, name = "PPE, Both")
#8. Vaccinate captive deer (start with r_captive = 1)
ranch_suburban.inits_alt8 <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1)
ranch_suburban.inits_alt8.steady <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1, steady = TRUE)
proj_ranch_suburban_alt8 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits_alt8, initial_compartments_steady = ranch_suburban.inits_alt8.steady,params = ranch_suburban.params, times = times, name = "Vaccinate captive deer")
#9. Vaccinate and boost captive deer (start with r_captive = 1)
ranch_suburban.inits_alt9 <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1)
ranch_suburban.inits_alt9.steady <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1, steady = TRUE)
boost_alt9 <- rep(0.9,nsamples)
ranch_suburban.params_alt9 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_suburban, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_suburban, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_ranch, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = boost_alt9)
proj_ranch_suburban_alt9 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits_alt9, initial_compartments_steady = ranch_suburban.inits_alt9.steady, params = ranch_suburban.params_alt9, times = times, name = "Vaccinate and boost captive deer")
#10. Double fence captive facilities
omega_cw_alt10 <- rep(0, nsamples)
ranch_suburban.params_alt10 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_suburban, omega_cw = omega_cw_alt10, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_suburban, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_ranch, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0,nsamples))
proj_ranch_suburban_alt10 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits.fall, initial_compartments_steady = ranch_suburban.inits.steady, params = ranch_suburban.params_alt10, times = times, name = "Double-fence captive facilities")
#11. Restrict human interaction (using proportional reduction from suburban to rural)
human_restrictions <- (get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Rural (per 120 days)") / 120)/(get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Suburban (per 120 days)") / 120)
ranch_suburban.params_alt11 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_suburban, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_suburban*human_restrictions, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_ranch, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0,nsamples))
proj_ranch_suburban_alt11 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits.fall, initial_compartments_steady = ranch_suburban.inits.steady, params = ranch_suburban.params_alt11, times = times, name = "Restrict human interactions")
#12. Improve indoor facility air exchange (does not apply to outdoor ranch facilities)
proj_ranch_suburban <- rbind(proj_ranch_suburban_null, proj_ranch_suburban_alt1, proj_ranch_suburban_alt2, proj_ranch_suburban_alt3, proj_ranch_suburban_alt4, proj_ranch_suburban_alt5, proj_ranch_suburban_alt6, proj_ranch_suburban_alt7, proj_ranch_suburban_alt8, proj_ranch_suburban_alt9, proj_ranch_suburban_alt10, proj_ranch_suburban_alt11)
Calcualte prevalence, cumulative infections, and persistence for interventions in the outdoor ranch : wild, suburban system
persist.threshold <- 0.001
proj_ranch_suburban %>%
mutate(steady_sir = map(steady_state, "y"),
steady_sir = map(steady_sir, as_tibble_row)) %>%
pull(steady_sir) %>%
list_rbind(names_to = "run_id") %>%
mutate(., persistWild = I_wild > persist.threshold, persistCaptive = I_captive > persist.threshold) %>%
select(., run_id, persistWild, persistCaptive) %>%
arrange(., run_id) -> Persist_ranch_suburban_df
proj_ranch_suburban %>%
mutate(ode_df = map(ode_proj, as.data.frame)) %>%
pull(ode_df) %>%
list_rbind(names_to = "run_id") %>%
mutate(., Complex = rep("Outdoor ranch and suburban deer",nrow(.)), Alternative = c(rep(unique(proj_ranch_suburban$Context)[1],length(times)*nsamples),rep(unique(proj_ranch_suburban$Context)[2],length(times)*nsamples), rep(unique(proj_ranch_suburban$Context)[3],length(times)*nsamples),rep(unique(proj_ranch_suburban$Context)[4],length(times)*nsamples), rep(unique(proj_ranch_suburban$Context)[5],length(times)*nsamples), rep(unique(proj_ranch_suburban$Context)[6],length(times)*nsamples), rep(unique(proj_ranch_suburban$Context)[7],length(times)*nsamples),rep(unique(proj_ranch_suburban$Context)[8],length(times)*nsamples),rep(unique(proj_ranch_suburban$Context)[9],length(times)*nsamples),rep(unique(proj_ranch_suburban$Context)[10],length(times)*nsamples),rep(unique(proj_ranch_suburban$Context)[11],length(times)*nsamples),rep(unique(proj_ranch_suburban$Context)[12],length(times)*nsamples))) %>%
group_by(Complex,Alternative, run_id) %>%
summarize(meanWild = mean(I_wild), meanCaptive = mean(I_captive),
cumulativeWild = last(I_wild_cumulative), cumulativeCaptive = last(I_captive_cumulative)) %>%
merge(.,Persist_ranch_suburban_df, by = "run_id")-> alternatives_for_ranch_suburban_complex
We can see the result of this processing, is a dataframe listing outbreak metrics for all iterations of each management alternative (including the default baseline). We will repeat this for the other systems, and then use summary functions to compare the effect of each alternative to the baseline, for each system.
head(alternatives_for_ranch_suburban_complex)
#> run_id Complex Alternative meanWild meanCaptive
#> 1 1 Outdoor ranch and suburban deer Baseline 6.676796e-02 0.0688039960
#> 2 2 Outdoor ranch and suburban deer Baseline 9.735172e-05 0.0000827675
#> 3 3 Outdoor ranch and suburban deer Baseline 5.937537e-02 0.0598487381
#> 4 4 Outdoor ranch and suburban deer Baseline 3.225770e-05 0.0002175466
#> 5 5 Outdoor ranch and suburban deer Baseline 4.790617e-02 0.0479719483
#> 6 6 Outdoor ranch and suburban deer Baseline 2.432780e-02 0.0266367681
#> cumulativeWild cumulativeCaptive persistWild persistCaptive
#> 1 1.397087341 1.435637721 TRUE TRUE
#> 2 0.002060018 0.001751416 FALSE FALSE
#> 3 1.236692239 1.246979083 TRUE TRUE
#> 4 0.000686242 0.004626696 FALSE FALSE
#> 5 0.967220064 0.968339379 TRUE TRUE
#> 6 0.537436553 0.573908352 TRUE TRUE
Combined system of outdoor ranch facilities and wild deer in a rural setting
ranch_rural.inits.fall <- initial_compartments(draws = nsamples)
ranch_rural.inits.steady <- initial_compartments(draws = nsamples, steady = TRUE)
ranch_rural.params <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_rural, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_rural, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_ranch, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_ranch_rural_null <- run(iter = nsamples, initial_compartments = ranch_rural.inits.fall, initial_compartments_steady = ranch_rural.inits.steady, params = ranch_rural.params, times = times, name = "Baseline")
#Alternatives
#1. 10% density reduction
omega_ww_alt1_density.red.10 <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.9)*(get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17)
ranch_rural.params_alt1 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_alt1_density.red.10, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_rural, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_ranch, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_ranch_rural_alt1 <- run(iter = nsamples, initial_compartments = ranch_rural.inits.fall, initial_compartments_steady = ranch_rural.inits.steady, params = ranch_rural.params_alt1, times = times, name = "Reduce wild density by 10%")
#2. 25% density reduction
omega_ww_alt2_density.red.25 <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.75)*(get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17)
ranch_rural.params_alt2 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_alt2_density.red.25, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_rural, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_ranch, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_ranch_rural_alt2 <- run(iter = nsamples, initial_compartments = ranch_rural.inits.fall, initial_compartments_steady = ranch_rural.inits.steady, params = ranch_rural.params_alt2, times = times, name = "Reduce wild density by 25%")
#3. Eliminate baiting
omega_ww_alt3_no.baiting <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.75)
ranch_rural.params_alt3 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_alt3_no.baiting, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_rural, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_ranch, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_ranch_rural_alt3 <- run(iter = nsamples, initial_compartments = ranch_rural.inits.fall, initial_compartments_steady = ranch_rural.inits.steady, params = ranch_rural.params_alt3, times = times, name = "Eliminate baiting")
#4. Pause research (assuming random 1-50% of human-wild deer interaction)
omega_hw_alt4_pause.research <- omega_hw_rural*runif(nsamples, min = 0.5, max = 0.99) #reduce human deer proximity rate in the wild up to 50%
ranch_rural.params_alt4 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_rural, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_alt4_pause.research, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_ranch, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_ranch_rural_alt4 <- run(iter = nsamples, initial_compartments = ranch_rural.inits.fall, initial_compartments_steady = ranch_rural.inits.steady, params = ranch_rural.params_alt4, times = times, name = "Pause Research")
#5. Require PPE for humans in wild settings
sigma_aero_deer_human_wild_alt5_PPE <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_rural/60, r = r_deer, V_d = 0.009*runif(nsamples, min = (1-0.994), max = (1-0.509)), nsamples = nsamples)
ranch_rural.params_alt5 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_rural, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_rural, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, sigma_aero_deer_human_capt = sigma_aero_deer_human_ranch, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_ranch_rural_alt5 <- run(iter = nsamples, initial_compartments = ranch_rural.inits.fall, initial_compartments_steady = ranch_rural.inits.steady, params = ranch_rural.params_alt5, times = times, name = "PPE, Wild")
#6. Require PPE for humans in captive settings
sigma_aero_deer_human_capt_alt6_PPE <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_ranch/60, r = r_deer, V_d = 0.009*runif(nsamples, min = (1-0.994), max = (1-0.509)), nsamples = nsamples)
ranch_rural.params_alt6 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_rural, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_rural, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_ranch_rural_alt6 <- run(iter = nsamples, initial_compartments = ranch_rural.inits.fall, initial_compartments_steady = ranch_rural.inits.steady, params = ranch_rural.params_alt6, times = times, name = "PPE, Captive")
#7. Require PPE for humans in wild and captive settings (use derived parameters from 5 and 6)
ranch_rural.params_alt7 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_rural, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_rural, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_ranch_rural_alt7 <- run(iter = nsamples, initial_compartments = ranch_rural.inits.fall, initial_compartments_steady = ranch_rural.inits.steady, params = ranch_rural.params_alt7, times = times, name = "PPE, Both")
#8. Vaccinate captive deer (start with r_captive = 1)
ranch_rural.inits_alt8 <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1)
ranch_rural.inits_alt8_steady <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1, steady = T)
proj_ranch_rural_alt8 <- run(iter = nsamples, initial_compartments = ranch_rural.inits_alt8, initial_compartments_steady = ranch_rural.inits_alt8_steady, params = ranch_rural.params, times = times, name = "Vaccinate captive deer")
#9. Vaccinate and boost captive deer (start with r_captive = 1)
ranch_rural.inits_alt9 <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1)
ranch_rural.inits_alt9_steady <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1, steady = T)
boost_alt9 <- rep(0.9,nsamples)
ranch_rural.params_alt9 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_rural, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_rural, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_ranch, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = boost_alt9)
proj_ranch_rural_alt9 <- run(iter = nsamples, initial_compartments = ranch_rural.inits_alt9, initial_compartments_steady = ranch_rural.inits_alt9_steady, params = ranch_rural.params_alt9, times = times, name = "Vaccinate and boost captive deer")
#10. Double fence captive facilities
omega_cw_alt10 <- rep(0, nsamples)
ranch_rural.params_alt10 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_rural, omega_cw = omega_cw_alt10, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_rural, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_ranch, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0,nsamples))
proj_ranch_rural_alt10 <- run(iter = nsamples, initial_compartments = ranch_rural.inits.fall, initial_compartments_steady = ranch_rural.inits.steady, params = ranch_rural.params_alt10, times = times, name = "Double-fence captive facilities")
#11. Restrict human interaction (using proportional reduction from suburban to rural)
human_restrictions <- (get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Rural (per 120 days)") / 120)/(get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Suburban (per 120 days)") / 120)
ranch_rural.params_alt11 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_rural, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_rural*human_restrictions, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_ranch, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0,nsamples))
proj_ranch_rural_alt11 <- run(iter = nsamples, initial_compartments = ranch_rural.inits.fall, initial_compartments_steady = ranch_rural.inits.steady, params = ranch_rural.params_alt11, times = times, name = "Restrict human interactions")
#12. Improve indoor facility air exchange (does not apply to outdoor ranch facilities)
proj_ranch_rural <- rbind(proj_ranch_rural_null, proj_ranch_rural_alt1, proj_ranch_rural_alt2, proj_ranch_rural_alt3, proj_ranch_rural_alt4, proj_ranch_rural_alt5, proj_ranch_rural_alt6, proj_ranch_rural_alt7, proj_ranch_rural_alt8, proj_ranch_rural_alt9, proj_ranch_rural_alt10, proj_ranch_rural_alt11)
Calculate prevalence, cumulative infections, and persistence for interventions in the outdoor ranch : wild, rural system
persist.threshold <- 0.001
proj_ranch_rural %>%
mutate(steady_sir = map(steady_state, "y"),
steady_sir = map(steady_sir, as_tibble_row)) %>%
pull(steady_sir) %>%
list_rbind(names_to = "run_id") %>%
mutate(., persistWild = I_wild > persist.threshold, persistCaptive = I_captive > persist.threshold) %>%
select(., run_id, persistWild, persistCaptive) %>%
arrange(., run_id) -> Persist_ranch_rural_df
proj_ranch_rural %>%
mutate(ode_df = map(ode_proj, as.data.frame)) %>%
pull(ode_df) %>%
list_rbind(names_to = "run_id") %>%
mutate(., Complex = rep("Outdoor ranch and rural deer",nrow(.)), Alternative = c(rep(unique(proj_ranch_rural$Context)[1],length(times)*nsamples),rep(unique(proj_ranch_rural$Context)[2],length(times)*nsamples), rep(unique(proj_ranch_rural$Context)[3],length(times)*nsamples),rep(unique(proj_ranch_rural$Context)[4],length(times)*nsamples), rep(unique(proj_ranch_rural$Context)[5],length(times)*nsamples), rep(unique(proj_ranch_rural$Context)[6],length(times)*nsamples), rep(unique(proj_ranch_rural$Context)[7],length(times)*nsamples),rep(unique(proj_ranch_rural$Context)[8],length(times)*nsamples),rep(unique(proj_ranch_rural$Context)[9],length(times)*nsamples),rep(unique(proj_ranch_rural$Context)[10],length(times)*nsamples),rep(unique(proj_ranch_rural$Context)[11],length(times)*nsamples),rep(unique(proj_ranch_rural$Context)[12],length(times)*nsamples))) %>%
group_by(Complex,Alternative, run_id) %>%
summarize(meanWild = mean(I_wild), meanCaptive = mean(I_captive),
cumulativeWild = last(I_wild_cumulative), cumulativeCaptive = last(I_captive_cumulative)) %>%
merge(.,Persist_ranch_rural_df, by = "run_id")-> alternatives_for_ranch_rural_complex
Combined system of intensive captive facilities and wild deer in a suburban setting
intensive_suburban.inits.fall <- initial_compartments(draws = nsamples)
intensive_suburban.inits.steady <- initial_compartments(draws = nsamples, steady = TRUE)
intensive_suburban.params <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_suburban, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_suburban, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_intensive_suburban_null <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params = intensive_suburban.params, times = times, name = "Baseline")
#Alternatives
#1. 10% density reduction
omega_ww_alt1_density.red.10 <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.9)*(get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17)
intensive_suburban.params_alt1 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_alt1_density.red.10, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_suburban, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_intensive_suburban_alt1 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params = intensive_suburban.params_alt1, times = times, name = "Reduce wild density by 10%")
#2. 25% density reduction
omega_ww_alt2_density.red.25 <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.75)*(get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17)
intensive_suburban.params_alt2 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_alt2_density.red.25, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_suburban, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_intensive_suburban_alt2 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params = intensive_suburban.params_alt2, times = times, name = "Reduce wild density by 25%")
#3. Eliminate baiting
omega_ww_alt3_no.baiting <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.75)
intensive_suburban.params_alt3 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_alt3_no.baiting, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_suburban, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_intensive_suburban_alt3 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params = intensive_suburban.params_alt3, times = times, name = "Eliminate baiting")
#4. Pause research (assuming random 1-50% of human-wild deer interaction)
omega_hw_alt4_pause.research <- omega_hw_suburban*runif(nsamples, min = 0.5, max = 0.99) #reduce human deer proximity rate in the wild up to 50%
intensive_suburban.params_alt4 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_suburban, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_alt4_pause.research, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_intensive_suburban_alt4 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params = intensive_suburban.params_alt4, times = times, name = "Pause Research")
#5. Require PPE for humans in wild settings
sigma_aero_deer_human_wild_alt5_PPE <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_suburban/60, r = r_deer, V_d = 0.009*runif(nsamples, min = (1-0.994), max = (1-0.509)), nsamples = nsamples)
intensive_suburban.params_alt5 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_suburban, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_suburban, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_intensive_suburban_alt5 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params = intensive_suburban.params_alt5, times = times, name = "PPE, Wild")
#6. Require PPE for humans in captive settings
sigma_aero_deer_human_capt_alt6_PPE <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_intensive/60, r = r_deer, V_d = 0.009*runif(nsamples, min = (1-0.994), max = (1-0.509)), nsamples = nsamples)
intensive_suburban.params_alt6 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_suburban, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_suburban, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_intensive_suburban_alt6 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params = intensive_suburban.params_alt6, times = times, name = "PPE, Captive")
#7. Require PPE for humans in wild and captive settings (use derived parameters from 5 and 6)
intensive_suburban.params_alt7 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_suburban, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_suburban, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_intensive_suburban_alt7 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params = intensive_suburban.params_alt7, times = times, name = "PPE, Both")
#8. Vaccinate captive deer (start with r_captive = 1)
intensive_suburban.inits_alt8.fall <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1)
intensive_suburban.inits_alt8.steady <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1, steady = TRUE)
proj_intensive_suburban_alt8 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits_alt8.fall, initial_compartments_steady = intensive_suburban.inits_alt8.steady, params = intensive_suburban.params, times = times, name = "Vaccinate captive deer")
#9. Vaccinate and boost captive deer (start with r_captive = 1)
intensive_suburban.inits_alt9.fall <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1)
intensive_suburban.inits_alt9.steady <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1, steady = TRUE)
boost_alt9 <- rep(0.9,nsamples)
intensive_suburban.params_alt9 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_suburban, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_suburban, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = boost_alt9)
proj_intensive_suburban_alt9 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits_alt9.fall, initial_compartments_steady = intensive_suburban.inits_alt9.steady, params = intensive_suburban.params_alt9, times = times, name = "Vaccinate and boost captive deer")
#10. Double fence captive facilities
omega_cw_alt10 <- rep(0, nsamples)
intensive_suburban.params_alt10 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_suburban, omega_cw = omega_cw_alt10, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_suburban, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0,nsamples))
proj_intensive_suburban_alt10 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params = intensive_suburban.params_alt10, times = times, name = "Double-fence captive facilities")
#11. Restrict human interaction (using proportional reduction from suburban to rural)
human_restrictions <- (get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Rural (per 120 days)") / 120)/(get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Suburban (per 120 days)") / 120)
intensive_suburban.params_alt11 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_suburban, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_suburban*human_restrictions, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0,nsamples))
proj_intensive_suburban_alt11 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params = intensive_suburban.params_alt11, times = times, name = "Restrict human interactions")
#12. Improve indoor facility air exchange
sigma_aero_deer_human_intensive_alt12 <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_intensive/60, r = r_deer, nsamples = nsamples)
sigma_aero_deer_deer_intensive_alt12 <- calc_sigma_aero(C_nu = C_nu_deer, t_contact = t_contact_deer_deer_null/60, r = r_deer, nsamples = nsamples)
intensive_suburban.params_alt12 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_suburban, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_suburban, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive_alt12, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive_alt12, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0,nsamples))
proj_intensive_suburban_alt12 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params = intensive_suburban.params_alt12, times = times, name = "Improve indoor facility air quality")
proj_intensive_suburban <- rbind(proj_intensive_suburban_null, proj_intensive_suburban_alt1, proj_intensive_suburban_alt2, proj_intensive_suburban_alt3, proj_intensive_suburban_alt4, proj_intensive_suburban_alt5, proj_intensive_suburban_alt6, proj_intensive_suburban_alt7, proj_intensive_suburban_alt8, proj_intensive_suburban_alt9, proj_intensive_suburban_alt10, proj_intensive_suburban_alt11, proj_intensive_suburban_alt12)
Calculate prevalence, cumulative infections, and persistence for interventions in the intensive facility : wild, suburban system
persist.threshold <- 0.001
proj_intensive_suburban %>%
mutate(steady_sir = map(steady_state, "y"),
steady_sir = map(steady_sir, as_tibble_row)) %>%
pull(steady_sir) %>%
list_rbind(names_to = "run_id") %>%
mutate(., persistWild = I_wild > persist.threshold, persistCaptive = I_captive > persist.threshold) %>%
select(., run_id, persistWild, persistCaptive) %>%
arrange(., run_id) -> Persist_intensive_suburban_df
proj_intensive_suburban %>%
mutate(ode_df = map(ode_proj, as.data.frame)) %>%
pull(ode_df) %>%
list_rbind(names_to = "run_id") %>%
mutate(., Complex = rep("Intensive facility and suburban deer",nrow(.)), Alternative = c(rep(unique(proj_intensive_suburban$Context)[1],length(times)*nsamples),rep(unique(proj_intensive_suburban$Context)[2],length(times)*nsamples), rep(unique(proj_intensive_suburban$Context)[3],length(times)*nsamples),rep(unique(proj_intensive_suburban$Context)[4],length(times)*nsamples), rep(unique(proj_intensive_suburban$Context)[5],length(times)*nsamples), rep(unique(proj_intensive_suburban$Context)[6],length(times)*nsamples), rep(unique(proj_intensive_suburban$Context)[7],length(times)*nsamples),rep(unique(proj_intensive_suburban$Context)[8],length(times)*nsamples),rep(unique(proj_intensive_suburban$Context)[9],length(times)*nsamples),rep(unique(proj_intensive_suburban$Context)[10],length(times)*nsamples),rep(unique(proj_intensive_suburban$Context)[11],length(times)*nsamples),rep(unique(proj_intensive_suburban$Context)[12],length(times)*nsamples), rep(unique(proj_intensive_suburban$Context)[13],length(times)*nsamples))) %>%
group_by(Complex,Alternative, run_id) %>%
summarize(meanWild = mean(I_wild), meanCaptive = mean(I_captive),
cumulativeWild = last(I_wild_cumulative), cumulativeCaptive = last(I_captive_cumulative)) %>%
merge(.,Persist_intensive_suburban_df, by = "run_id") -> alternatives_for_intensive_suburban_complex
Combined system of intensive captive facilities and wild deer in a rural setting
intensive_rural.inits.fall <- initial_compartments(draws = nsamples)
intensive_rural.inits.steady <- initial_compartments(draws = nsamples, steady = TRUE)
intensive_rural.params <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_rural, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_rural, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_intensive_rural_null <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params, times = times, name = "Baseline")
#Alternatives
#1. 10% density reduction
omega_ww_alt1_density.red.10 <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.9)*(get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17)
intensive_rural.params_alt1 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_alt1_density.red.10, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_rural, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_intensive_rural_alt1 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params_alt1, times = times, name = "Reduce wild density by 10%")
#2. 25% density reduction
omega_ww_alt2_density.red.25 <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.75)*(get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17)
intensive_rural.params_alt2 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_alt2_density.red.25, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_rural, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_intensive_rural_alt2 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params_alt2, times = times, name = "Reduce wild density by 25%")
#3. Eliminate baiting
omega_ww_alt3_no.baiting <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.75)
intensive_rural.params_alt3 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_alt3_no.baiting, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_rural, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_intensive_rural_alt3 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params_alt3, times = times, name = "Eliminate baiting")
#4. Pause research (assuming random 1-50% of human-wild deer interaction)
omega_hw_alt4_pause.research <- omega_hw_rural*runif(nsamples, min = 0.5, max = 0.99) #reduce human deer proximity rate in the wild up to 50%
intensive_rural.params_alt4 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_rural, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_alt4_pause.research, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_intensive_rural_alt4 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params_alt4, times = times, name = "Pause Research")
#5. Require PPE for humans in wild settings
sigma_aero_deer_human_wild_alt5_PPE <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_rural/60, r = r_deer, V_d = 0.009*runif(nsamples, min = (1-0.994), max = (1-0.509)), nsamples = nsamples)
intensive_rural.params_alt5 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_rural, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_rural, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_intensive_rural_alt5 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params_alt5, times = times, name = "PPE, Wild")
#6. Require PPE for humans in captive settings
sigma_aero_deer_human_capt_alt6_PPE <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_intensive/60, r = r_deer, V_d = 0.009*runif(nsamples, min = (1-0.994), max = (1-0.509)), nsamples = nsamples)
intensive_rural.params_alt6 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_rural, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_rural, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_intensive_rural_alt6 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params_alt6, times = times, name = "PPE, Captive")
#7. Require PPE for humans in wild and captive settings (use derived parameters from 5 and 6)
intensive_rural.params_alt7 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_rural, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_rural, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
proj_intensive_rural_alt7 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params_alt7, times = times, name = "PPE, Both")
#8. Vaccinate captive deer (start with r_captive = 1)
intensive_rural.inits.fall_alt8 <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1)
intensive_rural.inits.steady_alt8 <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1, steady = TRUE)
proj_intensive_rural_alt8 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall_alt8, initial_compartments_steady = intensive_rural.inits.steady_alt8, params = intensive_rural.params, times = times, name = "Vaccinate captive deer")
#9. Vaccinate and boost captive deer (start with r_captive = 1)
intensive_rural.inits.fall_alt9 <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1)
intensive_rural.inits.steady_alt9 <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1, steady = TRUE)
boost_alt9 <- rep(0.9,nsamples)
intensive_rural.params_alt9 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_rural, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_rural, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = boost_alt9)
proj_intensive_rural_alt9 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall_alt9, initial_compartments_steady = intensive_rural.inits.steady_alt9, params = intensive_rural.params_alt9, times = times, name = "Vaccinate and boost captive deer")
#10. Double fence captive facilities
omega_cw_alt10 <- rep(0, nsamples)
intensive_rural.params_alt10 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_rural, omega_cw = omega_cw_alt10, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_rural, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0,nsamples))
proj_intensive_rural_alt10 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params_alt10, times = times, name = "Double-fence captive facilities")
#11. Restrict human interaction (using proportional reduction from suburban to rural)
human_restrictions <- (get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Rural (per 120 days)") / 120)/(get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Suburban (per 120 days)") / 120)
intensive_rural.params_alt11 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_rural, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_rural*human_restrictions, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0,nsamples))
proj_intensive_rural_alt11 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params_alt11, times = times, name = "Restrict human interactions")
#12. Improve indoor facility air exchange
sigma_aero_deer_human_intensive_alt12 <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_intensive/60, r = r_deer, nsamples = nsamples)
sigma_aero_deer_deer_intensive_alt12 <- calc_sigma_aero(C_nu = C_nu_deer, t_contact = t_contact_deer_deer_null/60, r = r_deer, nsamples = nsamples)
intensive_rural.params_alt12 <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_rural, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_rural, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive_alt12, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive_alt12, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0,nsamples))
proj_intensive_rural_alt12 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params_alt12, times = times, name = "Improve indoor facility air quality")
proj_intensive_rural <- rbind(proj_intensive_rural_null, proj_intensive_rural_alt1, proj_intensive_rural_alt2, proj_intensive_rural_alt3, proj_intensive_rural_alt4, proj_intensive_rural_alt5, proj_intensive_rural_alt6, proj_intensive_rural_alt7, proj_intensive_rural_alt8, proj_intensive_rural_alt9, proj_intensive_rural_alt10, proj_intensive_rural_alt11, proj_intensive_rural_alt12)
Calculate prevalence, cumulative infections, and persistence for interventions in the intensive facility : wild, rural system
persist.threshold <- 0.001
proj_intensive_rural %>%
mutate(steady_sir = map(steady_state, "y"),
steady_sir = map(steady_sir, as_tibble_row)) %>%
pull(steady_sir) %>%
list_rbind(names_to = "run_id") %>%
mutate(., persistWild = I_wild > persist.threshold, persistCaptive = I_captive > persist.threshold) %>%
select(., run_id, persistWild, persistCaptive) %>%
arrange(., run_id) -> Persist_intensive_rural_df
proj_intensive_rural %>%
mutate(ode_df = map(ode_proj, as.data.frame)) %>%
pull(ode_df) %>%
list_rbind(names_to = "run_id") %>%
mutate(., Complex = rep("Intensive facility and rural deer",nrow(.)), Alternative = c(rep(unique(proj_intensive_rural$Context)[1],length(times)*nsamples),rep(unique(proj_intensive_rural$Context)[2],length(times)*nsamples), rep(unique(proj_intensive_rural$Context)[3],length(times)*nsamples),rep(unique(proj_intensive_rural$Context)[4],length(times)*nsamples), rep(unique(proj_intensive_rural$Context)[5],length(times)*nsamples), rep(unique(proj_intensive_rural$Context)[6],length(times)*nsamples), rep(unique(proj_intensive_rural$Context)[7],length(times)*nsamples),rep(unique(proj_intensive_rural$Context)[8],length(times)*nsamples),rep(unique(proj_intensive_rural$Context)[9],length(times)*nsamples),rep(unique(proj_intensive_rural$Context)[10],length(times)*nsamples),rep(unique(proj_intensive_rural$Context)[11],length(times)*nsamples),rep(unique(proj_intensive_rural$Context)[12],length(times)*nsamples), rep(unique(proj_intensive_rural$Context)[13],length(times)*nsamples))) %>%
group_by(Complex,Alternative, run_id) %>%
summarize(meanWild = mean(I_wild), meanCaptive = mean(I_captive),
cumulativeWild = last(I_wild_cumulative), cumulativeCaptive = last(I_captive_cumulative)) %>%
merge(.,Persist_intensive_rural_df, by = "run_id") -> alternatives_for_intensive_rural_complex
4.Combine and summarize results
Below, we combine the results for all alternatives, for all
combinations of captive and wild scenarios (systems). These results are
saved in the package’s data folder, and can be accessed through
whitetailedSIRS::alternatives_across_systems
.
alternatives_across_systems <- rbind(alternatives_for_ranch_suburban_complex, alternatives_for_ranch_rural_complex, alternatives_for_intensive_suburban_complex, alternatives_for_intensive_rural_complex)
alternatives_across_systems %>%
group_by(., Complex, Alternative) %>%
summarize(., Wild_Cumulative = median(cumulativeWild), Captive_Cumulative = median(cumulativeCaptive))
#> `summarise()` has grouped output by 'Complex'. You can override using the
#> `.groups` argument.
#> # A tibble: 50 × 4
#> # Groups: Complex [4]
#> Complex Alternative Wild_Cumulative Captive_Cumulative
#> <chr> <chr> <dbl> <dbl>
#> 1 Intensive facility and rural … Baseline 0.861 1.47
#> 2 Intensive facility and rural … Double-fen… 0.494 1.47
#> 3 Intensive facility and rural … Eliminate … 0.00270 1.47
#> 4 Intensive facility and rural … Improve in… 0.832 1.14
#> 5 Intensive facility and rural … PPE, Both 0.856 1.46
#> 6 Intensive facility and rural … PPE, Capti… 0.857 1.46
#> 7 Intensive facility and rural … PPE, Wild 0.861 1.47
#> 8 Intensive facility and rural … Pause Rese… 0.861 1.47
#> 9 Intensive facility and rural … Reduce wil… 0.783 1.47
#> 10 Intensive facility and rural … Reduce wil… 0.521 1.47
#> # ℹ 40 more rows
5.Impact of Cross-Sector Collaboration
Finally, we simulated SARS-CoV-2 outbreaks in each system with
multiple alternatives applied. Here, we implement the most effective
alternative for each sector (agriculture, public health, and wildlife)
in concert to understand the cumulative effects of coordinated actions
across sectors. The results of this final simulation are stored in
whitetailedSIRS::cross_sector_results
.
#Vaccinate..
cross_sector.inits.fall <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1)
cross_sector.inits.steady <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1, steady = TRUE)
#...and boost
boost_alt9 <- rep(0.9,nsamples)
#PPE for all interactions
sigma_aero_deer_human_capt_alt6_PPE <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_ranch/60, r = r_deer, V_d = 0.009*runif(nsamples, min = (1-0.994), max = (1-0.509)), nsamples = nsamples)
sigma_aero_deer_human_wild_alt5_PPE <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_rural/60, r = r_deer, V_d = 0.009*runif(nsamples, min = (1-0.994), max = (1-0.509)), nsamples = nsamples)
#Eliminate baiting
omega_ww_alt3_no.baiting <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.75)
#Ranch and suburban complex
ranch_suburban_cross_sector.params <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_alt3_no.baiting, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_suburban, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = boost_alt9)
proj_ranch_suburban_cross_sector <- run(iter = nsamples, initial_compartments = cross_sector.inits.fall, initial_compartments_steady = cross_sector.inits.steady, params = ranch_suburban_cross_sector.params, times = times, name = "Cross-Sector Collaboration: Vaccinate and Boost, PPE for all interactions, and eliminating baiting")
proj_ranch_suburban_cross_sector %>%
mutate(steady_sir = map(steady_state, "y"),
steady_sir = map(steady_sir, as_tibble_row)) %>%
pull(steady_sir) %>%
list_rbind(names_to = "run_id") %>%
mutate(., persistWild = I_wild > persist.threshold, persistCaptive = I_captive > persist.threshold) %>%
select(., run_id, persistWild, persistCaptive) %>%
arrange(., run_id) -> Persist_ranch_suburban_cross_sector_df
proj_ranch_suburban_cross_sector %>%
mutate(ode_df = map(ode_proj, as.data.frame)) %>%
pull(ode_df) %>%
list_rbind(names_to = "run_id") %>%
mutate(., Complex = rep("Ranch and suburban deer",nrow(.)), Alternative = c(rep(unique(proj_ranch_suburban_cross_sector$Context)[1],length(times)*nsamples))) %>%
group_by(Complex,Alternative, run_id) %>%
summarize(meanWild = mean(I_wild), meanCaptive = mean(I_captive),
cumulativeWild = last(I_wild_cumulative), cumulativeCaptive = last(I_captive_cumulative)) %>%
merge(.,Persist_ranch_suburban_cross_sector_df, by = "run_id") -> cross_sector_for_ranch_suburban_complex
#> `summarise()` has grouped output by 'Complex', 'Alternative'. You can override
#> using the `.groups` argument.
#Ranch and rural complex
ranch_rural_cross_sector.params <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_alt3_no.baiting, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
omega_hw = omega_hw_rural, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = boost_alt9)
proj_ranch_rural_cross_sector <- run(iter = nsamples, initial_compartments = cross_sector.inits.fall, initial_compartments_steady = cross_sector.inits.steady, params = ranch_rural_cross_sector.params, times = times, name = "Cross-Sector Collaboration: Vaccinate and Boost, PPE for all interactions, and eliminating baiting")
proj_ranch_rural_cross_sector %>%
mutate(steady_sir = map(steady_state, "y"),
steady_sir = map(steady_sir, as_tibble_row)) %>%
pull(steady_sir) %>%
list_rbind(names_to = "run_id") %>%
mutate(., persistWild = I_wild > persist.threshold, persistCaptive = I_captive > persist.threshold) %>%
select(., run_id, persistWild, persistCaptive) %>%
arrange(., run_id) -> Persist_ranch_rural_cross_sector_df
proj_ranch_rural_cross_sector %>%
mutate(ode_df = map(ode_proj, as.data.frame)) %>%
pull(ode_df) %>%
list_rbind(names_to = "run_id") %>%
mutate(., Complex = rep("Ranch and rural deer",nrow(.)), Alternative = c(rep(unique(proj_ranch_rural_cross_sector$Context)[1],length(times)*nsamples))) %>%
group_by(Complex,Alternative, run_id) %>%
summarize(meanWild = mean(I_wild), meanCaptive = mean(I_captive),
cumulativeWild = last(I_wild_cumulative), cumulativeCaptive = last(I_captive_cumulative)) %>%
merge(.,Persist_ranch_rural_cross_sector_df, by = "run_id")-> cross_sector_for_ranch_rural_complex
#> `summarise()` has grouped output by 'Complex', 'Alternative'. You can override
#> using the `.groups` argument.
#Intensive facility and suburban complex
intensive_suburban_cross_sector.params <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_alt3_no.baiting, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_suburban, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = boost_alt9)
proj_intensive_suburban_cross_sector <- run(iter = nsamples, initial_compartments = cross_sector.inits.fall, initial_compartments_steady = cross_sector.inits.steady, params = intensive_suburban_cross_sector.params, times = times, name = "Cross-Sector Collaboration: Vaccinate and Boost, PPE for all interactions, and eliminating baiting")
proj_intensive_suburban_cross_sector %>%
mutate(steady_sir = map(steady_state, "y"),
steady_sir = map(steady_sir, as_tibble_row)) %>%
pull(steady_sir) %>%
list_rbind(names_to = "run_id") %>%
mutate(., persistWild = I_wild > persist.threshold, persistCaptive = I_captive > persist.threshold) %>%
select(., run_id, persistWild, persistCaptive) %>%
arrange(., run_id) -> Persist_intensive_suburban_cross_sector_df
proj_intensive_suburban_cross_sector %>%
mutate(ode_df = map(ode_proj, as.data.frame)) %>%
pull(ode_df) %>%
list_rbind(names_to = "run_id") %>%
mutate(., Complex = rep("Intensive facility and suburban deer",nrow(.)), Alternative = c(rep(unique(proj_intensive_suburban_cross_sector$Context)[1],length(times)*nsamples))) %>%
group_by(Complex,Alternative, run_id) %>%
summarize(meanWild = mean(I_wild), meanCaptive = mean(I_captive),
cumulativeWild = last(I_wild_cumulative), cumulativeCaptive = last(I_captive_cumulative)) %>%
merge(.,Persist_intensive_suburban_cross_sector_df, by = "run_id")-> cross_sector_for_intensive_suburban_complex
#> `summarise()` has grouped output by 'Complex', 'Alternative'. You can override
#> using the `.groups` argument.
#Intensive facility and rural complex
intensive_rural_cross_sector.params <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_alt3_no.baiting, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_rural, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = boost_alt9)
proj_intensive_rural_cross_sector <- run(iter = nsamples, initial_compartments = cross_sector.inits.fall, initial_compartments_steady = cross_sector.inits.steady, params = intensive_rural_cross_sector.params, times = times, name = "Cross-Sector Collaboration: Vaccinate and Boost, PPE for all interactions, and eliminating baiting")
proj_intensive_rural_cross_sector %>%
mutate(steady_sir = map(steady_state, "y"),
steady_sir = map(steady_sir, as_tibble_row)) %>%
pull(steady_sir) %>%
list_rbind(names_to = "run_id") %>%
mutate(., persistWild = I_wild > persist.threshold, persistCaptive = I_captive > persist.threshold) %>%
select(., run_id, persistWild, persistCaptive) %>%
arrange(., run_id) -> Persist_intensive_rural_cross_sector_df
proj_intensive_rural_cross_sector %>%
mutate(ode_df = map(ode_proj, as.data.frame)) %>%
pull(ode_df) %>%
list_rbind(names_to = "run_id") %>%
mutate(., Complex = rep("Intensive facility and rural deer",nrow(.)), Alternative = c(rep(unique(proj_intensive_rural_cross_sector$Context)[1],length(times)*nsamples))) %>%
group_by(Complex,Alternative, run_id) %>%
summarize(meanWild = mean(I_wild), meanCaptive = mean(I_captive),
cumulativeWild = last(I_wild_cumulative), cumulativeCaptive = last(I_captive_cumulative)) %>%
merge(.,Persist_intensive_rural_cross_sector_df, by = "run_id")-> cross_sector_for_intensive_rural_complex
#> `summarise()` has grouped output by 'Complex', 'Alternative'. You can override
#> using the `.groups` argument.
cross_sector_results <- rbind(cross_sector_for_ranch_suburban_complex, cross_sector_for_ranch_rural_complex, cross_sector_for_intensive_suburban_complex, cross_sector_for_intensive_rural_complex)
cross_sector_results %>%
group_by(Complex) %>%
summarize(., prevalenceWild = median(meanWild), prevalenceCaptive = median(meanCaptive), cumulativeWild = median(cumulativeWild), cumulativeCaptive = median(cumulativeCaptive), meanPersistWild = (sum(persistWild)/nsamples), meanPersistCaptive= (sum(persistCaptive)/nsamples))
#> # A tibble: 4 × 7
#> Complex prevalenceWild prevalenceCaptive cumulativeWild cumulativeCaptive
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 Intensive f… 0.00000147 0.000000661 0.0000313 0.0000140
#> 2 Intensive f… 0.00000523 0.000000682 0.000113 0.0000143
#> 3 Ranch and r… 0.00000147 0.000000101 0.0000313 0.00000226
#> 4 Ranch and s… 0.00000522 0.000000107 0.000113 0.00000226
#> # ℹ 2 more variables: meanPersistWild <dbl>, meanPersistCaptive <dbl>
We can take these results and compare to baseline results, as well as single management alternatives.