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

Figures for paper 1 after change to using rr_restoration instead of selection coeff from sensi analysis

version 1 - in progress

  #library(devtools)
  #install_github('AndySouth/resistance')
  library(resistance)

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

  #default setting of the random seed in sensiAnPaperPart() should make the input files mostly the same
  #the exception is the exposure array a. The single exposure input will be the same, but the array changes
  #according to setExposure(  )

  #nScenarios <- 10  
  #nScenarios <- 100
  #nScenarios <- 500 
  #nScenarios <- 1000
  nScenarios <- 10000

  #to run 'extended' or not
  experiment <- 'extended' #'curtis'


  ## do model runs
  listOutMix <- sensiAnPaperPart( nScenarios, insecticideUsed = 'mixture', experiment=experiment )
  listOutI1 <- sensiAnPaperPart( nScenarios, insecticideUsed = 'insecticide1', experiment=experiment )
  listOutI2 <- sensiAnPaperPart( nScenarios, insecticideUsed = 'insecticide2', experiment=experiment )


  ## save results objects as rda for re-use
  save(listOutMix, file=paste0(outFolder,'listOutMix_rr_10000.rda'))
  save(listOutI1, file=paste0(outFolder,'listOutI1_rr_10000.rda'))
  save(listOutI2, file=paste0(outFolder,'listOutI2_rr_10000.rda'))  

  library(resistance)

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

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

  #doing; for all inputs
  if (experiment=='extended') {
    num_inputs <- 13
    } else {
    num_inputs <- 11  
    }

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

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


  #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 rr_restoration these not to be used anymore
  #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
  library(ggplot2)

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

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


  # #23/6/2016 add a check for which runs are being excluded (note that some are excluded in a different place too) 
  # #just those runs that didn't make threshold
  ggSubset <- ggInsOuts[didntReachThresh,]
  names_inputs <- names(ggInsOuts)[1:num_inputs]
  # 
  # ggSubset <- ggSubset[ !ggSubset$strategy %in% c("insecticide 1","insecticide 2"), ]
  # 
  # for(i in names_inputs)
  # {
  #   y <- 'time_to_resistance0.5'    
  #   
  #   print( ggplot(ggSubset, aes_string(x=i, y=y, colour="strategy")) + 
  #            
  #            geom_point(shape=3) #, show.legend=FALSE) + 
  #            #geom_smooth(linetype='dashed',size=1)
  #          )
  # }  

  # #seems that most failed runs are at low start freqs
  # hist(ggSubset$start_freq_allele1)
  # #compare to all runs
  # hist(ggInsOuts$start_freq_allele1)  
  # #and for start freqs combined
  # hist(ggSubset$start_freq_allele1+ggSubset$start_freq_allele2)
  # #compare to all runs
  # hist(ggInsOuts$start_freq_allele1+ggInsOuts$start_freq_allele2)
  # #also low exposure is a large contributor to failed runs
  # hist(ggSubset$exposure)  
  # #compared to all runs
  # hist(ggInsOuts$exposure) 

  #quick look at low start freq runs. Can they make it to thresholds ? Yes
  # for(i in names_inputs)
  # {
  #   y <- 'time_to_resistance0.5'    
  #   
  #   print( ggplot(ggInsOuts[ggInsOuts$start_freq_allele1 < 0.001 & ggInsOuts$start_freq_allele2 < 0.001,], aes_string(x=i, y=y, colour="strategy")) + 
  #            
  #            geom_point(shape=3) + #, show.legend=FALSE) + 
  #            geom_smooth(linetype='dashed',size=1)
  #          )
  # }   


  ####### end of new check

  #subset by runs not making threshold
  ggInsOuts <- ggInsOuts[-didntReachThresh,]


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

