In AndySouth/resistance: an insecticide resistance population genetics model with 2 loci and 2 insecticides

Figures for paper 1 : version 9 : produces Tiffs for submission (fails to produce pdf due to tiff options set in first chunk)

using rr_restoration instead of selection coefficient

library(knitr)

# options to create final publication quality figures
# causes the pdf document to fail, but the figs do get stored in pap1figs/
opts_chunk$set(dev="tiff",
               dev.args=list(compression="lzw"),
               dpi=300,
               cache=FALSE,
               fig.path='pap1figs/')
# opts_chunk$set(dev='postscript')

  library(resistance)

  outFolder <- "C:\\Dropbox\\resistanceResults\\"

  ## to load previously saved runs

  ## trying with the extended experiment
  experiment <- 'extended'

  ## to load previously saved runs  
  load(file=paste0(outFolder,'listOutMix_rr_10000.rda'))
  load(file=paste0(outFolder,'listOutI1_rr_10000.rda'))
  load(file=paste0(outFolder,'listOutI2_rr_10000.rda'))

  # load(file=paste0(outFolder,'listOutMix_ex2_10000.rda'))
  # load(file=paste0(outFolder,'listOutI1_ex2_10000.rda'))
  # load(file=paste0(outFolder,'listOutI2_ex2_10000.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'))

Fig1

#, fig.align='left'}

#PLOS fig dimensions
#Maximum width 7.5 or 5.2 to fit in text column.
#Maximum height 8.75

library(ggplot2)
library(cowplot)

#list to put plots in for cowplot
plotlist <- list(2)

#Fig 1A

dfd1 <- selection_simple( freq=0.01, selection_co = 1, dominance = 1, num_gens = 150, plot=FALSE)
dfd05 <- selection_simple( freq=0.01, selection_co = 1, dominance = 0.5, num_gens = 150, plot=FALSE)
dfd0 <- selection_simple( freq=0.01, selection_co = 1, dominance = 0, num_gens = 150, plot=FALSE)

df_to_plot <- rbind(dfd1, dfd05, dfd0)

#old way of getting in wide format  
# df_to_plot <- data.frame(generation=dfd1$generation,
#                          dominance1=dfd1$freq,
#                          dominance0.5=dfd05$freq,
#                          dominance0=dfd0$freq)

plotlist[[1]] <- ggplot(df_to_plot, aes(x=generation, y=freq, 
                              colour=factor(dominance,labels=c("0   recessive",
                                                               "0.5 semi-dominant",
                                                               "1   dominant")))) + 
       theme_bw() +  
       #theme(legend.position = "bottom", legend.key = element_blank()) + 
       #guides(colour = guide_legend(reverse=TRUE)) +  
       labs(colour = "dominance") +
       ylab("Resistance Allele Frequency") +   
       #coord_trans(y = "log10") +    
       geom_line() +
       theme(axis.text.x = element_blank(),
             axis.ticks = element_blank(),
             panel.grid = element_blank())  
       #annotate('text',label='A)',x=-10, y=1.1)

#Fig 1B

#22/11/16 adding extra section on left where all killed following request from Ian
#first try adding on left @-10 and modify x extent so I don't have to fiddle with all coords
df_labs <- data.frame(lab=rep(c('RR','RS','SS'),5),
                      x=c(-10,-10,-10,10,10,10,30,30,30,50,50,50,70,70,70),
                      y=seq(62,6,-4),
                      #y=c(50,46,42,38,34,30,26,22,18,14,10,6),
                      col=c(0,0,0,1,0,0,1,0.5,0,1,1,0,1,1,1))

df_boxes <- data.frame(lab=c('All killed\n(no spread)',
                             'R recessive\n(slow spread)',
                             'R semi-dominant\n(rapid spread)',
                             'R dominant\n(rapid spread)',
                             'All survive\n(no spread)'),
                      x=c(-10,10,30,50,70),
                      y=c(105,90,75,60,45))

