In ian-hastings/rotations: Simulating the Evolution of Insecticide Resistance under Rotations
library(rotations)
library(ggplot2)
library(gridExtra)
outFolder <- "C:\\Dropbox\\resistanceResults\\rotations\\"
10,000 runs comparing rotations and sequences.
#eval=FALSE just to temporarily hide this text
# Starting question : is there any benefit to rotations when there are no costs and no migration to an untreated area ?
#
# Test 1. No refugia (i.e. Coverage=100%) or fitness costs. Would expect all policies to be equivalent.
#
# Test 2. Dominance of resistance and fitness costs both set to 0.5.
# I suspect one of the advantages of rotations (and the presence of refugia) is that they maintain lower resistance allele frequencies which slows the spread of recessive resistance alleles. I suspect/guess/intuit that if we make them semi-dominate these effects should disappear and all policies should be the same. I’m probably wrong but it’s a nice test…..
#
#\pagebreak
#results='hide',
#comment=NA removes ## before text output
#start small & build it up
#maybe create an equiv. of resistance::sensiAnPaperPart.r later
# example for 1 run
# linputs <- list(max_gen = 50, n_insecticides = 2)
# do.call(run_rot, linputs)
# dfinputs <- as.data.frame(linputs)
# So one way of doing is to create an linputs list for each run & pass to run_rot()
# so run_rot() only ever runs for 1 scenario
# in contrast resistance::runModel2() accepted an input matrix & would run for multiple scenarios
# not sure which is better
# set random seed so that results can be reproduced
#set.seed(1)
set.seed(2)
max_gen <- 500
plot_scenarios <- FALSE
n_scenarios <- 10 #000
already_ran <- FALSE
#already_ran <- TRUE
# if already done runs, just load data for plots
if (already_ran)
{
#C:/rsprojects/rotations/inst/documents/
#load("df_res_all_rot10000.rda")
load("df_res_all_rot_mig1000.rda")
}else
{
df_res_all <- data.frame()
for(i in 1:n_scenarios)
{
#message should output to the R markdown console
if (i%%200 == 1) message("scenario ",i," of ",n_scenarios," ",Sys.time())
# inputs (could also be saved to a dataframe first & extracted one row at a time)
linputs <- list()
## constant inputs ##
linputs$coverage <- 1
linputs$migration <- 0
#linputs$cost <- 0 #variable cost below
linputs$max_gen <- max_gen
#beware min_rwr_interval, it's to stop repeated small rotations
linputs$min_rwr_interval <- 5
linputs$min_gens_switch_back <- 1
linputs$no_r_below_start <- FALSE #BEWARE this has been set to TRUE in all earlier runs
linputs$no_r_below_mut <- FALSE #BEWARE trying out new mutation-selection balance option
#set plot to F if using grid.arrange as in phaps-rotations-better.Rmd
linputs$plot <- FALSE
## variable inputs ##
linputs$cost <- runif(1, min=0, max=0.1)
#linputs$cost <- 0
linputs$expo_hi <- runif(1, min=0.4, max=0.9) #min=0.1, max=0.9
linputs$male_expo_prop <- runif(1, min=0, max=1)
linputs$eff <- runif(1, min=0.5, max=1) # min=0.3, max=1
linputs$rr <- runif(1, min=0.1, max=0.9)
linputs$dom_sel <- runif(1, min=0, max=1)
linputs$dom_cos <- runif(1, min=0, max=1)
#TODO make this lognormal to get lower frequencies too
linputs$start_freqs <- runif(1, min=0.005, max=0.1) #note not log yet
linputs$n_insecticides <- sample( 2:5, 1 ) #beware set replace=T if more than 1
# dispersal (initially with cost set to 0)
#linputs$coverage <- runif(1, min=0.1, max=0.9)
#linputs$migration <- runif(1, min=0.1, max=0.9)
# create null objects used later for plotting
ggrot <- ggseq <- NULL
# to compare sequence to rotation run all the scenarios above for both seq & rot
for(rot_or_not in 0:1)
{
#to get rot_int of 0 every other run
if (rot_or_not == 0) linputs$rot_interval <- 0
else #random rot interval integer
linputs$rot_interval <- sample( 5:50, 1 ) #beware set replace=T if more than 1
##################
# run one scenario
dfres <- do.call(run_rot, linputs)
# summarise results per scenario
# note this is also currently calc at end of run_rot() here is probably better
res <- gens_under_thresh(dfres)
# to plot results for checking
if ( linputs$rot_interval == 0 )
{
gens_seq <- as.numeric(res)
} else
{
gens_rot <- as.numeric(res)
}
if (plot_scenarios)
{
# make plot
ggres <- rot_plot_resistance(dfres, plot_refuge=FALSE, logy=FALSE, plot=FALSE)
# to plot results for checking
if ( linputs$rot_interval == 0 )
{
ggseq <- ggres
gens_seq <- as.numeric(res)
} else
{
ggrot <- ggres
gens_rot <- as.numeric(res)
}
}
# TODO I may want to concat all dfres into one object (adda column for scenario num) to save for later analysis
# e.g. if I wanted to calculate average resistance frequencies over time in addition to gens below threshold
} #end rot_or_not loop
# make an output row for each scenario with inputs, gens_rot, gens_seq
# need to make sure I retain the rot_interval for rot
# add scenario id column onto results
df_in_out <- as.data.frame(c(linputs, gens_rot=gens_rot, gens_seq=gens_seq, id=i))
df_in_out$rot_minus_seq <- df_in_out$gens_rot - df_in_out$gens_seq
# previously I used rbind this scenario inputs & results onto previous
# not a good speed or memory way of doing
#df_res_all <- rbind(df_res_all, df_in_out)
# now replaced with this better way
# set a dataframe for the outputs to the required size in scenario1
if (i==1)
{
df_res_all <- data.frame(matrix(NA, ncol=length(names(df_in_out)), nrow=n_scenarios))
names(df_res_all) <- names(df_in_out)
}
#putting results into one row of overall results dataframe
df_res_all[i,] <- df_in_out
if (plot_scenarios)
{
# plot input vals
# extract out those between 0 & 1
dfi <- data.frame(t(as.data.frame(linputs)))
names(dfi) <- 'value'
dfi$inputs <- rownames(dfi)
#take out columns not between 0 & 1
dfi <- filter(dfi, !inputs %in% c('max_gen','rot_interval','n_insecticides','plot','min_rwr_interval'))
#also remove coverage & migration both set to 0 for now
dfi <- filter(dfi, !inputs %in% c('coverage','migration'))
ggins <- ggplot(dfi, aes(x=inputs,y=value)) +
geom_point() +
theme_minimal() +
ylim(0,1) +
theme(axis.text.x = element_text(angle = 20,hjust = 1, vjust = 1))
grid.arrange(ggrot, ggseq, ggins, ncol=1)
cat('scenario',i,'sum generations deployed insecticide under thresh rotation:', gens_rot,' sequence:',gens_seq,"\n\n")
#can I give the code to rerun this scenario
cat(paste0("run_rot(n_insecticides=",linputs$n_insecticides,", ",
"cost=",signif(linputs$cost,2),", ",
"start_freqs=",signif(linputs$start_freqs,2),", ",
"rot_interval=",linputs$rot_interval,", ",
"eff=",signif(linputs$eff,2),",\n",
"dom_sel=",signif(linputs$dom_sel,2),", ",
"dom_cos=",signif(linputs$dom_cos,2),", ",
"rr=",signif(linputs$rr,2),", ",
"expo_hi=",signif(linputs$expo_hi,2),", ",
"coverage=",signif(linputs$coverage,2),",\n",
"migration=",signif(linputs$migration,2),", ",
"max_gen=",linputs$max_gen,", ",
"min_rwr_interval=",linputs$min_rwr_interval,
")\n") )
}
#I could save partial outputs every 200 scenarios
#just in case doesn't make to end
if (i%%1000 == 0)
{
#df_res_all$rot_minus_seq <- df_res_all$gens_rot - df_res_all$gens_seq
# save object containing inputs and results as rda for analysis
save(df_res_all, file='df_res_all_rot.rda') #paste0(outFolder,'*.rda'))
}
} #end scenario loop
}
proportion of scenarios not reaching resistance thresholds
# check proportion of scenarios not reaching thresholds by max_runs for both seq & rot
nruns_fail <- df_res_all %>%
filter(gens_seq == max_gen & gens_rot == max_gen) %>%
summarise(n=n()) %>%
unlist()
cat('proportion runs not reaching thresholds : ', nruns_fail/nrow(df_res_all) )
# remove runs where thresholds not reached for seq & rot
df_res_all <- df_res_all %>%
filter(gens_seq != max_gen | gens_rot != max_gen)
nruns_same <- df_res_all %>%
filter(gens_seq == gens_rot) %>%
summarise(n=n()) %>%
unlist()
cat('proportion runs (that did reach threshold) where rotations and sequence the same : ', nruns_same/nrow(df_res_all) )
nruns_rot_20pc_better <- df_res_all %>%
filter(rot_minus_seq/gens_seq >= 0.2 ) %>%
summarise(n=n()) %>%
unlist()
cat('proportion runs (that did reach threshold) where rotations >20 % better than sequence : ', nruns_rot_20pc_better/nrow(df_res_all) )
difference between rotations and sequences across all runs
mostly (~67%) no difference but infrequent where rotations better
~ 7% of runs, rotation at least 20% better than sequence
print( ggplot(df_res_all, aes(x=rot_minus_seq)) +
geom_histogram(binwidth=10) +
xlab('rot_minus_seq : rotation generations minus sequence') +
ylab('count (sqrt scale to show low values)') +
scale_y_sqrt(breaks=c(0,10,20,50,200,500,2000,5000)) + #can't use log scale because anchored at 0
theme_minimal() +
theme(panel.grid.minor.y=element_blank()))
\pagebreak
Examples of runs where rotations > 20% better than sequences
#eval false because smaller runs don't have these
scenarios_to_plot <- c(624,2167,602,4028,5771)
for(scenario in scenarios_to_plot)
{
cat(scenario,"\n")
linputs <- as.list( df_res_all[scenario, 1:15] )
# run one scenario
dfres <- do.call(run_rot, linputs)
rot_plot_resistance(dfres, plot_refuge=FALSE, logy=FALSE, plot=TRUE)
#now run as sequence
linputs$rot_interval <- 0
dfres <- do.call(run_rot, linputs)
rot_plot_resistance(dfres, plot_refuge=FALSE, logy=FALSE, plot=TRUE)
dfi <- data.frame(t(as.data.frame(linputs)))
names(dfi) <- 'value'
dfi$inputs <- rownames(dfi)
#take out columns not between 0 & 1
dfi <- filter(dfi, !inputs %in% c('max_gen','rot_interval','n_insecticides','plot','min_rwr_interval'))
#also remove coverage & migration both set to 0 for now
dfi <- filter(dfi, !inputs %in% c('coverage','migration'))
ggins <- ggplot(dfi, aes(x=inputs,y=value)) +
geom_point() +
theme_minimal() +
ylim(0,1) +
theme(axis.text.x = element_text(angle = 20,hjust = 1, vjust = 1))
plot(ggins)
}
How does the difference between rotations and sequences respond to model inputs ?
# go through selected inputs and plot them on x with rot_minus_seq on y
# initially with 200 reps doesn't show anything obvious, try it with 2000 overnight runs
in_cols <- c('cost','dom_cos','dom_sel','rr','expo_hi','male_expo_prop','eff','start_freqs')
y <- 'rot_minus_seq'
for(i in in_cols)
{
#plotnum <- plotnum + 1
#print( ggplot(ggInsOuts, aes_string(x=i, y=y, colour="strategy")) +
#plotlist[[plotnum]] <-
print( ggplot(df_res_all, aes_string(x=i, y=y)) +
geom_point(shape=1, size=1, alpha=0.5, show.legend=FALSE) +
theme_minimal() +
geom_smooth(linetype='dashed',size=0.5) )
}
PRCC
library(sensitivity)
# pcc(X, y, rank = FALSE, nboot = 0, conf = 0.95)
#
# Arguments
# X a data frame containing model input variables
# y a vector containing the responses
# rank logical. If TRUE, the analysis is done on the ranks.
# nboot the number of bootstrap replicates.
# conf the confidence level of the bootstrap confidence intervals.
#same as above
in_cols <- c('cost','dom_cos','dom_sel','rr','expo_hi','male_expo_prop','eff','start_freqs')
x <- df_res_all[,in_cols]
y <- df_res_all['rot_minus_seq']
pcc_res <- pcc(x, y, rank=TRUE)
#if you add bootstrap, error bars are added to the default plot & extra columns to the PRCC object
#pcc_res <- pcc(x, y, rank=TRUE, nboot=100)
#plot(pcc_res)
#results are here I can probably rbind them together into a df that I can ggplot
#pcc_res$PRCC
to_plot <- pcc_res$PRCC
#rename column 1 from 'original to PRCC
names(to_plot)[1] <- 'PRCC'
to_plot$inputs <- rownames(to_plot)
axis_text_x <- element_text(angle = 45,hjust = 1, vjust = 1)
height <- 0.9
gg <- ggplot( to_plot, aes_string(x='inputs',y='PRCC') ) +
geom_point(shape=1, colour='red') +
theme(axis.text.x = axis_text_x) +
geom_hline(yintercept = 0, linetype=3) +
#theme(panel.grid.major.y = element_line(color = "grey", size = 1)) +
ylim(-1,1) +
#ggtitle(paste(strat_name[strategy_num])) +
xlab(NULL)
plot(gg)
#to put back in some ref gridlines that are removed by cowplot
#theme_set(theme_gray())
#theme_set(theme_bw())
#plot_grid( plotlist[[1]],plotlist[[2]],plotlist[[3]],plotlist[[4]],ncol=1, rel_heights=c(1,1,1,1.8))
#plot_grid( plotlist[[1]],plotlist[[2]],plotlist[[3]],ncol=1, rel_heights=c(1,1,1.6), labels='AUTO')
Runs where rotations are >20% better (in red) are associated with higher dominance of cost and lower dominance of selection
but that does not guarantee, must be interaction with other inputs too.
#'cost','dom_cos','dom_sel','rr','expo_hi','male_expo_prop','eff','start_freqs'
#'cost','dom_cos','dom_sel'
y <- 'rot_minus_seq'
ggplot(df_res_all, aes(x=dom_sel, y=dom_cos)) +
geom_point(shape=3, size=1, alpha=0.2, show.legend=FALSE) +
theme_minimal() +
geom_point(data=filter(df_res_all, rot_minus_seq/gens_seq >= 0.3 ), shape=1, size=1, col='red', alpha=0.8)
# looking at modelling it
# library(mgcv)
# mod_lm2 <- gam(rot_minus_seq ~ cost + dom_cos + dom_sel, data=df_res_all)
# summary(mod_lm2)
#
# mod_gam2 <- gam(rot_minus_seq ~ s(cost) + s(dom_cos) + s(dom_sel) + s(rr) + s(expo_hi) + s(male_expo_prop) + s(eff) + s(start_freqs), data=df_res_all)
#
# summary(mod_gam2)
# #Deviance explained = 19.5%
#
# mod_gam3 <- gam(rot_minus_seq ~ te(cost,dom_cos,dom_sel), data=df_res_all)
# summary(mod_gam3)
#actually decision trees may be useful
#eval=FALSE tmp because changes above will break this
# summary plot of all runs, would show any differences between sequence & rotations as non horizontal lines
# want to plot line between sequence and rotation value
# any deviation from horizontal will be obvious
# could use rotation interval as the x axis, so that longer intervals will be clearer
# maybe not for publication
# but useful for quick view
#it was useful for me but now that I've moved to one row per scenario it won't work !!
gg <- ggplot( df_res_all, aes_string(x='rot_interval',
y='tot_gens_dep_under50',
group='scenario',
colour='scenario') ) +
geom_line( alpha=0.5, lwd=1 )
plot(gg)
# plot gens againts exposure and effectiveness
# as a quick test to see if I can set the starting values of these higher
gg <- ggplot( df_res_all, aes_string(x='eff',
y='tot_gens_dep_under50') ) +
geom_point( alpha=0.5 )
plot(gg)
gg <- ggplot( df_res_all, aes_string(x='expo_hi',
y='tot_gens_dep_under50') ) +
geom_point( alpha=0.5 )
plot(gg)
gg <- ggplot( df_res_all, aes_string(x='expo_hi*eff*male_expo_prop',
y='tot_gens_dep_under50') ) +
geom_point( alpha=0.5 )
plot(gg)
Have a go at running for insecticides with different parameters
starting to work ...
# set random seed so that results can be reproduced
set.seed(2)
max_gen <- 500
plot_scenarios <- FALSE
n_scenarios <- 100 #1000 #000
#df_res_all <- data.frame()
for(i in 1:n_scenarios)
{
#message should output to the R markdown console
if (i%%200 == 1) message("scenario ",i," of ",n_scenarios," ",Sys.time())
# inputs (could also be saved to a dataframe first & extracted one row at a time)
linputs <- list()
## constant inputs ##
linputs$coverage <- 1
linputs$migration <- 0
#linputs$cost <- 0 #variable cost below
linputs$max_gen <- max_gen
#beware min_rwr_interval, it's to stop repeated small rotations
linputs$min_rwr_interval <- 5
linputs$min_gens_switch_back <- 1
linputs$no_r_below_start <- FALSE #BEWARE this has been set to TRUE in all earlier runs
linputs$no_r_below_mut <- TRUE #BEWARE trying out new mutation-selection balance option
#set plot to F if using grid.arrange as in phaps-rotations-better.Rmd
linputs$plot <- FALSE
## variable inputs ##
linputs$n_insecticides <- sample( 2:5, 1 ) #beware set replace=T if more than 1
n_ins <- linputs$n_insecticides
## new those below can be different for each insecticide
linputs$cost <- runif(n_ins, min=0, max=0.1)
#linputs$cost <- 0
linputs$expo_hi <- runif(n_ins, min=0.4, max=0.9) #min=0.1, max=0.9
linputs$male_expo_prop <- runif(1, min=0, max=1)
linputs$eff <- runif(n_ins, min=0.5, max=1) # min=0.3, max=1
linputs$rr <- runif(n_ins, min=0.1, max=0.9)
linputs$dom_sel <- runif(n_ins, min=0, max=1)
linputs$dom_cos <- runif(n_ins, min=0, max=1)
#DIFFERENT
linputs$start_freqs <- runif(n_ins, min=0.0001, max=0.01) #note not log yet
# dispersal (initially with cost set to 0)
#linputs$coverage <- runif(1, min=0.1, max=0.9)
#linputs$migration <- runif(1, min=0.1, max=0.9)
# create null objects used later for plotting
ggrot <- ggseq <- NULL
# to compare sequence to rotation run all the scenarios above for both seq & rot
for(rot_or_not in 0:1)
{
#to get rot_int of 0 every other run
if (rot_or_not == 0) linputs$rot_interval <- 0
else #random rot interval integer
linputs$rot_interval <- sample( 5:50, 1 ) #beware set replace=T if more than 1
##################
# run one scenario
dfres <- do.call(run_rot, linputs)
# summarise results per scenario
# note this is also currently calc at end of run_rot() here is probably better
res <- gens_under_thresh(dfres)
# to plot results for checking
if ( linputs$rot_interval == 0 )
{
gens_seq <- as.numeric(res)
} else
{
gens_rot <- as.numeric(res)
}
if (plot_scenarios)
{
# make plot
ggres <- rot_plot_resistance(dfres, plot_refuge=FALSE, logy=FALSE, plot=FALSE)
# to plot results for checking
if ( linputs$rot_interval == 0 )
{
ggseq <- ggres
gens_seq <- as.numeric(res)
} else
{
ggrot <- ggres
gens_rot <- as.numeric(res)
}
}
} #end rot_or_not loop
# make an output row for each scenario with inputs, gens_rot, gens_seq
# need to make sure I retain the rot_interval for rot
# to cope with multiple values per scenario convert them to strings
linputs2 <- lapply(linputs, function(x) ifelse(length(x)>1,toString(x),x))
#added scenario id column
df_in_out <- as.data.frame(c(linputs2, gens_rot=gens_rot, gens_seq=gens_seq, id=i), stringsAsFactors=FALSE)
df_in_out$rot_minus_seq <- df_in_out$gens_rot - df_in_out$gens_seq
#instead set up the dataframe to the required size in scenario1
if (i==1)
{
df_res_all <- data.frame(matrix(NA, ncol=length(names(df_in_out)), nrow=n_scenarios))
names(df_res_all) <- names(df_in_out)
}
#putting results into one row of overall results dataframe
df_res_all[i,] <- df_in_out
if (plot_scenarios)
{
# plot input vals
# extract out those between 0 & 1
dfi <- data.frame(t(as.data.frame(linputs)))
names(dfi) <- 'value'
dfi$inputs <- rownames(dfi)
#take out columns not between 0 & 1
dfi <- filter(dfi, !inputs %in% c('max_gen','rot_interval','n_insecticides','plot','min_rwr_interval'))
#also remove coverage & migration both set to 0 for now
dfi <- filter(dfi, !inputs %in% c('coverage','migration'))
ggins <- ggplot(dfi, aes(x=inputs,y=value)) +
geom_point() +
theme_minimal() +
ylim(0,1) +
theme(axis.text.x = element_text(angle = 20,hjust = 1, vjust = 1))
grid.arrange(ggrot, ggseq, ggins, ncol=1)
cat('scenario',i,'sum generations deployed insecticide under thresh rotation:', gens_rot,' sequence:',gens_seq,"\n\n")
#can I give the code to rerun this scenario
cat(paste0("run_rot(n_insecticides=",linputs$n_insecticides,", ",
"cost=",signif(linputs$cost,2),", ",
"start_freqs=",signif(linputs$start_freqs,2),", ",
"rot_interval=",linputs$rot_interval,", ",
"eff=",signif(linputs$eff,2),",\n",
"dom_sel=",signif(linputs$dom_sel,2),", ",
"dom_cos=",signif(linputs$dom_cos,2),", ",
"rr=",signif(linputs$rr,2),", ",
"expo_hi=",signif(linputs$expo_hi,2),", ",
"coverage=",signif(linputs$coverage,2),",\n",
"migration=",signif(linputs$migration,2),", ",
"max_gen=",linputs$max_gen,", ",
"min_rwr_interval=",linputs$min_rwr_interval,
")\n") )
}
#save partial outputs just in case doesn't make to end
if (i%%1000 == 0)
{
#df_res_all$rot_minus_seq <- df_res_all$gens_rot - df_res_all$gens_seq
# save object containing inputs and results as rda for analysis
save(df_res_all, file='df_res_all_rot.rda') #paste0(outFolder,'*.rda'))
}
} #end scenario loop
#to allow plotting of chosen scenarios by re-running them
#initially didn't work for diff insecticides because the inputs are saved as strings
#now the strings get converted back to vectors below
scenarios_to_plot <- c(21,23,32) #big difference ones all have insecticide1 as the good one
scenarios_to_plot <- c(29,38,46)
scenarios_to_plot <- c(13)
for(scenario in scenarios_to_plot)
{
#cat(scenario,"\n")
linputs <- as.list( df_res_all[scenario, 1:15] )
#converting multiple insecticide inputs back from strings
#have to ignore the ones that aren't strings
#note the space after the comma in split=", "
linputs <- lapply(linputs, function(x) if (!grepl(",",x)) x else
as.numeric(unlist(strsplit(as.character(x),split=", "))))
#TESTING changing min_rwr_interval from 5
#if min_rwr_interval is set to 0 you get frequent switch backs to other insecticides for 1 generation
#but when it's set to 5 it pushes frequences above the resistance freq threshold
linputs$min_rwr_interval <- 5
#min_gens_switch_back has to be 0 to allow one insecticide to continue when all others dead
linputs$min_gens_switch_back = 0
# run one scenario (rotation)
dfres <- do.call(run_rot, linputs)
ggrot <- rot_plot_resistance(dfres, plot_refuge=FALSE, logy=TRUE, plot=FALSE)
# now run as sequence
linputs$rot_interval <- 0
dfres <- do.call(run_rot, linputs)
ggseq <- rot_plot_resistance(dfres, plot_refuge=FALSE, logy=TRUE, plot=FALSE)
grid.arrange(ggrot, ggseq, ncol=1)
}
ian-hastings/rotations documentation built on Dec. 14, 2020, 11:42 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(rotations) library(ggplot2) library(gridExtra) outFolder <- "C:\\Dropbox\\resistanceResults\\rotations\\"
10,000 runs comparing rotations and sequences.
#eval=FALSE just to temporarily hide this text # Starting question : is there any benefit to rotations when there are no costs and no migration to an untreated area ? # # Test 1. No refugia (i.e. Coverage=100%) or fitness costs. Would expect all policies to be equivalent. # # Test 2. Dominance of resistance and fitness costs both set to 0.5. # I suspect one of the advantages of rotations (and the presence of refugia) is that they maintain lower resistance allele frequencies which slows the spread of recessive resistance alleles. I suspect/guess/intuit that if we make them semi-dominate these effects should disappear and all policies should be the same. I’m probably wrong but it’s a nice test….. # #\pagebreak
#results='hide', #comment=NA removes ## before text output #start small & build it up #maybe create an equiv. of resistance::sensiAnPaperPart.r later # example for 1 run # linputs <- list(max_gen = 50, n_insecticides = 2) # do.call(run_rot, linputs) # dfinputs <- as.data.frame(linputs) # So one way of doing is to create an linputs list for each run & pass to run_rot() # so run_rot() only ever runs for 1 scenario # in contrast resistance::runModel2() accepted an input matrix & would run for multiple scenarios # not sure which is better # set random seed so that results can be reproduced #set.seed(1) set.seed(2) max_gen <- 500 plot_scenarios <- FALSE n_scenarios <- 10 #000 already_ran <- FALSE #already_ran <- TRUE # if already done runs, just load data for plots if (already_ran) { #C:/rsprojects/rotations/inst/documents/ #load("df_res_all_rot10000.rda") load("df_res_all_rot_mig1000.rda") }else { df_res_all <- data.frame() for(i in 1:n_scenarios) { #message should output to the R markdown console if (i%%200 == 1) message("scenario ",i," of ",n_scenarios," ",Sys.time()) # inputs (could also be saved to a dataframe first & extracted one row at a time) linputs <- list() ## constant inputs ## linputs$coverage <- 1 linputs$migration <- 0 #linputs$cost <- 0 #variable cost below linputs$max_gen <- max_gen #beware min_rwr_interval, it's to stop repeated small rotations linputs$min_rwr_interval <- 5 linputs$min_gens_switch_back <- 1 linputs$no_r_below_start <- FALSE #BEWARE this has been set to TRUE in all earlier runs linputs$no_r_below_mut <- FALSE #BEWARE trying out new mutation-selection balance option #set plot to F if using grid.arrange as in phaps-rotations-better.Rmd linputs$plot <- FALSE ## variable inputs ## linputs$cost <- runif(1, min=0, max=0.1) #linputs$cost <- 0 linputs$expo_hi <- runif(1, min=0.4, max=0.9) #min=0.1, max=0.9 linputs$male_expo_prop <- runif(1, min=0, max=1) linputs$eff <- runif(1, min=0.5, max=1) # min=0.3, max=1 linputs$rr <- runif(1, min=0.1, max=0.9) linputs$dom_sel <- runif(1, min=0, max=1) linputs$dom_cos <- runif(1, min=0, max=1) #TODO make this lognormal to get lower frequencies too linputs$start_freqs <- runif(1, min=0.005, max=0.1) #note not log yet linputs$n_insecticides <- sample( 2:5, 1 ) #beware set replace=T if more than 1 # dispersal (initially with cost set to 0) #linputs$coverage <- runif(1, min=0.1, max=0.9) #linputs$migration <- runif(1, min=0.1, max=0.9) # create null objects used later for plotting ggrot <- ggseq <- NULL # to compare sequence to rotation run all the scenarios above for both seq & rot for(rot_or_not in 0:1) { #to get rot_int of 0 every other run if (rot_or_not == 0) linputs$rot_interval <- 0 else #random rot interval integer linputs$rot_interval <- sample( 5:50, 1 ) #beware set replace=T if more than 1 ################## # run one scenario dfres <- do.call(run_rot, linputs) # summarise results per scenario # note this is also currently calc at end of run_rot() here is probably better res <- gens_under_thresh(dfres) # to plot results for checking if ( linputs$rot_interval == 0 ) { gens_seq <- as.numeric(res) } else { gens_rot <- as.numeric(res) } if (plot_scenarios) { # make plot ggres <- rot_plot_resistance(dfres, plot_refuge=FALSE, logy=FALSE, plot=FALSE) # to plot results for checking if ( linputs$rot_interval == 0 ) { ggseq <- ggres gens_seq <- as.numeric(res) } else { ggrot <- ggres gens_rot <- as.numeric(res) } } # TODO I may want to concat all dfres into one object (adda column for scenario num) to save for later analysis # e.g. if I wanted to calculate average resistance frequencies over time in addition to gens below threshold } #end rot_or_not loop # make an output row for each scenario with inputs, gens_rot, gens_seq # need to make sure I retain the rot_interval for rot # add scenario id column onto results df_in_out <- as.data.frame(c(linputs, gens_rot=gens_rot, gens_seq=gens_seq, id=i)) df_in_out$rot_minus_seq <- df_in_out$gens_rot - df_in_out$gens_seq # previously I used rbind this scenario inputs & results onto previous # not a good speed or memory way of doing #df_res_all <- rbind(df_res_all, df_in_out) # now replaced with this better way # set a dataframe for the outputs to the required size in scenario1 if (i==1) { df_res_all <- data.frame(matrix(NA, ncol=length(names(df_in_out)), nrow=n_scenarios)) names(df_res_all) <- names(df_in_out) } #putting results into one row of overall results dataframe df_res_all[i,] <- df_in_out if (plot_scenarios) { # plot input vals # extract out those between 0 & 1 dfi <- data.frame(t(as.data.frame(linputs))) names(dfi) <- 'value' dfi$inputs <- rownames(dfi) #take out columns not between 0 & 1 dfi <- filter(dfi, !inputs %in% c('max_gen','rot_interval','n_insecticides','plot','min_rwr_interval')) #also remove coverage & migration both set to 0 for now dfi <- filter(dfi, !inputs %in% c('coverage','migration')) ggins <- ggplot(dfi, aes(x=inputs,y=value)) + geom_point() + theme_minimal() + ylim(0,1) + theme(axis.text.x = element_text(angle = 20,hjust = 1, vjust = 1)) grid.arrange(ggrot, ggseq, ggins, ncol=1) cat('scenario',i,'sum generations deployed insecticide under thresh rotation:', gens_rot,' sequence:',gens_seq,"\n\n") #can I give the code to rerun this scenario cat(paste0("run_rot(n_insecticides=",linputs$n_insecticides,", ", "cost=",signif(linputs$cost,2),", ", "start_freqs=",signif(linputs$start_freqs,2),", ", "rot_interval=",linputs$rot_interval,", ", "eff=",signif(linputs$eff,2),",\n", "dom_sel=",signif(linputs$dom_sel,2),", ", "dom_cos=",signif(linputs$dom_cos,2),", ", "rr=",signif(linputs$rr,2),", ", "expo_hi=",signif(linputs$expo_hi,2),", ", "coverage=",signif(linputs$coverage,2),",\n", "migration=",signif(linputs$migration,2),", ", "max_gen=",linputs$max_gen,", ", "min_rwr_interval=",linputs$min_rwr_interval, ")\n") ) } #I could save partial outputs every 200 scenarios #just in case doesn't make to end if (i%%1000 == 0) { #df_res_all$rot_minus_seq <- df_res_all$gens_rot - df_res_all$gens_seq # save object containing inputs and results as rda for analysis save(df_res_all, file='df_res_all_rot.rda') #paste0(outFolder,'*.rda')) } } #end scenario loop }
proportion of scenarios not reaching resistance thresholds
# check proportion of scenarios not reaching thresholds by max_runs for both seq & rot nruns_fail <- df_res_all %>% filter(gens_seq == max_gen & gens_rot == max_gen) %>% summarise(n=n()) %>% unlist() cat('proportion runs not reaching thresholds : ', nruns_fail/nrow(df_res_all) ) # remove runs where thresholds not reached for seq & rot df_res_all <- df_res_all %>% filter(gens_seq != max_gen | gens_rot != max_gen) nruns_same <- df_res_all %>% filter(gens_seq == gens_rot) %>% summarise(n=n()) %>% unlist() cat('proportion runs (that did reach threshold) where rotations and sequence the same : ', nruns_same/nrow(df_res_all) ) nruns_rot_20pc_better <- df_res_all %>% filter(rot_minus_seq/gens_seq >= 0.2 ) %>% summarise(n=n()) %>% unlist() cat('proportion runs (that did reach threshold) where rotations >20 % better than sequence : ', nruns_rot_20pc_better/nrow(df_res_all) )
difference between rotations and sequences across all runs
mostly (~67%) no difference but infrequent where rotations better
~ 7% of runs, rotation at least 20% better than sequence
print( ggplot(df_res_all, aes(x=rot_minus_seq)) + geom_histogram(binwidth=10) + xlab('rot_minus_seq : rotation generations minus sequence') + ylab('count (sqrt scale to show low values)') + scale_y_sqrt(breaks=c(0,10,20,50,200,500,2000,5000)) + #can't use log scale because anchored at 0 theme_minimal() + theme(panel.grid.minor.y=element_blank()))
\pagebreak
Examples of runs where rotations > 20% better than sequences
#eval false because smaller runs don't have these scenarios_to_plot <- c(624,2167,602,4028,5771) for(scenario in scenarios_to_plot) { cat(scenario,"\n") linputs <- as.list( df_res_all[scenario, 1:15] ) # run one scenario dfres <- do.call(run_rot, linputs) rot_plot_resistance(dfres, plot_refuge=FALSE, logy=FALSE, plot=TRUE) #now run as sequence linputs$rot_interval <- 0 dfres <- do.call(run_rot, linputs) rot_plot_resistance(dfres, plot_refuge=FALSE, logy=FALSE, plot=TRUE) dfi <- data.frame(t(as.data.frame(linputs))) names(dfi) <- 'value' dfi$inputs <- rownames(dfi) #take out columns not between 0 & 1 dfi <- filter(dfi, !inputs %in% c('max_gen','rot_interval','n_insecticides','plot','min_rwr_interval')) #also remove coverage & migration both set to 0 for now dfi <- filter(dfi, !inputs %in% c('coverage','migration')) ggins <- ggplot(dfi, aes(x=inputs,y=value)) + geom_point() + theme_minimal() + ylim(0,1) + theme(axis.text.x = element_text(angle = 20,hjust = 1, vjust = 1)) plot(ggins) }
How does the difference between rotations and sequences respond to model inputs ?
# go through selected inputs and plot them on x with rot_minus_seq on y # initially with 200 reps doesn't show anything obvious, try it with 2000 overnight runs in_cols <- c('cost','dom_cos','dom_sel','rr','expo_hi','male_expo_prop','eff','start_freqs') y <- 'rot_minus_seq' for(i in in_cols) { #plotnum <- plotnum + 1 #print( ggplot(ggInsOuts, aes_string(x=i, y=y, colour="strategy")) + #plotlist[[plotnum]] <- print( ggplot(df_res_all, aes_string(x=i, y=y)) + geom_point(shape=1, size=1, alpha=0.5, show.legend=FALSE) + theme_minimal() + geom_smooth(linetype='dashed',size=0.5) ) }
PRCC
library(sensitivity) # pcc(X, y, rank = FALSE, nboot = 0, conf = 0.95) # # Arguments # X a data frame containing model input variables # y a vector containing the responses # rank logical. If TRUE, the analysis is done on the ranks. # nboot the number of bootstrap replicates. # conf the confidence level of the bootstrap confidence intervals. #same as above in_cols <- c('cost','dom_cos','dom_sel','rr','expo_hi','male_expo_prop','eff','start_freqs') x <- df_res_all[,in_cols] y <- df_res_all['rot_minus_seq'] pcc_res <- pcc(x, y, rank=TRUE) #if you add bootstrap, error bars are added to the default plot & extra columns to the PRCC object #pcc_res <- pcc(x, y, rank=TRUE, nboot=100) #plot(pcc_res) #results are here I can probably rbind them together into a df that I can ggplot #pcc_res$PRCC to_plot <- pcc_res$PRCC #rename column 1 from 'original to PRCC names(to_plot)[1] <- 'PRCC' to_plot$inputs <- rownames(to_plot) axis_text_x <- element_text(angle = 45,hjust = 1, vjust = 1) height <- 0.9 gg <- ggplot( to_plot, aes_string(x='inputs',y='PRCC') ) + geom_point(shape=1, colour='red') + theme(axis.text.x = axis_text_x) + geom_hline(yintercept = 0, linetype=3) + #theme(panel.grid.major.y = element_line(color = "grey", size = 1)) + ylim(-1,1) + #ggtitle(paste(strat_name[strategy_num])) + xlab(NULL) plot(gg) #to put back in some ref gridlines that are removed by cowplot #theme_set(theme_gray()) #theme_set(theme_bw()) #plot_grid( plotlist[[1]],plotlist[[2]],plotlist[[3]],plotlist[[4]],ncol=1, rel_heights=c(1,1,1,1.8)) #plot_grid( plotlist[[1]],plotlist[[2]],plotlist[[3]],ncol=1, rel_heights=c(1,1,1.6), labels='AUTO')
Runs where rotations are >20% better (in red) are associated with higher dominance of cost and lower dominance of selection
but that does not guarantee, must be interaction with other inputs too.
#'cost','dom_cos','dom_sel','rr','expo_hi','male_expo_prop','eff','start_freqs' #'cost','dom_cos','dom_sel' y <- 'rot_minus_seq' ggplot(df_res_all, aes(x=dom_sel, y=dom_cos)) + geom_point(shape=3, size=1, alpha=0.2, show.legend=FALSE) + theme_minimal() + geom_point(data=filter(df_res_all, rot_minus_seq/gens_seq >= 0.3 ), shape=1, size=1, col='red', alpha=0.8) # looking at modelling it # library(mgcv) # mod_lm2 <- gam(rot_minus_seq ~ cost + dom_cos + dom_sel, data=df_res_all) # summary(mod_lm2) # # mod_gam2 <- gam(rot_minus_seq ~ s(cost) + s(dom_cos) + s(dom_sel) + s(rr) + s(expo_hi) + s(male_expo_prop) + s(eff) + s(start_freqs), data=df_res_all) # # summary(mod_gam2) # #Deviance explained = 19.5% # # mod_gam3 <- gam(rot_minus_seq ~ te(cost,dom_cos,dom_sel), data=df_res_all) # summary(mod_gam3) #actually decision trees may be useful
#eval=FALSE tmp because changes above will break this # summary plot of all runs, would show any differences between sequence & rotations as non horizontal lines # want to plot line between sequence and rotation value # any deviation from horizontal will be obvious # could use rotation interval as the x axis, so that longer intervals will be clearer # maybe not for publication # but useful for quick view #it was useful for me but now that I've moved to one row per scenario it won't work !! gg <- ggplot( df_res_all, aes_string(x='rot_interval', y='tot_gens_dep_under50', group='scenario', colour='scenario') ) + geom_line( alpha=0.5, lwd=1 ) plot(gg) # plot gens againts exposure and effectiveness # as a quick test to see if I can set the starting values of these higher gg <- ggplot( df_res_all, aes_string(x='eff', y='tot_gens_dep_under50') ) + geom_point( alpha=0.5 ) plot(gg) gg <- ggplot( df_res_all, aes_string(x='expo_hi', y='tot_gens_dep_under50') ) + geom_point( alpha=0.5 ) plot(gg) gg <- ggplot( df_res_all, aes_string(x='expo_hi*eff*male_expo_prop', y='tot_gens_dep_under50') ) + geom_point( alpha=0.5 ) plot(gg)
Have a go at running for insecticides with different parameters
starting to work ...
# set random seed so that results can be reproduced set.seed(2) max_gen <- 500 plot_scenarios <- FALSE n_scenarios <- 100 #1000 #000 #df_res_all <- data.frame() for(i in 1:n_scenarios) { #message should output to the R markdown console if (i%%200 == 1) message("scenario ",i," of ",n_scenarios," ",Sys.time()) # inputs (could also be saved to a dataframe first & extracted one row at a time) linputs <- list() ## constant inputs ## linputs$coverage <- 1 linputs$migration <- 0 #linputs$cost <- 0 #variable cost below linputs$max_gen <- max_gen #beware min_rwr_interval, it's to stop repeated small rotations linputs$min_rwr_interval <- 5 linputs$min_gens_switch_back <- 1 linputs$no_r_below_start <- FALSE #BEWARE this has been set to TRUE in all earlier runs linputs$no_r_below_mut <- TRUE #BEWARE trying out new mutation-selection balance option #set plot to F if using grid.arrange as in phaps-rotations-better.Rmd linputs$plot <- FALSE ## variable inputs ## linputs$n_insecticides <- sample( 2:5, 1 ) #beware set replace=T if more than 1 n_ins <- linputs$n_insecticides ## new those below can be different for each insecticide linputs$cost <- runif(n_ins, min=0, max=0.1) #linputs$cost <- 0 linputs$expo_hi <- runif(n_ins, min=0.4, max=0.9) #min=0.1, max=0.9 linputs$male_expo_prop <- runif(1, min=0, max=1) linputs$eff <- runif(n_ins, min=0.5, max=1) # min=0.3, max=1 linputs$rr <- runif(n_ins, min=0.1, max=0.9) linputs$dom_sel <- runif(n_ins, min=0, max=1) linputs$dom_cos <- runif(n_ins, min=0, max=1) #DIFFERENT linputs$start_freqs <- runif(n_ins, min=0.0001, max=0.01) #note not log yet # dispersal (initially with cost set to 0) #linputs$coverage <- runif(1, min=0.1, max=0.9) #linputs$migration <- runif(1, min=0.1, max=0.9) # create null objects used later for plotting ggrot <- ggseq <- NULL # to compare sequence to rotation run all the scenarios above for both seq & rot for(rot_or_not in 0:1) { #to get rot_int of 0 every other run if (rot_or_not == 0) linputs$rot_interval <- 0 else #random rot interval integer linputs$rot_interval <- sample( 5:50, 1 ) #beware set replace=T if more than 1 ################## # run one scenario dfres <- do.call(run_rot, linputs) # summarise results per scenario # note this is also currently calc at end of run_rot() here is probably better res <- gens_under_thresh(dfres) # to plot results for checking if ( linputs$rot_interval == 0 ) { gens_seq <- as.numeric(res) } else { gens_rot <- as.numeric(res) } if (plot_scenarios) { # make plot ggres <- rot_plot_resistance(dfres, plot_refuge=FALSE, logy=FALSE, plot=FALSE) # to plot results for checking if ( linputs$rot_interval == 0 ) { ggseq <- ggres gens_seq <- as.numeric(res) } else { ggrot <- ggres gens_rot <- as.numeric(res) } } } #end rot_or_not loop # make an output row for each scenario with inputs, gens_rot, gens_seq # need to make sure I retain the rot_interval for rot # to cope with multiple values per scenario convert them to strings linputs2 <- lapply(linputs, function(x) ifelse(length(x)>1,toString(x),x)) #added scenario id column df_in_out <- as.data.frame(c(linputs2, gens_rot=gens_rot, gens_seq=gens_seq, id=i), stringsAsFactors=FALSE) df_in_out$rot_minus_seq <- df_in_out$gens_rot - df_in_out$gens_seq #instead set up the dataframe to the required size in scenario1 if (i==1) { df_res_all <- data.frame(matrix(NA, ncol=length(names(df_in_out)), nrow=n_scenarios)) names(df_res_all) <- names(df_in_out) } #putting results into one row of overall results dataframe df_res_all[i,] <- df_in_out if (plot_scenarios) { # plot input vals # extract out those between 0 & 1 dfi <- data.frame(t(as.data.frame(linputs))) names(dfi) <- 'value' dfi$inputs <- rownames(dfi) #take out columns not between 0 & 1 dfi <- filter(dfi, !inputs %in% c('max_gen','rot_interval','n_insecticides','plot','min_rwr_interval')) #also remove coverage & migration both set to 0 for now dfi <- filter(dfi, !inputs %in% c('coverage','migration')) ggins <- ggplot(dfi, aes(x=inputs,y=value)) + geom_point() + theme_minimal() + ylim(0,1) + theme(axis.text.x = element_text(angle = 20,hjust = 1, vjust = 1)) grid.arrange(ggrot, ggseq, ggins, ncol=1) cat('scenario',i,'sum generations deployed insecticide under thresh rotation:', gens_rot,' sequence:',gens_seq,"\n\n") #can I give the code to rerun this scenario cat(paste0("run_rot(n_insecticides=",linputs$n_insecticides,", ", "cost=",signif(linputs$cost,2),", ", "start_freqs=",signif(linputs$start_freqs,2),", ", "rot_interval=",linputs$rot_interval,", ", "eff=",signif(linputs$eff,2),",\n", "dom_sel=",signif(linputs$dom_sel,2),", ", "dom_cos=",signif(linputs$dom_cos,2),", ", "rr=",signif(linputs$rr,2),", ", "expo_hi=",signif(linputs$expo_hi,2),", ", "coverage=",signif(linputs$coverage,2),",\n", "migration=",signif(linputs$migration,2),", ", "max_gen=",linputs$max_gen,", ", "min_rwr_interval=",linputs$min_rwr_interval, ")\n") ) } #save partial outputs just in case doesn't make to end if (i%%1000 == 0) { #df_res_all$rot_minus_seq <- df_res_all$gens_rot - df_res_all$gens_seq # save object containing inputs and results as rda for analysis save(df_res_all, file='df_res_all_rot.rda') #paste0(outFolder,'*.rda')) } } #end scenario loop #to allow plotting of chosen scenarios by re-running them #initially didn't work for diff insecticides because the inputs are saved as strings #now the strings get converted back to vectors below scenarios_to_plot <- c(21,23,32) #big difference ones all have insecticide1 as the good one scenarios_to_plot <- c(29,38,46) scenarios_to_plot <- c(13) for(scenario in scenarios_to_plot) { #cat(scenario,"\n") linputs <- as.list( df_res_all[scenario, 1:15] ) #converting multiple insecticide inputs back from strings #have to ignore the ones that aren't strings #note the space after the comma in split=", " linputs <- lapply(linputs, function(x) if (!grepl(",",x)) x else as.numeric(unlist(strsplit(as.character(x),split=", ")))) #TESTING changing min_rwr_interval from 5 #if min_rwr_interval is set to 0 you get frequent switch backs to other insecticides for 1 generation #but when it's set to 5 it pushes frequences above the resistance freq threshold linputs$min_rwr_interval <- 5 #min_gens_switch_back has to be 0 to allow one insecticide to continue when all others dead linputs$min_gens_switch_back = 0 # run one scenario (rotation) dfres <- do.call(run_rot, linputs) ggrot <- rot_plot_resistance(dfres, plot_refuge=FALSE, logy=TRUE, plot=FALSE) # now run as sequence linputs$rot_interval <- 0 dfres <- do.call(run_rot, linputs) ggseq <- rot_plot_resistance(dfres, plot_refuge=FALSE, logy=TRUE, plot=FALSE) grid.arrange(ggrot, ggseq, ncol=1) }
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.