### 4/2016 NEW part to analyse what happens when a new insecticide is added to an existing one
### either to replace it or to continue in a mixture

  # find the resistance points for I1 & I2 in a mixture 
  #(diff than Mix_1 & 2 which is for 1st & 2nd threshold to be reached)
  resistPointsMixI1 <- findResistancePoints(listOutMix, locus=1)
  resistPointsMixI2 <- findResistancePoints(listOutMix, locus=2)

  resistPointsMixI1_T <- addStrategyColumn(resistPointsMixI1,"I1inMix",ggInput) 
  resistPointsMixI2_T <- addStrategyColumn(resistPointsMixI2,"I2inMix",ggInput)     

  #and use these from above
  #resistPointsI1_T
  #resistPointsI2_T

  #OldPlusNew 
  # find whether I1 or I2 has lower starting freq this is the 'new' insecticide
  # if I1 < I2 use resistPointsMixI1_T
  # if I2 < I1 use resistPointsMixI2_T                                               

  #NewOnly
  # if I1 < I2 use resistPointsI1_T
  # if I2 < I1 use resistPointsI2_T     

  # find whether I1 or I2 has lower starting freq this is the 'new' insecticide
  # rember the inputs are copied onto the start of all resistPoints objects

  indicesI2_lessthan_I1 <- which( resistPointsI1_T$start_freq_allele2 < resistPointsI1_T$start_freq_allele1 )
  #this shows about half, I also checked no equals which is good
  #length(indicesI2_lessthan_I2)
  #[1] 4983
  #check whether this works
  #check <- resistPointsI2_T[indicesI2_lessthan_I1, c('gen_cP0.1','gen_cP0.25','gen_cP0.5')]
  #aha! problem I was having is that I1 & I2 inputs are not swapped around
  #if I wanted to swap the inputs for I1 & I2 around, would need to do other coefficients etc.
  #much potential confusion
  #might be better just to restrict to the runs in which I2 < I1
  resistPointsNewOnly <- resistPointsI2_T[indicesI2_lessthan_I1, ]
  resistPointsOldPlusNew <- resistPointsMixI2_T[indicesI2_lessthan_I1, ]  

  # #this bit now replaced by the 2 previous lines
  # #NewOnly first make a copy from existing
  # resistPointsNewOnly <- resistPointsI1_T
  # #replace those values where I2 is the 'new' insecticide with lower resistance
  # resistPointsNewOnly[indicesI2_lessthan_I1, c('gen_cP0.1','gen_cP0.25','gen_cP0.5')] <-
  #    resistPointsI2_T[indicesI2_lessthan_I1, c('gen_cP0.1','gen_cP0.25','gen_cP0.5')]  
  # 
  # #OldPlusNew first make a copy from existing
  # resistPointsOldPlusNew <- resistPointsMixI1_T
  # #replace those values where I2 is the 'new' insecticide with lower resistance
  # resistPointsNewOnly[indicesI2_lessthan_I1, c('gen_cP0.1','gen_cP0.25','gen_cP0.5')] <-
  # resistPointsMixI2_T[indicesI2_lessthan_I1, c('gen_cP0.1','gen_cP0.25','gen_cP0.5')] 


  resistPointsNewOnly$strategy <- "new only"
  resistPointsOldPlusNew$strategy <- "old plus new"

  #BE CAREFUL

  #format data to enable PRCC on difference between the 2 strategies
  #follow previous naming where
  #OldPlusNew="gen_cP0.5", NewOnly="gen_cP0.5seq"


  #columns to remove from first
  indices1 <- which(names(resistPointsOldPlusNew) %in% c("gen_cP0.1","gen_cP0.25","strategy"))
  #columns to add to 2nd  
  indices2 <- which(names(resistPointsNewOnly) %in% c("gen_cP0.5")) 
  #rename column in 2nd
  tmp <- resistPointsNewOnly[indices2]
  names(tmp) <- "gen_cP0.5seq"

  dif_oldnew_newonly <- cbind(resistPointsOldPlusNew[-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_oldnew_newonly$gen_cP0.5 > 500 | dif_oldnew_newonly$gen_cP0.5seq > 500 )
  #subset by runs not making threshold
  dif_oldnew_newonly <- dif_oldnew_newonly[-didntReachThresh,]  

  dif_oldnew_newonly["mix_minus_sole0.5"] <- dif_oldnew_newonly["gen_cP0.5"] - dif_oldnew_newonly["gen_cP0.5seq"]


  #*****************
  # to check
  # are there any cases where startfreqI1 > startfreqI2
  # by my method above i don't think there should be
  # i think min should be 1
  min(dif_oldnew_newonly$resist_start_1_div_2)
  max(dif_oldnew_newonly$resist_start_1_div_2)  

  min( dif_oldnew_newonly$start_freq_allele1 - dif_oldnew_newonly$start_freq_allele2)
  max( dif_oldnew_newonly$start_freq_allele1 - dif_oldnew_newonly$start_freq_allele2)

  #resistPointsNewOnly
  min(resistPointsNewOnly$resist_start_1_div_2)
  max(resistPointsNewOnly$resist_start_1_div_2) 

  # because I now just use runs where startfreq12 < I1
  # resist_start_1_div_2 should be the inverse of resist_start_hi_div_lo
  # or at least directly correlated
  # yes they are all on a single curved line
  # plot( resistPointsNewOnly$resist_start_1_div_2 ~ resistPointsNewOnly$resist_start_lo_div_hi)

  #not sure I need to do this because all results already show mixture better
  #now I want to subset just those runs where start freq old is > 10 * new 


  #proportion of runs in which mix better than sole
  num_mix_better <- length(which(dif_oldnew_newonly["mix_minus_sole0.5"] > 0))
  num_all <- nrow(dif_oldnew_newonly)
  num_mix_better/num_all
  #1
  #?now mix always seems to be better
  #which kind of makes sense, any insecticide better than none


  #28/4/16 testing Ians statement
  #does resistance always arise faster in sole use than in a mixture (irrespective of starting frequencies)
  #resistPointsI1_T #"insecticide 1" 
  #resistPointsI2_T #"insecticide 2"
  #resistPointsMixI1_T #"I1inMix" 
  #resistPointsMixI2_T #"I2inMix"

  mix_minus_soleI1 <- resistPointsMixI1_T["gen_cP0.5"] - resistPointsI1_T["gen_cP0.5"]
  which( mix_minus_soleI1 < 0 ) 
  #0
  mix_minus_soleI2 <- resistPointsMixI2_T["gen_cP0.5"] - resistPointsI2_T["gen_cP0.5"]
  which( mix_minus_soleI2 < 0 ) 
  #0

  #so yes in our param space adding an extra insecticide always increases resistance time for the first

  #i could plot it against exposure & effectiveness

---

Fig. x1 Time-to-resistance across all input values

  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 ) +
      #scale_x_discrete("Resistance to which insecticide", limits=c("insecticide 1","Mix either","Sequential","Mix adaptive","Mix both"), labels = c("insecticide 1"="Sole use", "Mix either"="1st in mix","Sequential"="2nd in sequence", "Mix adaptive"="2nd in mix\nswitch to sole", "Mix both"="2nd in mix" ) )
      scale_x_discrete("Resistance to which insecticide", limits=c("insecticide 1","Mix either","Sequential","Mix adaptive","Mix both"), labels = c("insecticide 1"="Sole use", "Mix either"="1st in mix","Sequential"="2nd in sequence", "Mix adaptive"="adaptive mix", "Mix both"="2nd in mix" ) )

  ) #end print
}  

