### ## Data figure A Seasonal net use [Koenker et al 2019] ## Data Figure B Individual variability in biting risk [Guglielmo et al 2021] ## A model comparison of ## - screening where insecticide wanes over 4 years but the barrier effect remains, ## - nets with waning efficacy and adherence to use over 3 years, ## - IRS withdrawn, then housing turned on devtools::load_all() library(malariasimulation) ########################################## ## ## malsim test to simulate trial arms (n = 4) ## later we will simulate every village (n = 53) year <- 365 month <- 30 sim_length <- 12 * year ## Jan 2020 - Dec 2031 ## background = ITN campaign every 3 years ## 60% treatment of people 45% ACT, 15% non-ACT ## IRS campaigns annually ## output children and adults ## HOUSE screening implemented in 2025 ## year 6 run_mod_f = function(eir, housing_coverages, ## vector length 2 itn_cov, ## vector length 4 # spray_cover, ## either 0 or 0.8 # end_irs, ## value of year for withdrawing IRS DH_IN, RH_IN, GAMMAH_IN, DN_IN, RN_IN, GAMMAN_IN ){ human_population <- 30000 starting_EIR <- eir ## to be estimated later simparams <- get_parameters( list( human_population = human_population, # irs_correlation = prevalence_rendering_min_ages = c(0.5, 5, 15, 0) * 365, ## Prev in 6 months to 5 years measured, prevalence_rendering_max_ages = c(5, 15, 49, 100) * 365, clinical_incidence_rendering_min_ages = c(0.5, 5, 15, 0) * 365, ## as prev clinical_incidence_rendering_max_ages = c(5, 15, 49, 100) * 365, ## Try for a Burkina / Tanzania seasonal scenario model_seasonality = TRUE, ## ## These are from rainfall # OLD seasonal_a0 seasonal_a1 seasonal_b1 seasonal_a2 seasonal_b2 seasonal_a3 seasonal_b3 g0 = 0.2852297, g = c(-0.2952712, -0.03408224, 0.07596435), h = c(-0.1126063, 0.07789561, -0.007051094), individual_mosquitoes = FALSE, ## Update next carrying_capacity = FALSE ) ) # set species simparams <- set_species(simparams, species=list(arab_params, fun_params), proportions = c(0.5,0.5)) ## default gambiae # set treatment # set treatment simparams <- set_drugs(simparams, list(AL_params, ## whichever is ACT drug SP_AQ_params))## whichever is non-ACT drug simparams <- set_clinical_treatment(simparams, drug=1, time=1, coverage=0.45) # simparams <- set_clinical_treatment(simparams, drug=2, time=1, coverage=0.15) # ## housing housingtimesteps <- c(year*6, ## year*9, year*12) housingparams <- set_housing( simparams, timesteps = housingtimesteps, coverages = housing_coverages, phi_housing = c(1,1,1), dh0 = matrix(rep(DH_IN,each=3), nrow=3, ncol=2), rh = matrix(rep(RH_IN,each=3), nrow=3, ncol=2), rhm = matrix(rep(0.25,6), nrow=3, ncol=2), ## This stays high (gets higher for survived...) gammah = matrix(rep(mean(GAMMAH_IN),each=3), nrow=3, ncol=1) ) ## Set up bed nets bednetparams <- housingparams ## as done bednet_events = data.frame( timestep = seq(1,12,by=3)*year - 364, ## first jan 2020 name=c("2020", "2023", "2026",## housing screens implemented simultaneously "2029") ) simparams <- set_bednets( bednetparams, timesteps = bednet_events$timestep, coverages = c(itn_cov), ## each 3 years, retention = 3 * year, ## observed during RCT at high retention levels thanks to workshops ## each row needs to show the efficacy parameter across years (and cols are diff mosquito) ## resistance 46% dn0 = matrix(rep(DN_IN,4), nrow=4, ncol=2), rn = matrix(rep(RN_IN,4), nrow=4, ncol=2), rnm = matrix(rep(.24, 8), nrow=4, ncol=2), gamman = as.numeric(rep(mean(GAMMAN_IN),4)) ) # # ## Set IRS # peak <- peak_season_offset(simparams) # # sprayingtimesteps <- seq(1, end_irs) * year + peak - 3 * month # A round of IRS is implemented in the 1st and second year 3 months prior to peak transmission. # # sprayparams <- set_spraying( # simparams, # timesteps = sprayingtimesteps, # coverages = rep(spray_cover,length(sprayingtimesteps)), # # The first round covers 30% of the population and the second covers 80%. # ls_theta = matrix(2.025, nrow=length(sprayingtimesteps), ncol=1), # Matrix of mortality parameters; nrows=length(timesteps), ncols=length(species) # ls_gamma = matrix(-0.009, nrow=length(sprayingtimesteps), ncol=1), # Matrix of mortality parameters per round of IRS and per species # ks_theta = matrix(-2.222, nrow=length(sprayingtimesteps), ncol=1), # Matrix of feeding success parameters per round of IRS and per species # ks_gamma = matrix(0.008, nrow=length(sprayingtimesteps), ncol=1), # Matrix of feeding success parameters per round of IRS and per species # ms_theta = matrix(-1.232, nrow=length(sprayingtimesteps), ncol=1), # Matrix of deterrence parameters per round of IRS and per species # ms_gamma = matrix(-0.009, nrow=length(sprayingtimesteps), ncol=1) # Matrix of deterrence parameters per round of IRS and per species # ) # # full_itn_params <- sprayparams ## assume the same people are getting nets each round correlations <- get_correlation_parameters(simparams) correlations$inter_round_rho('bednets', 1) simparams <- set_equilibrium(simparams, starting_EIR) ## Run the simulations output <- run_simulation(sim_length, simparams,correlations) output$prev_182.5_1825 = output$p_detect_182.5_1825/output$n_182.5_1825 output$prev_1825_5475 = output$p_detect_1825_5475/output$n_1825_5475 output$prev_0_36500 = output$p_detect_0_36500/output$n_0_36500 output$cases_182.5_1825 = output$n_inc_clinical_182.5_1825/output$n_182.5_1825 output$cases_1825_5475 = output$n_inc_clinical_1825_5475/output$n_1825_5475 output$cases_0_36500 = output$n_inc_clinical_0_36500/output$n_0_36500 return(list(output, simparams)) } ## Exploring the effects of the different assumptions on how ## housing improvement might alter the pathway for mosquitoes ##control ## Standard of care with all interventions on and no house screening ## Standard of care without IRS or house screening ## WANT TO LOOP THROUGH THE RANGES IN PARAMETER VALUES ## AND PRODUCE OUTPUT FIGURE... ## Generate input values dt = expand.grid(GAMMAH_IN = seq(1,5,length=11)*365/log(2), housing_coverages = seq(0,0.8,length=11)) dt2 = data.frame(id = c(501:516,1:121), info = c("counterfactual", rep("comparisons: meta_nets",3), rep("comparisons: sfs_nets",3), rep("comparisons: meta_nets_screens",3), rep("comparisons: sfs_nets_screens",3), rep("comparisons: screens",3), rep("sens",121)), housing_cov_yr6 = c(rep(0,7),rep(0.6,9),dt$housing_coverages), housing_cov_yr9 = c(rep(0,7),rep(0.6,9),dt$housing_coverages), itn_cov_yr6 = c(0,rep(0.6,12),0,0,0,rep(0.6,121)), itn_cov_yr9 = c(0,rep(0.6,12),0,0,0,rep(0.6,121)), DH_IN_ara = c(0.72, ## off: 0.72 (0.59 – 0.82) rep(0.72,3), ## off rep(0.72,3), ## off 0.59,0.72,0.82, ## on 0.59,0.72,0.82, ## on 0.59,0.72,0.82, ## on rep(0.72,121)), ## on sensitivity analysis DH_IN_fun = c(0.91, ## off: 0.91 (0.80 – 0.95) rep(0.91,3), ## off rep(0.91,3), ## off 0.80,0.91,0.95, ## on 0.80,0.91,0.95, ## on 0.80,0.91,0.95, ## on rep(0.91,121)), ## on sensitivity analysis RH_IN_ara = c(0.13, ## off:0.13 (0.04 – 0.25) rep(0.13,3), ## off rep(0.13,3), ## off 0.04,0.13,0.25, ## on 0.04,0.13,0.25, ## on 0.04,0.13,0.25, ## on rep(0,121)), ## on sensitivity analysis RH_IN_fun = c(0.04, ## off:0.04 (0.00 – 0.13) rep(0.04,3), ## off rep(0.04,3), ## off 0,0.04,0.13, ## on 0,0.04,0.13, ## on 0,0.04,0.13, ## on rep(0,121)), ## on sensitivity analysis GAMMAH_IN_ara= c(774, ## off rep(774,3), ## off rep(774,3), ## off 467,774,1795, ## on 467,774,1795, ## on 467,774,1795, ## on dt$GAMMAH_IN), ## on sensitivity analysis GAMMAH_IN_fun= c(1284, ## off rep(1284,3), ## off rep(1284,3), ## off 774,1284,2101, ## on 774,1284,2101, ## on 774,1284,2101, ## on dt$GAMMAH_IN), ## on sensitivity analysis DN_IN_ara = c(0.453905, ## off 0.261,0.453905,0.504,## meta 0.48,0.57,0.57, ## sfs: 0.57 (0.48 – 0.57) rep(0.453905,3), ## meta median rep(0.57,3), ## sfs median rep(0.453905,3), ## off rep(0.453905,121)), ## meta median DN_IN_fun = c(0.50096, 0.365,0.50096,0.541, 0.63, 0.72, 0.75,##:0.72 (0.63 – 0.75) rep(0.50096,3), rep(0.72,3), rep(0.50096,3), rep(0.50096,121)), RN_IN_ara = c(0.500, ## off 0.567,0.500,0.471,## meta 0.46,0.28,0.24, ## sfs:0.28 (0.46 – 0.24) rep(0.500,3), ## meta median rep(0.28,3), ## sfs median rep(0.500,3), ## off rep(0.500,121)), ## meta median RN_IN_fun = c(0.467, ## off 0.551,0.467,0.444,## meta 0.34,0.27,0.24, ## sfs: 0.27(0.34 – 0.24) rep(0.467,3), ## meta median rep(0.27,3), ## sfs median rep(0.467,3), ## off rep(0.467,121)), ## meta median GAMMAN_IN = c(902, 378,902,1277, 1136,1405,1472, ## 1392 (1016 – 1472) days ; 1418 (1257 – 1472) days rep(902,3), ## meta median rep(1405,3), ## sfs median rep(902,3), ## off rep(902,121)) ) ## Loop through the simulations screening = readRDS("C:/Users/esherrar/Documents/HOUSING GRANT 2023/odufuwa publication 2024/model_simulations") for(i in c(5,6,7,11,12,13)){ screening[[i]] = run_mod_f(eir = 180, housing_coverages = c(dt2$housing_cov_yr6[i],dt2$housing_cov_yr9[i],0), itn_cov = c(0.6,0.6,dt2$itn_cov_yr6[i],dt2$itn_cov_yr9[i]), # spray_cover = 0.8, # end_irs = 12, DH_IN = c(dt2$DH_IN_ara[i],dt2$DH_IN_fun[i]), ## arabiensis resistant and new 0.5856466 0.7070226 0.7989045 RH_IN = c(dt2$RH_IN_ara[i],dt2$RH_IN_fun[i]), ## Has to be summing to 1 with DH_IN and fed GAMMAH_IN = c(dt2$GAMMAH_IN_ara[i],dt2$GAMMAH_IN_fun[i]), ## how to get at this? Differs for killing and repellence DN_IN = c(dt2$DN_IN_ara[i],dt2$DN_IN_fun[i]), ## from SFS closest to bioassay survival of 51% so using these PBO ITN parameters from bet binomial RN_IN = c(dt2$RN_IN_ara[i],dt2$RN_IN_fun[i]), ## GAMMAN_IN = dt2$GAMMAN_IN[i]) print(i) } saveRDS(screening,"C:/Users/esherrar/Documents/HOUSING GRANT 2023/odufuwa publication 2024/model_simulations") ## screening = readRDS("C:/Users/esherrar/Documents/HOUSING GRANT 2023/odufuwa publication 2024/model_simulations") year = 365 layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE)) par(mar=c(4,4,1,1)) plot(screening[[1]][[1]]$prev_182.5_1825 ~ screening[[1]][[1]]$timestep, pch="", ylim = c(0,1), xlab = "Years",xlim = c(3*365,12*365), ylab = "Prevalence in 6 to 59 months (%)", yaxt="n",xaxt="n" ) axis(1, at = c(1:12)*year,labels = c(1:12)-6) axis(2, las = 2, at = seq(0,1,0.2),labels = seq(0,100,20)) housingtimesteps <- c(year*6, year*9, year*12) abline(v = housingtimesteps,lty = 1,col = "yellow", lwd=3) bednet_events = data.frame( timestep = seq(1,12,by=3)*year - 364 ## first jan 2020 ) abline(v = bednet_events$timestep,lty=2,col="grey",lwd=2) cols_compare = c("blue","gray50","gray80","orange","darkred","aquamarine3") lty_ITN = c(1,1,2,1,2,1) legend("topright",legend = c("counterfactual (no additional ITN campaign post year 0)", "Standard-of-care (60% ITN use & treatment), A. meta-analysis", "Standard-of-care (60% ITN use & treatment), B. SFS", "A with house screens 60% cover (at year 0)", "B with house screens 60% cover (at year 0)", "with ITN halted, and house screens instead"), col = cols_compare, lwd=3, lty = lty_ITN) lines(screening[[1]][[1]]$prev_182.5_1825 ~ screening[[1]][[1]]$timestep,lwd=2,col="blue") ## No housing but IRS and nets seq_lowers = c(2,5,8,11,14) seq_uppers = c(4,7,10,13,16) seq_median = c(3,6,9,12,15) for(i in 3:5){ polygon(c(screening[[seq_median[i]]][[1]]$timestep, rev(screening[[seq_median[i]]][[1]]$timestep)), c(screening[[seq_lowers[i]]][[1]]$prev_182.5_1825, rev(screening[[seq_uppers[i]]][[1]]$prev_182.5_1825)), col = adegenet::transp(cols_compare[i+1],0.4),border=NA) lines(screening[[seq_median[i]]][[1]]$prev_182.5_1825 ~ screening[[1]][[1]]$timestep, lwd=2,col=cols_compare[i+1], lty = lty_ITN[i+1]) ## No housing but IRS and nets } # for(i in 1:3){ lines(screening[[seq_median[i]]][[1]]$prev_182.5_1825 ~ screening[[1]][[1]]$timestep, lwd=2,col=cols_compare[i+1], lty = lty_ITN[i+1]) ## No housing but IRS and nets } ######################### ## estimates for reduction in prevalence ## estimates for cases averted ## counterfactual mean prev prev_u5_cfl = mean(screening[[1]][[1]]$prev_182.5_1825[c(6*365):c(9*365)]) prev_5_15_cfl = mean(screening[[1]][[1]]$prev_1825_5475[c(6*365):c(9*365)]) prev_all_Age_cfl = mean(screening[[1]][[1]]$prev_0_36500[c(6*365):c(9*365)]) ## counterfactual cases per 1000 people over 3 years case_u5_cfl = 1000*sum(screening[[1]][[1]]$cases_182.5_1825[c(6*365):c(9*365)]) case_5_15_cfl = 1000*sum(screening[[1]][[1]]$cases_1825_5475[c(6*365):c(9*365)]) case_all_Age_cfl = 1000*sum(screening[[1]][[1]]$cases_0_36500[c(6*365):c(9*365)]) ## treatment PBO meta-analysis prev_u5_pboMA_MIN = mean(screening[[2]][[1]]$prev_182.5_1825[c(6*365):c(9*365)]) prev_5_15_pboMA_MIN = mean(screening[[2]][[1]]$prev_1825_5475[c(6*365):c(9*365)]) prev_all_Age_pboMA_MIN = mean(screening[[2]][[1]]$prev_0_36500[c(6*365):c(9*365)]) prev_u5_pboMA = mean(screening[[3]][[1]]$prev_182.5_1825[c(6*365):c(9*365)]) prev_5_15_pboMA = mean(screening[[3]][[1]]$prev_1825_5475[c(6*365):c(9*365)]) prev_all_Age_pboMA = mean(screening[[3]][[1]]$prev_0_36500[c(6*365):c(9*365)]) prev_u5_pboMA_MAX = mean(screening[[4]][[1]]$prev_182.5_1825[c(6*365):c(9*365)]) prev_5_15_pboMA_MAX = mean(screening[[4]][[1]]$prev_1825_5475[c(6*365):c(9*365)]) prev_all_Age_pboMA_MAX = mean(screening[[4]][[1]]$prev_0_36500[c(6*365):c(9*365)]) (prev_u5_cfl - prev_u5_pboMA)/prev_u5_cfl (prev_u5_cfl - prev_u5_pboMA_MIN)/prev_u5_cfl (prev_u5_cfl - prev_u5_pboMA_MAX)/prev_u5_cfl ## CASES case_u5_pboMA_MIN = mean(screening[[2]][[1]]$cases_182.5_1825[c(6*365):c(9*365)]) case_5_15_pboMA_MIN = mean(screening[[2]][[1]]$cases_1825_5475[c(6*365):c(9*365)]) case_all_Age_pboMA_MIN = mean(screening[[2]][[1]]$cases_0_36500[c(6*365):c(9*365)]) case_u5_pboMA = mean(screening[[3]][[1]]$cases_182.5_1825[c(6*365):c(9*365)]) case_5_15_pboMA = mean(screening[[3]][[1]]$cases_1825_5475[c(6*365):c(9*365)]) case_all_Age_pboMA = mean(screening[[3]][[1]]$cases_0_36500[c(6*365):c(9*365)]) case_u5_pboMA_MAX = mean(screening[[4]][[1]]$cases_182.5_1825[c(6*365):c(9*365)]) case_5_15_pboMA_MAX = mean(screening[[4]][[1]]$cases_1825_5475[c(6*365):c(9*365)]) case_all_Age_pboMA_MAX = mean(screening[[4]][[1]]$cases_0_36500[c(6*365):c(9*365)]) ## treatment PBO SFS prev_u5_pboSFS_MIN = mean(screening[[5]][[1]]$prev_182.5_1825[c(6*365):c(9*365)]) prev_5_15_pboSFS_MIN = mean(screening[[5]][[1]]$prev_1825_5475[c(6*365):c(9*365)]) prev_all_Age_pboSFS_MIN = mean(screening[[5]][[1]]$prev_0_36500[c(6*365):c(9*365)]) prev_u5_pboSFS = mean(screening[[6]][[1]]$prev_182.5_1825[c(6*365):c(9*365)]) prev_5_15_pboSFS = mean(screening[[6]][[1]]$prev_1825_5475[c(6*365):c(9*365)]) prev_all_Age_pboSFS = mean(screening[[6]][[1]]$prev_0_36500[c(6*365):c(9*365)]) prev_u5_pboSFS_MAX = mean(screening[[7]][[1]]$prev_182.5_1825[c(6*365):c(9*365)]) prev_5_15_pboSFS_MAX = mean(screening[[7]][[1]]$prev_1825_5475[c(6*365):c(9*365)]) prev_all_Age_pboSFS_MAX = mean(screening[[7]][[1]]$prev_0_36500[c(6*365):c(9*365)]) (prev_u5_cfl - prev_u5_pboSFS)/prev_u5_cfl (prev_u5_cfl - prev_u5_pboSFS_MIN)/prev_u5_cfl (prev_u5_cfl - prev_u5_pboSFS_MAX)/prev_u5_cfl ## treatment PBO-MA + Screens prev_u5_pboMAsc_MIN = mean(screening[[8]][[1]]$prev_182.5_1825[c(6*365):c(9*365)]) prev_5_15_pboMAsc_MIN = mean(screening[[8]][[1]]$prev_1825_5475[c(6*365):c(9*365)]) prev_all_Age_pboMAsc_MIN = mean(screening[[8]][[1]]$prev_0_36500[c(6*365):c(9*365)]) prev_u5_pboMAsc = mean(screening[[9]][[1]]$prev_182.5_1825[c(6*365):c(9*365)]) prev_5_15_pboMAsc = mean(screening[[9]][[1]]$prev_1825_5475[c(6*365):c(9*365)]) prev_all_Age_pboMAsc = mean(screening[[9]][[1]]$prev_0_36500[c(6*365):c(9*365)]) prev_u5_pboMAsc_MAX = mean(screening[[10]][[1]]$prev_182.5_1825[c(6*365):c(9*365)]) prev_5_15_pboMAsc_MAX = mean(screening[[10]][[1]]$prev_1825_5475[c(6*365):c(9*365)]) prev_all_Age_pboMAsc_MAX = mean(screening[[10]][[1]]$prev_0_36500[c(6*365):c(9*365)]) (prev_u5_cfl - prev_u5_pboMAsc)/prev_u5_cfl (prev_u5_cfl - prev_u5_pboMAsc_MIN)/prev_u5_cfl (prev_u5_cfl - prev_u5_pboMAsc_MAX)/prev_u5_cfl ## treatment PBO-SFS + Screens prev_u5_pboSFSSC_MIN = mean(screening[[11]][[1]]$prev_182.5_1825[c(6*365):c(9*365)]) prev_5_15_pboSFSSC_MIN = mean(screening[[11]][[1]]$prev_1825_5475[c(6*365):c(9*365)]) prev_all_Age_pboSFSSC_MIN = mean(screening[[11]][[1]]$prev_0_36500[c(6*365):c(9*365)]) prev_u5_pboSFSSC = mean(screening[[12]][[1]]$prev_182.5_1825[c(6*365):c(9*365)]) prev_5_15_pboSFSSC = mean(screening[[12]][[1]]$prev_1825_5475[c(6*365):c(9*365)]) prev_all_Age_pboSFSSC = mean(screening[[12]][[1]]$prev_0_36500[c(6*365):c(9*365)]) prev_u5_pboSFSSC_MAX = mean(screening[[13]][[1]]$prev_182.5_1825[c(6*365):c(9*365)]) prev_5_15_pboSFSSC_MAX = mean(screening[[13]][[1]]$prev_1825_5475[c(6*365):c(9*365)]) prev_all_Age_pboSFSSC_MAX = mean(screening[[13]][[1]]$prev_0_36500[c(6*365):c(9*365)]) (prev_u5_cfl - prev_u5_pboSFSSC)/prev_u5_cfl (prev_u5_cfl - prev_u5_pboSFSSC_MIN)/prev_u5_cfl (prev_u5_cfl - prev_u5_pboSFSSC_MAX)/prev_u5_cfl ## treatment Screens prev_u5_SC_MIN = mean(screening[[14]][[1]]$prev_182.5_1825[c(6*365):c(9*365)]) prev_5_15_SC_MIN = mean(screening[[14]][[1]]$prev_1825_5475[c(6*365):c(9*365)]) prev_all_Age_SC_MIN = mean(screening[[14]][[1]]$prev_0_36500[c(6*365):c(9*365)]) prev_u5_SC = mean(screening[[15]][[1]]$prev_182.5_1825[c(6*365):c(9*365)]) prev_5_15_SC = mean(screening[[15]][[1]]$prev_1825_5475[c(6*365):c(9*365)]) prev_all_Age_SC = mean(screening[[15]][[1]]$prev_0_36500[c(6*365):c(9*365)]) prev_u5_SC_MAX = mean(screening[[16]][[1]]$prev_182.5_1825[c(6*365):c(9*365)]) prev_5_15_SC_MAX = mean(screening[[16]][[1]]$prev_1825_5475[c(6*365):c(9*365)]) prev_all_Age_SC_MAX = mean(screening[[16]][[1]]$prev_0_36500[c(6*365):c(9*365)]) (prev_u5_cfl - prev_u5_SC)/prev_u5_cfl (prev_u5_cfl - prev_u5_SC_MIN)/prev_u5_cfl (prev_u5_cfl - prev_u5_SC_MAX)/prev_u5_cfl (prev_u5_pboSFS - prev_u5_SC)/prev_u5_pboSFS (prev_u5_pboSFS - prev_u5_SC_MIN)/prev_u5_pboSFS (prev_u5_pboSFS - prev_u5_SC_MAX)/prev_u5_pboSFS ########################################## ## ## Add the heatmap plot ## treatment sensitivity runs sa_prev_u5 = sa_prev_5_15 = sa_prev_all_Age = sa_case_u5 = sa_case_5_15 = sa_case_all = numeric(121) for(i in 1: 121){ sa_prev_u5[i] = mean(screening[[i+16]][[1]]$prev_182.5_1825[c(6*365):c(9*365)]) sa_prev_5_15[i] = mean(screening[[i+16]][[1]]$prev_1825_5475[c(6*365):c(9*365)]) sa_prev_all_Age[i] = mean(screening[[i+16]][[1]]$prev_0_36500[c(6*365):c(9*365)]) sa_case_u5[i] = 1000*sum(screening[[i+16]][[1]]$cases_182.5_1825[c(6*365):c(9*365)]) sa_case_5_15[i] = 1000*sum(screening[[i+16]][[1]]$cases_1825_5475[c(6*365):c(9*365)]) sa_case_all[i] = 1000*sum(screening[[i+16]][[1]]$cases_0_36500[c(6*365):c(9*365)]) } # library library(latticeExtra) # create data reduction_prev_u5 = reduction_prev_5_15 = reduction_prev_all_age = numeric(121) cases_averted_u5 = cases_averted_5_15 = cases_averted_all_age = numeric(121) for(i in 1:121){ reduction_prev_u5[i] = 100*c(prev_u5_cfl - sa_prev_u5[i])/prev_u5_cfl reduction_prev_5_15[i] = 100*c(prev_5_15_cfl - sa_prev_5_15[i])/prev_5_15_cfl reduction_prev_all_age[i] = 100*c(prev_all_Age_cfl - sa_prev_all_Age[i])/prev_all_Age_cfl cases_averted_u5[i] = 100*(case_u5_cfl - sa_case_u5[i])/case_u5_cfl cases_averted_5_15[i] = 100*(case_5_15_cfl - sa_case_5_15[i])/case_5_15_cfl cases_averted_all_age[i] = 100*(case_all_Age_cfl - sa_case_all[i])/case_all_Age_cfl } set.seed(1) data <- data.frame(x = dt$housing_coverages, ## coverage y = dt$GAMMAH_IN, z = reduction_prev_u5, z2 = cases_averted_u5, z3 = cases_averted_5_15, z4 = cases_averted_all_age) ## durability # showing data points on the same color scale aa = levelplot(z ~ x * y, data, panel = panel.levelplot.points, cex = 1.2, ylab = "Mean duration of mortality impact in days", xlab = "Intervention cover for population (%)", # scales = list(x = seq(0,80,20), y = data$y), pch="" ) + layer_(panel.2dsmoother(..., n = 200)) # showing data points on the same color scale ## All age cases averyred bb = levelplot(z4 ~ x * y, data, panel = panel.levelplot.points, cex = 1.2, ylab = "Mean duration of mortality impact in days", xlab = "Intervention cover for population (%)", col.regions = hcl.colors(110, palette = "Red-blue", rev = T), contour = FALSE, pch="" ) + layer_(panel.2dsmoother(..., n = 200)) library(gridExtra) print(aa, split=c(1,1,2,1), more=TRUE) print(bb, split=c(2,1,2,1)) data2 = matrix(data$z2,nrow=11, ncol = 11) heatmap(data2, Colv = NA, Rowv = NA, xlab="ITENS & ITWS percentage cover", ylab="Durability") # ## OR # # install.packages('plot.matrix') # # numeric matrix # x <- matrix(runif(35), ncol=5) # create a numeric matrix object # class(x) # #> [1] "matrix" "array" # par(mar=c(5.1, 4.1, 4.1, 4.1)) # adapt margins # plot.matrix::plot(x)