# df_labs <- data.frame(lab=rep(c('RR','RS','SS'),4),
#                       x=c(10,10,10,30,30,30,50,50,50,70,70,70),
#                       y=c(50,46,42,38,34,30,26,22,18,14,10,6),
#                       #col=c(2,1,1,2,1.5,1,2,2,1,2,2,2))
#                       col=c(1,0,0,1,0.5,0,1,1,0,1,1,1))
# 
# df_boxes <- data.frame(lab=c('R recessive\n(slow spread)',
#                              'R semi-dominant\n(rapid spread)',
#                              'R dominant\n(rapid spread)',
#                              'No selection for R'),
#                       x=c(10,30,50,80),
#                       y=c(90,75,60,35))

plotlist[[2]] <- ggplot(data.frame(x=c(-20,85),y=c(90,15)), aes(x=x, y=y)) + 
       theme_bw() + 
       ylim(0,115) +
       xlim(-20,85) +
       coord_cartesian(expand = FALSE) + #so that axes fit data exactly
       ylab("Insecticide concentration") +
       xlab("Time since treatment") +         
       geom_line(size=2, colour=rgb(0.5,0,0.5)) +
       geom_vline(xintercept = c(0,20,40,60), linetype=2) +
       geom_text(data=df_labs,aes(x=x, y=y, label=lab, colour=factor(col)),size=4) +
       scale_colour_manual(values = c("0" = "red","0.5" = "orange","1" = "green"), guide=FALSE) +  
       geom_label(data=df_boxes,aes(x=x, y=y, label=lab),size=3) +
       theme(axis.text = element_blank(),
             axis.ticks = element_blank(),
             panel.grid = element_blank())
       #draw_plot_label('B)') +



plot_grid( plotlist[[1]],plotlist[[2]],ncol=1, rel_heights=c(0.7,1), labels='AUTO') 

#draw_figure_label('B)', pos='top.left')

Fig2 replicating Curtis, and now adding in copies of Curtis plots

  #Beths code for these figures taken from resistanceMaster

  #set layout for 4 plots, A&C Curtis, B&D new ones
  layout(matrix(1:4,2,2,byrow=FALSE))

  #add in images of curtis plots here
  #to get into the layout I need to do from R
  #may be ugly
  #par(mar=c(0,0,0,0)) #bltr default is c(5, 4, 4, 2)
  #pty="s" #to create a square plot region
  par(pty="s",mar=c(0.4,2.9,0,0.4)) #bltr default is c(5, 4, 4, 2)

  library(png)
  f1png <- readPNG(system.file("documents","pap1_curtis1.png", package="resistance"))
  #both the ld and the rfreq plots (ld on its own means little with the way axes are labelled)
  #f1png <- readPNG(system.file("documents","pap1_curtis1_both.png", package="resistance"))
  plot( 0, 0, type="n", axes=F, xlim=c(0,1), ylim=c(0,1), xlab="", ylab="")
  rasterImage(f1png, xleft=0, ybottom=0, xright=1, ytop=1, interpolate = FALSE)
  mtext('A)',side=3, adj=0, line=1)


  f2png <- readPNG(system.file("documents","pap1_curtis2.png", package="resistance"))
  plot( 0, 0, type="n", axes=F, xlim=c(0,1), ylim=c(0,1), xlab="", ylab="")
  rasterImage(f2png, xleft=0, ybottom=0, xright=1, ytop=1, interpolate = FALSE)
  mtext('C)',side=3, adj=0, line=1)



  par(mar=c(0.4,2.9,0,0.4)) #bltr default is c(5, 4, 4, 2)

  #2B
  calibration <- 1011
  input <- createInputMatrix( params.csv=FALSE, calibration=calibration )
  listOut <- runModel2( input=input )
  #we submitted the LD one but should have done the Fig1
  #plotcurtis_ld( listOut$results[[1]], listOut$results[[2]], 1, 4 )
  plotcurtis_f1( listOut$results[[1]], listOut$results[[2]], 1, 2 )
  #side=1=bottom, 2=left, 3=top, 4=right, adj=0 left justify, lines move outwards
  mtext('B)',side=3, adj=0, line=1)

  #2D
  calibration <- 1012
  input <- createInputMatrix( params.csv=FALSE, calibration=calibration )
  listOut <- runModel2( input=input )  
  plotcurtis_f2( listOut$results[[3]], listOut$results[[1]], listOut$results[[2]], 1, 2, 3 )
  #side=1=bottom, 2=left, 3=top, 4=right, adj=0 left justify, lines move outwards
  mtext('D)',side=3, adj=0, line=1)

  #Beths code for these figures taken from resistanceMaster

  #set layout for 2 plots
  layout(t(1:2))
  par(mar=c(0.4,2.9,0,0.4)) #bltr default is c(5, 4, 4, 2)

  #2A
  calibration <- 1011
  input <- createInputMatrix( params.csv=FALSE, calibration=calibration )
  listOut <- runModel2( input=input )
  #we submitted the LD one but should have done the Fig1
  #plotcurtis_ld( listOut$results[[1]], listOut$results[[2]], 1, 4 )
  plotcurtis_f1( listOut$results[[1]], listOut$results[[2]], 1, 2 )
  #side=1=bottom, 2=left, 3=top, 4=right, adj=0 left justify, lines move outwards
  mtext('A)',side=3, adj=0, line=1)

  #2B
  calibration <- 1012
  input <- createInputMatrix( params.csv=FALSE, calibration=calibration )
  listOut <- runModel2( input=input )  
  plotcurtis_f2( listOut$results[[3]], listOut$results[[1]], listOut$results[[2]], 1, 2, 3 )
  #side=1=bottom, 2=left, 3=top, 4=right, adj=0 left justify, lines move outwards
  mtext('B)',side=3, adj=0, line=1)