Fig. x2 How time-to-resistance is effected by input values for each insecticide use strategy

  names_inputs <- names(ggInsOuts)[1:num_inputs]

  ggSubset <- ggInsOuts[ !ggInsOuts$strategy %in% c("insecticide 1","insecticide 2"), ]

  ggSubset$effectiveness_added <- ggSubset$effectiveness_ins1 + ggSubset$effectiveness_ins2
  ggSubset$dominance_added <- ggSubset$dominance_allele1 + ggSubset$dominance_allele2
  names_inputs <- c(names_inputs, "effectiveness_added", "dominance_added")  

  for(i in names_inputs)
  {
    y <- 'time_to_resistance0.5'    

    print( ggplot(ggSubset, aes_string(x=i, y=y, colour="strategy")) + 

             geom_point(shape=3, size=0.2, alpha=0.2, show.legend=FALSE) + 
             geom_smooth(linetype='dashed',size=1) +

             #labs(title = i) +
             #add line for Curtis inputs
             geom_vline(data = curtisInputs, aes_string(xintercept = i), linetype='dotted', colour='red') +
             geom_text(data = curtisInputs, inherit.aes = FALSE, aes_string(x = i, y=210, label='"Curtis Fig.2"'), size=3, colour='red')
           )
  }

Fig. xx2 Quick check whether difference between effectiveness 1 & 2 dissapears when minimising excluded runs

  #names_inputs <- names(ggInsOuts)[1:num_inputs]
  names_inputs <- c("effectiveness_ins1", "effectiveness_ins2", "effectiveness_added")  

  ggSubset <- ggInsOuts[ !ggInsOuts$strategy %in% c("insecticide 1","insecticide 2","Mix either","Sequential"), ]

  ggSubset$effectiveness_added <- ggSubset$effectiveness_ins1 + ggSubset$effectiveness_ins2

  # ggSubset <- ggSubset[ ggSubset$start_freq_allele1 > 0.005 & 
  #                       ggSubset$start_freq_allele2 > 0.005 &
  #                       ggSubset$dominance_allele1 > 0.1 &
  #                       ggSubset$dominance_allele2 > 0.1 &
  #                       ggSubset$exposure > 0.3  , ]

  for(i in names_inputs)
  {
    y <- 'time_to_resistance0.5'    

    print( ggplot(ggSubset, aes_string(x=i, y=y, colour="strategy")) + 

             geom_point(shape=1, size=1, alpha=0.2, show.legend=FALSE) + 
             geom_smooth(linetype='dashed',size=1)

             #labs(title = i) +
             #add line for Curtis inputs
             #geom_vline(data = curtisInputs, aes_string(xintercept = i), linetype='dotted', colour='red') +
             #geom_text(data = curtisInputs, inherit.aes = FALSE, aes_string(x = i, y=210, label='"Curtis Fig.2"'), size=3, colour='red')
           )
  }

