#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
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 }
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') ) }
#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') ) }
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") ) }
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") ) }
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) )
#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) )
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) )
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) ) }
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 )
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.