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

Main Figures for paper 1 : version 3 - in progress

This file is automagically generated by R from paper1_results_figs_slimmed.Rmd

Notable results

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.
3. Exposure and dominance have the greatest on time-to-resistance as shown in PRCC analyses.
4. 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.
5. Exposure and effectiveness come out as the most important effect on the difference between the sequential and adpative mixture strategies.
6. 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.
7. The implication of additional operational costs of adaptive mixtures can be investigated. If a requirement of a 20% improvement in resistance is imposed then more runs favour sequential than mixture.

  library(resistance)

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

  ## to load previously saved runs
#   load(file=paste0(outFolder,'listOutMix_1000.rda'))
#   load(file=paste0(outFolder,'listOutI1_1000.rda'))
#   load(file=paste0(outFolder,'listOutI2_1000.rda'))

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

  load(file=paste0(outFolder,'listOutMix_ex2_10000.rda'))
  load(file=paste0(outFolder,'listOutI1_ex2_10000.rda'))
  load(file=paste0(outFolder,'listOutI2_ex2_10000.rda'))

  # load(file=paste0(outFolder,'listOutMix_ex2_1000.rda'))
  # load(file=paste0(outFolder,'listOutI1_ex2_1000.rda'))
  # load(file=paste0(outFolder,'listOutI2_ex2_1000.rda'))

# very quick test data
  # load(file=paste0(outFolder,'listOutMix_ex2_3.rda'))
  # load(file=paste0(outFolder,'listOutI1_ex2_3.rda'))
  # load(file=paste0(outFolder,'listOutI2_ex2_3.rda'))

### calculate times to reach critical points and add them on to the input file 
### for different insecticide strategies, sequential, mix1 and mix2  


  #1) sequential : time to resistance for each insecticide in isolation
  #inputs : inAndOutI1, inAndOutI2
  #find time to criticalPoint for insecticide1
  #find time to criticalPoint for insecticide2
  #add together
  resistPointsI1 <- findResistancePoints(listOutI1, locus=1)
  resistPointsI2 <- findResistancePoints(listOutI2, locus=2)  
  resistPointsSeq <- resistPointsI1 + resistPointsI2


  #2) mixture1 : time to resistance for either insecticide when used in a mixture
  #inputs : inAndOutMix
  #find time to criticalPoint for EITHER insecticide in mixture  
  resistPointsMix_1 <- findResistancePoints(listOutMix, locus='either')  
  #todo - to be comparable I think this should be for when resistance to BOTH insecticides is reached

  #3) mixture2 : when resistance to one insecticide in the mixture reached, switch to sole use of the 
  #   other until that too reaches the critical point. Record total time.
  # what I actually need to do is start with mixture find the first critical point
  # (need to know which of the insecticides it is)
  # then I need to go to the single run for the other insecticide & starting at 
  # it's current resistance point find out how many more generations to go
  #inputs : inAndOutI1, inAndOutI2, inAndOutMix
  resistPointsMix_A <- findResistancePointsMixResponsive(listOutMix, listOutI1, listOutI2)

  #4) mixture3 : time to resistance for both insecticides when used in a mixture
  #inputs : inAndOutMix
  #find time to criticalPoint for BOTH insecticide in mixture  
  resistPointsMix_2 <- findResistancePoints(listOutMix, locus='both')