Fig. x3 How difference in time-to-resistance between sequential and adaptive mixture is effected by input values

  ggSubset <- dif_mixA_seq

  names_inputs <- names(ggSubset)[1:num_inputs]


  for(i in names_inputs)
  {
    #x11()

    y <- "mixA_minus_seq0.5"    

    print( ggplot(ggSubset, aes_string(x=i, y=y)) + #, colour="strategy")) + 
            #points not wanted if 10000
             geom_point(shape=1, colour="darkblue", alpha=0.2, show.legend=FALSE) + 
             #geom_smooth(colour='red', linetype='dashed',size=0.5) +
             geom_smooth(linetype='dashed',size=1, colour="red") +
             #facet_wrap( ~ strategy) +
             #geom_smooth(aes_string(x=i, y=y, color=NULL)) )
             labs(title = i) +
             geom_hline(yintercept = 0, colour="white")
           )
  }

Fig. xx3 Using a new insecticide on its own or in combination with old

  ggSubset <- dif_oldnew_newonly

  names_inputs <- names(ggSubset)[1:num_inputs]


  for(i in names_inputs)
  {
    #x11()

    y <- "mix_minus_sole0.5"    

    print( ggplot(ggSubset, aes_string(x=i, y=y)) + #, colour="strategy")) + 
             geom_point(shape=1, colour="darkblue", alpha=0.2, show.legend=FALSE) + 
             #geom_smooth(colour='red', linetype='dashed',size=0.5) +
             geom_smooth(linetype='dashed',size=1, colour="red") +
             #facet_wrap( ~ strategy) +
             #geom_smooth(aes_string(x=i, y=y, color=NULL)) )
             labs(title = i) +
             geom_hline(yintercept = 0, colour="white")
           )
  }