### 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)
  #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"

  #27/6/16 to convert any varnames in old output files
  rownames(treeInput)[ "maleExposureProp" == rownames(treeInput) ] <- "male_exposure_prop"
  rownames(treeInput)[ "correctMixDeployProp" == rownames(treeInput) ] <- "correct_mix_deploy"
  rownames(treeInput)[ "rr_advantage_I1" == rownames(treeInput) ] <- "rr_restoration_ins1"
  rownames(treeInput)[ "rr_advantage_I2" == rownames(treeInput) ] <- "rr_restoration_ins2"  

  #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')
  #14/6/16 change to rr_restoration
  treePredictors <- c('P_1','P_2','exposure','phi.SS1_A0','phi.SS2_0B','h.RS1_A0','h.RS2_0B','rr_restoration_ins1','rr_restoration_ins2')

  #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"

  # 14/6/16
  # 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,
                             "effectiveness_added"=1.73,                             
                             "start_freq_allele1"=0.01,
                             "start_freq_allele2"=0.01,
                             "dominance_allele1"=0.17,
                             "dominance_allele2"=0.0016,
                             "dominance_added"=0.1716,  
                             # 14/6/16 rr_restoration = s / effectiveness
                             #"selection_coef_allele1"=0.23,
                             #"selection_coef_allele2"=0.43,
                             "rr_restoration_ins1"=0.23/0.73,
                             "rr_restoration_ins2"=0.43/1, 
                             "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 )


  #22/6/2016 to check % of runs excluded 
  #its complicated given that there are runs for I1,I2 & mix 
  #this is an approximation
  #TODO, are any others excluded later ?
  percentExcluded <- 100*length(didntReachThresh)/nrow(ggInsOuts)
  #16.5

  #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 )

  #TODO 22/6/16 RECHECK that nothing dodgy going on here, seems that more runs being excluded
  #but this is because its going back to the raw files
  percentExcluded2 <- 100*length(didntReachThresh)/nrow(dif_mixA_seq)
  #24%

  #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 )

  percentExcluded3 <- 100*length(didntReachThresh)/nrow(dif_mix2_seq)
  #30%  

  ####################
  #23/6/2016 experimenting to see if can reduce the % excluded by tweaking input ranges
  #20%
  #dif_mix2_seq2 <- dif_mix2_seq[ dif_mix2_seq$start_freq_allele1 > 0.005, ]
  #16%
  #dif_mix2_seq2 <- dif_mix2_seq[ dif_mix2_seq$start_freq_allele1 > 0.005 & dif_mix2_seq$start_freq_allele2 > 0.005, ]
  #18%
  #dif_mix2_seq2 <- dif_mix2_seq[ dif_mix2_seq$exposure > 0.25, ]
  #16%
  #dif_mix2_seq2 <- dif_mix2_seq[ dif_mix2_seq$exposure > 0.3, ]
  #7%
  #dif_mix2_seq2 <- dif_mix2_seq[ dif_mix2_seq$start_freq_allele1 > 0.005 & dif_mix2_seq$start_freq_allele2 > 0.005 & dif_mix2_seq$exposure > 0.25, ]   
  #5%
  #dif_mix2_seq2 <- dif_mix2_seq[ dif_mix2_seq$start_freq_allele1 > 0.005 & dif_mix2_seq$start_freq_allele2 > 0.005 & dif_mix2_seq$exposure > 0.3, ] 
  #26%
  #dif_mix2_seq2 <- dif_mix2_seq[ dif_mix2_seq$effectiveness_ins1 > 0.4 & dif_mix2_seq$effectiveness_ins2 > 0.4, ]
  #13%
  #dif_mix2_seq2 <- dif_mix2_seq[ dif_mix2_seq$effectiveness_ins1 > 0.4 & dif_mix2_seq$effectiveness_ins2 > 0.4 & dif_mix2_seq$exposure > 0.3, ]
  #10%
  #dif_mix2_seq2 <- dif_mix2_seq[ dif_mix2_seq$start_freq_allele1 > 0.001 & dif_mix2_seq$start_freq_allele2 > 0.001 & dif_mix2_seq$exposure > 0.3, ] 
  #8%
  #dif_mix2_seq2 <- dif_mix2_seq[ dif_mix2_seq$effectiveness_ins1 > 0.4 & dif_mix2_seq$effectiveness_ins2 > 0.4 & dif_mix2_seq$start_freq_allele1 > 0.001 & dif_mix2_seq$start_freq_allele2 > 0.001 & dif_mix2_seq$exposure > 0.3, ] 
  #22%
  #dif_mix2_seq2 <- dif_mix2_seq[ dif_mix2_seq$dominance_allele1 > 0.1 & dif_mix2_seq$dominance_allele2 > 0.1, ] 
  #5%
  #dif_mix2_seq2 <- dif_mix2_seq[ dif_mix2_seq$start_freq_allele1 > 0.001 & dif_mix2_seq$start_freq_allele2 > 0.001 & dif_mix2_seq$exposure > 0.3 & dif_mix2_seq$dominance_allele1 > 0.1 & dif_mix2_seq$dominance_allele2 > 0.1, ] 

  # didntReachThresh <- which( dif_mix2_seq2$gen_cP0.5 > 500 | dif_mix2_seq2$gen_cP0.5seq > 500 )
  # 
  # percentExcluded4 <- 100*length(didntReachThresh)/nrow(dif_mix2_seq2)
  # 
  # percentExcluded4  

  #### end of expt 23/6/2016  

  #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 3