### chunk copied from part of one in sensiAnPaper1All
### bind results onto input file

  treeInput <- listOutMix$input

  #input files in listOutMix, listOutIn1 & listOutI2 are the same if the runs are done with default randomSeed 
  #except that exposure will be in a.f_AB, a.f_A0 and a.f_B0 respectively (& a.m* too)
  #I could just rename one to exposure
  #BEWARE risk if future changes
  #1/2/16 don't need this now because I've saved exposure from the single original random value
  #rownames(treeInput)[rownames(treeInput)=="a.f_AB"] <- "exposure"


  #hardcode which variables to include in analysis to keep it simple and transparent
  treePredictors <- c('P_1','P_2','exposure','phi.SS1_A0','phi.SS2_0B','h.RS1_A0','h.RS2_0B','s.RR1_A0','s.RR2_0B')

  #add these for extended analysis
  if (experiment=='extended')
    treePredictors <- c(treePredictors,'male_exposure_prop','correct_mix_deploy')

  treeInput <- treeInput[ treePredictors, ]

  #add an extra predictor, the lower starting freq of resistance divided by the larger
  resist_start_lo_div_hi <- ifelse( treeInput['P_1',] < treeInput['P_2',], treeInput['P_1',]/treeInput['P_2',], treeInput['P_2',]/treeInput['P_1',])
  treeInput <- rbind(treeInput,resist_start_lo_div_hi)    

  #20160122 add test for Ian of resistance1/resistance2
  resist_start_1_div_2 <- treeInput['P_1',]/treeInput['P_2',]
  treeInput <- rbind(treeInput,resist_start_1_div_2)   

  #renaming other rownames to make nicer plots
  rownames(treeInput)[rownames(treeInput)=="phi.SS1_A0"] <- "effectiveness_ins1"
  rownames(treeInput)[rownames(treeInput)=="phi.SS2_0B"] <- "effectiveness_ins2"

  rownames(treeInput)[rownames(treeInput)=="P_1"] <- "start_freq_allele1"
  rownames(treeInput)[rownames(treeInput)=="P_2"] <- "start_freq_allele2"

  rownames(treeInput)[rownames(treeInput)=="h.RS1_A0"] <- "dominance_allele1"
  rownames(treeInput)[rownames(treeInput)=="h.RS2_0B"] <- "dominance_allele2"

  rownames(treeInput)[rownames(treeInput)=="s.RR1_A0"] <- "selection_coef_allele1"
  rownames(treeInput)[rownames(treeInput)=="s.RR2_0B"] <- "selection_coef_allele2"


  #get the inputs here to use in ggplot investigation of model responses to inputs
  #used in a later chunk
  ggInput <- treeInput

  #create a curtisInputs dataframe for use later
  curtisInputs <- data.frame("exposure"=0.9,
                             "effectiveness_ins1"=0.73,
                             "effectiveness_ins2"=1,
                             "start_freq_allele1"=0.01,
                             "start_freq_allele2"=0.01,
                             "dominance_allele1"=0.17,
                             "dominance_allele2"=0.0016,
                             "selection_coef_allele1"=0.23,
                             "selection_coef_allele2"=0.43,
                             "male_exposure_prop"=1,
                             "correct_mix_deploy"=1,
                             "resist_start_lo_div_hi"=1,
                             "resist_start_1_div_2"=1)

# getting data into format for ggplot


  #uses ggInput calculated in earlier chunk

  #first needs to transpose rows to cols
  #ggInput_T <- t(ggInput)

  #function to transpose and add strategy column with the passed value
  #also add on the input columns  
  addStrategyColumn <- function(x, value, inputs){
    x <- as.data.frame( t(x) )
    x$strategy <- value

    #inputs <- as.numeric(inputs)
    #transpose inputs
    inputs <- as.data.frame( t(inputs), stringsAsFactors=FALSE )
    #cbind onto the outputs
    x <- cbind(inputs,x)
    x
  }    

  resistPointsI1_T <- addStrategyColumn(resistPointsI1,"insecticide 1",ggInput) 
  resistPointsI2_T <- addStrategyColumn(resistPointsI2,"insecticide 2",ggInput)   
  resistPointsSeq_T <- addStrategyColumn(resistPointsSeq,"Sequential",ggInput) 
  resistPointsMix_1_T <- addStrategyColumn(resistPointsMix_1,"Mix either",ggInput) 
  resistPointsMix_A_T <- addStrategyColumn(resistPointsMix_A,"Mix adaptive",ggInput) 
  resistPointsMix_2_T <- addStrategyColumn(resistPointsMix_2,"Mix both",ggInput)     

  ggInsOuts <- rbind( resistPointsI1_T, resistPointsI2_T, resistPointsSeq_T, 
                      resistPointsMix_1_T, resistPointsMix_A_T, resistPointsMix_2_T)    

  #remove runs that didn't reach a resistance threshold (999), even if just for one strategy
  #>1000 excludes sequential strategy that had a 999 in
  #BEWARE would need to increase 1000 if I increase max generations in the runs
  didntReachThresh <- which(ggInsOuts$gen_cP0.5 == 999 | ggInsOuts$gen_cP0.5 > 500 |
                            ggInsOuts$gen_cP0.25 == 999 | ggInsOuts$gen_cP0.25 > 500 |
                            ggInsOuts$gen_cP0.1 == 999 | ggInsOuts$gen_cP0.1 > 500 )

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


  #prettify output names
  names(ggInsOuts)[names(ggInsOuts)=="gen_cP0.5"] <- "time_to_resistance0.5"
  names(ggInsOuts)[names(ggInsOuts)=="gen_cP0.25"] <- "time_to_resistance0.25"
  names(ggInsOuts)[names(ggInsOuts)=="gen_cP0.1"] <- "time_to_resistance0.1"

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

  #format data in a different way to enable PRCC on difference between sequential & mix2
  #resistPointsMix_A_T
  #resistPointsSeq_T

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

  dif_mixA_seq <- cbind(resistPointsMix_A_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.25 > 500 | dif_mixA_seq$gen_cP0.25seq > 500 )

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


  #dif_mixA_seq$mixA_minus_seq0.25 <- dif_mixA_seq["gen_cP0.25"] - dif_mixA_seq["gen_cP0.25seq"]
  dif_mixA_seq["mixA_minus_seq0.25"] <- dif_mixA_seq["gen_cP0.25"] - dif_mixA_seq["gen_cP0.25seq"]

  #TODO : later do PRCC on dif_mixA_seq["mixA_minus_seq0.25"]