Fig xxx PRCC difference in time-to-resistance between sole new and mix old+new

  library(sensitivity)

    x <- dif_oldnew_newonly[,1:num_inputs]
    y <- dif_oldnew_newonly["mix_minus_sole0.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)  

    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) +
             ylim(-1,1) +
             ggtitle(paste("PRCC benefit of mix old+new over sole new (old=1, new=2)")) +
             xlab(NULL)
          )

  #can I do a 2D hexbin plot with x=exposure, y=effectiveness
  #but do I just want to use 1 effectiveness, or maybe add both effectivenesses together ?
  #bin colour : proportion of runs in which sequential better & could have diverging scale about 0.5
  #I probably still need to tweak data first
  #currently I have differences in times to resistance : mixA_minus_seq0.5
  #what I think I want is the proportion of runs in which mix is better
  #I can get that by setting to 0 & 1 and calc average
  # tst_runs <- c(1,0,0,0)
  # tst_runs <- c(1,1,1,0)
  # mean(tst_runs)


  ##BEWARE toogle this to look at plots for mixboth vs seq
  ggSubset <- dif_mixA_seq
  #ggSubset <- dif_mix2_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_mixA_minus_seq0.5 <- 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$mixA_minus_seq0.5==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, probably not needed ??
  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$mixA_not20pc_better_seq0.5 <- ifelse( ggSubset$mixA_minus_seq0.5 / ggSubset$gen_cP0.5seq <= 0.2, 1, 0 )   


    #failed hex plot just gave me counts
    # print( ggplot(ggSubset, aes_string(x='exposure', y='effectiveness_ins1', colour='mixA_better_seq0.5'))+
    #          geom_hex()
    #        )

    #interesting plot showing where sequence is better is restricted to higher values of exposure and effectiveness 
    print( ggplot(ggSubset, aes_string(x='exposure', y='effectiveness_ins1', colour='mixA_better_seq0.5'))+
             #geom_hex()
             geom_point(alpha=0.2)  
           ) 

    #can I multiply the two effectivenesses together ?
    #aha! yes that's even better
    print( ggplot(ggSubset, aes(x=exposure, y=effectiveness_ins1 * effectiveness_ins2, colour=mixA_better_seq0.5))+
             #geom_hex()
             geom_point(alpha=0.2)  
           )  

    #i think adding the effectivenesses may give an even better plot
    #and makes it easier for the reader to interpret
    #aha! yes that's even better
    print( ggplot(ggSubset, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2, colour=mixA_better_seq0.5))+

             geom_point(alpha=0.2)  
           )  

    #facetted 0/1, labelled mix better, seq better
    #if I get this sorted I could try to have +ve & -ve on the same scale
    print( ggplot(ggSubset, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2, colour=seq_better_mixA0.5))+
             facet_grid( ~ seq_better_mixA0.5, labeller=as_labeller( c('0'='adaptive mix better','1'='sequential better'))) +

             theme(legend.position="none") +
             geom_point(alpha=1)  
           )     


    #what if I do the difference rather than just 0/1
    #aha! this is close to being the killer plot, just need to get the colour scale right 
    #but the previous one may be clearer
    print( ggplot(ggSubset, aes(x=exposure, y=effectiveness_ins1 * effectiveness_ins2, colour=mixA_minus_seq0.5))+
             #geom_hex()
             geom_point(alpha=1) +
             #no adding grey outlines obscures the pattern
             #geom_point(shape=21, colour='grey') +
             #this one does look good but slight potential to see negatives above even though they aren't there
             #could I modify the shapes to be different
             #would a repeated, faceted plot solve that ?
             #scale looks good but colours could do with being denser
             #and I'm not sure that making the points that are close to zero white is a good idea ?
             #although it does mean points that are only a few generations different do not obscure the pattern
             #maybe I could explicitly do Ians 20% thing & colour those white
             scale_colour_gradient2()

             #didn't manage to get the colorBrewer scale to diverge from 0
             #scale_colour_distiller("RdYlGn",type='div')
            )     

    #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_mixA_minus_seq0.5, size="Curtis"))+
    print( ggplot(ggSubset, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2, colour=abs_mixA_minus_seq0.5))+

             facet_grid( ~ seq_better_mixA0.5, labeller=as_labeller( c('0'='adaptive mix better','1'='sequential better'))) +

             #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")) + 
             theme_bw() +
             #theme_minimal() +
             #add point for Curtis inputs
             geom_point(data = curtisInputs, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2), colour='red', 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) + 
             #add text for Curtis point
             geom_text(data = curtisInputs, inherit.aes = FALSE, aes_string(x = 'exposure', y='effectiveness_ins1 + effectiveness_ins2', label='"Curtis Fig.2"'), nudge_x=-0.05, size=3, colour='red')
           )       

    #NOW what if I subset the 0 values & fit a linear model to it
    #could that be used to decide between the 2 strategies ?
    #I might even be able to do the lm plot within ggplot

    #trickiness I have is that I have previously put 0 values in the adptive mix better plot 
    #which is maybe a little misleading anyway
    #and also means that I only get the lm line in the one plot
    #maybe what I can do is add the zeros into both plots then  will get lines in both
    #might even be good to overplot the zero values after
    #or even remove the zero values ?

    print( ggplot(ggSubset, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2, colour=abs_mixA_minus_seq0.5))+

             facet_grid( ~ seq_better_mixA0.5, labeller=as_labeller( c('0'='adaptive mix better','1'='sequential better'))) +

             scale_color_viridis(begin=1, end=0, option="viridis", 
                                 guide=guide_colourbar(title="difference in \ngenerations to reach\nresistance threshold")) + 
             theme_bw() +
             #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) + 
             #!!!NOW can I add a lm fitted to the 0 points mixA_minus_seq0.5
             geom_smooth(data=ggSubset[which(ggSubset$mixA_minus_seq0.5==0),], method = "lm", se = FALSE) +
             #can I subtley highlight the zero points which are repeated in each plot ?
             geom_point(data=ggSubset[which(ggSubset$mixA_minus_seq0.5==0),], size=0.8, colour='lightgrey', shape=2, fill='white') +            
             #add text for Curtis point
             geom_text(data = curtisInputs, inherit.aes = FALSE, aes_string(x = 'exposure', y='effectiveness_ins1 + effectiveness_ins2', label='"Curtis Fig.2"'), nudge_x=-0.05, size=3, colour='red')
           )      

  # doing the lm outside of ggplot    
  df_same = ggSubset[which(ggSubset$mixA_minus_seq0.5==0),]
  x1=df_same$exposure 
  y1=df_same$effectiveness_ins1 + df_same$effectiveness_ins2

  l_mod <- lm(y1 ~ x1)
  l_mod$coefficients
  #(Intercept)           x 
  #  0.9790948   0.6095132 

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


  #modify the killer plot to show where mixtures are 20% better ?
  print( ggplot(ggSubset, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2))+ #, colour=abs_mixA_minus_seq0.5))+

             facet_grid( ~ mixA_not20pc_better_seq0.5, labeller=as_labeller( c('0'='adaptive mix\nresistance >20% slower','1'='adaptive mix\nresistance <20% slower'))) +

             theme_bw() +
             #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
             geom_smooth(data=df_same, method = "lm", se = FALSE) 
           )  