#get input values from Ian, then I can re-generate, might want to change mix2 to mix
#90% exposure, equal for female and male
#Starting freq. at locus A =0.01 and at locus B is 0.001
#Fitness =1 for RR genotypes at 0 for SR and SS at both loci

runcurtis_f2( P_1 = 0.01 , P_2 = 0.001, h.RS1_A0 = 1, h.RS2_0B = 1, exposure = 0.9,                        phi.SS1_A0 = 1 , phi.SS2_0B = 1 , s.RR1_A0 = 1 , s.RR2_0B = 1,
              addCombinedStrategy=FALSE, addStrategyLabels=FALSE )

Fig4 Time-to-resistance across all input values

1. 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.
2. 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"), ]  

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, 350)) +
      geom_violin(draw_quantiles = c(0.25, 0.75), show.legend=FALSE, adjust = .5 ) +
      ylab('Time to resistance (gens)') +
      scale_x_discrete("", limits=c("insecticide 1","Mix either","Sequential","Mix adaptive","Mix both"), labels = c("insecticide 1"="Seq 1st", "Mix either"="Mix 1st","Sequential"="Seq both", "Mix adaptive"="Adaptive both", "Mix both"="Mix both" ) )

  ) #end print
}  

\pagebreak

Fig5 PRCC analysis for time-to-resistance 0.5 for sequential and mixture strategies