\pagebreak

Fig. x1 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.25")

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
}  

\pagebreak

Fig x2 PRCC analysis for time-to-resistance 0.25 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 <- strategies[-1:-2]

  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.25']

    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

    #draw_plot(p,x=0,y=0,width=.5,height = .5)

    # print( ggdraw() + #panel_border(remove = TRUE) +
    #          draw_plot( ggplot( to_plot, aes_string(x='inputs',y='PRCC') ) +     
    # #print( ggplot( to_plot, aes_string(x='inputs',y='PRCC') ) + 
    #          geom_point(shape=1, colour='red') +
    #          theme(axis.text.x = axis_text_x) +
    #          #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.25")) +
    #          xlab(NULL)
    #       , height= height))

    plotlist[[strategy_num]] <- ggplot( to_plot, aes_string(x='inputs',y='PRCC') ) +     
             geom_point(shape=1, colour='red') +
             theme(axis.text.x = axis_text_x) +
             geom_hline(yintercept = 0, linetype=3) +
             #theme_bw() + # to add in gridline but messes up other bits of plot
             #theme(panel.grid.major.y = element_line(color = "grey", size = 1)) +
             ylim(-1,1) +
             ggtitle(paste(strategy)) +
             xlab(NULL) 
  }

    plot_grid( plotlist[[1]],plotlist[[2]],plotlist[[3]],plotlist[[4]],ncol=1, rel_heights=c(1,1,1,1.8))  

\pagebreak

Fig. x3 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.

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

  #remove insecticide1 & 2 strategies
  ggSubset <- ggInsOuts[ !ggInsOuts$strategy %in% c("insecticide 1","insecticide 2"), ]
  #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")
  out_cols <- c("time_to_resistance0.25","strategy")
  ggSubset <- ggSubset[, c(in_cols,out_cols) ]

  #num_inputs <- length(in_cols)
  #names_inputs <- in_cols

  #to rearrange order of legend entries
  ggSubset$strategy <- factor(ggSubset$strategy, levels=c("Mix both", "Mix adaptive", "Sequential", "Mix either")) #, labels=c("MBB", "MAA", "MCC"))

  #TODO could I facet the plots below rather than using the loop
  #probably need to put all inputs into a single column and then facet by that

  #for(i in names_inputs)
  for(i in in_cols)
    {
    #x11()

    y <- 'time_to_resistance0.25'    

    #print( ggplot(ggInsOuts, aes_string(x=i, y=y, colour="strategy")) + 
    print( ggplot(ggSubset, aes_string(x=i, y=y, colour="strategy")) + 
            #points not wanted if 10000
             #geom_point(shape=3, show.legend=FALSE) + 
             #geom_smooth(colour='red', linetype='dashed',size=0.5) +
             geom_smooth(linetype='dashed',size=1) +
             #facet_wrap( ~ strategy) +
             #geom_smooth(aes_string(x=i, y=y, color=NULL)) )
             #labs(title = i) +
             #add line for Curtis inputs
             geom_vline(data = curtisInputs, aes_string(xintercept = i), linetype='dotted', colour='red') +
             #tryting to get the vline into the legend fails
             #geom_vline(data = curtisInputs, aes_string(xintercept = i, colour='"Curtis Fig.2"' ), linetype='dotted', show.legend=TRUE) +
             geom_text(data = curtisInputs, inherit.aes = FALSE, aes_string(x = i, y=210, label='"Curtis Fig.2"'), size=3, colour='red')

           )
  }

