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

Investigation into effects of individual inputs on time to reach resistance thresholds.

version 2 - in progress

Based on 1000 runs with all inputs varying at once which explains the variability at particular input values.

Within each figure the plots are divided into 6 sub-plots according to the insecticide use strategy and when the resistance threshold is reached

1. insecticide 1 : sole use
2. insecticide 2 : sole use
3. Mixture 1 : threshold reached for either insecticide in mixture
4. Mixture 2 : once threshold reached for either insecticide in mixture, switch to sole use of other until it too reaches threshold
5. Mixture 3 : threshold reached for both insecticides in mixture
6. Sequential : sole use of one insecticide, switch to other when threshold reached until it too reaches threshold

These plots show gen_cP0.2, which is the number of generations to reach 20% resistance.

Red dashed lines are a smoothed mean.

Runs where resistance thresholds have not been reached have been removed.

Previous plots showed a set of points at 1000 generations. The model is run for 500 generations, any runs which have not reached the resistance threshold by this time are given a value of 1000. This has little effect on questions of whether mixtures or sequential use is better.

  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'

  #remember these don't have new exposure column
  load(file=paste0(outFolder,'listOutMix_ex2_10000.rda'))
  load(file=paste0(outFolder,'listOutI1_ex2_10000.rda'))
  load(file=paste0(outFolder,'listOutI2_ex2_10000.rda'))  

# ex10000 runs showed effect of selection    
#   load(file=paste0(outFolder,'listOutMix_ex10000.rda'))
#   load(file=paste0(outFolder,'listOutI1_ex10000.rda'))
#   load(file=paste0(outFolder,'listOutI2_ex10000.rda'))

### chunk copied from sensiAnPaper1All
### now want to calculate the 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  
  resistPointsMix1 <- 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
  resistPointsMix2 <- 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  
  resistPointsMix3 <- 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

