In AndySouth/resistance: an insecticide resistance population genetics model with 2 loci and 2 insecticides
ALL FIGS FOR UNEQUAL SCENARIOS WHERE RESISTANCE FOR INSECTICIDE STARTS AT 0.5
initially just 500 test runs
ALL FIGS AT 50% resistance threshold
library(resistance)
outFolder <- "C:\\Dropbox\\resistanceResults\\"
## to load previously saved runs
## trying with the extended experiment
experiment <- 'extended'
# load(file=paste0(outFolder,'listOutMix_ex2_10000.rda'))
# load(file=paste0(outFolder,'listOutI1_ex2_10000.rda'))
# load(file=paste0(outFolder,'listOutI2_ex2_10000.rda'))
load(file=paste0(outFolder,'listOutMix_unequal.rda'))
load(file=paste0(outFolder,'listOutI1_unequal.rda'))
load(file=paste0(outFolder,'listOutI2_unequal.rda'))
# very quick test data
# load(file=paste0(outFolder,'listOutMix_ex2_3.rda'))
# load(file=paste0(outFolder,'listOutI1_ex2_3.rda'))
# load(file=paste0(outFolder,'listOutI2_ex2_3.rda'))
### calculate times to reach critical points and add them on to the input file
### for different insecticide strategies, sequential, mix1 and mix2
#1) sequential : time to resistance for each insecticide in isolation
#inputs : inAndOutI1, inAndOutI2
#find time to criticalPoint for insecticide1
#find time to criticalPoint for insecticide2
#add together
resistPointsI1 <- findResistancePoints(listOutI1, locus=1)
resistPointsI2 <- findResistancePoints(listOutI2, locus=2)
resistPointsSeq <- resistPointsI1 + resistPointsI2
#2) mixture1 : time to resistance for either insecticide when used in a mixture
#inputs : inAndOutMix
#find time to criticalPoint for EITHER insecticide in mixture
resistPointsMix_1 <- findResistancePoints(listOutMix, locus='either')
#todo - to be comparable I think this should be for when resistance to BOTH insecticides is reached
#3) mixture2 : when resistance to one insecticide in the mixture reached, switch to sole use of the
# other until that too reaches the critical point. Record total time.
# what I actually need to do is start with mixture find the first critical point
# (need to know which of the insecticides it is)
# then I need to go to the single run for the other insecticide & starting at
# it's current resistance point find out how many more generations to go
#inputs : inAndOutI1, inAndOutI2, inAndOutMix
resistPointsMix_A <- findResistancePointsMixResponsive(listOutMix, listOutI1, listOutI2)
#4) mixture3 : time to resistance for both insecticides when used in a mixture
#inputs : inAndOutMix
#find time to criticalPoint for BOTH insecticide in mixture
resistPointsMix_2 <- findResistancePoints(listOutMix, locus='both')
### chunk copied from part of one in sensiAnPaper1All
### bind results onto input file
treeInput <- listOutMix$input
#input files in listOutMix, listOutIn1 & listOutI2 are the same if the runs are done with default randomSeed
#except that exposure will be in a.f_AB, a.f_A0 and a.f_B0 respectively (& a.m* too)
#I could just rename one to exposure
#BEWARE risk if future changes
#1/2/16 don't need this now because I've saved exposure from the single original random value
#rownames(treeInput)[rownames(treeInput)=="a.f_AB"] <- "exposure"
#hardcode which variables to include in analysis to keep it simple and transparent
treePredictors <- c('P_1','P_2','exposure','phi.SS1_A0','phi.SS2_0B','h.RS1_A0','h.RS2_0B','s.RR1_A0','s.RR2_0B')
#add these for extended analysis
if (experiment=='extended')
treePredictors <- c(treePredictors,'male_exposure_prop','correct_mix_deploy')
treeInput <- treeInput[ treePredictors, ]
#add an extra predictor, the lower starting freq of resistance divided by the larger
resist_start_lo_div_hi <- ifelse( treeInput['P_1',] < treeInput['P_2',], treeInput['P_1',]/treeInput['P_2',], treeInput['P_2',]/treeInput['P_1',])
treeInput <- rbind(treeInput,resist_start_lo_div_hi)
#20160122 add test for Ian of resistance1/resistance2
resist_start_1_div_2 <- treeInput['P_1',]/treeInput['P_2',]
treeInput <- rbind(treeInput,resist_start_1_div_2)
#renaming other rownames to make nicer plots
rownames(treeInput)[rownames(treeInput)=="phi.SS1_A0"] <- "effectiveness_ins1"
rownames(treeInput)[rownames(treeInput)=="phi.SS2_0B"] <- "effectiveness_ins2"
rownames(treeInput)[rownames(treeInput)=="P_1"] <- "start_freq_allele1"
rownames(treeInput)[rownames(treeInput)=="P_2"] <- "start_freq_allele2"
rownames(treeInput)[rownames(treeInput)=="h.RS1_A0"] <- "dominance_allele1"
rownames(treeInput)[rownames(treeInput)=="h.RS2_0B"] <- "dominance_allele2"
rownames(treeInput)[rownames(treeInput)=="s.RR1_A0"] <- "selection_coef_allele1"
rownames(treeInput)[rownames(treeInput)=="s.RR2_0B"] <- "selection_coef_allele2"
#get the inputs here to use in ggplot investigation of model responses to inputs
#used in a later chunk
ggInput <- treeInput
#create a curtisInputs dataframe for use later
curtisInputs <- data.frame("exposure"=0.9,
"effectiveness_ins1"=0.73,
"effectiveness_ins2"=1,
"start_freq_allele1"=0.01,
"start_freq_allele2"=0.01,
"dominance_allele1"=0.17,
"dominance_allele2"=0.0016,
"selection_coef_allele1"=0.23,
"selection_coef_allele2"=0.43,
"male_exposure_prop"=1,
"correct_mix_deploy"=1,
"resist_start_lo_div_hi"=1,
"resist_start_1_div_2"=1)
# getting data into format for ggplot
#uses ggInput calculated in earlier chunk
#first needs to transpose rows to cols
#ggInput_T <- t(ggInput)
#function to transpose and add strategy column with the passed value
#also add on the input columns
addStrategyColumn <- function(x, value, inputs){
x <- as.data.frame( t(x) )
x$strategy <- value
#inputs <- as.numeric(inputs)
#transpose inputs
inputs <- as.data.frame( t(inputs), stringsAsFactors=FALSE )
#cbind onto the outputs
x <- cbind(inputs,x)
x
}
resistPointsI1_T <- addStrategyColumn(resistPointsI1,"insecticide 1",ggInput)
resistPointsI2_T <- addStrategyColumn(resistPointsI2,"insecticide 2",ggInput)
resistPointsSeq_T <- addStrategyColumn(resistPointsSeq,"Sequential",ggInput)
resistPointsMix_1_T <- addStrategyColumn(resistPointsMix_1,"Mix either",ggInput)
resistPointsMix_A_T <- addStrategyColumn(resistPointsMix_A,"Mix adaptive",ggInput)
resistPointsMix_2_T <- addStrategyColumn(resistPointsMix_2,"Mix both",ggInput)
ggInsOuts <- rbind( resistPointsI1_T, resistPointsI2_T, resistPointsSeq_T,
resistPointsMix_1_T, resistPointsMix_A_T, resistPointsMix_2_T)
#remove runs that didn't reach a resistance threshold (999), even if just for one strategy
#>1000 excludes sequential strategy that had a 999 in
#BEWARE would need to increase 1000 if I increase max generations in the runs
didntReachThresh <- which(ggInsOuts$gen_cP0.5 == 999 | ggInsOuts$gen_cP0.5 > 500 |
ggInsOuts$gen_cP0.25 == 999 | ggInsOuts$gen_cP0.25 > 500 |
ggInsOuts$gen_cP0.1 == 999 | ggInsOuts$gen_cP0.1 > 500 )
#subset by runs not making threshold
ggInsOuts <- ggInsOuts[-didntReachThresh,]
#prettify output names
names(ggInsOuts)[names(ggInsOuts)=="gen_cP0.5"] <- "time_to_resistance0.5"
names(ggInsOuts)[names(ggInsOuts)=="gen_cP0.25"] <- "time_to_resistance0.25"
names(ggInsOuts)[names(ggInsOuts)=="gen_cP0.1"] <- "time_to_resistance0.1"
#doing; for all inputs
if (experiment=='extended') {
num_inputs <- 13
} else {
num_inputs <- 11
}
#format data in a different way to enable PRCC on difference between sequential & mix2
#resistPointsMix_A_T
#resistPointsSeq_T
#columns to remove from first
indices1 <- which(names(resistPointsMix_A_T) %in% c("gen_cP0.1","gen_cP0.25","strategy"))
#columns to add to 2nd
indices2 <- which(names(resistPointsMix_A_T) %in% c("gen_cP0.5"))
#rename column in 2nd
tmp <- resistPointsSeq_T[indices2]
names(tmp) <- "gen_cP0.5seq"
dif_mixA_seq <- cbind(resistPointsMix_A_T[-indices1], tmp)
dif_mix2_seq <- cbind(resistPointsMix_2_T[-indices1], tmp)
#remove runs that didn't reach a resistance threshold (999), even if just for one strategy
#>1000 excludes sequential strategy that had a 999 in
#BEWARE would need to increase 1000 if I increase max generations in the runs
didntReachThresh <- which( dif_mixA_seq$gen_cP0.5 > 500 | dif_mixA_seq$gen_cP0.5seq > 500 )
#subset by runs not making threshold
dif_mixA_seq <- dif_mixA_seq[-didntReachThresh,]
dif_mixA_seq["mixA_minus_seq0.5"] <- dif_mixA_seq["gen_cP0.5"] - dif_mixA_seq["gen_cP0.5seq"]
#for dif_mix2_seq
didntReachThresh <- which( dif_mix2_seq$gen_cP0.5 > 500 | dif_mix2_seq$gen_cP0.5seq > 500 )
#subset by runs not making threshold
dif_mix2_seq <- dif_mix2_seq[-didntReachThresh,]
dif_mix2_seq["mix2_minus_seq0.5"] <- dif_mix2_seq["gen_cP0.5"] - dif_mix2_seq["gen_cP0.5seq"]
\pagebreak
Fig. x1 Time-to-resistance across all input values
- Time-to-resistance is lowest for single use and for the first insecticide in a mixture. Time-to-resistance is highest for the 2nd insecticide in a mixture.
- In between these extremes, there is less difference between the alternative strategies of a) sequential and b) adaptive mixture, switching to sole use. The latter takes slightly longer to reach the threshold.
library(ggplot2)
#names_results <- c("time_to_resistance0.5","time_to_resistance0.25","time_to_resistance0.1")
names_results <- c("time_to_resistance0.5")
#ggSubset <- ggInsOuts[ !ggInsOuts$strategy %in% c("insecticide 2"), ]
#4unequal also remove Seq 1st & Mix 1st which are meaningless
#unique(ggInsOuts$strategy)
#"insecticide 1" "insecticide 2" "Sequential" "Mix either" "Mix adaptive" "Mix both"
ggSubset <- ggInsOuts[ ggInsOuts$strategy %in% c("Sequential", "Mix adaptive", "Mix both"), ]
for(i in names_results)
{
#x11()
print( ggplot(ggSubset, aes_string(x='strategy',y=i, color='strategy')) +
#ylim(0,450) +
coord_cartesian( ylim=c(0, 450)) +
geom_violin(draw_quantiles = c(0.25, 0.75), show.legend=FALSE, adjust = .5 ) +
ylab('Time to resistance (gens)') +
#4unequal removed 2 strategies
scale_x_discrete("", limits=c("Sequential","Mix adaptive","Mix both"), labels = c("Sequential"="Seq both", "Mix adaptive"="Adaptive", "Mix both"="Mix both" ) )
) #end print
}
\pagebreak
Fig x2 PRCC analysis for time-to-resistance 0.5 for sequential and mixture strategies
- Exposure and dominance have the greatest on time-to-resistance as shown in PRCC analyses.
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.
#strategies <- unique(ggInsOuts$strategy)
#in main paper I want to just display sequential and mixture strategies
#strategies <- strategies[-1:-2]
#4unequal restrict to these strategies
strategies <- c("Sequential", "Mix adaptive", "Mix both")
plotlist <- list(length(strategies))
for(strategy_num in 1:length(strategies))
{
strategy <- strategies[strategy_num]
by_strategy <- ggInsOuts[ggInsOuts$strategy==strategy,]
x <- by_strategy[,1:num_inputs]
y <- by_strategy['time_to_resistance0.5']
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)
#just labelling axis for final plot
if (strategy_num == length(strategies))
{
axis_text_x <- element_text(angle = 45,hjust = 1, vjust = 1)
height <- 0.9
} else
{
axis_text_x <- element_blank()
height = 0.5
}
library(cowplot) # to enable setting relative heights of plots
plotlist[[strategy_num]] <- 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_bw() +
#theme(panel.grid.major.y = element_line(color = "grey", size = 1)) +
ylim(-1,1) +
ggtitle(paste(strategy)) +
xlab(NULL)
}
#plot_grid( plotlist[[1]],plotlist[[2]],plotlist[[3]],plotlist[[4]],ncol=1, rel_heights=c(1,1,1,1.8))
#4unequal only 3 strategies now
plot_grid( plotlist[[1]],plotlist[[2]],plotlist[[3]],ncol=1, rel_heights=c(1,1,1.8))
\pagebreak
Fig. x3 How time-to-resistance is effected by input values for each insecticide use strategy
- Effectiveness of the insecticides is the one input that can change the average relative positioning of the strategies a) sequential and b) adaptive mixture. Below an effectiveness of ~0.55 to 0.6 sequential takes longer to reach the resistance threshold, and above this the adaptive mixture takes longer.
#remove insecticide1 & 2 strategies
#ggSubset <- ggInsOuts[ !ggInsOuts$strategy %in% c("insecticide 1","insecticide 2"), ]
#4unequal
ggSubset <- ggInsOuts[ ggInsOuts$strategy %in% c("Sequential", "Mix adaptive", "Mix both"), ]
#in_cols <- c("exposure","effectiveness_ins1","dominance_allele1")
#4unequal
in_cols <- c("exposure","dominance_allele1","selection_coef_allele1","effectiveness_ins1")
out_cols <- c("time_to_resistance0.5","strategy")
ggSubset <- ggSubset[, c(in_cols,out_cols) ]
#num_inputs <- length(in_cols)
#names_inputs <- in_cols
#to rearrange order of legend entries
#ggSubset$strategy <- factor(ggSubset$strategy, levels=c("Mix both", "Mix adaptive", "Sequential", "Mix either"))
#4unequal
ggSubset$strategy <- factor(ggSubset$strategy, levels=c("Mix both", "Mix adaptive", "Sequential"))
#TODO could I facet the plots below rather than using the loop
#probably need to put all inputs into a single column and then facet by that
#for(i in names_inputs)
for(i in in_cols)
{
#x11()
y <- 'time_to_resistance0.5'
#print( ggplot(ggInsOuts, aes_string(x=i, y=y, colour="strategy")) +
print( ggplot(ggSubset, aes_string(x=i, y=y, colour="strategy")) +
#points not wanted if 10000
#geom_point(shape=3, show.legend=FALSE) +
#geom_smooth(colour='red', linetype='dashed',size=0.5) +
geom_smooth(linetype='dashed',size=1) +
#facet_wrap( ~ strategy) +
#geom_smooth(aes_string(x=i, y=y, color=NULL)) )
#labs(title = i) +
#add line for Curtis inputs
geom_vline(data = curtisInputs, aes_string(xintercept = i), linetype='dotted', colour='red') +
#tryting to get the vline into the legend fails
#geom_vline(data = curtisInputs, aes_string(xintercept = i, colour='"Curtis Fig.2"' ), linetype='dotted', show.legend=TRUE) +
geom_text(data = curtisInputs, inherit.aes = FALSE, aes_string(x = i, y=210, label='"Curtis Fig.2"'), size=3, colour='red')
)
}
\pagebreak
Fig. x4 PRCC difference in time-to-resistance 0.5 between sequential and adaptive mixture strategies
- Exposure and effectiveness come out as the most important effect on the difference between the sequential and adpative mixture strategies.
library(sensitivity)
#mixture2
x <- dif_mix2_seq[,1:num_inputs]
y <- dif_mix2_seq["mix2_minus_seq0.5"]
pcc_res <- pcc(x, y, rank=TRUE)
#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)
print( ggplot( to_plot, aes_string(x='inputs',y='PRCC') ) +
geom_point(shape=1, colour='red') +
theme(axis.text.x = element_text(angle = 45,hjust = 1, vjust = 1)) +
geom_hline(yintercept = 0, linetype=3) +
#theme_bw() +
ylim(-1,1) +
xlab(NULL)
)
## repeat code for adaptive
#4unequal comment out adaptive stuff
# x <- dif_mixA_seq[,1:num_inputs]
# y <- dif_mixA_seq["mixA_minus_seq0.5"]
#
# pcc_res <- pcc(x, y, rank=TRUE)
#
# #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)
#
# print( ggplot( to_plot, aes_string(x='inputs',y='PRCC') ) +
# geom_point(shape=1, colour='red') +
# theme(axis.text.x = element_text(angle = 45,hjust = 1, vjust = 1)) +
# geom_hline(yintercept = 0, linetype=3) +
# #theme_bw() +
# ylim(-1,1) +
# xlab(NULL)
# )
\pagebreak
Fig. x5 with unequal starting resistance this plot works less well, because resistance is always slower for mixture versuis sequential (and maybe other inputs have greater influence)
(previous text : Effectiveness and exposure together determine most of whether resistance arises slower under mixture or sequential strategies with unequal starting resistance)
- Exposure and effectiveness can be used to predict most of whether resistance will arise slower under the sequential or adaptive mixture strategy, either by a simple plot or using classification trees.
#4unequal change mixA to mix2 (mixA not relevant for unequal scenarios)
#in the test 500 runs resistance was always slower in mixtures vs sequential
#although not always >20% slower
ggSubset <- dif_mix2_seq
#if mix-seq >1 resistance took longer for mix so mix is better
ggSubset$mix2_better_seq0.5 <- ifelse( ggSubset$mix2_minus_seq0.5 > 0, 1, 0 )
ggSubset$seq_better_mix20.5 <- ifelse( ggSubset$mix2_minus_seq0.5 < 0, 1, 0 )
ggSubset$abs_diff <- abs(ggSubset$mix2_minus_seq0.5)
#4unequal commented out code below
# #tricky
# #be a bit careful what to do with zeros where results are same for both strategies
# #I use these better vars to facet the plots
# #may want the zeroes in both facets
# #to get that I can rbind the zeros from one onto the other BEWARE this may screw up existing plots
# df_same <- ggSubset[which(ggSubset$abs_diff==0),]
# #set these to 1 (they were 0)
# df_same$mix2_better_seq0.5 <- 1
# df_same$seq_better_mix20.5 <- 1
# #rbind them back onto other dataset
# ggSubset <- rbind( ggSubset, df_same)
#add an option for 20% better, put the not in to get facets in correct order in ggplot later
ggSubset$mix2_not20pc_better_seq0.5 <- ifelse( ggSubset$mix2_minus_seq0.5 / ggSubset$gen_cP0.5seq <= 0.2, 1, 0 )
#COOLIO !! I really like this figure.
library(viridis)
#facet showing generation differences, labelled mix better, seq better
#+ve & -ve on the same colour scale
#failed attempt at getting Curtis point in legend
#print( ggplot(ggSubset, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2, colour=abs_diff, size="Curtis"))+
print( ggplot(ggSubset, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2, colour=abs_diff))+
facet_grid( ~ seq_better_mix20.5, labeller=as_labeller( c('0'='mixture\nresistance slower','1'='sequential\nresistance slower'))) +
#YlGnBu looks really good except that the yellow points are very faint
#scale_colour_distiller(palette='YlGnBu', direction=1) +
#option="plasma" looks good too
scale_color_viridis(begin=1, end=0, option="viridis",
guide=guide_colourbar(title="difference in \ngenerations to reach\nresistance threshold", direction='horizontal')) +
theme_bw() +
theme(legend.position = "bottom") +
#theme_minimal() +
#add point for Curtis inputs
geom_point(data = curtisInputs, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2), colour='black', shape=3) +
#failing to add legend for curtis point
#scale_fill_manual(name = "Curtis Fig.2", values = 'red', labels = 'test') +
#scale_size_manual(name = "Curtis Fig.2", values="test") +
#scale_shape_manual(name = "Curtis Fig.2", values="test") +
ylab('effectiveness of insecticides added') +
xlab('exposure to both insecticides') +
#main points for the plot
geom_point(alpha=1, size=0.5) +
#4unequal copmment out below
#NOW can I add a lm fitted to the 0 points abs_diff
# geom_smooth(data=ggSubset[which(ggSubset$abs_diff==0),], method = "lm", se = FALSE, colour='red', lwd=0.5) +
# #can I subtley highlight the zero points which are repeated in each plot ?
# geom_point(data=ggSubset[which(ggSubset$abs_diff==0),], size=0.8, colour='red', shape=2, fill='white') +
#add text for Curtis point
geom_label(data = curtisInputs, inherit.aes = FALSE, aes_string(x = 'exposure', y='effectiveness_ins1 + effectiveness_ins2', label='"Curtis Fig.2"', hjust='"right"'), nudge_x=-0.02, size=3, colour='black')
)
# doing the lm outside of ggplot
# x1=df_same$exposure
# y1=df_same$effectiveness_ins1 + df_same$effectiveness_ins2
#
# l_mod <- lm(y1 ~ x1)
# l_mod$coefficients
\pagebreak
Fig. x6 Modifying the previous plot to show where resistance arises over 20% slower in a mixture strategy versus sequential.
- The line fitted to runs where resistance is reached at the same time for mixture and sequential is retained. Less of the parameter space favours the mixture under these constraints.
#4unequal changing from mixA to mix2
df_20pc <- ggSubset[which(ggSubset$mix2_minus_seq0.5 / ggSubset$gen_cP0.5seq > 0.19 &
ggSubset$mix2_minus_seq0.5 / ggSubset$gen_cP0.5seq < 0.21 ),]
#remove column it is faceted on to get line in both facets
col_indices <- colnames(df_20pc) != 'mix2_not20pc_better_seq0.5'
df_20pc <- df_20pc[,col_indices]
#modify the killer plot to show where mixtures are 20% better ?
print( ggplot(ggSubset, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2, colour=abs_diff))+
facet_grid( ~ mix2_not20pc_better_seq0.5, labeller=as_labeller( c('0'='mixture\nresistance >20% slower','1'='mixture\nresistance <20% slower'))) +
theme_bw() +
theme(legend.position = "bottom") +
scale_color_viridis(begin=1, end=0, option="viridis",
guide=guide_colourbar(title="difference in \ngenerations to reach\nresistance threshold", direction='horizontal')) +
#add point for Curtis inputs
#geom_point(data = curtisInputs, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2), colour='red', shape=3) +
ylab('effectiveness of insecticides added') +
xlab('exposure to both insecticides') +
#main points for the plot
geom_point(alpha=1, size=0.5)
#add a lm fitted to where exactly 20% diff
#geom_smooth(data=df_20pc, method = "lm", se = FALSE,colour='red', lwd=0.5)+
#can I subtley highlight the zero points which are repeated in each plot ?
#geom_point(data=df_20pc, size=0.8, colour='red', shape=2) #, fill='white')
)
# doing the lm outside of ggplot
# x1=df_20pc$exposure
# y1=df_20pc$effectiveness_ins1 + df_20pc$effectiveness_ins2
#
# l_mod <- lm(y1 ~ x1)
# l_mod$coefficients
# (Intercept) x1
# 1.3914547 0.2904911
Fig. x7 Classification Trees for whether resistance arises slower with a mixture or sequential strategy. The upper plot shows when mixture is slower than sequential, the lower plot is for when the mixture is more than 20% slower.
#code copied from sensiAnPaper1All
#4unequal changed mix_A to mix_2
require('plyr')
## is mix adaptive better than sequential ?
#T/F
resistBetterMix_2SeqBoolean <- resistPointsMix_2 > resistPointsSeq
#convert to 0/1
resistBetterMix_2Seq <- plyr::aaply(resistBetterMix_2SeqBoolean,.margins=c(1,2),.fun=as.numeric)
#rename rows ready for binding on to inputs
rownames(resistBetterMix_2Seq) <- gsub('gen','betterMix_2Seq', rownames(resistBetterMix_2Seq))
#20% better
resistBetterMix_2SeqBoolean20 <- resistPointsMix_2 > (1.2 * resistPointsSeq)
#convert to 0/1
resistBetterMix_2Seq20 <- plyr::aaply(resistBetterMix_2SeqBoolean20,.margins=c(1,2),.fun=as.numeric)
#rename rows ready for binding on to inputs
rownames(resistBetterMix_2Seq20) <- gsub('gen','betterMix_2Seq20', rownames(resistBetterMix_2Seq20))
#replace 0s & 1s with mixture/sequence to make tree plot clearer
repMix <- function(x){
x[x==0] <- 'sequence'
x[x==1] <- 'mixture'
x
}
resistBetterMix_2Seq <- repMix(resistBetterMix_2Seq)
resistBetterMix_2Seq20 <- repMix(resistBetterMix_2Seq20)
treeInput <- rbind( treeInput, resistBetterMix_2Seq, resistBetterMix_2Seq20 )
#transpose
treeInput <- t(treeInput)
#convert to a dataframe
treeInput <- data.frame(treeInput, stringsAsFactors = FALSE)
require(rpart.plot)
#convert predictor columns back to numeric
#following this above : replace 0s & 1s with mixture/sequence to make plot clearer
results_columns <- substr(names(treeInput),1,6) %in% 'better'
input_columns <- !results_columns
treeInput[,input_columns] <- lapply(treeInput[,input_columns],as.numeric)
treeInput[,results_columns] <- lapply(treeInput[,results_columns],as.factor)
#create string with predictor names in
treePredictorString <- paste(colnames(treeInput)[input_columns], collapse="+")
#just doing for a few strategies, sensiAnPpaer1All.Rmd does for a bunch
treeResponses <- c( rownames(resistBetterMix_2Seq)[2], #2 is 25%
rownames(resistBetterMix_2Seq20)[2]
)
#to do trees & plots for all response variables
for( treeResponse in treeResponses )
{
tree <- rpart::rpart(as.formula(paste(treeResponse,"~",treePredictorString)), data = treeInput, method = 'class')
#http://stackoverflow.com/questions/29197213/what-is-the-difference-between-rel-error-and-x-error-in-a-rpart-decision-tree
#A rule of thumb is to choose the lowest level where the rel_error + xstd < xerror.
cpOptimal <- tree$cptable[ tree$cptable[,"rel error"] + tree$cptable[,"xstd"] < tree$cptable[,"xerror"],"CP"][1]
treePruned <- prune(tree, cp=cpOptimal)
#to set box colours
cols <- ifelse(treePruned$frame$yval == 1, "red", "green3")
#add fig label
if (treeResponse == treeResponses[1]) main <- "A)"
else main <- "B)"
#adds in no. scenarios correctly and incorrectly classified
#prp(tree, extra=1, varlen = 10, main=treeResponse)
#prp(treePruned, extra=1, varlen = 10, main=treeResponse, col=cols, border.col=cols)
#prp(treePruned, extra=1, varlen = 0, main=treeResponse, box.col=cols, under=TRUE)
prp(treePruned, extra=1, varlen = 0, main=main, box.col=cols, under=TRUE)
}
Fig. x8 Introducing a new insecticide when there is already resistance to an existing one. This takes Curtis Fig 2 and sets the starting frequency for each allele in turn to a higher level to represent an 'old' insecticide to which there is resistance. In both cases resistance to the new insecticide develops slower when the 'old' insecticide continues to be used (mixture) than when it doesn't (sequence).
#curtis_f2 with P_1 changed from 0.01 to 0.45 : mix slower than seq
runcurtis_f2( P_1 = 0.49 , P_2 = 0.01 , h.RS1_A0 = 0.17 , h.RS2_0B = 0.0016 , exposure = 0.9 , phi.SS1_A0 = 0.73 , phi.SS2_0B = 1 , s.RR1_A0 = 0.23 , s.RR2_0B = 0.43, addCombinedStrategy=FALSE )
#curtis_f2 with P_2 changed from 0.01 to 0.45 : mix slower than seq
runcurtis_f2( P_1 = 0.01 , P_2 = 0.49 , h.RS1_A0 = 0.17 , h.RS2_0B = 0.0016 , exposure = 0.9 , phi.SS1_A0 = 0.73 , phi.SS2_0B = 1 , s.RR1_A0 = 0.23 , s.RR2_0B = 0.43, addCombinedStrategy=FALSE )
AndySouth/resistance documentation built on Nov. 12, 2020, 3:39 a.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.
ALL FIGS FOR UNEQUAL SCENARIOS WHERE RESISTANCE FOR INSECTICIDE STARTS AT 0.5
initially just 500 test runs
ALL FIGS AT 50% resistance threshold
library(resistance) outFolder <- "C:\\Dropbox\\resistanceResults\\" ## to load previously saved runs ## trying with the extended experiment experiment <- 'extended' # load(file=paste0(outFolder,'listOutMix_ex2_10000.rda')) # load(file=paste0(outFolder,'listOutI1_ex2_10000.rda')) # load(file=paste0(outFolder,'listOutI2_ex2_10000.rda')) load(file=paste0(outFolder,'listOutMix_unequal.rda')) load(file=paste0(outFolder,'listOutI1_unequal.rda')) load(file=paste0(outFolder,'listOutI2_unequal.rda')) # very quick test data # load(file=paste0(outFolder,'listOutMix_ex2_3.rda')) # load(file=paste0(outFolder,'listOutI1_ex2_3.rda')) # load(file=paste0(outFolder,'listOutI2_ex2_3.rda'))
### calculate times to reach critical points and add them on to the input file ### for different insecticide strategies, sequential, mix1 and mix2 #1) sequential : time to resistance for each insecticide in isolation #inputs : inAndOutI1, inAndOutI2 #find time to criticalPoint for insecticide1 #find time to criticalPoint for insecticide2 #add together resistPointsI1 <- findResistancePoints(listOutI1, locus=1) resistPointsI2 <- findResistancePoints(listOutI2, locus=2) resistPointsSeq <- resistPointsI1 + resistPointsI2 #2) mixture1 : time to resistance for either insecticide when used in a mixture #inputs : inAndOutMix #find time to criticalPoint for EITHER insecticide in mixture resistPointsMix_1 <- findResistancePoints(listOutMix, locus='either') #todo - to be comparable I think this should be for when resistance to BOTH insecticides is reached #3) mixture2 : when resistance to one insecticide in the mixture reached, switch to sole use of the # other until that too reaches the critical point. Record total time. # what I actually need to do is start with mixture find the first critical point # (need to know which of the insecticides it is) # then I need to go to the single run for the other insecticide & starting at # it's current resistance point find out how many more generations to go #inputs : inAndOutI1, inAndOutI2, inAndOutMix resistPointsMix_A <- findResistancePointsMixResponsive(listOutMix, listOutI1, listOutI2) #4) mixture3 : time to resistance for both insecticides when used in a mixture #inputs : inAndOutMix #find time to criticalPoint for BOTH insecticide in mixture resistPointsMix_2 <- findResistancePoints(listOutMix, locus='both')
### chunk copied from part of one in sensiAnPaper1All ### bind results onto input file treeInput <- listOutMix$input #input files in listOutMix, listOutIn1 & listOutI2 are the same if the runs are done with default randomSeed #except that exposure will be in a.f_AB, a.f_A0 and a.f_B0 respectively (& a.m* too) #I could just rename one to exposure #BEWARE risk if future changes #1/2/16 don't need this now because I've saved exposure from the single original random value #rownames(treeInput)[rownames(treeInput)=="a.f_AB"] <- "exposure" #hardcode which variables to include in analysis to keep it simple and transparent treePredictors <- c('P_1','P_2','exposure','phi.SS1_A0','phi.SS2_0B','h.RS1_A0','h.RS2_0B','s.RR1_A0','s.RR2_0B') #add these for extended analysis if (experiment=='extended') treePredictors <- c(treePredictors,'male_exposure_prop','correct_mix_deploy') treeInput <- treeInput[ treePredictors, ] #add an extra predictor, the lower starting freq of resistance divided by the larger resist_start_lo_div_hi <- ifelse( treeInput['P_1',] < treeInput['P_2',], treeInput['P_1',]/treeInput['P_2',], treeInput['P_2',]/treeInput['P_1',]) treeInput <- rbind(treeInput,resist_start_lo_div_hi) #20160122 add test for Ian of resistance1/resistance2 resist_start_1_div_2 <- treeInput['P_1',]/treeInput['P_2',] treeInput <- rbind(treeInput,resist_start_1_div_2) #renaming other rownames to make nicer plots rownames(treeInput)[rownames(treeInput)=="phi.SS1_A0"] <- "effectiveness_ins1" rownames(treeInput)[rownames(treeInput)=="phi.SS2_0B"] <- "effectiveness_ins2" rownames(treeInput)[rownames(treeInput)=="P_1"] <- "start_freq_allele1" rownames(treeInput)[rownames(treeInput)=="P_2"] <- "start_freq_allele2" rownames(treeInput)[rownames(treeInput)=="h.RS1_A0"] <- "dominance_allele1" rownames(treeInput)[rownames(treeInput)=="h.RS2_0B"] <- "dominance_allele2" rownames(treeInput)[rownames(treeInput)=="s.RR1_A0"] <- "selection_coef_allele1" rownames(treeInput)[rownames(treeInput)=="s.RR2_0B"] <- "selection_coef_allele2" #get the inputs here to use in ggplot investigation of model responses to inputs #used in a later chunk ggInput <- treeInput #create a curtisInputs dataframe for use later curtisInputs <- data.frame("exposure"=0.9, "effectiveness_ins1"=0.73, "effectiveness_ins2"=1, "start_freq_allele1"=0.01, "start_freq_allele2"=0.01, "dominance_allele1"=0.17, "dominance_allele2"=0.0016, "selection_coef_allele1"=0.23, "selection_coef_allele2"=0.43, "male_exposure_prop"=1, "correct_mix_deploy"=1, "resist_start_lo_div_hi"=1, "resist_start_1_div_2"=1)
# getting data into format for ggplot #uses ggInput calculated in earlier chunk #first needs to transpose rows to cols #ggInput_T <- t(ggInput) #function to transpose and add strategy column with the passed value #also add on the input columns addStrategyColumn <- function(x, value, inputs){ x <- as.data.frame( t(x) ) x$strategy <- value #inputs <- as.numeric(inputs) #transpose inputs inputs <- as.data.frame( t(inputs), stringsAsFactors=FALSE ) #cbind onto the outputs x <- cbind(inputs,x) x } resistPointsI1_T <- addStrategyColumn(resistPointsI1,"insecticide 1",ggInput) resistPointsI2_T <- addStrategyColumn(resistPointsI2,"insecticide 2",ggInput) resistPointsSeq_T <- addStrategyColumn(resistPointsSeq,"Sequential",ggInput) resistPointsMix_1_T <- addStrategyColumn(resistPointsMix_1,"Mix either",ggInput) resistPointsMix_A_T <- addStrategyColumn(resistPointsMix_A,"Mix adaptive",ggInput) resistPointsMix_2_T <- addStrategyColumn(resistPointsMix_2,"Mix both",ggInput) ggInsOuts <- rbind( resistPointsI1_T, resistPointsI2_T, resistPointsSeq_T, resistPointsMix_1_T, resistPointsMix_A_T, resistPointsMix_2_T) #remove runs that didn't reach a resistance threshold (999), even if just for one strategy #>1000 excludes sequential strategy that had a 999 in #BEWARE would need to increase 1000 if I increase max generations in the runs didntReachThresh <- which(ggInsOuts$gen_cP0.5 == 999 | ggInsOuts$gen_cP0.5 > 500 | ggInsOuts$gen_cP0.25 == 999 | ggInsOuts$gen_cP0.25 > 500 | ggInsOuts$gen_cP0.1 == 999 | ggInsOuts$gen_cP0.1 > 500 ) #subset by runs not making threshold ggInsOuts <- ggInsOuts[-didntReachThresh,] #prettify output names names(ggInsOuts)[names(ggInsOuts)=="gen_cP0.5"] <- "time_to_resistance0.5" names(ggInsOuts)[names(ggInsOuts)=="gen_cP0.25"] <- "time_to_resistance0.25" names(ggInsOuts)[names(ggInsOuts)=="gen_cP0.1"] <- "time_to_resistance0.1" #doing; for all inputs if (experiment=='extended') { num_inputs <- 13 } else { num_inputs <- 11 } #format data in a different way to enable PRCC on difference between sequential & mix2 #resistPointsMix_A_T #resistPointsSeq_T #columns to remove from first indices1 <- which(names(resistPointsMix_A_T) %in% c("gen_cP0.1","gen_cP0.25","strategy")) #columns to add to 2nd indices2 <- which(names(resistPointsMix_A_T) %in% c("gen_cP0.5")) #rename column in 2nd tmp <- resistPointsSeq_T[indices2] names(tmp) <- "gen_cP0.5seq" dif_mixA_seq <- cbind(resistPointsMix_A_T[-indices1], tmp) dif_mix2_seq <- cbind(resistPointsMix_2_T[-indices1], tmp) #remove runs that didn't reach a resistance threshold (999), even if just for one strategy #>1000 excludes sequential strategy that had a 999 in #BEWARE would need to increase 1000 if I increase max generations in the runs didntReachThresh <- which( dif_mixA_seq$gen_cP0.5 > 500 | dif_mixA_seq$gen_cP0.5seq > 500 ) #subset by runs not making threshold dif_mixA_seq <- dif_mixA_seq[-didntReachThresh,] dif_mixA_seq["mixA_minus_seq0.5"] <- dif_mixA_seq["gen_cP0.5"] - dif_mixA_seq["gen_cP0.5seq"] #for dif_mix2_seq didntReachThresh <- which( dif_mix2_seq$gen_cP0.5 > 500 | dif_mix2_seq$gen_cP0.5seq > 500 ) #subset by runs not making threshold dif_mix2_seq <- dif_mix2_seq[-didntReachThresh,] dif_mix2_seq["mix2_minus_seq0.5"] <- dif_mix2_seq["gen_cP0.5"] - dif_mix2_seq["gen_cP0.5seq"]
\pagebreak
Fig. x1 Time-to-resistance across all input values
- Time-to-resistance is lowest for single use and for the first insecticide in a mixture. Time-to-resistance is highest for the 2nd insecticide in a mixture.
- In between these extremes, there is less difference between the alternative strategies of a) sequential and b) adaptive mixture, switching to sole use. The latter takes slightly longer to reach the threshold.
library(ggplot2) #names_results <- c("time_to_resistance0.5","time_to_resistance0.25","time_to_resistance0.1") names_results <- c("time_to_resistance0.5") #ggSubset <- ggInsOuts[ !ggInsOuts$strategy %in% c("insecticide 2"), ] #4unequal also remove Seq 1st & Mix 1st which are meaningless #unique(ggInsOuts$strategy) #"insecticide 1" "insecticide 2" "Sequential" "Mix either" "Mix adaptive" "Mix both" ggSubset <- ggInsOuts[ ggInsOuts$strategy %in% c("Sequential", "Mix adaptive", "Mix both"), ] for(i in names_results) { #x11() print( ggplot(ggSubset, aes_string(x='strategy',y=i, color='strategy')) + #ylim(0,450) + coord_cartesian( ylim=c(0, 450)) + geom_violin(draw_quantiles = c(0.25, 0.75), show.legend=FALSE, adjust = .5 ) + ylab('Time to resistance (gens)') + #4unequal removed 2 strategies scale_x_discrete("", limits=c("Sequential","Mix adaptive","Mix both"), labels = c("Sequential"="Seq both", "Mix adaptive"="Adaptive", "Mix both"="Mix both" ) ) ) #end print }
\pagebreak
Fig x2 PRCC analysis for time-to-resistance 0.5 for sequential and mixture strategies
- Exposure and dominance have the greatest on time-to-resistance as shown in PRCC analyses.
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. #strategies <- unique(ggInsOuts$strategy) #in main paper I want to just display sequential and mixture strategies #strategies <- strategies[-1:-2] #4unequal restrict to these strategies strategies <- c("Sequential", "Mix adaptive", "Mix both") plotlist <- list(length(strategies)) for(strategy_num in 1:length(strategies)) { strategy <- strategies[strategy_num] by_strategy <- ggInsOuts[ggInsOuts$strategy==strategy,] x <- by_strategy[,1:num_inputs] y <- by_strategy['time_to_resistance0.5'] 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) #just labelling axis for final plot if (strategy_num == length(strategies)) { axis_text_x <- element_text(angle = 45,hjust = 1, vjust = 1) height <- 0.9 } else { axis_text_x <- element_blank() height = 0.5 } library(cowplot) # to enable setting relative heights of plots plotlist[[strategy_num]] <- 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_bw() + #theme(panel.grid.major.y = element_line(color = "grey", size = 1)) + ylim(-1,1) + ggtitle(paste(strategy)) + xlab(NULL) } #plot_grid( plotlist[[1]],plotlist[[2]],plotlist[[3]],plotlist[[4]],ncol=1, rel_heights=c(1,1,1,1.8)) #4unequal only 3 strategies now plot_grid( plotlist[[1]],plotlist[[2]],plotlist[[3]],ncol=1, rel_heights=c(1,1,1.8))
\pagebreak
Fig. x3 How time-to-resistance is effected by input values for each insecticide use strategy
- Effectiveness of the insecticides is the one input that can change the average relative positioning of the strategies a) sequential and b) adaptive mixture. Below an effectiveness of ~0.55 to 0.6 sequential takes longer to reach the resistance threshold, and above this the adaptive mixture takes longer.
#remove insecticide1 & 2 strategies #ggSubset <- ggInsOuts[ !ggInsOuts$strategy %in% c("insecticide 1","insecticide 2"), ] #4unequal ggSubset <- ggInsOuts[ ggInsOuts$strategy %in% c("Sequential", "Mix adaptive", "Mix both"), ] #in_cols <- c("exposure","effectiveness_ins1","dominance_allele1") #4unequal in_cols <- c("exposure","dominance_allele1","selection_coef_allele1","effectiveness_ins1") out_cols <- c("time_to_resistance0.5","strategy") ggSubset <- ggSubset[, c(in_cols,out_cols) ] #num_inputs <- length(in_cols) #names_inputs <- in_cols #to rearrange order of legend entries #ggSubset$strategy <- factor(ggSubset$strategy, levels=c("Mix both", "Mix adaptive", "Sequential", "Mix either")) #4unequal ggSubset$strategy <- factor(ggSubset$strategy, levels=c("Mix both", "Mix adaptive", "Sequential")) #TODO could I facet the plots below rather than using the loop #probably need to put all inputs into a single column and then facet by that #for(i in names_inputs) for(i in in_cols) { #x11() y <- 'time_to_resistance0.5' #print( ggplot(ggInsOuts, aes_string(x=i, y=y, colour="strategy")) + print( ggplot(ggSubset, aes_string(x=i, y=y, colour="strategy")) + #points not wanted if 10000 #geom_point(shape=3, show.legend=FALSE) + #geom_smooth(colour='red', linetype='dashed',size=0.5) + geom_smooth(linetype='dashed',size=1) + #facet_wrap( ~ strategy) + #geom_smooth(aes_string(x=i, y=y, color=NULL)) ) #labs(title = i) + #add line for Curtis inputs geom_vline(data = curtisInputs, aes_string(xintercept = i), linetype='dotted', colour='red') + #tryting to get the vline into the legend fails #geom_vline(data = curtisInputs, aes_string(xintercept = i, colour='"Curtis Fig.2"' ), linetype='dotted', show.legend=TRUE) + geom_text(data = curtisInputs, inherit.aes = FALSE, aes_string(x = i, y=210, label='"Curtis Fig.2"'), size=3, colour='red') ) }
\pagebreak
Fig. x4 PRCC difference in time-to-resistance 0.5 between sequential and adaptive mixture strategies
- Exposure and effectiveness come out as the most important effect on the difference between the sequential and adpative mixture strategies.
library(sensitivity) #mixture2 x <- dif_mix2_seq[,1:num_inputs] y <- dif_mix2_seq["mix2_minus_seq0.5"] pcc_res <- pcc(x, y, rank=TRUE) #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) print( ggplot( to_plot, aes_string(x='inputs',y='PRCC') ) + geom_point(shape=1, colour='red') + theme(axis.text.x = element_text(angle = 45,hjust = 1, vjust = 1)) + geom_hline(yintercept = 0, linetype=3) + #theme_bw() + ylim(-1,1) + xlab(NULL) ) ## repeat code for adaptive #4unequal comment out adaptive stuff # x <- dif_mixA_seq[,1:num_inputs] # y <- dif_mixA_seq["mixA_minus_seq0.5"] # # pcc_res <- pcc(x, y, rank=TRUE) # # #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) # # print( ggplot( to_plot, aes_string(x='inputs',y='PRCC') ) + # geom_point(shape=1, colour='red') + # theme(axis.text.x = element_text(angle = 45,hjust = 1, vjust = 1)) + # geom_hline(yintercept = 0, linetype=3) + # #theme_bw() + # ylim(-1,1) + # xlab(NULL) # )
\pagebreak
Fig. x5 with unequal starting resistance this plot works less well, because resistance is always slower for mixture versuis sequential (and maybe other inputs have greater influence)
(previous text : Effectiveness and exposure together determine most of whether resistance arises slower under mixture or sequential strategies with unequal starting resistance)
- Exposure and effectiveness can be used to predict most of whether resistance will arise slower under the sequential or adaptive mixture strategy, either by a simple plot or using classification trees.
#4unequal change mixA to mix2 (mixA not relevant for unequal scenarios) #in the test 500 runs resistance was always slower in mixtures vs sequential #although not always >20% slower ggSubset <- dif_mix2_seq #if mix-seq >1 resistance took longer for mix so mix is better ggSubset$mix2_better_seq0.5 <- ifelse( ggSubset$mix2_minus_seq0.5 > 0, 1, 0 ) ggSubset$seq_better_mix20.5 <- ifelse( ggSubset$mix2_minus_seq0.5 < 0, 1, 0 ) ggSubset$abs_diff <- abs(ggSubset$mix2_minus_seq0.5) #4unequal commented out code below # #tricky # #be a bit careful what to do with zeros where results are same for both strategies # #I use these better vars to facet the plots # #may want the zeroes in both facets # #to get that I can rbind the zeros from one onto the other BEWARE this may screw up existing plots # df_same <- ggSubset[which(ggSubset$abs_diff==0),] # #set these to 1 (they were 0) # df_same$mix2_better_seq0.5 <- 1 # df_same$seq_better_mix20.5 <- 1 # #rbind them back onto other dataset # ggSubset <- rbind( ggSubset, df_same) #add an option for 20% better, put the not in to get facets in correct order in ggplot later ggSubset$mix2_not20pc_better_seq0.5 <- ifelse( ggSubset$mix2_minus_seq0.5 / ggSubset$gen_cP0.5seq <= 0.2, 1, 0 ) #COOLIO !! I really like this figure. library(viridis) #facet showing generation differences, labelled mix better, seq better #+ve & -ve on the same colour scale #failed attempt at getting Curtis point in legend #print( ggplot(ggSubset, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2, colour=abs_diff, size="Curtis"))+ print( ggplot(ggSubset, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2, colour=abs_diff))+ facet_grid( ~ seq_better_mix20.5, labeller=as_labeller( c('0'='mixture\nresistance slower','1'='sequential\nresistance slower'))) + #YlGnBu looks really good except that the yellow points are very faint #scale_colour_distiller(palette='YlGnBu', direction=1) + #option="plasma" looks good too scale_color_viridis(begin=1, end=0, option="viridis", guide=guide_colourbar(title="difference in \ngenerations to reach\nresistance threshold", direction='horizontal')) + theme_bw() + theme(legend.position = "bottom") + #theme_minimal() + #add point for Curtis inputs geom_point(data = curtisInputs, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2), colour='black', shape=3) + #failing to add legend for curtis point #scale_fill_manual(name = "Curtis Fig.2", values = 'red', labels = 'test') + #scale_size_manual(name = "Curtis Fig.2", values="test") + #scale_shape_manual(name = "Curtis Fig.2", values="test") + ylab('effectiveness of insecticides added') + xlab('exposure to both insecticides') + #main points for the plot geom_point(alpha=1, size=0.5) + #4unequal copmment out below #NOW can I add a lm fitted to the 0 points abs_diff # geom_smooth(data=ggSubset[which(ggSubset$abs_diff==0),], method = "lm", se = FALSE, colour='red', lwd=0.5) + # #can I subtley highlight the zero points which are repeated in each plot ? # geom_point(data=ggSubset[which(ggSubset$abs_diff==0),], size=0.8, colour='red', shape=2, fill='white') + #add text for Curtis point geom_label(data = curtisInputs, inherit.aes = FALSE, aes_string(x = 'exposure', y='effectiveness_ins1 + effectiveness_ins2', label='"Curtis Fig.2"', hjust='"right"'), nudge_x=-0.02, size=3, colour='black') ) # doing the lm outside of ggplot # x1=df_same$exposure # y1=df_same$effectiveness_ins1 + df_same$effectiveness_ins2 # # l_mod <- lm(y1 ~ x1) # l_mod$coefficients
\pagebreak
Fig. x6 Modifying the previous plot to show where resistance arises over 20% slower in a mixture strategy versus sequential.
- The line fitted to runs where resistance is reached at the same time for mixture and sequential is retained. Less of the parameter space favours the mixture under these constraints.
#4unequal changing from mixA to mix2 df_20pc <- ggSubset[which(ggSubset$mix2_minus_seq0.5 / ggSubset$gen_cP0.5seq > 0.19 & ggSubset$mix2_minus_seq0.5 / ggSubset$gen_cP0.5seq < 0.21 ),] #remove column it is faceted on to get line in both facets col_indices <- colnames(df_20pc) != 'mix2_not20pc_better_seq0.5' df_20pc <- df_20pc[,col_indices] #modify the killer plot to show where mixtures are 20% better ? print( ggplot(ggSubset, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2, colour=abs_diff))+ facet_grid( ~ mix2_not20pc_better_seq0.5, labeller=as_labeller( c('0'='mixture\nresistance >20% slower','1'='mixture\nresistance <20% slower'))) + theme_bw() + theme(legend.position = "bottom") + scale_color_viridis(begin=1, end=0, option="viridis", guide=guide_colourbar(title="difference in \ngenerations to reach\nresistance threshold", direction='horizontal')) + #add point for Curtis inputs #geom_point(data = curtisInputs, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2), colour='red', shape=3) + ylab('effectiveness of insecticides added') + xlab('exposure to both insecticides') + #main points for the plot geom_point(alpha=1, size=0.5) #add a lm fitted to where exactly 20% diff #geom_smooth(data=df_20pc, method = "lm", se = FALSE,colour='red', lwd=0.5)+ #can I subtley highlight the zero points which are repeated in each plot ? #geom_point(data=df_20pc, size=0.8, colour='red', shape=2) #, fill='white') ) # doing the lm outside of ggplot # x1=df_20pc$exposure # y1=df_20pc$effectiveness_ins1 + df_20pc$effectiveness_ins2 # # l_mod <- lm(y1 ~ x1) # l_mod$coefficients # (Intercept) x1 # 1.3914547 0.2904911
Fig. x7 Classification Trees for whether resistance arises slower with a mixture or sequential strategy. The upper plot shows when mixture is slower than sequential, the lower plot is for when the mixture is more than 20% slower.
#code copied from sensiAnPaper1All #4unequal changed mix_A to mix_2 require('plyr') ## is mix adaptive better than sequential ? #T/F resistBetterMix_2SeqBoolean <- resistPointsMix_2 > resistPointsSeq #convert to 0/1 resistBetterMix_2Seq <- plyr::aaply(resistBetterMix_2SeqBoolean,.margins=c(1,2),.fun=as.numeric) #rename rows ready for binding on to inputs rownames(resistBetterMix_2Seq) <- gsub('gen','betterMix_2Seq', rownames(resistBetterMix_2Seq)) #20% better resistBetterMix_2SeqBoolean20 <- resistPointsMix_2 > (1.2 * resistPointsSeq) #convert to 0/1 resistBetterMix_2Seq20 <- plyr::aaply(resistBetterMix_2SeqBoolean20,.margins=c(1,2),.fun=as.numeric) #rename rows ready for binding on to inputs rownames(resistBetterMix_2Seq20) <- gsub('gen','betterMix_2Seq20', rownames(resistBetterMix_2Seq20)) #replace 0s & 1s with mixture/sequence to make tree plot clearer repMix <- function(x){ x[x==0] <- 'sequence' x[x==1] <- 'mixture' x } resistBetterMix_2Seq <- repMix(resistBetterMix_2Seq) resistBetterMix_2Seq20 <- repMix(resistBetterMix_2Seq20) treeInput <- rbind( treeInput, resistBetterMix_2Seq, resistBetterMix_2Seq20 ) #transpose treeInput <- t(treeInput) #convert to a dataframe treeInput <- data.frame(treeInput, stringsAsFactors = FALSE) require(rpart.plot) #convert predictor columns back to numeric #following this above : replace 0s & 1s with mixture/sequence to make plot clearer results_columns <- substr(names(treeInput),1,6) %in% 'better' input_columns <- !results_columns treeInput[,input_columns] <- lapply(treeInput[,input_columns],as.numeric) treeInput[,results_columns] <- lapply(treeInput[,results_columns],as.factor) #create string with predictor names in treePredictorString <- paste(colnames(treeInput)[input_columns], collapse="+") #just doing for a few strategies, sensiAnPpaer1All.Rmd does for a bunch treeResponses <- c( rownames(resistBetterMix_2Seq)[2], #2 is 25% rownames(resistBetterMix_2Seq20)[2] ) #to do trees & plots for all response variables for( treeResponse in treeResponses ) { tree <- rpart::rpart(as.formula(paste(treeResponse,"~",treePredictorString)), data = treeInput, method = 'class') #http://stackoverflow.com/questions/29197213/what-is-the-difference-between-rel-error-and-x-error-in-a-rpart-decision-tree #A rule of thumb is to choose the lowest level where the rel_error + xstd < xerror. cpOptimal <- tree$cptable[ tree$cptable[,"rel error"] + tree$cptable[,"xstd"] < tree$cptable[,"xerror"],"CP"][1] treePruned <- prune(tree, cp=cpOptimal) #to set box colours cols <- ifelse(treePruned$frame$yval == 1, "red", "green3") #add fig label if (treeResponse == treeResponses[1]) main <- "A)" else main <- "B)" #adds in no. scenarios correctly and incorrectly classified #prp(tree, extra=1, varlen = 10, main=treeResponse) #prp(treePruned, extra=1, varlen = 10, main=treeResponse, col=cols, border.col=cols) #prp(treePruned, extra=1, varlen = 0, main=treeResponse, box.col=cols, under=TRUE) prp(treePruned, extra=1, varlen = 0, main=main, box.col=cols, under=TRUE) }
Fig. x8 Introducing a new insecticide when there is already resistance to an existing one. This takes Curtis Fig 2 and sets the starting frequency for each allele in turn to a higher level to represent an 'old' insecticide to which there is resistance. In both cases resistance to the new insecticide develops slower when the 'old' insecticide continues to be used (mixture) than when it doesn't (sequence).
#curtis_f2 with P_1 changed from 0.01 to 0.45 : mix slower than seq runcurtis_f2( P_1 = 0.49 , P_2 = 0.01 , h.RS1_A0 = 0.17 , h.RS2_0B = 0.0016 , exposure = 0.9 , phi.SS1_A0 = 0.73 , phi.SS2_0B = 1 , s.RR1_A0 = 0.23 , s.RR2_0B = 0.43, addCombinedStrategy=FALSE ) #curtis_f2 with P_2 changed from 0.01 to 0.45 : mix slower than seq runcurtis_f2( P_1 = 0.01 , P_2 = 0.49 , h.RS1_A0 = 0.17 , h.RS2_0B = 0.0016 , exposure = 0.9 , phi.SS1_A0 = 0.73 , phi.SS2_0B = 1 , s.RR1_A0 = 0.23 , s.RR2_0B = 0.43, addCombinedStrategy=FALSE )
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.