In AndySouth/resistance: an insecticide resistance population genetics model with 2 loci and 2 insecticides
code to start the sensitivity analysis for the first paper
main analyses of results is in : paper1_results_figs_slimmed_50_rr
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'))
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
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')
calculate the difference between sequential & mixture scenarios
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'))
finding when mixtures are >20% better
repeating calculation of the difference between sequential & mixture scenarios
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
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')
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)
do tree analysis
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')
}
testing tree analysis for paper
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)
AndySouth/resistance documentation built on Nov. 12, 2020, 3:39 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
code to start the sensitivity analysis for the first paper
main analyses of results is in : paper1_results_figs_slimmed_50_rr
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'))
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 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')
calculate the difference between sequential & mixture scenarios
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'))
finding when mixtures are >20% better
repeating calculation of the difference between sequential & mixture scenarios
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
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') 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)
do tree analysis
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') }
testing tree analysis for paper
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)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.