The main question is under which scenarios are mixtures better than sequential ?
#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 with exception of exposure nScenarios <- 10 #nScenarios <- 100 #nScenarios <- 500 #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 ) ## 23/3/16 expt with an 'old' insecticide # listOutMix <- sensiAnPaperPart( nScenarios, insecticideUsed = 'mixture', experiment=experiment, old_insecticide=TRUE ) # listOutI1 <- sensiAnPaperPart( nScenarios, insecticideUsed = 'insecticide1', experiment=experiment, old_insecticide=TRUE ) # listOutI2 <- sensiAnPaperPart( nScenarios, insecticideUsed = 'insecticide2', experiment=experiment, old_insecticide=TRUE ) ## save results objects as rda for re-use save(listOutMix, file=paste0(outFolder,'listOutMix_rr_10.rda')) save(listOutI1, file=paste0(outFolder,'listOutI1_rr_10.rda')) save(listOutI2, file=paste0(outFolder,'listOutI2_rr_10.rda'))
## 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')) load(file=paste0(outFolder,'listOutMix_ex2_10000.rda')) load(file=paste0(outFolder,'listOutI1_ex2_10000.rda')) load(file=paste0(outFolder,'listOutI2_ex2_10000.rda'))
#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') #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')
require('plyr') ## is mix1 better than sequential ? #T/F resistBetterMix_1SeqBoolean <- resistPointsMix_1 > resistPointsSeq #convert to 0/1 resistBetterMix_1Seq <- plyr::aaply(resistBetterMix_1SeqBoolean,.margins=c(1,2),.fun=as.numeric) #rename rows ready for binding on to inputs rownames(resistBetterMix_1Seq) <- gsub('gen','betterMix_1Seq', rownames(resistBetterMix_1Seq)) rowSums(resistBetterMix_1Seq)/ncol(resistBetterMix_1Seq) #new1000 #betterMix_1Seq_cP0.1 betterMix_1Seq_cP0.25 betterMix_1Seq_cP0.5 # 0.223 0.245 0.243 ## 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)) rowSums(resistBetterMix_ASeq)/ncol(resistBetterMix_ASeq) #new1000 #betterMix_ASeq_cP0.1 betterMix_ASeq_cP0.25 betterMix_ASeq_cP0.5 # 0.558 0.563 0.534 ## is mix2 better than sequential ? #T/F resistBetterMix_2SeqBoolean <- resistPointsMix_2 > resistPointsSeq #convert to 0/1 resistBetterMix_2Seq <- plyr::aaply(resistBetterMix_2SeqBoolean,.margins=c(1,2),.fun=as.numeric) #rename rows ready for binding on to inputs rownames(resistBetterMix_2Seq) <- gsub('gen','betterMix_2Seq', rownames(resistBetterMix_2Seq)) #tally runs where mix3 better than seq rowSums(resistBetterMix_2Seq)/ncol(resistBetterMix_2Seq) #new1000 #betterMix_2Seq_cP0.1 betterMix_2Seq_cP0.25 betterMix_2Seq_cP0.5 # 0.820 0.747 0.677 #check1 #plot the 2 individual results & the mixture, just to check what is going on #x11() # plotallele.freq.andy(listOutMix) # #x11() # plotallele.freq.andy(listOutI1) # #x11() # plotallele.freq.andy(listOutI2) #I should be able to compare single scenarios #yes this works, but not too useful # plotallele.freq.andy( list(results=list(listOutMix$results[[1]], listOutI1$results[[1]] )) ) #this is a modifdeied version of Beths function #allowing curtis fig2 to be applied to any scenario #plotcurtis_f2_generic( listOutMix$results[[1]], listOutI2$results[[1]], listOutI1$results[[1]] ) #check2 # #what is the absolute difference between mix2 & sequential ? # resistDifferenceMix_2Seq <- resistPointsMix_2 - resistPointsSeq # # #check3 is their a difference between the mix & the individual strategies # #looks mostly plausible, occasionally there are 0s # resistDifferenceMix_1I1 <- resistPointsMix_1 - resistPointsI1 # resistDifferenceMix_1I2 <- resistPointsMix_1 - resistPointsI2 # comparing sensitivity inputs to those used by Beth to reproduce curtis fig2 # inputCurtis <- createInputMatrix(FALSE) # inputSensiMix <- listOutMix$input # # input <- cbind(inputCurtis, inputSensiMix) # write.csv(input, file=paste0(outFolder,'inputCurtisSensiCompare.csv'))
require('plyr') ## is mix1 better than sequential ? #T/F resistBetterMix_1SeqBoolean20 <- resistPointsMix_1 > (1.2 * resistPointsSeq) #convert to 0/1 resistBetterMix_1Seq20 <- plyr::aaply(resistBetterMix_1SeqBoolean20,.margins=c(1,2),.fun=as.numeric) #rename rows ready for binding on to inputs rownames(resistBetterMix_1Seq20) <- gsub('gen','betterMix_1Seq20', rownames(resistBetterMix_1Seq20)) rowSums(resistBetterMix_1Seq20)/ncol(resistBetterMix_1Seq20) #1000runs #betterMix_1Seq20_cP0.1 betterMix_1Seq20_cP0.25 betterMix_1Seq20_cP0.5 # 0.136 0.175 0.173 #new1000 #betterMix_1Seq20_cP0.1 betterMix_1Seq20_cP0.25 betterMix_1Seq20_cP0.5 # 0.163 0.178 0.182 ## is mix2 better than sequential ? #T/F 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)) rowSums(resistBetterMix_ASeq20)/ncol(resistBetterMix_ASeq20) #1000runs # betterMix_ASeq20_cP0.1 betterMix_ASeq20_cP0.25 betterMix_ASeq20_cP0.5 # 0.241 0.326 0.313 #new1000 # betterMix_ASeq20_cP0.1 betterMix_ASeq20_cP0.25 betterMix_ASeq20_cP0.5 # 0.283 0.302 0.283 ## is mix3 better than sequential ? #T/F resistBetterMix_2SeqBoolean20 <- resistPointsMix_2 > (1.2 * resistPointsSeq) #convert to 0/1 resistBetterMix_2Seq20 <- plyr::aaply(resistBetterMix_2SeqBoolean20,.margins=c(1,2),.fun=as.numeric) #rename rows ready for binding on to inputs rownames(resistBetterMix_2Seq20) <- gsub('gen','betterMix_2Seq20', rownames(resistBetterMix_2Seq20)) #tally runs where mix3 better than seq rowSums(resistBetterMix_2Seq20)/ncol(resistBetterMix_2Seq20) #1000runs # betterMix_2Seq20_cP0.1 betterMix_2Seq20_cP0.25 betterMix_2Seq20_cP0.5 # 0.672 0.560 0.488 #new1000 # betterMix_2Seq20_cP0.1 betterMix_2Seq20_cP0.25 betterMix_2Seq20_cP0.5 # 0.622 0.536 0.458
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') 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" #subsetting now done above before names prettified #treePredictors <- c('P_1','P_2','startResistanceRatio','exposure','effectiveness_ins1','effectiveness_ins2','h.RS1_A0','h.RS2_0B','s.RR1_A0','s.RR2_0B') #treeInput <- treeInput[ treePredictors, ] #get the inputs here to use in ggplot investigation of model responses to inputs #used in a later chunk ggInput <- treeInput #replace 0s & 1s with mixture/sequence to make tree plot clearer repMix <- function(x){ x[x==0] <- 'sequence' x[x==1] <- 'mixture' x } resistBetterMix_1Seq <- repMix(resistBetterMix_1Seq) resistBetterMix_ASeq <- repMix(resistBetterMix_ASeq) resistBetterMix_2Seq <- repMix(resistBetterMix_2Seq) resistBetterMix_1Seq20 <- repMix(resistBetterMix_1Seq20) resistBetterMix_ASeq20 <- repMix(resistBetterMix_ASeq20) resistBetterMix_2Seq20 <- repMix(resistBetterMix_2Seq20) #bind results onto predictors treeInput <- rbind( treeInput, resistBetterMix_1Seq ) treeInput <- rbind( treeInput, resistBetterMix_ASeq ) treeInput <- rbind( treeInput, resistBetterMix_2Seq ) #add >20% results treeInput <- rbind( treeInput, resistBetterMix_1Seq20 ) treeInput <- rbind( treeInput, resistBetterMix_ASeq20 ) treeInput <- rbind( treeInput, resistBetterMix_2Seq20 ) #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(treePredictors, collapse="+") treePredictorString <- paste(colnames(treeInput)[input_columns], collapse="+") treeResponses <- c( rownames(resistBetterMix_1Seq), rownames(resistBetterMix_ASeq), rownames(resistBetterMix_2Seq), rownames(resistBetterMix_1Seq20), rownames(resistBetterMix_ASeq20), rownames(resistBetterMix_2Seq20) ) #to do trees & plots for all response variables for( treeResponse in treeResponses ) { x11() #png(paste0(outFolder,"tree",treeResponse,".png")) tree <- rpart::rpart(as.formula(paste(treeResponse,"~",treePredictorString)), data = treeInput, method = 'class') #treePruned <- prune(tree, cp=tree$cptable[which.min(tree$cptable[,"xerror"]),"CP"]) #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, col=cols, border.col=cols) prp(treePruned, extra=1, varlen = 0, main=treeResponse, box.col=cols, under=TRUE) #dev.off() #savePlot(paste0(outFolder,"tree",treeResponse,".png"), type='png') savePlot(paste0(outFolder,"tree",treeResponse,".jpg"), type='jpg') #trying to limit num levels a different way using maxdepth #tree <- rpart::rpart(as.formula(paste(treeResponse,"~",treePredictorString)), data = treeInput, method = 'class', maxdepth=3) #adds in no. scenarios correctly and incorrectly classified #prp(tree, extra=1, varlen = 10, main=treeResponse) #savePlot(paste0(outFolder,"tree",treeResponse,"depth3.jpg"), type='jpg') }
x11() treeResponse <- "betterMix_ASeq_cP0.5" tree <- rpart::rpart(as.formula(paste(treeResponse,"~",treePredictorString)), data = treeInput, method = 'class') treePruned <- prune(tree, cp=tree$cptable[which.min(tree$cptable[,"xerror"]),"CP"]) #to set box colours cols <- ifelse(treePruned$frame$yval == 1, "red", "green3") #adds in no. scenarios correctly and incorrectly classified #prp(treePruned, extra=1, varlen = 0, main=treeResponse, box.col=cols, under=TRUE) cpOptimal <- tree$cptable[ tree$cptable[,"rel error"] + tree$cptable[,"xstd"] < tree$cptable[,"xerror"],"CP"][1] treePruned <- prune(tree, cp=cpOptimal) prp(treePruned, extra=1, varlen = 0, main=treeResponse, box.col=cols, under=TRUE) #comparing it with 0/1 gives the same big tree treeInput$betterMix_ASeq_cP0.5numeric <- as.numeric(treeInput$betterMix_ASeq_cP0.5) treeResponse <- "betterMix_ASeq_cP0.5numeric" tree <- rpart::rpart(as.formula(paste(treeResponse,"~",treePredictorString)), data = treeInput, method = 'class') treePruned <- prune(tree, cp=tree$cptable[which.min(tree$cptable[,"xerror"]),"CP"]) #to set box colours cols <- ifelse(treePruned$frame$yval == 1, "red", "green3") prp(treePruned, extra=1, varlen = 0, main=treeResponse, box.col=cols, under=TRUE) #testing alternative pruning #A rule of thumb is to choose the lowest level where the rel_error + xstd < xerror. treePruned <- prune(tree, cp=tree$cptable[which.min(tree$cptable[,"xerror"]),"CP"])
# treeResponse <- 'betterMix_1Seq_cP0.5' # tree <- rpart::rpart(as.formula(paste(treeResponse,"~",treePredictorString)), data = treeInput, method = 'class') # # #if get error number of rows of matrices must match (see arg 2) # #is because all vals are 0 # # plot(tree) #plots tree # text(tree) #labels tree # #branch predictors 2 phis & P_2 # # #for mix3 (resistance to both) # x11() # treeResponse <- 'betterMix_2Seq_cP0.5' # tree <- rpart::rpart(as.formula(paste(treeResponse,"~",treePredictorString)), data = treeInput, method = 'class') # plot(tree) #plots tree # text(tree) #labels tree # # #branch predictors 2 phis & P_1 # # #or from rpart.plot # #prp(tree) # # #adds in no. scenarios correctly and incorrectly classified # prp(tree, extra=1, varlen = 10, main=treeResponse) # savePlot(paste0(outFolder,"tree",treeResponse,".png"), type='png') #prp(tree, extra=1,box.col=c("pink","palegreen3")[tree$frame$yval],cex=0.5,yesno.yshift=-2) #susanas code #pruning tree # treePruned <- prune(tree, cp=tree$cptable[which.min(tree$cptable[,"xerror"]),"CP"]) # plot(treePruned) #plots tree # text(treePruned) #labels tree # prp(prunned_tree, extra=1,box.col=c("pink", "palegreen3")[prunned_tree$frame$yval],cex=0.5,yesno.yshift=-2)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.