repeating killer plot for mixtures both versus sequential

  #try repeating the killer plot for mixtures both versus sequential


  ##BEWARE toogle this to look at plots for mixboth vs seq
  #ggSubset <- dif_mixA_seq
  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_mix2_minus_seq0.5 <- abs(ggSubset$mix2_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$mix2_minus_seq0.5==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, probably not needed ??
  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 )   

    #i think adding the effectivenesses may give an even better plot
    #and makes it easier for the reader to interpret
    #aha! yes that's even better
    print( ggplot(ggSubset, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2, colour=mix2_better_seq0.5))+

             geom_point(alpha=0.2)  
           )  

    #facetted 0/1, labelled mix better, seq better
    #if I get this sorted I could try to have +ve & -ve on the same scale
    print( ggplot(ggSubset, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2, colour=seq_better_mix20.5))+
             facet_grid( ~ seq_better_mix20.5, labeller=as_labeller( c('0'='mix better','1'='sequential better'))) +

             theme(legend.position="none") +
             geom_point(alpha=1)  
           )     


    #killer plot mix2 vs seq
    library(viridis)
    print( ggplot(ggSubset, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2, colour=abs_mix2_minus_seq0.5))+

             facet_grid( ~ seq_better_mix20.5, labeller=as_labeller( c('0'='mix better','1'='sequential better'))) +

             #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")) + 
             theme_bw() +
             #theme_minimal() +
             #add point for Curtis inputs
             geom_point(data = curtisInputs, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2), colour='red', 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) + 
             #add text for Curtis point
             geom_text(data = curtisInputs, inherit.aes = FALSE, aes_string(x = 'exposure', y='effectiveness_ins1 + effectiveness_ins2', label='"Curtis Fig.2"'), nudge_x=-0.05, size=3, colour='red')
           )       


    #modify killer plot mix2 vs seq to add in male exposure prop on x axis
    library(viridis)
    print( ggplot(ggSubset, aes(x=exposure+exposure*male_exposure_prop, y=effectiveness_ins1 + effectiveness_ins2, colour=abs_mix2_minus_seq0.5))+

             facet_grid( ~ seq_better_mix20.5, labeller=as_labeller( c('0'='mix better','1'='sequential better'))) +
             scale_color_viridis(begin=1, end=0, option="viridis", 
                                 guide=guide_colourbar(title="difference in \ngenerations to reach\nresistance threshold")) + 
             theme_bw() +
             #add point for Curtis inputs
             geom_point(data = curtisInputs, aes(x=exposure+exposure*male_exposure_prop, y=effectiveness_ins1 + effectiveness_ins2), colour='red', shape=3) +
             ylab('effectiveness of insecticides added') +
             xlab('exposure to both insecticides for both sexes') + 
             #main points for the plot
             geom_point(alpha=1, size=0.5) + 
             #add text for Curtis point
             geom_text(data = curtisInputs, inherit.aes = FALSE, aes_string(x = 'exposure', y='effectiveness_ins1 + effectiveness_ins2', label='"Curtis Fig.2"'), nudge_x=-0.05, size=3, colour='red')
           )       

  #modify the killer plot to show where mixtures are 20% better ?
  print( ggplot(ggSubset, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2))+ #, colour=abs_mixA_minus_seq0.5))+

             facet_grid( ~ mix2_not20pc_better_seq0.5, labeller=as_labeller( c('0'='mix\nresistance >20% slower','1'='mix\nresistance <20% slower'))) +

             theme_bw() +
             #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
             geom_smooth(data=df_same, method = "lm", se = FALSE) 
           )  

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

  library(sensitivity)

    x <- dif_mixA_seq[,1:num_inputs]
    y <- dif_mixA_seq["mixA_minus_seq0.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)  

    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) +
             ylim(-1,1) +
             ggtitle(paste("PRCC difference between sequential and mixture adaptive")) +
             xlab(NULL)
          )