1. 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 <- c("Sequential", "Mix adaptive", "Mix both")
  strat_name <- c("Sequential both","Adaptive both","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(panel.grid.major.y = element_line(color = "grey", size = 1)) +
             ylim(-1,1) +
             ggtitle(paste(strat_name[strategy_num])) +
             xlab(NULL) 
  }

  #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')    

\pagebreak

Fig6 How time-to-resistance is effected by input values for each insecticide use strategy

1. 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.

  #layout for 3 panel plot (t=horiz, c=vertical)
  #doesn't work for ggplot
  #layout(c(1:3))

  #list to put plots in for cowplot
  plotlist <- list(3)

  #remove some strategies for plots
  ggSubset <- ggInsOuts[ !ggInsOuts$strategy %in% c("insecticide 1","insecticide 2","Mix either"), ]

  # create added columns for effectiveness and dominance
  ggSubset$effectiveness_added <- ggSubset$effectiveness_ins1 + ggSubset$effectiveness_ins2
  ggSubset$dominance_added <- ggSubset$dominance_allele1 + ggSubset$dominance_allele2

  #subset which output plots to show in the main paper
  # in_cols <- c("start_freq_allele1", "exposure","effectiveness_ins1",
  #             "dominance_allele1","selection_coef_allele1","male_exposure_prop")
  #in_cols <- c("exposure","effectiveness_ins1","dominance_allele1")  
  in_cols <- c("exposure","effectiveness_added","dominance_added")  

  out_cols <- c("time_to_resistance0.5","strategy")
  ggSubset <- ggSubset[, c(in_cols,out_cols) ]

  #to rearrange order of legend entries
  ggSubset$strategy <- factor(ggSubset$strategy, levels=c("Mix both", "Mix adaptive", "Sequential", "Mix either")) 
  #20/6/16 try removing "Mix either" 
  #ggSubset$strategy <- factor(ggSubset$strategy, levels=c("Mix both", "Mix adaptive", "Sequential"))   

  #renaming levels just for fig labels
  levels(ggSubset$strategy) <- c("Mix both","Adaptive both","Sequential both","Mix either")

  #for(i in names_inputs)

  plotnum <- 0

  for(i in in_cols)
    {

    y <- 'time_to_resistance0.5'    

    plotnum <- plotnum + 1

    #print( ggplot(ggInsOuts, aes_string(x=i, y=y, colour="strategy")) + 
    plotlist[[plotnum]] <- ggplot(ggSubset, aes_string(x=i, y=y, colour="strategy")) + 

             #ideally would include points but they obscure means
             #geom_point(shape=1, size=1, alpha=0.2, show.legend=FALSE) + 
             geom_smooth(linetype='dashed',size=0.5) +
             #facet_wrap( ~ strategy) +
             #labs(title = i) +
             #add line for Curtis inputs
             geom_vline(data = curtisInputs, aes_string(xintercept = i), linetype='dotted', colour='red') +
             #trying 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_label(data = curtisInputs, inherit.aes=FALSE, aes_string(x=i, y=210, label='"Curtis\nFig.2"'), size=3, colour='red', nudge_x = -0.05)

  }

  #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]],ncol=1, rel_heights=c(1,1,1), labels='AUTO', vjust=1)   

\pagebreak

Fig7 PRCC difference in time-to-resistance 0.5 between sequential and adaptive mixture strategies

