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