# 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) 
  resistPointsMix1_T <- addStrategyColumn(resistPointsMix1,"Mixture 1",ggInput) 
  resistPointsMix2_T <- addStrategyColumn(resistPointsMix2,"Mixture 2",ggInput) 
  resistPointsMix3_T <- addStrategyColumn(resistPointsMix3,"Mixture 3",ggInput)     

  ggInsOuts <- rbind( resistPointsI1_T, resistPointsI2_T, resistPointsSeq_T, 
                      resistPointsMix1_T, resistPointsMix2_T, resistPointsMix3_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"

Plot time to resistance for each strategy using all sensitivity analysis values

names_results <- c("time_to_resistance0.5","time_to_resistance0.25","time_to_resistance0.1")

for(i in names_results)
{
  #x11()

  print( ggplot(ggInsOuts, aes_string(x='strategy',y=i, color='strategy')) + 
      #ylim(0,450) +
      coord_cartesian(xlim=c(1, 6), ylim=c(0, 350)) +
      geom_violin(draw_quantiles = c(0.25, 0.75), show.legend=FALSE)
  )#end print
}  

  #can I facet by exposure ?
  #exposure <- runif(1, min=0.1, max=0.9)
  #why does it not look uniform ? must be because runs have been removed if too high
  hist(ggInsOuts$exposure)
  #create exposure bins
  ggInsOuts$exposure_cat <- cut(ggInsOuts$exposure,breaks=c(0.1,0.3,0.5,0.7,0.9), labels=c("exposure 0.1-0.3","exposure 0.3-0.5","exposure 0.5-0.7","exposure 0.7-0.9"))
  #or just 2
  #ggInsOuts$exposure_cat <- cut(ggInsOuts$exposure,breaks=c(0.1,0.5,0.9), labels=c("exposure 0.1-0.5","exposure 0.5-0.9"))
  #actually might be better to show just the 2 extremes
  #ggSubset <- subset(ggInsOuts, exposure_cat==c("exposure 0.1-0.3","exposure 0.7-0.9"))

  print( ggplot(ggInsOuts, aes_string(x='strategy',y=i, color='strategy')) + 
      #ylim(0,450) +
      #scale_y_continuous( limits = c( 0,450 ) ) +
      #to zoom in without changing data  
      coord_cartesian(xlim=c(1, 6), ylim=c(0, 350)) +
      facet_wrap( ~ exposure_cat, scales='fixed') +  
      geom_violin(draw_quantiles = c(0.25, 0.75), show.legend=FALSE, adjust = .5 ) #, scale='count')
  )#end print


# or trying as a dotplot, with these options it's very similar to violin
#   print( ggplot(ggInsOuts, aes_string(x='strategy',y=i, color='strategy')) + 
#       #ylim(0,450) +
#       scale_y_continuous( limits = c( 0,450 ) ) +
#       facet_wrap( ~ exposure_cat, scales='fixed') +  
#       geom_dotplot(stackdir = 'center', binaxis='y', binwidth = 5 ,dotsize = 0.1, show.legend=FALSE) #, scale='count')
#   )#end print  


#now try the same faceting by effectiveness_ins1
#doesn't really show much ...
#so why does effectiveness come out so important in the trees ??  
  #hist(ggInsOuts$effectiveness_ins1)  
  #create effect1 bins
  ggInsOuts$effect1_cat <- cut(ggInsOuts$effectiveness_ins1,breaks=c(0.4,0.6,0.8,1), labels=c("effect1 0.4-0.6","effect1 0.6-0.8","effect1 0.8-1"))    

  print( ggplot(ggInsOuts, aes_string(x='strategy',y=i, color='strategy')) + 
      #ylim(0,450) +
      #scale_y_continuous( limits = c( 0,450 ) ) +
      #to zoom in without changing data  
      coord_cartesian(xlim=c(1, 6), ylim=c(0, 250)) +
      facet_wrap( ~ effect1_cat, scales='fixed', nrow=1) +  
      geom_violin(draw_quantiles = c(0.25, 0.75), show.legend=FALSE, adjust = .5 ) #, scale='count')
  )#end print  

## BEWARE what am I doing with putting mixture better than seq on each repeated input ....  

## !!try facetting by whether mixture or sequence better
## initially just with the time to resistance
## relies on ggInsOuts being prepared slightly differently from sensiAnPaper1All
## aha! this is weird it actually shows much shorter times for mixture1, & mix2 but sequential fairly similar  
#   print( ggplot(ggInsOuts, aes_string(x='strategy',y=i, color='strategy')) + 
#       #ylim(0,450) +
#       #scale_y_continuous( limits = c( 0,450 ) ) +
#       #to zoom in without changing data  
#       coord_cartesian(xlim=c(1, 6), ylim=c(0, 250)) +
#       facet_wrap( ~ betterMix2Seq20_cP0.1, scales='fixed', nrow=1) +  
#       geom_violin(draw_quantiles = c(0.25, 0.75), show.legend=FALSE, adjust = .5 ) #, scale='count')
#   )#end print    

# or I could look at the values for exposure on the y axis
#   print( ggplot(ggInsOuts, aes_string(x='strategy',y='exposure', color='strategy')) + 
#       #ylim(0,450) +
#       #scale_y_continuous( limits = c( 0,450 ) ) +
#       #to zoom in without changing data  
#       #coord_cartesian(xlim=c(1, 6), ylim=c(0, 250)) +
#       facet_wrap( ~ betterMix2Seq20_cP0.1, scales='fixed', nrow=1) +  
#       geom_violin(draw_quantiles = c(0.25, 0.75), show.legend=FALSE, adjust = .5 ) #, scale='count')
#   )#end print  

trying a PRCC analysis using the package sensitivity

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

  library(sensitivity)

# pcc(X, y, rank = FALSE, nboot = 0, conf = 0.95)
# 
# ## S3 method for class 'pcc'
# plot(x, ylim = c(-1,1), ...)
# 
# 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.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) +
             ylim(-1,1) +
             ggtitle(paste(strategy,"PRCC")) +
             xlab(NULL)
          )
  }

trying to highlight differences between sequential and mixture strategies in response of time-to-resistance to inputs

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

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

  for(i in names_inputs)
  {
    #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.5) +
             #facet_wrap( ~ strategy) +
             #geom_smooth(aes_string(x=i, y=y, color=NULL)) )
             labs(title = i))
  }

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.25'    
    #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))
  }

Looking more closely at ratio between the level of resistance to the first insecticide divided by that of the second one. Restricting the x axis to lower values.

?not very useful

  i <- 'resist_start_1_div_2'
  color <- "start_freq_allele2"

  print( ggplot(ggInsOuts, aes_string(x=i, y=y, color=color)) + 
           geom_point(shape=3) + 
           facet_wrap( ~ strategy) + 
           #geom_smooth(aes_string(x=i, y=y, color=NULL)) )
           geom_smooth(colour='red', linetype='dashed',size=0.5) +
           labs(title = i) +
           xlim(0,100)
         )

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