\pagebreak

Fig. x4 PRCC difference in time-to-resistance 0.25 between sequential and adaptive mixture strategies

1. Exposure and effectiveness come out as the most important effect on the difference between the sequential and adpative mixture strategies.

  library(sensitivity)


    x <- dif_mixA_seq[,1:num_inputs]
    y <- dif_mixA_seq["mixA_minus_seq0.25"]

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

\pagebreak

Fig. x5 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.25 <- ifelse( ggSubset$mixA_minus_seq0.25 > 0, 1, 0 )
  ggSubset$seq_better_mixA0.25 <- ifelse( ggSubset$mixA_minus_seq0.25 < 0, 1, 0 )
  ggSubset$abs_diff <- abs(ggSubset$mixA_minus_seq0.25)

  #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.25 <- 1
  df_same$seq_better_mixA0.25 <- 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$mixA_not20pc_better_seq0.25 <- ifelse( ggSubset$mixA_minus_seq0.25 / ggSubset$gen_cP0.25seq <= 0.2, 1, 0 )   

    #COOLIO !! I really like this figure.
    library(viridis)
    #facet showing generation differences, labelled mix better, seq better
    #+ve & -ve on the same colour scale
    #failed attempt at getting Curtis point in legend
    #print( ggplot(ggSubset, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2, colour=abs_diff, size="Curtis"))+
    print( ggplot(ggSubset, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2, colour=abs_diff))+

             facet_grid( ~ seq_better_mixA0.25, labeller=as_labeller( c('0'='adaptive mix\nresistance slower','1'='sequential\nresistance slower'))) +

             #YlGnBu looks really good except that the yellow points are very faint
             #scale_colour_distiller(palette='YlGnBu', direction=1) +
             #option="plasma" looks good too
             scale_color_viridis(begin=1, end=0, option="viridis", 
                                 guide=guide_colourbar(title="difference in \ngenerations to reach\nresistance threshold", direction='horizontal')) + 
             theme_bw() +
             theme(legend.position = "bottom") +
             #theme_minimal() +
             #add point for Curtis inputs
             geom_point(data = curtisInputs, aes(x=exposure, y=effectiveness_ins1 + effectiveness_ins2), colour='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) + 
             #!!!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='lightgrey', 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='lightgrey', 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='red')
           )       

Grey 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 grey line is a linear model through these same points.

\pagebreak

Fig. x6 Modifying the previous plot to show where resistance arises over 20% slower in an adaptive mixture strategy 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.

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

Fig. x7 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.

  #code copied from sensiAnPaper1All

  require('plyr')

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

  #just doing for a few strategies, sensiAnPpaer1All.Rmd does for a bunch
  treeResponses <- c( rownames(resistBetterMix_ASeq)[2], #2 is 25%
                      rownames(resistBetterMix_ASeq20)[2]
                    )


  #to do trees & plots for all response variables
  for( treeResponse in treeResponses )
  {
    tree <- rpart::rpart(as.formula(paste(treeResponse,"~",treePredictorString)), data = treeInput, method = 'class') 
    #http://stackoverflow.com/questions/29197213/what-is-the-difference-between-rel-error-and-x-error-in-a-rpart-decision-tree
    #A rule of thumb is to choose the lowest level where the rel_error + xstd < xerror.

    cpOptimal <- tree$cptable[ tree$cptable[,"rel error"] + tree$cptable[,"xstd"]  < tree$cptable[,"xerror"],"CP"][1]

    treePruned <- prune(tree, cp=cpOptimal)

    #to set box colours
    cols <- ifelse(treePruned$frame$yval == 1, "red", "green3")

    #adds in no. scenarios correctly and incorrectly classified
    #prp(tree, extra=1, varlen = 10, main=treeResponse)
    #prp(treePruned, extra=1, varlen = 10, main=treeResponse, col=cols, border.col=cols)
    #prp(treePruned, extra=1, varlen = 0, main=treeResponse, box.col=cols, under=TRUE)
    prp(treePruned, extra=1, varlen = 0, main='', box.col=cols, under=TRUE)
    }

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