1. Exposure and effectiveness come out as the most important effect on the difference between A) sequential and mixture, and B) difference between sequential and adaptive mixture strategies.

  library(sensitivity)

  #list to put plots in for cowplot
  plotlist <- list(2)

    #mixture2
    x <- dif_mix2_seq[,1:num_inputs]
    y <- dif_mix2_seq["mix2_minus_seq0.5"]

    pcc_res <- pcc(x, y, rank=TRUE)


    to_plot <- pcc_res$PRCC
    #rename column 1 from 'original to PRCC
    names(to_plot)[1] <- 'PRCC'
    to_plot$inputs <- rownames(to_plot)  

    #cat("A)")

    plotlist[[1]] <- ggplot( to_plot, aes_string(x='inputs',y='PRCC') ) + 
             geom_point(shape=1, colour='red') +
             theme(axis.text.x = element_blank()) +
             geom_hline(yintercept = 0, linetype=3) +
             ggtitle(paste("A) Mixture minus sequential")) +
             ylim(-1,1) +
             xlab(NULL)

    ## repeat code for adaptive
    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)  

    #cat("B)")

    plotlist[[2]] <- 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) +
             ggtitle(paste("B) Adaptive minus sequential")) +
             ylim(-1,1) +
             xlab(NULL)



  #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]],ncol=1, rel_heights=c(1,1.6), labels='AUTO') 

\pagebreak

Fig9 Effectiveness and exposure together determine most of whether resistance arises slower under adaptive mixture or sequential strategies

1. 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.

  ggSubset <- dif_mixA_seq
  #if mix-seq >1 resistance took longer for mix so mix is better
  ggSubset$mixA_better_seq0.5 <- ifelse( ggSubset$mixA_minus_seq0.5 > 0, 1, 0 )
  ggSubset$seq_better_mixA0.5 <- ifelse( ggSubset$mixA_minus_seq0.5 < 0, 1, 0 )
  ggSubset$abs_diff <- abs(ggSubset$mixA_minus_seq0.5)

  #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$mixA_better_seq0.5 <- 1
  df_same$seq_better_mixA0.5 <- 1
  #rbind them back onto other dataset
  ggSubset <- rbind( ggSubset, df_same)


    #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_mixA0.5, labeller=as_labeller( c('0'='adaptive mix\nresistance slower','1'='sequential\nresistance slower'))) +

             facet_grid( ~ seq_better_mixA0.5, labeller=as_labeller( c('0'='adaptive mix\ntime to resistance\n longer','1'='sequential\ntime to resistance\n longer'))) +

             #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) + 
             #!!!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
  #(Intercept)          x1 
  ## 1.0746184 0.5590287

  #20/6/16 calculated T/F positives & negatives, according to the linear model calculation
  #this is a very heath-robinson way of doing, must be a better way
  #so I want to find whether a point is above or below the line (and probably ignore zeros where mix & sequence the same)
  #the eq is approx y=0.56*x+1.07
  #so I could create a new column that has the expected lm val of added effectivenesses
  #where +ve is mix better
  #true+ve : effectiveness1+2 > lm & mixA_minus_seq0.5>0
  #true-ve : effectiveness1+2 < lm & mixA_minus_seq0.5<0
  #false+ve : effectiveness1+2 > lm & mixA_minus_seq0.5<0
  #false-ve : effectiveness1+2 < lm & mixA_minus_seq0.5>0  
  ggSubset$l_mod_eff <- l_mod$coefficients[2] * ggSubset$exposure + l_mod$coefficients[1]

  #BEWARE that this sets 0s (where mix & sequence are equal) to NA
  ggSubset$l_mod_confusion <- 
    ifelse ( (ggSubset$effectiveness_ins1+ggSubset$effectiveness_ins2) > ggSubset$l_mod_eff,
               ifelse(ggSubset$mixA_minus_seq0.5>0,"TP",ifelse(ggSubset$mixA_minus_seq0.5<0,"FP",NA))

            ,ifelse(ggSubset$mixA_minus_seq0.5<0,"TN",ifelse(ggSubset$mixA_minus_seq0.5>0,"FN",NA))
            )

  table(ggSubset$l_mod_confusion)

  #  FN   FP   TN   TP 
  #602  219 3834 2787 

  #sensitivity or true positive rate TPR=TP/(TP+FN)
  TPR = length(which(ggSubset$l_mod_confusion=="TP")) / (length(which(ggSubset$l_mod_confusion=="TP"))+length(which(ggSubset$l_mod_confusion=="FN")))
  #specificity (SPC) or true negative rate SPC=TN/(TN+FP)
  SPC = length(which(ggSubset$l_mod_confusion=="TN")) / (length(which(ggSubset$l_mod_confusion=="TN"))+length(which(ggSubset$l_mod_confusion=="FP")))

  cat("sensitivity(TPR) : ",TPR)#,"\n")
  cat("specificity(SPC) : ",SPC)#,"\n")

  #approximation
  #x2 <- seq(0,1,0.1)
  #y2 <- 0.6*x2+0.98
  #plot(x2,y2)    