PRCC analysis for time-to-resistance subsetted by each insecticide use strategy

  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)

  for(strategy in strategies)
  {  


    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)  

    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) +
             ylim(-1,1) +
             #ggtitle(paste(strategy,"PRCC")) +
             ggtitle(paste(strategy,"time-to-resistance 0.5")) +
             xlab(NULL)
          )
  }

response of resistance thresholds to individual inputs faceted by strategy

  names_inputs <- names(ggInsOuts)[1:num_inputs]

  #todo find out to take sample this doesn't work
  #ggSample <- sample(ggInsOuts,500)

  for(i in names_inputs)
  {
    #x11()

    y <- 'time_to_resistance0.5'    
    #y <- 'time_to_resistance0.5'
    color <- 'strategy'
    #coloring by exposure shows its effect even when other params are varying
    #color <- 'exposure'    

    print( ggplot(ggInsOuts, aes_string(x=i, y=y, color=color)) + 
             #points not wanted if 10000
             #geom_point(shape=3, show.legend=FALSE) + 
             geom_smooth(colour='red', linetype='dashed',size=0.5, show.legend=FALSE) +
             facet_wrap( ~ strategy) + #, show_guide = FALSE) + 
             #geom_smooth(aes_string(x=i, y=y, color=NULL)) )
             labs(title = i))
  }

#eval=FALSE so this one not displayed

## 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 , rr_restoration_ins1=0.23/0.73, rr_restoration_ins2 = 0.43/1, 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 , rr_restoration_ins1=0.23/0.73, rr_restoration_ins2 = 0.43/1, addCombinedStrategy=FALSE )

Fig. x5.1 modified for mix2 vs seq

NOT USED IN PAPER1 fewer runs where sequential gives slower resistance

  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)

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

    #find unequal runs
    unequal <- ggSubset[ggSubset$resist_start_lo_div_hi < 0.01,]
    #nrow(ggSubset)
    #nrow(unequal)
    #2721 out of the 6189 remaining runs

    #in how many runs was seq better than mix
    #sum(unequal$seq_better_mix20.5)
    #311

    library(viridis)

    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\ntime to resistance\nlonger','1'='sequential\ntime to resistance\nlonger'))) +

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

             #now add the unequal points
             #showed that they were on both sides
             #geom_point(data = unequal, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2), colour='red', shape=2) +

             #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.019330    0.560508

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

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