Red triangles are those runs in which resistance reached the threshold at the same time in the two strategies (and are repeated in both panels). The red line is a linear model through these same points. \pagebreak

Fig10 Modifying previous plot to show where resistance arises over 20% slower for adaptive mixture versus sequential.

1. 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 adaptive mixture under these constraints.

  #add an option for 20% better, put the not in to get facets in correct order in ggplot later
  ggSubset$mixA_not20pc_better_seq0.5 <- ifelse( ggSubset$mixA_minus_seq0.5 / ggSubset$gen_cP0.5seq <= 0.2, 1, 0 )   


  df_20pc <- ggSubset[which(ggSubset$mixA_minus_seq0.5 / ggSubset$gen_cP0.5seq > 0.19 & 
                            ggSubset$mixA_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) != 'mixA_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( ~ mixA_not20pc_better_seq0.5, labeller=as_labeller( c('0'='adaptive mix\ntime to resistance\n >=20% longer','1'='adaptive mix\ntime to resistance\n <=20% longer'))) +

             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 the 0 points abs_diff
             #geom_smooth(data=df_same, method = "lm", se = FALSE)
             #todo change 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 

#20/6/16 calculated T/F positives & negatives, according to the linear model calculation
  #this is a very heath-robinson way of doing, must be a better way
  #so I want to find whether a point is above or below the line (and probably ignore zeros where mix & sequence the same)
  #the eq is approx y=0.56*x+1.07
  #so I could create a new column that has the expected lm val of added effectivenesses
  #where +ve is mix better
  #true+ve : effectiveness1+2 > lm & mixA_minus_seq0.5>0
  #true-ve : effectiveness1+2 < lm & mixA_minus_seq0.5<0
  #false+ve : effectiveness1+2 > lm & mixA_minus_seq0.5<0
  #false-ve : effectiveness1+2 < lm & mixA_minus_seq0.5>0  
  ggSubset$l_mod_eff <- l_mod$coefficients[2] * ggSubset$exposure + l_mod$coefficients[1]

  #BEWARE that this sets 0s (where mix & sequence are equal) to NA
  ggSubset$l_mod_confusion <- 
    ifelse ( (ggSubset$effectiveness_ins1+ggSubset$effectiveness_ins2) > ggSubset$l_mod_eff,
               ifelse(ggSubset$mixA_not20pc_better_seq0.5==0,"TP",ifelse(ggSubset$mixA_not20pc_better_seq0.5==1,"FP",NA))

            ,ifelse(ggSubset$mixA_not20pc_better_seq0.5==1,"TN",ifelse(ggSubset$mixA_not20pc_better_seq0.5==0,"FN",NA))
            )

  table(ggSubset$l_mod_confusion)

  #  FN   FP   TN   TP 
  ## 251 468 6178 915

  #sensitivity or true positive rate TPR=TP/(TP+FN)
  TPR = length(which(ggSubset$l_mod_confusion=="TP")) / (length(which(ggSubset$l_mod_confusion=="TP"))+length(which(ggSubset$l_mod_confusion=="FN")))
  #specificity (SPC) or true negative rate SPC=TN/(TN+FP)
  SPC = length(which(ggSubset$l_mod_confusion=="TN")) / (length(which(ggSubset$l_mod_confusion=="TN"))+length(which(ggSubset$l_mod_confusion=="FP")))

  cat("sensitivity(TPR) : ",TPR)#,"\n")
  cat("specificity(SPC) : ",SPC)#,"\n")  

\pagebreak

Fig8 Classification Trees for whether resistance arises slower with an adaptive 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.

  #layout for 2 panel plot (t=horiz, c=vertical)
  layout(c(1:2))

  require('plyr')

  #4/7/16 need to remove runs where gens >= 500
  length(which(resistPointsSeq['gen_cP0.5',] >= 500 | resistPointsMix_A['gen_cP0.5',] >= 500))
  #2376
  length(which(resistPointsSeq['gen_cP0.5',] >= 500 ))
  #2151
  length(which(resistPointsMix_A['gen_cP0.5',] >= 500))
  #2306

  #BEWARE1 of subtracting if diff num in each ...
  #this removes runs that fail either

  failed_run_indices <- (resistPointsSeq['gen_cP0.5',] >= 500 | resistPointsMix_A['gen_cP0.5',] >= 500)

  resistPointsMix_A <- resistPointsMix_A[,!failed_run_indices]
  resistPointsSeq <- resistPointsSeq[,!failed_run_indices]
  #BEWARE2 I need to remove these indices from the predictors too
  treeInput <- treeInput[,!failed_run_indices]

  ## is mix adaptive better than sequential ?
  #T/F
  resistBetterMix_ASeqBoolean <- resistPointsMix_A > resistPointsSeq
  #convert to 0/1
  resistBetterMix_ASeq <- plyr::aaply(resistBetterMix_ASeqBoolean,.margins=c(1,2),.fun=as.numeric)
  #rename rows ready for binding on to inputs
  rownames(resistBetterMix_ASeq) <- gsub('gen','betterMix_ASeq', rownames(resistBetterMix_ASeq))

  #20% better
  resistBetterMix_ASeqBoolean20 <- resistPointsMix_A > (1.2 * resistPointsSeq)
  #convert to 0/1
  resistBetterMix_ASeq20 <- plyr::aaply(resistBetterMix_ASeqBoolean20,.margins=c(1,2),.fun=as.numeric)
  #rename rows ready for binding on to inputs
  rownames(resistBetterMix_ASeq20) <- gsub('gen','betterMix_ASeq20', rownames(resistBetterMix_ASeq20))


  #replace 0s & 1s with mixture/sequence to make tree plot clearer
  repMix <- function(x){
    x[x==0] <- 'sequence'
    x[x==1] <- 'mixture'
    x
  }

  resistBetterMix_ASeq <- repMix(resistBetterMix_ASeq)
  resistBetterMix_ASeq20 <- repMix(resistBetterMix_ASeq20)

  treeInput <- rbind( treeInput, resistBetterMix_ASeq, resistBetterMix_ASeq20 )

  #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="+") 

  treeResponses <- c( "betterMix_ASeq_cP0.5","betterMix_ASeq20_cP0.5")


  #pruned trees commented out, not used in paper1
  # 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)
  #   
  #   #test tree pruning, seemed to suggest we are pruning more than needed
  #   # plotcp(tree)  
  #   # cat("optimal cp used in pruning=",cpOptimal,"\n")
  #   
  #   #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(treePruned, extra=1, varlen = 0, main=main, box.col=cols, under=TRUE)
  # }


  #4/7/2016 unpruned trees, repeats previous loop with pruning removed
  for( treeResponse in treeResponses )
  {
    tree <- rpart::rpart(as.formula(paste(treeResponse,"~",treePredictorString)), data = treeInput, method = 'class')

    #to set box colours
    cols <- ifelse(tree$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
    #cex to set textsize, good that A&B are now same size
    #round=0 to make boxes square
    #split.font=1 to stop bold text in labels
    prp(tree, extra=1, varlen = 0, main=main, box.col=cols, under=TRUE, 
        cex=0.6, round=0, split.font=1, ycompress = TRUE,
        eq="\n= ", lt="\n< ", ge="\n>= " )
  }

Fig11 How sequential can be better than mixture. When the time taken for resistance to the first insecticide to be reached is greater than the difference in time taken for the second resistance to be reached in the presence and absence of the first insecticide.

#sequential best. effectivenesses bothg reduced by 0.2 from curtis fig2
runcurtis_f2( P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.17 , h.RS2_0B = 0.0016 , exposure = 0.9 , phi.SS1_A0 = 0.53 , phi.SS2_0B = 0.8 , s.RR1_A0 = 0.23 , s.RR2_0B = 0.43, addCombinedStrategy=FALSE )

AndySouth/resistance documentation built on Nov. 12, 2020, 3:39 a.m.