library(gridExtra) require(ggplot2) require(plyr) ptm <- proc.time() # PLOT PlotFigure1 = F PlotFigure2 = T PlotFigure3 = F PlotFigure4 = F # Plot PlotFigure1 PlotFigure1.alpha<- c(1,1,1) # Dirichlet priors PlotFigure1.f1<- 5 # costs and benefits of attacking the three prey types PlotFigure1.f2<- -0.5 PlotFigure1.f3<- -1 PlotFigure1.mutantTypeAdded = 2 PlotFigure1.Nmax = 50 PlotFigure1.Nmin = 10 PlotFigure1.stepNtot = 3 #10 PlotFigure1.Ntot = c(25,50,100) PlotFigure1.stepQHigh = 5 PlotFigure1.nbRepForward = 2000 # number of forward iteration PlotFigure1.propCommunityExperienced = 1 PlotFigure1.colorLow = "#ECECEC" PlotFigure1.colorHigh = "#000000" PlotFigure1.topLimit = 60 PlotFigure1.topLimitAttackRate = NA # Plot PlotFigure2 PlotFigure2.alpha<- c(1,1,1) # Dirichlet priors (alpha's) good (r1), bad (r2), very bad (r3) PlotFigure2.f1<- 5 # Costs and benefits of attacking the three prey types PlotFigure2.f2<- -0.5 PlotFigure2.f3<- -1 PlotFigure2.mutantTypeAdded = 2 PlotFigure2.n1 = 50 # prey mutant in total PlotFigure2.n2 = 10 # prey f3 PlotFigure2.n2bis = NA # prey f3 (second simulation) PlotFigure2.stepNbMutantNb = 15 PlotFigure2.nbRepForward = 2000 # number of forward iteration PlotFigure2.propCommunityExperienced = 1 PlotFigure2.yLim = c(0,1) PlotFigure2.yLimNbSampled = c(0,100) PlotFigure2.yLimNbSampledModel = c(0,50) # qB-Mimicry Arrow1 = F Arrow2 = F xArrow1 = c(0,10,20,29,39,48) yArrow1I = c(0.28,0.33,0.42,0.53,0.79,0.98) yArrow1E = c(0.2,0.19,0.18,0.17,0.16,0.15) colorArrow1 = "#A1A1A1" xArrow2 = c(0) yArrow2I = c(0.57) yArrow2E = c(0.33) colorArrow2 = "black" # B-Mimicry # Arrow1 = T # Arrow2 = T # xArrow1 = c(0,10) # yArrow1I = c(0.98,0.98) # yArrow1E = c(0.21,0.77) # colorArrow1 = "#A1A1A1" # xArrow2 = c(0) # yArrow2I = c(0.98) # yArrow2E = c(0.72) # colorArrow2 = "#A1A1A1" # qB-Mimicry : reverse # Arrow1 = T # Arrow2 = F # xArrow1 = c(0,3,6,9) # yArrow1I = c(0.57,0.68,0.98,0.98) # yArrow1E = c(0.33,0.3,0.29,0.27) # colorArrow1 = "#A1A1A1" # xArrow2 = c(1) # yArrow2I = c(0.98) # yArrow2E = c(0.80) # colorArrow2 = "#A1A1A1" # qB-Mimicry pessimistic predator # PlotFigure2.yLimNbSampled = c(0,8) # PlotFigure2.yLimNbSampledModel = c(0,4) # PlotFigure2.yLim = c(0.0,0.12) # # Arrow1 = T # Arrow2 = T # xArrow1 = c(0,10,20,29,39,48) # yArrow1I = c(0.097,0.097,0.097,0.093,0.79,0.98) # yArrow1E = c(0.065,0.069,0.065,0.067,0.16,0.15) # colorArrow1 = "#A1A1A1" # xArrow2 = c(0) # yArrow2I = c(0.57) # yArrow2E = c(0.33) # colorArrow2 = "black" # qB-Mimicry f1 = -5 f2 = -0.5 # PlotFigure2.yLimNbSampledModel = c(0,3) # PlotFigure2.yLim = c(0,0.5) # Arrow1 = T # Arrow2 = T # xArrow1 = c(0,10,20,29,39,48) # yArrow1I = c(0.18,0.21,0.24,0.28,0.34,0.46) # yArrow1E = c(0.06,0.06,0.07,0.07,0.07,0.07) # colorArrow1 = "#A1A1A1" # xArrow2 = c(0) # yArrow2I = c(0.57) # yArrow2E = c(0.33) # colorArrow2 = "black" # qB-Mimicry f1 = -5 f2 = -2 # PlotFigure2.yLimNbSampledModel = c(0,3) # PlotFigure2.yLim = c(0,0.5) # Arrow1 = T # Arrow2 = T # xArrow1 = c(0,10,20,29,39) # yArrow1I = c(0.073,0.072,0.095,0.091,0.15) # yArrow1E = c(0.046,0.044,0.041,0.039,0.039) # colorArrow1 = "#A1A1A1" # xArrow2 = c(0) # yArrow2I = c(0.57) # yArrow2E = c(0.33) # colorArrow2 = "black" # Plot PlotFigure3 PlotFigure3.alpha<- c(1,1,1) # Dirichlet priors (alpha's) good (r1), bad (r2), very bad (r3) PlotFigure3.f1<- 5 # costs and benefits of attacking the three prey types PlotFigure3.f2<- 0.05 PlotFigure3.f3<- -1 PlotFigure3.mutantTypeAdded = 2 PlotFigure3.n1 = 50 # prey mutant in total PlotFigure3.n2 = 50 # prey f3 PlotFigure3.stepNbMutantNb = 10 PlotFigure3.nbRepForward = 500 # number of forward iteration PlotFigure3.propCommunityExperienced = 1 # PlotFigure3.yLim = c(-0.05,1) PlotFigure3.yLim = c(-0.03,0.5) PlotFigure3.yLimNbSampled = c(0,NA) ## qB-Mimicry # LineFig3=0 # ArrowFig3 = T # ArrowFig3Bis = F # xArrowFig3F = c(-0.05,-0.05) # yArrowFig31I = c(0.48,0.52) # yArrowFig31E = c(0,1) # colorArrowFig31 = "black" # PointFig3 = F # XPoint = 0.5 # YPoint = -0.05 ## B-Mimicry LineFig3=0 ArrowFig3 = T xArrowFig3F = c(-0.05,-0.05) yArrowFig31E = c(0.48,0.52) yArrowFig31I = c(0,1) colorArrowFig31 = "black" ArrowFig3Bis = T xArrowFig3FBis = c(-0,0.2) yArrowFig31IBis = c(0.97,0.97) yArrowFig31EBis = c(0.25,0.75) colorArrowFig32 = "#A1A1A1" PointFig3 = T XPoint = 0.5 YPoint = -0.05 ##qB-Mimicry Pessimistic # LineFig3=0.055 # PlotFigure3.yLim = c(0.054,0.07) # # ArrowFig3 = T # xArrowFig3F = c(0.0541,0.0541,0.0541,0.0541,0.0541,0.0541) # yArrowFig31E = c(0.115,0.155,0.48,0.52,0.85,0.89) # yArrowFig31I = c(0,0.34,0.37,0.65,0.67,1) # colorArrowFig31 = "black" # # ArrowFig3Bis = T # xArrowFig3FBis = c(0,0.08) # yArrowFig31IBis = c(0.0612,0.0614) # yArrowFig31EBis = c(0.0602,0.0596) # colorArrowFig32 = "#A1A1A1" # # PointFig3 = T # XPoint = c(0.135,0.5,0.87) # YPoint = c(0.0541,0.0541,0.0541) ##qB-Mimicry f1 = -5 f2 = -0.5 # LineFig3=0.03 # PlotFigure3.yLim = c(0.028,0.06) # # # ArrowFig3 = T # xArrowFig3F = c(0.028,0.028,0.028,0.028) # yArrowFig31E = c(0.29,0.33,0.69,0.73) # yArrowFig31I = c(0,0.48,0.52,1) # colorArrowFig31 = "black" # # ArrowFig3Bis = T # xArrowFig3FBis = c(0,0.08,0.16,0.24) # yArrowFig31IBis = c(0.0445,0.0455,0.0465,0.0475) # yArrowFig31EBis = c(0.041,0.041,0.0406,0.0445) # colorArrowFig32 = "#A1A1A1" # # PointFig3 = T # XPoint = c(0.31,0.71) # YPoint = c(0.028,0.028) ##qB-Mimicry f1 = -5 f2 = -2 # LineFig3=0.03 # PlotFigure3.yLim = c(0.028,0.06) # # # ArrowFig3 = T # xArrowFig3F = c(0.028,0.028) # yArrowFig31E = c(0,1) # yArrowFig31I = c(0.48,0.52) # colorArrowFig31 = "black" # # ArrowFig3Bis = F # xArrowFig3FBis = c(0,0.08,0.16,0.24) # yArrowFig31IBis = c(0.0445,0.0455,0.0465,0.0475) # yArrowFig31EBis = c(0.041,0.041,0.0406,0.0445) # colorArrowFig32 = "#A1A1A1" # # PointFig3 = F # XPoint = c(0.31,0.71) # YPoint = c(0.028,0.028) # Plot PlotFigure4 PlotFigure4.rePlot = F PlotFigure4.saveTable = T PlotFigure4.loadTable = F # PlotFigure4.fileName = "QuasiBatesianBenefit5Priors111Densities30_30b" # PlotFigure4.fileName = "QuasiBatesianBenefit10Priors111Densities30_30b" # PlotFigure4.fileName = "QuasiBatesianBenefit5Priors111Densities30_30Var0" # PlotFigure4.fileName = "QuasiBatesianBenefit5Priors111Densities30_30Var1" # PlotFigure4.fileName = "QuasiBatesianBenefit5Priors111Densities30_30Var2" # PlotFigure4.fileName = "QuasiBatesianBenefit5Priors111Densities30_30Var05" # PlotFigure4.fileName = "QuasiBatesianBenefit5Priors111Densities30_30Var05bis" # with more strategies (30 instead of 10) # PlotFigure4.fileName = "BatesianBenefit5Priors111Densities30_30b" # PlotFigure4.fileName = "BatesianBenefit10Priors111Densities30_30b" # PlotFigure4.fileName = "QuasiBatesianBenefit5Priors111Densities50_50Var0" # PlotFigure4.fileName = "QuasiBatesianBenefit5Priors111Densities50_50Var2" # # PlotFigure4.fileName = c("QuasiBatesianBenefit5Priors111Densities50_50Var2BIS", # "QuasiBatesianBenefit5Priors111Densities50_50Var2BIS2") # # # # PlotFigure4.fileName = c("QuasiBatesianBenefit5Priors111Densities50_50Var1BISBIS", # "QuasiBatesianBenefit5Priors111Densities50_50Var1BISBIS2", # "QuasiBatesianBenefit5Priors111Densities50_50Var1BISBIS3", # "QuasiBatesianBenefit5Priors111Densities50_50Var1BISBIS4", # "QuasiBatesianBenefit5Priors111Densities50_50Var1BISBIS5" # # "QuasiBatesianBenefit5Priors111Densities50_50Var1BISBIS6", # # "QuasiBatesianBenefit5Priors111Densities50_50Var1BISBIS7", # # "QuasiBatesianBenefit5Priors111Densities50_50Var1BISBIS8", # # "QuasiBatesianBenefit5Priors111Densities50_50Var1BISBIS9", # # "QuasiBatesianBenefit5Priors111Densities50_50Var1BISBIS10" # ) # PlotFigure4.fileName = "QuasiBatesianBenefit5Priors111Densities50_50Var1bis" PlotFigure4.fileName = "New2Var0NegativeValuesRep8" # get only one portion CHOOSE limf2f3Min = -0.23 #Rep9 limf2f3Max = -0.20 limf2f3Min = -0.50 #Rep8 limf2f3Max = -0.22 # limf2f3Min = -0.50 #Rep1 Var2 # limf2f3Max = -0.2 # limf2f3Min = -0.08 #Rep4 Var2 # limf2f3Max = -0.04 # limf2f3Min = -0.11 #Rep5 Var2 # limf2f3Max = -0.08 # limf2f3Min = 0.0 #Rep2 Var2 # limf2f3Max = 0.04 # limf2f3Min = -0.04 #Rep3 Var2 # limf2f3Max = 0.0 PlotFigure4.alpha<- c(1,1,1) # Dirichlet priors (alpha's) good (r1), bad (r2), very bad (r3) PlotFigure4.f1<- 5 # benefit of attacking the three prey types PlotFigure4.mutantTypeAdded = 2 PlotFigure4.nbStepsGraph = 40 # 40 PlotFigure4.nbStepsGraphY = 61 # 80 PlotFigure4.f3Min = -5 PlotFigure4.f3Max = -0.1 PlotFigure4.f2f3Min = -0.5 PlotFigure4.f2f3Max = 1 PlotFigure4.fSd = 0 PlotFigure4.nbCombinationf = 1 PlotFigure4.n1 = 50 # prey f2 in total PlotFigure4.n2 = 50 # prey f3 in each pop PlotFigure4.stepNbMutant = 2 PlotFigure4.nbRepForward = 2000 # number of forward iteration 500 PlotFigure4.propCommunityExperienced = 1 PlotFigure4.yLim = c(0,NA) PlotFigure4.colorLow = "#ECECEC" PlotFigure4.colorHigh = "#000000" PlotFigure4.topLimit = 0.5 ########################################### g_legend<-function(a.gplot){ tmp <- ggplot_gtable(ggplot_build(a.gplot)) leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box") legend <- tmp$grobs[[leg]] return(legend)} # Function to create the optimal sampling strategy matrixes SamplingStrategyMatrix <- function(alpha,bigN,b,c1,c2){ S <- array(0, dim = c(bigN+1,bigN+1,bigN+1)) # first dimension r1, second dimension r2, third dimension r3 A <- array(-1, dim = c(bigN+1,bigN+1,bigN+1)) # first dimension r1, second dimension r2, third dimension r3 rownames(S)<- 0:bigN; colnames(S)<- 0:bigN; dimnames(S)[[3]] <- 0:bigN rownames(A)<- 0:bigN; colnames(A)<- 0:bigN; dimnames(A)[[3]] <- 0:bigN probe<- numeric(3) # probability of encountering each prey type for (n in (bigN-1):0){ for (r1i in (n):0){ for (r2i in (n-r1i):0) { ## r1, r2, r3 corresponds to the coordinates of r1i, r2i and r3i in the matrix r1<- r1i+1 r2<- r2i+1 r3i<- n-r1i-r2i r3<- r3i + 1 v1 = S[r1, r2, r3] sumv<- alpha[1]+r1i + alpha[2]+r2i + alpha[3]+r3i probe[1]<- (alpha[1]+r1i) / sumv probe[2]<- (alpha[2]+r2i) / sumv probe[3]<- (alpha[3]+r3i) / sumv v2<-probe[1]*(S[r1+1,r2,r3] + b) + probe[2]*(S[r1,r2+1,r3] - c1) + probe[3]*(S[r1,r2,r3+1] - c2) if (v1 >= v2) { A[r1, r2, r3]<-0 S[r1, r2, r3]<- v1 } else { A[r1, r2, r3]<-1 S[r1, r2, r3]<- v2 } } } } return(A) } if(PlotFigure1==TRUE){ qVec <- c() NtotVec <- c() attackRateVec <- c() attackRateVec1 <- c() attackRateVec2 <- c() attackRateVec3 <- c() NbPreyEatenVec <- c() NbPreyEatenVec1 <- c() NbPreyEatenVec2 <- c() NbPreyEatenVec3 <- c() popIndexVec <- c() # if 1 : bad prey, 2: very bad prey pop PlotFigure1.stepQHigh2 = 1/PlotFigure1.stepQHigh PlotFigure1.stepNtot2 = (PlotFigure1.Nmax - PlotFigure1.Nmin)/PlotFigure1.stepNtot PlotFigure1.qHigh = seq(0,1,by=PlotFigure1.stepQHigh2) # PlotFigure1.Ntot = seq(PlotFigure1.Nmin,PlotFigure1.Nmax,by=PlotFigure1.stepNtot2) for(Ntot in PlotFigure1.Ntot){ for(q in PlotFigure1.qHigh){ pop=1 PlotFigure1.n2 = floor(q * Ntot) if(PlotFigure1.n2==0){ PlotFigure1.n2 = 1 } if(PlotFigure1.n2==Ntot){ PlotFigure1.n2 = Ntot-1 } PlotFigure1.n1 = Ntot - PlotFigure1.n2 # community of prey # nbBadPrey = PlotFigure1.n1-n1 preyCommunity <- c(rep(2,PlotFigure1.n2),rep(PlotFigure1.mutantTypeAdded-1,PlotFigure1.n1)) # 0 : profitable prey, nbGoodPrey = length(which(preyCommunity==0)) nbBadPrey = length(which(preyCommunity==1)) nbVeryBadPrey = length(which(preyCommunity==2)) # preyCommunity <- c(rep(1,PlotFigure1.n1-n1)) # 0 : profitable prey, # nbVeryBadPrey = 0 # 1 : toxic prey (c1), 2 : toxic prey (c2) # forward iteration print(floor(PlotFigure1.propCommunityExperienced*length(preyCommunity))) PlotFigure1.bigN = length(preyCommunity) A = SamplingStrategyMatrix(PlotFigure1.alpha, floor(PlotFigure1.propCommunityExperienced*length(preyCommunity)), PlotFigure1.f1, -PlotFigure1.f2, -PlotFigure1.f3) for(rep in 1:PlotFigure1.nbRepForward){ preyCommunityOnePred <- sample(preyCommunity,size=floor(PlotFigure1.propCommunityExperienced*length(preyCommunity))) # encounter in random order predKnowledge.ri1 = 0 # profitable prey predKnowledge.ri2 = 0 # toxic prey (c1) predKnowledge.ri3 = 0 # toxic prey (c2) nbPreyEaten = 0 nbPreyEaten1 = 0 nbPreyEaten2 = 0 nbPreyEaten3 = 0 for(prey in preyCommunityOnePred){ optBehavior = A[predKnowledge.ri1+1, predKnowledge.ri2+1, predKnowledge.ri3+1] if(optBehavior==1){ # attack nbPreyEaten = nbPreyEaten + 1 if(prey==0){ predKnowledge.ri1 = predKnowledge.ri1 + 1 nbPreyEaten1 = nbPreyEaten1 + 1 }else if(prey==1){ predKnowledge.ri2 = predKnowledge.ri2 + 1 nbPreyEaten2 = nbPreyEaten2 + 1 }else if(prey==2){ predKnowledge.ri3 = predKnowledge.ri3 + 1 nbPreyEaten3 = nbPreyEaten3 + 1 } } } popIndexVec <- c(popIndexVec, pop) qVec <- c(qVec, 1-q) NtotVec <- c(NtotVec, Ntot) attackRateVec <- c(attackRateVec, nbPreyEaten/PlotFigure1.bigN) attackRateVec1 <- c(attackRateVec1, nbPreyEaten1/nbGoodPrey) attackRateVec2 <- c(attackRateVec2, nbPreyEaten2/nbBadPrey) attackRateVec3 <- c(attackRateVec3, nbPreyEaten3/nbVeryBadPrey) NbPreyEatenVec <- c(NbPreyEatenVec,nbPreyEaten) NbPreyEatenVec1 <- c(NbPreyEatenVec1,nbPreyEaten1) NbPreyEatenVec2 <- c(NbPreyEatenVec2,nbPreyEaten2) NbPreyEatenVec3 <- c(NbPreyEatenVec3,nbPreyEaten3) } } } print('ok') tableAttackRate = data.frame(q=qVec,Ntot=NtotVec, attackRate=attackRateVec,attackRate1=attackRateVec1,attackRate2=attackRateVec2,attackRate3=attackRateVec3,NbPreyEaten=NbPreyEatenVec,NbPreyEaten1=NbPreyEatenVec1,NbPreyEaten2=NbPreyEatenVec2,NbPreyEaten3=NbPreyEatenVec3,popIndex=popIndexVec) # tableAttackRate <- tableAttackRate[tableAttackRate$popIndex==1,] tableAttackRateMean = ddply(tableAttackRate, .(q,popIndex,Ntot), summarize, mean = mean(attackRate,na.rm=T), sd = sd(attackRate,na.rm=T), mean1 = mean(attackRate1,na.rm=T), sd1 = sd(attackRate1,na.rm=T), mean2 = mean(attackRate2,na.rm=T), sd2 = sd(attackRate2,na.rm=T), mean3 = mean(attackRate3,na.rm=T), sd3 = sd(attackRate3,na.rm=T), meanNb = mean(NbPreyEaten,na.rm=T), sdNb = sd(NbPreyEaten,na.rm=T), meanNb1 = mean(NbPreyEaten1,na.rm=T), sdNb1 = sd(NbPreyEaten1,na.rm=T), meanNb2 = mean(NbPreyEaten2,na.rm=T), sdNb2 = sd(NbPreyEaten2,na.rm=T), meanNb3 = mean(NbPreyEaten3,na.rm=T), sdNb3 = sd(NbPreyEaten3,na.rm=T)) tableAttackRateMean$sampleNb = NA for(i in 1:dim(tableAttackRateMean)[1]){ tableAttackRateMean$sampleNb[i] = length(which(tableAttackRate$q==tableAttackRateMean$q[i] & tableAttackRate$Ntot==tableAttackRateMean$Ntot[i])) } tableAttackRateMean$confidenceInterval <- qt(0.975,df=tableAttackRateMean$sampleNb-1)*tableAttackRateMean$sd/sqrt(tableAttackRateMean$sampleNb) tableAttackRateMean$confidenceInterval2 <- qt(0.975,df=tableAttackRateMean$sampleNb-1)*tableAttackRateMean$sd2/sqrt(tableAttackRateMean$sampleNb) tableAttackRateMean$confidenceInterval3 <- qt(0.975,df=tableAttackRateMean$sampleNb-1)*tableAttackRateMean$sd3/sqrt(tableAttackRateMean$sampleNb) tableAttackRateMean$confidenceIntervalNb <- qt(0.975,df=tableAttackRateMean$sampleNb-1)*tableAttackRateMean$sdNb/sqrt(tableAttackRateMean$sampleNb) tableAttackRateMean$confidenceIntervalNb2 <- qt(0.975,df=tableAttackRateMean$sampleNb-1)*tableAttackRateMean$sdNb2/sqrt(tableAttackRateMean$sampleNb) tableAttackRateMean$confidenceIntervalNb3 <- qt(0.975,df=tableAttackRateMean$sampleNb-1)*tableAttackRateMean$sdNb3/sqrt(tableAttackRateMean$sampleNb) # nb Prey of\nprey sampled\n with confidence interval # nbsampledPlot1 <- ggplot(tableAttackRateMean)+ # geom_rect(mapping=aes(xmin=q-PlotFigure1.stepQHigh2/2, xmax=q+PlotFigure1.stepQHigh2/2, ymin=Ntot-PlotFigure1.stepNtot2/2, ymax=Ntot+PlotFigure1.stepNtot2/2, # fill=meanNb, color=meanNb)) + # scale_fill_gradient(name="Number of\nprey sampled\n", low=PlotFigure1.colorLow, high=PlotFigure1.colorHigh, limits=c(0, PlotFigure1.topLimit))+ # scale_color_gradient(name="Number of\nprey sampled\n", low=PlotFigure1.colorLow, high=PlotFigure1.colorHigh, limits=c(0, PlotFigure1.topLimit))+ # xlab("Proportion of mimics")+ # ylab("Number of prey (N)")+ # theme(panel.background = element_blank(), # title = element_text(size=12), # axis.title.x = element_text(size=12), # axis.title.y = element_text(size=12), # legend.position="bottom")+ # ggtitle("All prey") # # nbsampledPlot2 <- ggplot(tableAttackRateMean)+ # geom_rect(mapping=aes(xmin=q-PlotFigure1.stepQHigh2/2, xmax=q+PlotFigure1.stepQHigh2/2, ymin=Ntot-PlotFigure1.stepNtot2/2, ymax=Ntot+PlotFigure1.stepNtot2/2, # fill=meanNb3, color=meanNb3)) + # scale_fill_gradient(name="Number of\nprey sampled\n(class 3)", low=PlotFigure1.colorLow, high=PlotFigure1.colorHigh, limits=c(0, PlotFigure1.topLimit))+ # scale_color_gradient(name="Number of\nprey sampled\n(class 3)", low=PlotFigure1.colorLow, high=PlotFigure1.colorHigh, limits=c(0, PlotFigure1.topLimit))+ # xlab("Proportion of mimics")+ # ylab("Number of prey (N)")+ # theme(panel.background = element_blank(), # title = element_text(size=12), # axis.title.x = element_text(size=12), # axis.title.y = element_text(size=12))+ # ggtitle("Highly unprofitable models") # # nbsampledPlot3 <- ggplot(tableAttackRateMean)+ # geom_rect(mapping=aes(xmin=q-PlotFigure1.stepQHigh2/2, xmax=q+PlotFigure1.stepQHigh2/2, ymin=Ntot-PlotFigure1.stepNtot2/2, ymax=Ntot+PlotFigure1.stepNtot2/2, # fill=meanNb2, color=meanNb2)) + # scale_fill_gradient(name="Number of\nprey sampled\n(class 2)", low=PlotFigure1.colorLow, high=PlotFigure1.colorHigh, limits=c(0, PlotFigure1.topLimit),na.value="red")+ # scale_color_gradient(name="Number of\nprey sampled\n(class 2)", low=PlotFigure1.colorLow, high=PlotFigure1.colorHigh, limits=c(0, PlotFigure1.topLimit),na.value="red")+ # xlab("Proportion of mimics")+ # ylab("Number of prey (N)")+ # theme(panel.background = element_blank(), # title = element_text(size=12), # axis.title.x = element_text(size=12), # axis.title.y = element_text(size=12))+ # ggtitle("Moderately unprofitable mimics") # # nbsampledPlot4 <- ggplot(tableAttackRateMean)+ # geom_rect(mapping=aes(xmin=q-PlotFigure1.stepQHigh2/2, xmax=q+PlotFigure1.stepQHigh2/2, ymin=Ntot-PlotFigure1.stepNtot2/2, ymax=Ntot+PlotFigure1.stepNtot2/2, # fill=meanNb1, color=meanNb1)) + # scale_fill_gradient(name="Number of\nprey sampled\n(class 2)", low=PlotFigure1.colorLow, high=PlotFigure1.colorHigh, limits=c(0, PlotFigure1.topLimit),na.value="red")+ # scale_color_gradient(name="Number of\nprey sampled\n(class 2)", low=PlotFigure1.colorLow, high=PlotFigure1.colorHigh, limits=c(0, PlotFigure1.topLimit),na.value="red")+ # xlab("Proportion of mimics")+ # ylab("Number of prey (N)")+ # theme(panel.background = element_blank(), # title = element_text(size=12), # axis.title.x = element_text(size=12), # axis.title.y = element_text(size=12))+ # ggtitle("Profitable mimics") XLAB = "Proportion of moderately unprofitable prey" # XLAB = "Proportion of profitable prey" # HERE # attack rate with confidence interval # # attPlot1 <- ggplot(tableAttackRateMean)+ # geom_rect(mapping=aes(xmin=q-PlotFigure1.stepQHigh2/2, xmax=q+PlotFigure1.stepQHigh2/2, ymin=Ntot-PlotFigure1.stepNtot2/2, ymax=Ntot+PlotFigure1.stepNtot2/2, # fill=mean, color=mean)) + # scale_fill_gradient(name="Per capita\npredation rate", low=PlotFigure1.colorLow, high=PlotFigure1.colorHigh, limits=c(0, PlotFigure1.topLimitAttackRate))+ # scale_color_gradient(name="Per capita\npredation rate", low=PlotFigure1.colorLow, high=PlotFigure1.colorHigh, limits=c(0, PlotFigure1.topLimitAttackRate))+ # xlab("Proportion of mimics")+ # ylab("Number of prey (N)")+ # theme(panel.background = element_blank(), # title = element_text(size=12), # axis.title.x = element_text(size=12), # axis.title.y = element_text(size=12), # legend.position="bottom")+ # ggtitle("All prey") # # attPlot2 <- ggplot(tableAttackRateMean)+ # geom_rect(mapping=aes(xmin=q-PlotFigure1.stepQHigh2/2, xmax=q+PlotFigure1.stepQHigh2/2, ymin=Ntot-PlotFigure1.stepNtot2/2, ymax=Ntot+PlotFigure1.stepNtot2/2, # fill=mean3, color=mean3)) + # scale_fill_gradient(name="Per capita\npredation rate\n(class 3)", low=PlotFigure1.colorLow, high=PlotFigure1.colorHigh, limits=c(0, PlotFigure1.topLimitAttackRate))+ # scale_color_gradient(name="Per capita\npredation rate\n(class 3)", low=PlotFigure1.colorLow, high=PlotFigure1.colorHigh, limits=c(0, PlotFigure1.topLimitAttackRate))+ # xlab("Proportion of mimics")+ # ylab("Number of prey (N)")+ # theme(panel.background = element_blank(), # title = element_text(size=12), # axis.title.x = element_text(size=12), # axis.title.y = element_text(size=12))+ # ggtitle("Highly unprofitable models") # # attPlot3 <- ggplot(tableAttackRateMean)+ # geom_rect(mapping=aes(xmin=q-PlotFigure1.stepQHigh2/2, xmax=q+PlotFigure1.stepQHigh2/2, ymin=Ntot-PlotFigure1.stepNtot2/2, ymax=Ntot+PlotFigure1.stepNtot2/2, # fill=mean2, color=mean2)) + # scale_fill_gradient(name="Per capita\npredation rate\n(class 2)", low=PlotFigure1.colorLow, high=PlotFigure1.colorHigh, limits=c(0, PlotFigure1.topLimitAttackRate))+ # scale_color_gradient(name="Per capita\npredation rate\n(class 2)", low=PlotFigure1.colorLow, high=PlotFigure1.colorHigh, limits=c(0, PlotFigure1.topLimitAttackRate))+ # xlab("Proportion of mimics")+ # ylab("Number of prey (N)")+ # theme(panel.background = element_blank(), # title = element_text(size=12), # axis.title.x = element_text(size=12), # axis.title.y = element_text(size=12))+ # ggtitle("Moderately unprofitable mimics") # # attPlot4 <- ggplot(tableAttackRateMean)+ # geom_rect(mapping=aes(xmin=q-PlotFigure1.stepQHigh2/2, xmax=q+PlotFigure1.stepQHigh2/2, ymin=Ntot-PlotFigure1.stepNtot2/2, ymax=Ntot+PlotFigure1.stepNtot2/2, # fill=mean1, color=mean1)) + # scale_fill_gradient(name="Per capita\npredation rate\n(class 2)", low=PlotFigure1.colorLow, high=PlotFigure1.colorHigh, limits=c(0, PlotFigure1.topLimitAttackRate))+ # scale_color_gradient(name="Per capita\npredation rate\n(class 2)", low=PlotFigure1.colorLow, high=PlotFigure1.colorHigh, limits=c(0, PlotFigure1.topLimitAttackRate))+ # xlab("Proportion of mimics")+ # ylab("Number of prey (N)")+ # theme(panel.background = element_blank(), # title = element_text(size=12), # axis.title.x = element_text(size=12), # axis.title.y = element_text(size=12))+ # ggtitle("Profitable mimics") nbsampledPlot1 <- ggplot(tableAttackRateMean)+ geom_point(aes(q,meanNb,color=as.factor(Ntot),shape=as.factor(Ntot)),size=3)+ geom_line(aes(q,meanNb,color=as.factor(Ntot)),size=0.5)+ scale_colour_manual(name="Number of prey (N):", breaks=PlotFigure1.Ntot, values=c("grey","#737373","black"))+ scale_shape_manual(name="Number of prey (N):", breaks=PlotFigure1.Ntot, values=c(15,16,17))+ xlab(XLAB)+ ylab("Number of prey sampled")+ ylim(0,PlotFigure1.topLimit)+ theme(panel.background = element_rect(fill='white', colour='black'), axis.title.x = element_text(size=20), axis.title.y = element_text(size=20), legend.text=element_text(size=15), axis.text = element_text(size=13), panel.grid.major = element_line(colour = "#DFDFDF",size=0.003), panel.grid.minor = element_blank(),legend.key = element_blank(), legend.title=element_text(size=15), legend.text.align=0, legend.key.width=unit(2,"line"), legend.position="bottom") nbsampledPlot2 <- ggplot(tableAttackRateMean)+ geom_point(aes(q,meanNb3,color=as.factor(Ntot),shape=as.factor(Ntot)),size=3)+ geom_line(aes(q,meanNb3,color=as.factor(Ntot)),size=0.5)+ scale_colour_manual(name="Number of prey (N) :", breaks=PlotFigure1.Ntot, values=c("grey","#737373","black"))+ scale_shape_manual(name="Number of prey :", breaks=PlotFigure1.Ntot, values=c(15,16,17))+ xlab(XLAB)+ ylab("Number of prey sampled")+ ylim(0,PlotFigure1.topLimit)+ theme(panel.background = element_rect(fill='white', colour='black'), axis.title.x = element_text(size=20), axis.title.y = element_text(size=20), legend.text=element_text(size=15), axis.text = element_text(size=13), panel.grid.major = element_line(colour = "#DFDFDF",size=0.003), panel.grid.minor = element_blank(),legend.key = element_blank(), legend.title=element_text(size=15), legend.text.align=0, legend.key.width=unit(2,"line"), legend.position="bottom") nbsampledPlot3 <- ggplot(tableAttackRateMean)+ geom_point(aes(q,meanNb2,color=as.factor(Ntot),shape=as.factor(Ntot)),size=3)+ geom_line(aes(q,meanNb2,color=as.factor(Ntot)),size=0.5)+ scale_colour_manual(name="Number of prey :", breaks=PlotFigure1.Ntot, values=c("grey","#737373","black"))+ scale_shape_manual(name="Number of prey :", breaks=PlotFigure1.Ntot, values=c(15,16,17))+ xlab(XLAB)+ ylab("Number of prey sampled")+ ylim(0,PlotFigure1.topLimit)+ theme(panel.background = element_rect(fill='white', colour='black'), axis.title.x = element_text(size=20), axis.title.y = element_text(size=20), legend.text=element_text(size=15), axis.text = element_text(size=13), panel.grid.major = element_line(colour = "#DFDFDF",size=0.003), panel.grid.minor = element_blank(),legend.key = element_blank(), legend.title=element_text(size=15), legend.text.align=0, legend.key.width=unit(2,"line"), legend.position="bottom") nbsampledPlot4 <- ggplot(tableAttackRateMean)+ geom_point(aes(q,meanNb1,color=as.factor(Ntot),shape=as.factor(Ntot)),size=3)+ geom_line(aes(q,meanNb1,color=as.factor(Ntot)),size=0.5)+ scale_colour_manual(name="Number of prey :", breaks=PlotFigure1.Ntot, values=c("grey","#737373","black"))+ scale_shape_manual(name="Number of prey :", breaks=PlotFigure1.Ntot, values=c(15,16,17))+ xlab(XLAB)+ ylab("Number of prey sampled")+ ylim(0,PlotFigure1.topLimit)+ theme(panel.background = element_rect(fill='white', colour='black'), axis.title.x = element_text(size=20), axis.title.y = element_text(size=20), legend.text=element_text(size=15), axis.text = element_text(size=13), panel.grid.major = element_line(colour = "#DFDFDF",size=0.003), panel.grid.minor = element_blank(),legend.key = element_blank(), legend.title=element_text(size=15), legend.text.align=0, legend.key.width=unit(2,"line"), legend.position="bottom") # grid.arrange(attPlot1,attPlot2,attPlot3,ncol=3) # grid.arrange(nbsampledPlot1,nbsampledPlot2,nbsampledPlot3,ncol=3) attPlot1 <- ggplot(tableAttackRateMean)+ geom_point(aes(q,mean,color=as.factor(Ntot),shape=as.factor(Ntot)),size=3)+ geom_line(aes(q,mean,color=as.factor(Ntot)),size=0.5)+ scale_colour_manual(name="Number of prey :", breaks=PlotFigure1.Ntot, values=c("grey","#737373","black"))+ scale_shape_manual(name="Number of prey :", breaks=PlotFigure1.Ntot, values=c(15,16,17))+ xlab(XLAB)+ ylab("Per capita predation rate")+ ylim(0,1)+ theme(panel.background = element_rect(fill='white', colour='black'), axis.title.x = element_text(size=20), axis.title.y = element_text(size=20), legend.text=element_text(size=15), axis.text = element_text(size=13), panel.grid.major = element_line(colour = "#DFDFDF",size=0.003), panel.grid.minor = element_blank(),legend.key = element_blank(), legend.title=element_text(size=15), legend.text.align=0, legend.key.width=unit(2,"line"), legend.position="bottom") attPlot4 <- ggplot(tableAttackRateMean)+ geom_point(aes(q,mean1,color=as.factor(Ntot),shape=as.factor(Ntot)),size=3)+ geom_line(aes(q,mean1,color=as.factor(Ntot)),size=0.5)+ scale_colour_manual(name="Number of prey :", breaks=PlotFigure1.Ntot, values=c("grey","#737373","black"))+ scale_shape_manual(name="Number of prey :", breaks=PlotFigure1.Ntot, values=c(15,16,17))+ xlab(XLAB)+ ylab("Per capita predation rate")+ ylim(0,1)+ theme(panel.background = element_rect(fill='white', colour='black'), axis.title.x = element_text(size=20), axis.title.y = element_text(size=20), legend.text=element_text(size=15), axis.text = element_text(size=13), panel.grid.major = element_line(colour = "#DFDFDF",size=0.003), panel.grid.minor = element_blank(),legend.key = element_blank(), legend.title=element_text(size=15), legend.text.align=0, legend.key.width=unit(2,"line"), legend.position="bottom") mylegend <- g_legend(nbsampledPlot1) if(PlotFigure1.mutantTypeAdded==1){ grid.arrange(arrangeGrob(nbsampledPlot1 + theme(legend.position="none"), nbsampledPlot2 + theme(legend.position="none"), nbsampledPlot4 + theme(legend.position="none"), nrow=1), mylegend, nrow=2,heights=c(6, 1)) }else if(PlotFigure1.mutantTypeAdded==2){ grid.arrange(arrangeGrob(nbsampledPlot1 + theme(legend.position="none"), nbsampledPlot2 + theme(legend.position="none"), nbsampledPlot3 + theme(legend.position="none"), nrow=1), mylegend, nrow=2,heights=c(6, 1)) } mylegend <- g_legend(nbsampledPlot1) if(PlotFigure1.mutantTypeAdded==1){ grid.arrange(arrangeGrob(attPlot1 + theme(legend.position="none"), attPlot2 + theme(legend.position="none"), attPlot4 + theme(legend.position="none"), nrow=1), nrow=2,heights=c(6, 1)) }else if(PlotFigure1.mutantTypeAdded==2){ grid.arrange(arrangeGrob(attPlot1 + theme(legend.position="none"), attPlot2 + theme(legend.position="none"), attPlot3 + theme(legend.position="none"), nrow=1), mylegend, nrow=2,heights=c(6, 1)) } } if(PlotFigure2==TRUE){ if(PlotFigure2.stepNbMutantNb>PlotFigure2.n1){ PlotFigure2.stepNbMutantNb = PlotFigure2.n1 } PlotFigure2.stepNbMutant = floor(PlotFigure2.n1/PlotFigure2.stepNbMutantNb) n1Vec <- c() attackRateVec <- c() attackRateVec2 <- c() attackRateVec3 <- c() NbPreyEatenVec <- c() NbPreyEatenVec2 <- c() NbPreyEatenVec3 <- c() popIndexVec <- c() # if 1 : bad prey, 2: very bad prey pop PlotFigure2.nMutant = seq(0,PlotFigure2.n1-1,by=PlotFigure2.stepNbMutant) # number mutants bad prey's phenotype to very bad prey's phenotype for(n1 in PlotFigure2.nMutant){ for(pop in 1:3){ if(pop==1){ # community of prey nbBadPrey = n1 preyCommunity <- c(rep(2,PlotFigure2.n2),rep(PlotFigure2.mutantTypeAdded-1,n1)) # 0 : profitable prey, nbGoodPrey = length(which(preyCommunity==0)) nbBadPrey = length(which(preyCommunity==1)) nbVeryBadPrey = length(which(preyCommunity==2)) # preyCommunity <- c(rep(1,PlotFigure2.n1-n1)) # 0 : profitable prey, # nbVeryBadPrey = 0 }else if(pop==2){ nbBadPrey = n1 # preyCommunity <- c(rep(2,PlotFigure2.n2bis),rep(PlotFigure2.mutantTypeAdded-1,n1)) # 0 : profitable prey, preyCommunity <- c(1) # No simulation nbGoodPrey = length(which(preyCommunity==0)) nbBadPrey = length(which(preyCommunity==1)) nbVeryBadPrey = length(which(preyCommunity==2)) }else if(pop==3){ # community of prey ### Rowland 2010: preyCommunity <- c(rep(PlotFigure2.mutantTypeAdded-1,PlotFigure2.n1-n1)) # 0 : profitable prey, ### HERE Rowland 2007: preyCommunity <- c(rep(PlotFigure2.mutantTypeAdded-1,PlotFigure2.n1)) # 0 : profitable prey, nbGoodPrey = length(which(preyCommunity==0)) nbBadPrey = length(which(preyCommunity==1)) nbVeryBadPrey = length(which(preyCommunity==2)) } # 1 : toxic prey (c1), 2 : toxic prey (c2) # forward iteration print(floor(PlotFigure2.propCommunityExperienced*length(preyCommunity))) PlotFigure2.bigN = length(preyCommunity) A = SamplingStrategyMatrix(PlotFigure2.alpha, floor(PlotFigure2.propCommunityExperienced*length(preyCommunity)), PlotFigure2.f1, -PlotFigure2.f2, -PlotFigure2.f3) for(rep in 1:PlotFigure2.nbRepForward){ preyCommunityOnePred <- sample(preyCommunity,size=floor(PlotFigure2.propCommunityExperienced*length(preyCommunity))) # encounter in random order predKnowledge.ri1 = 0 # profitable prey predKnowledge.ri2 = 0 # toxic prey (c1) predKnowledge.ri3 = 0 # toxic prey (c2) nbPreyEaten = 0 nbPreyEaten1 = 0 nbPreyEaten2 = 0 nbPreyEaten3 = 0 for(prey in preyCommunityOnePred){ optBehavior = A[predKnowledge.ri1+1, predKnowledge.ri2+1, predKnowledge.ri3+1] if(optBehavior==1){ # attack nbPreyEaten = nbPreyEaten + 1 if(prey==0){ predKnowledge.ri1 = predKnowledge.ri1 + 1 nbPreyEaten1 = nbPreyEaten1 + 1 }else if(prey==1){ predKnowledge.ri2 = predKnowledge.ri2 + 1 nbPreyEaten2 = nbPreyEaten2 + 1 }else if(prey==2){ predKnowledge.ri3 = predKnowledge.ri3 + 1 nbPreyEaten3 = nbPreyEaten3 + 1 } } } popIndexVec <- c(popIndexVec, pop) n1Vec <- c(n1Vec, n1) attackRateVec <- c(attackRateVec, nbPreyEaten/PlotFigure2.bigN) if(PlotFigure2.mutantTypeAdded==1){ attackRateVec2 <- c(attackRateVec2, nbPreyEaten1/nbGoodPrey) }else if(PlotFigure2.mutantTypeAdded==2){ attackRateVec2 <- c(attackRateVec2, nbPreyEaten2/nbBadPrey) } attackRateVec3 <- c(attackRateVec3, nbPreyEaten3/nbVeryBadPrey) NbPreyEatenVec <- c(NbPreyEatenVec,nbPreyEaten) NbPreyEatenVec2 <- c(NbPreyEatenVec2,nbPreyEaten2) NbPreyEatenVec3 <- c(NbPreyEatenVec3,nbPreyEaten3) } } } tableAttackRate = data.frame(n1=n1Vec,attackRate=attackRateVec,attackRate2=attackRateVec2,attackRate3=attackRateVec3,NbPreyEaten=NbPreyEatenVec,NbPreyEaten2=NbPreyEatenVec2,NbPreyEaten3=NbPreyEatenVec3,popIndex=popIndexVec) # tableAttackRate <- tableAttackRate[tableAttackRate$popIndex==1,] tableAttackRateMean = ddply(tableAttackRate, .(n1,popIndex), summarize, mean = mean(attackRate,na.rm=T), sd = sd(attackRate,na.rm=T), mean2 = mean(attackRate2,na.rm=T), sd2 = sd(attackRate2,na.rm=T), mean3 = mean(attackRate3,na.rm=T), sd3 = sd(attackRate3,na.rm=T), meanNb = mean(NbPreyEaten,na.rm=T), sdNb = sd(NbPreyEaten,na.rm=T), meanNb2 = mean(NbPreyEaten2,na.rm=T), sdNb2 = sd(NbPreyEaten2,na.rm=T), meanNb3 = mean(NbPreyEaten3,na.rm=T), sdNb3 = sd(NbPreyEaten3,na.rm=T)) tableAttackRateMean$sampleNb = NA for(i in 1:dim(tableAttackRateMean)[1]){ tableAttackRateMean$sampleNb[i] = length(which(tableAttackRate$n1==tableAttackRateMean$n1[i])) } tableAttackRateMean$confidenceInterval <- qt(0.975,df=tableAttackRateMean$sampleNb-1)*tableAttackRateMean$sd/sqrt(tableAttackRateMean$sampleNb) tableAttackRateMean$confidenceInterval2 <- qt(0.975,df=tableAttackRateMean$sampleNb-1)*tableAttackRateMean$sd2/sqrt(tableAttackRateMean$sampleNb) tableAttackRateMean$confidenceInterval3 <- qt(0.975,df=tableAttackRateMean$sampleNb-1)*tableAttackRateMean$sd3/sqrt(tableAttackRateMean$sampleNb) tableAttackRateMean$confidenceIntervalNb <- qt(0.975,df=tableAttackRateMean$sampleNb-1)*tableAttackRateMean$sdNb/sqrt(tableAttackRateMean$sampleNb) tableAttackRateMean$confidenceIntervalNb2 <- qt(0.975,df=tableAttackRateMean$sampleNb-1)*tableAttackRateMean$sdNb2/sqrt(tableAttackRateMean$sampleNb) tableAttackRateMean$confidenceIntervalNb3 <- qt(0.975,df=tableAttackRateMean$sampleNb-1)*tableAttackRateMean$sdNb3/sqrt(tableAttackRateMean$sampleNb) # tableAttackRateMean$n1 = tableAttackRateMean$n1/PlotFigure2.n1 # # xFigBreak = PlotFigure2.nMutant[seq(1, length(PlotFigure2.nMutant), 4)] # # attack rate with confidence interval # labelPhenotype = c(expression(paste(C[1]," (",N[1],"=50",")",sep="")), # expression(paste(C[1]," (",N[1],"=10",")",sep="")), # expression(paste(C[2]," (",N[d],"=50",")",sep=""))) # attPlot1 <- # ggplot(tableAttackRateMean)+ # geom_point(aes(n1,mean,color=as.factor(popIndex),shape=as.factor(popIndex)),size=3)+ # geom_line(aes(n1,mean,color=as.factor(popIndex),lty=as.factor(popIndex)),size=0.5)+ # xlab(expression(N[mimic]))+ # ylab("Per capita predation rate")+ # ylim(PlotFigure2.yLim)+ # scale_colour_manual(name="Phenotype :", # breaks=c(1, 2, 3), # labels=labelPhenotype, # values=c("#A1A1A1","#A1A1A1","black"))+ # scale_shape_manual(name="Phenotype :", # breaks=c(1, 2, 3), # labels=labelPhenotype, # values=c(16,17,16))+ # scale_linetype_manual(name="Phenotype :", # breaks=c(1, 2, 3), # labels=labelPhenotype, # values=c(1,2,1))+ # theme(panel.background = element_rect(fill='white', colour='black'), # axis.title.x = element_text(size=20), # axis.title.y = element_text(size=20), # legend.text=element_text(size=15), # panel.grid.major = element_line(colour = "#DFDFDF",size=0.003), # panel.grid.minor = element_blank(),legend.key = element_blank(), # legend.title=element_text(size=15), # legend.text.align=0, # legend.key.width=unit(2,"line"))+ # scale_x_continuous(breaks=xFigBreak) # # nbSampledPlot1 <- # ggplot(tableAttackRateMean)+ # geom_point(aes(n1,meanNb,color=as.factor(popIndex),shape=as.factor(popIndex)),size=3)+ # geom_line(aes(n1,meanNb,color=as.factor(popIndex),lty=as.factor(popIndex)),size=0.5)+ # xlab(expression(N[2]))+ # ylab("Number of prey sampled")+ # ylim(PlotFigure2.yLimNbSampled)+ # scale_colour_manual(name="Phenotype :", # breaks=c(1, 2, 3), # labels=labelPhenotype, # values=c("#A1A1A1","#A1A1A1","black"))+ # scale_shape_manual(name="Phenotype :", # breaks=c(1, 2, 3), # labels=labelPhenotype, # values=c(16,17,16))+ # scale_linetype_manual(name="Phenotype :", # breaks=c(1, 2, 3), # labels=labelPhenotype, # values=c(1,2,1))+ # theme(panel.background = element_rect(fill='white', colour='black'), # axis.title.x = element_text(size=20), # axis.title.y = element_text(size=20), # legend.text=element_text(size=15), # panel.grid.major = element_line(colour = "#DFDFDF",size=0.003), # panel.grid.minor = element_blank(),legend.key = element_blank(), # legend.title=element_text(size=15), # legend.text.align=0, # legend.key.width=unit(2,"line"))+ # scale_x_continuous(breaks=xFigBreak) # # # grid.arrange(attPlot1 + theme(legend.position="none")) # # grid.arrange(nbSampledPlot1 + theme(legend.position="none")) # # grid.arrange(attPlot1) # # grid.arrange(nbSampledPlot1) # # # mylegend <- g_legend(attPlot1) # # grid.arrange(arrangeGrob(attPlot1 + theme(legend.position="none"), # # nrow=1), # # mylegend, nrow=2,heights=c(3, 1)) xFigBreak = PlotFigure2.nMutant[seq(1, length(PlotFigure2.nMutant), 4)] # attack rate with confidence interval labelPhenotype = c("Non mimics", "Models and\nmimics") tableAttackRateMeanPlot1 <- tableAttackRateMean tableAttackRateMeanPlot1 <- tableAttackRateMeanPlot1[tableAttackRateMeanPlot1$popIndex%in%c(3,1),] popNb=1 # attPlot1 <- # ggplot(tableAttackRateMeanPlot1)+ # geom_point(aes(n1,mean,color=as.factor(popIndex),shape=as.factor(popIndex)),size=3)+ # geom_line(aes(n1,mean,color=as.factor(popIndex),lty=as.factor(popIndex)),size=0.5)+ # xlab(expression(N[mimic]))+ # ylab("Per capita predation rate")+ # ylim(PlotFigure2.yLim)+ # scale_colour_manual(name="Appearance :", # breaks=c(3, popNb), # labels=labelPhenotype, # values=c("black","#A1A1A1"))+ # scale_shape_manual(name="Appearance :", # breaks=c(3, popNb), # labels=labelPhenotype, # values=c(16,16))+ # scale_linetype_manual(name="Appearance :", # breaks=c(3, popNb), # labels=labelPhenotype, # values=c(1,1))+ # theme(panel.background = element_rect(fill='white', colour='black'), # axis.title.x = element_text(size=20), # axis.title.y = element_text(size=20), # legend.text=element_text(size=15), # axis.text = element_text(size=13), # panel.grid.major = element_line(colour = "#DFDFDF",size=0.003), # panel.grid.minor = element_blank(),legend.key = element_blank(), # legend.title=element_text(size=15), # legend.text.align=0, # legend.key.width=unit(2,"line"))+ # scale_x_continuous(breaks=xFigBreak) # # nbSampledPlot1 <- # ggplot(tableAttackRateMeanPlot1)+ # geom_point(aes(n1,meanNb,color=as.factor(popIndex),shape=as.factor(popIndex)),size=3)+ # geom_line(aes(n1,meanNb,color=as.factor(popIndex),lty=as.factor(popIndex)),size=0.5)+ # xlab(expression(N[mimic]))+ # ylab("Number of prey sampled")+ # ylim(PlotFigure2.yLimNbSampled)+ # scale_colour_manual(name="Appearance :", # breaks=c(3, popNb), # labels=labelPhenotype, # values=c("black","#A1A1A1"))+ # scale_shape_manual(name="Appearance :", # breaks=c(3, popNb), # labels=labelPhenotype, # values=c(16,16))+ # scale_linetype_manual(name="Appearance :", # breaks=c(3, popNb), # labels=labelPhenotype, # values=c(1,1))+ # theme(panel.background = element_rect(fill='white', colour='black'), # axis.title.x = element_text(size=20), # axis.title.y = element_text(size=20), # legend.text=element_text(size=15), # axis.text = element_text(size=13), # panel.grid.major = element_line(colour = "#DFDFDF",size=0.003), # panel.grid.minor = element_blank(),legend.key = element_blank(), # legend.title=element_text(size=15), # legend.text.align=0, # legend.key.width=unit(2,"line"))+ # scale_x_continuous(breaks=xFigBreak) # # nbSampledModelPlot1 <- # ggplot(tableAttackRateMeanPlot1[tableAttackRateMeanPlot1$popIndex==1,])+ # geom_point(aes(n1,meanNb3,color=as.factor(popIndex),shape=as.factor(popIndex)),size=3)+ # geom_line(aes(n1,meanNb3,color=as.factor(popIndex),lty=as.factor(popIndex)),size=0.5)+ # xlab(expression(N[mimic]))+ # ylab("Number of models sampled")+ # ylim(PlotFigure2.yLimNbSampledModel)+ # scale_colour_manual(name="Appearance :", # breaks=c(3, popNb), # labels=labelPhenotype, # values=c("black","#A1A1A1"))+ # scale_shape_manual(name="Appearance :", # breaks=c(3, popNb), # labels=labelPhenotype, # values=c(16,16))+ # scale_linetype_manual(name="Appearance :", # breaks=c(3, popNb), # labels=labelPhenotype, # values=c(1,1))+ # theme(legend.position="none", # panel.background = element_rect(fill='white', colour='black'), # axis.title.x = element_text(size=20), # axis.title.y = element_text(size=20), # legend.text=element_text(size=15), # axis.text = element_text(size=13), # panel.grid.major = element_line(colour = "#DFDFDF",size=0.003), # panel.grid.minor = element_blank(),legend.key = element_blank(), # legend.title=element_text(size=15), # legend.text.align=0, # legend.key.width=unit(2,"line"))+ # scale_x_continuous(breaks=xFigBreak) # # # tableAttackRateMeanPlot2 <- tableAttackRateMean # tableAttackRateMeanPlot2 <- tableAttackRateMeanPlot2[tableAttackRateMeanPlot2$popIndex%in%c(3,2),] # popNb = 2 # attPlot2 <- # ggplot(tableAttackRateMeanPlot2)+ # geom_point(aes(n1,mean,color=as.factor(popIndex),shape=as.factor(popIndex)),size=3)+ # geom_line(aes(n1,mean,color=as.factor(popIndex),lty=as.factor(popIndex)),size=0.5)+ # xlab(expression(N[mimic]))+ # ylab("Per capita predation rate")+ # ylim(PlotFigure2.yLim)+ # scale_colour_manual(name="", # breaks=c(3, popNb), # labels=labelPhenotype, # values=c("black","#A1A1A1"))+ # scale_shape_manual(name="", # breaks=c(3, popNb), # labels=labelPhenotype, # values=c(16,16))+ # scale_linetype_manual(name="", # breaks=c(3, popNb), # labels=labelPhenotype, # values=c(1,1))+ # theme(panel.background = element_rect(fill='white', colour='black'), # axis.title.x = element_text(size=20), # axis.text = element_text(size=13), # axis.title.y = element_text(size=20), # legend.text=element_text(size=15), # panel.grid.major = element_line(colour = "#DFDFDF",size=0.003), # panel.grid.minor = element_blank(),legend.key = element_blank(), # legend.title=element_text(size=15), # legend.text.align=0, # legend.key.width=unit(2,"line"))+ # scale_x_continuous(breaks=xFigBreak) # # nbSampledPlot2 <- # ggplot(tableAttackRateMeanPlot2)+ # geom_point(aes(n1,meanNb,color=as.factor(popIndex),shape=as.factor(popIndex)),size=3)+ # geom_line(aes(n1,meanNb,color=as.factor(popIndex),lty=as.factor(popIndex)),size=0.5)+ # xlab(expression(N[mimic]))+ # ylab("Number of prey sampled")+ # ylim(PlotFigure2.yLimNbSampled)+ # scale_colour_manual(name="", # breaks=c(3, popNb), # labels=labelPhenotype, # values=c("black","#A1A1A1"))+ # scale_shape_manual(name="", # breaks=c(3, popNb), # labels=labelPhenotype, # values=c(16,16))+ # scale_linetype_manual(name="", # breaks=c(3, popNb), # labels=labelPhenotype, # values=c(1,1))+ # theme(panel.background = element_rect(fill='white', colour='black'), # axis.title.x = element_text(size=20), # axis.title.y = element_text(size=20), # axis.text = element_text(size=13), # legend.text=element_text(size=15), # panel.grid.major = element_line(colour = "#DFDFDF",size=0.003), # panel.grid.minor = element_blank(),legend.key = element_blank(), # legend.title=element_text(size=15), # legend.text.align=0, # legend.key.width=unit(2,"line"))+ # scale_x_continuous(breaks=xFigBreak) # # if(Arrow1==TRUE){ # dataArrow <- data.frame(xArrow1,yArrow1I,yArrow1E) # attPlot1 <- attPlot1 + # geom_segment(data=dataArrow, aes(x = xArrow1, y = yArrow1I, xend = xArrow1, yend = yArrow1E), color=colorArrow1, arrow = arrow(length = unit(0.4, "cm"))) # } # if(Arrow2==TRUE){ # dataArrow <- data.frame(xArrow2,yArrow2I,yArrow2E) # attPlot2 <- attPlot2 + # geom_segment(data=dataArrow,aes(x = xArrow2, y = yArrow2I, xend = xArrow2, yend = yArrow2E), color=colorArrow2, arrow = arrow(length = unit(0.4, "cm"))) # } # # # grid.arrange(nbSampledPlot1 + theme(legend.position="none")) # # # grid.arrange(attPlot2) # # grid.arrange(attPlot1) # grid.arrange(attPlot1 + theme(legend.position="none")) # # grid.arrange(nbSampledPlot1 + theme(legend.position="none")) # grid.arrange(nbSampledModelPlot1) # # grid.arrange(nbSampledPlot1) # # # mylegend <- g_legend(attPlot1) # # grid.arrange(arrangeGrob(attPlot1 + theme(legend.position="none"), # # nrow=1), # # mylegend, nrow=2,heights=c(3, 1)) XLAB = expression(paste("Proportion of mimics (",P[mimic],")",sep="")) xFigBreak = c(0,PlotFigure2.n1/2,PlotFigure2.n1) xFigLabel = c(0,0.5,1) attPlot1 <- ggplot(tableAttackRateMeanPlot1)+ geom_point(aes(n1,mean,color=as.factor(popIndex),shape=as.factor(popIndex)),size=3)+ geom_line(aes(n1,mean,color=as.factor(popIndex),lty=as.factor(popIndex)),size=0.5)+ xlab(XLAB)+ ylab("Per capita predation rate")+ ylim(PlotFigure2.yLim)+ scale_colour_manual(name="Appearance :", breaks=c(3, popNb), labels=labelPhenotype, values=c("black","#A1A1A1"))+ scale_shape_manual(name="Appearance :", breaks=c(3, popNb), labels=labelPhenotype, values=c(16,16))+ scale_linetype_manual(name="Appearance :", breaks=c(3, popNb), labels=labelPhenotype, values=c(1,1))+ theme(panel.background = element_rect(fill='white', colour='black'), axis.title.x = element_text(size=20), axis.title.y = element_text(size=20), legend.text=element_text(size=15), axis.text = element_text(size=13), panel.grid.major = element_line(colour = "#DFDFDF",size=0.003), panel.grid.minor = element_blank(),legend.key = element_blank(), legend.title=element_text(size=15), legend.text.align=0, legend.key.width=unit(2,"line"), legend.key.height=unit(3,"line"))+ scale_x_continuous(breaks=xFigBreak,labels = xFigLabel,limits=c(0,PlotFigure2.n1)) nbSampledPlot1 <- ggplot(tableAttackRateMeanPlot1)+ geom_point(aes(n1,meanNb,color=as.factor(popIndex),shape=as.factor(popIndex)),size=3)+ geom_line(aes(n1,meanNb,color=as.factor(popIndex),lty=as.factor(popIndex)),size=0.5)+ xlab(XLAB)+ ylab("Number of prey sampled")+ ylim(PlotFigure2.yLimNbSampled)+ scale_colour_manual(name="Appearance :", breaks=c(3, popNb), labels=labelPhenotype, values=c("black","#A1A1A1"))+ scale_shape_manual(name="Appearance :", breaks=c(3, popNb), labels=labelPhenotype, values=c(16,16))+ scale_linetype_manual(name="Appearance :", breaks=c(3, popNb), labels=labelPhenotype, values=c(1,1))+ theme(panel.background = element_rect(fill='white', colour='black'), axis.title.x = element_text(size=20), axis.title.y = element_text(size=20), legend.text=element_text(size=15), axis.text = element_text(size=13), panel.grid.major = element_line(colour = "#DFDFDF",size=0.003), panel.grid.minor = element_blank(),legend.key = element_blank(), legend.title=element_text(size=15), legend.text.align=0, legend.key.width=unit(2,"line"))+ scale_x_continuous(breaks=xFigBreak,labels = xFigLabel,limits=c(0,PlotFigure2.n1)) nbSampledModelPlot1 <- ggplot(tableAttackRateMeanPlot1[tableAttackRateMeanPlot1$popIndex==1,])+ geom_point(aes(n1,meanNb3,color=as.factor(popIndex),shape=as.factor(popIndex)),size=3)+ geom_line(aes(n1,meanNb3,color=as.factor(popIndex),lty=as.factor(popIndex)),size=0.5)+ xlab(XLAB)+ ylab("Number of models sampled")+ ylim(PlotFigure2.yLimNbSampledModel)+ scale_colour_manual(name="Appearance :", breaks=c(3, popNb), labels=labelPhenotype, values=c("black","#A1A1A1"))+ scale_shape_manual(name="Appearance :", breaks=c(3, popNb), labels=labelPhenotype, values=c(16,16))+ scale_linetype_manual(name="Appearance :", breaks=c(3, popNb), labels=labelPhenotype, values=c(1,1))+ theme(legend.position="none", panel.background = element_rect(fill='white', colour='black'), axis.title.x = element_text(size=20), axis.title.y = element_text(size=20), legend.text=element_text(size=15), axis.text = element_text(size=13), panel.grid.major = element_line(colour = "#DFDFDF",size=0.003), panel.grid.minor = element_blank(),legend.key = element_blank(), legend.title=element_text(size=15), legend.text.align=0, legend.key.width=unit(2,"line"))+ scale_x_continuous(breaks=xFigBreak,labels = xFigLabel,limits=c(0,PlotFigure2.n1)) tableAttackRateMeanPlot2 <- tableAttackRateMean tableAttackRateMeanPlot2 <- tableAttackRateMeanPlot2[tableAttackRateMeanPlot2$popIndex%in%c(3,2),] popNb = 2 attPlot2 <- ggplot(tableAttackRateMeanPlot2)+ geom_point(aes(n1,mean,color=as.factor(popIndex),shape=as.factor(popIndex)),size=3)+ geom_line(aes(n1,mean,color=as.factor(popIndex),lty=as.factor(popIndex)),size=0.5)+ xlab(XLAB)+ ylab("Per capita predation rate")+ ylim(PlotFigure2.yLim)+ scale_colour_manual(name="", breaks=c(3, popNb), labels=labelPhenotype, values=c("black","#A1A1A1"))+ scale_shape_manual(name="", breaks=c(3, popNb), labels=labelPhenotype, values=c(16,16))+ scale_linetype_manual(name="", breaks=c(3, popNb), labels=labelPhenotype, values=c(1,1))+ theme(panel.background = element_rect(fill='white', colour='black'), axis.title.x = element_text(size=20), axis.text = element_text(size=13), axis.title.y = element_text(size=20), legend.text=element_text(size=15), panel.grid.major = element_line(colour = "#DFDFDF",size=0.003), panel.grid.minor = element_blank(),legend.key = element_blank(), legend.title=element_text(size=15), legend.text.align=0, legend.key.width=unit(2,"line"))+ scale_x_continuous(breaks=xFigBreak,labels = xFigLabel,limits=c(0,PlotFigure2.n1)) nbSampledPlot2 <- ggplot(tableAttackRateMeanPlot2)+ geom_point(aes(n1,meanNb,color=as.factor(popIndex),shape=as.factor(popIndex)),size=3)+ geom_line(aes(n1,meanNb,color=as.factor(popIndex),lty=as.factor(popIndex)),size=0.5)+ xlab(XLAB)+ ylab("Number of prey sampled")+ ylim(PlotFigure2.yLimNbSampled)+ scale_colour_manual(name="", breaks=c(3, popNb), labels=labelPhenotype, values=c("black","#A1A1A1"))+ scale_shape_manual(name="", breaks=c(3, popNb), labels=labelPhenotype, values=c(16,16))+ scale_linetype_manual(name="", breaks=c(3, popNb), labels=labelPhenotype, values=c(1,1))+ theme(panel.background = element_rect(fill='white', colour='black'), axis.title.x = element_text(size=20), axis.title.y = element_text(size=20), axis.text = element_text(size=13), legend.text=element_text(size=15), panel.grid.major = element_line(colour = "#DFDFDF",size=0.003), panel.grid.minor = element_blank(),legend.key = element_blank(), legend.title=element_text(size=15), legend.text.align=0, legend.key.width=unit(2,"line"))+ scale_x_continuous(breaks=xFigBreak,labels = xFigLabel,limits=c(0, PlotFigure2.n1)) if(Arrow1==TRUE){ dataArrow <- data.frame(xArrow1,yArrow1I,yArrow1E) attPlot1 <- attPlot1 + geom_segment(data=dataArrow, aes(x = xArrow1, y = yArrow1I, xend = xArrow1, yend = yArrow1E), color=colorArrow1, arrow = arrow(length = unit(0.4, "cm"))) } if(Arrow2==TRUE){ dataArrow <- data.frame(xArrow2,yArrow2I,yArrow2E) attPlot2 <- attPlot2 + geom_segment(data=dataArrow,aes(x = xArrow2, y = yArrow2I, xend = xArrow2, yend = yArrow2E), color=colorArrow2, arrow = arrow(length = unit(0.4, "cm"))) } grid.arrange(nbSampledPlot1 + theme(legend.position="none")) grid.arrange(nbSampledModelPlot1 + theme(legend.position="none")) # grid.arrange(attPlot2) # grid.arrange(attPlot1) grid.arrange(attPlot1 + theme(legend.position="none")) } if(PlotFigure3==TRUE){ PlotFigure3.stepNbMutant = floor(PlotFigure3.n1/PlotFigure3.stepNbMutantNb) n1Vec <- c() attackRateVec <- c() attackRateVec2 <- c() attackRateVec3 <- c() NbPreyEatenVec <- c() NbPreyEatenVec2 <- c() NbPreyEatenVec3 <- c() popIndexVec <- c() # if 1 : bad prey, 2: very bad prey pop PlotFigure3.nMutant = seq(0,PlotFigure3.n1,by=PlotFigure3.stepNbMutant) # number mutants bad prey's phenotype to very bad prey's phenotype PlotFigure3.qMutant = PlotFigure3.nMutant/PlotFigure3.n1 for(n1 in PlotFigure3.nMutant){ for(pop in 1:2){ if(pop==1){ # community of prey nbBadPrey = PlotFigure3.n1-n1 preyCommunity <- c(rep(2,PlotFigure3.n2),rep(PlotFigure3.mutantTypeAdded-1,PlotFigure3.n1-n1)) # 0 : profitable prey, nbGoodPrey = length(which(preyCommunity==0)) nbBadPrey = length(which(preyCommunity==1)) nbVeryBadPrey = length(which(preyCommunity==2)) # preyCommunity <- c(rep(1,PlotFigure3.n1-n1)) # 0 : profitable prey, # nbVeryBadPrey = 0 }else if(pop==2){ # community of prey preyCommunity <- c(rep(2,PlotFigure3.n2),rep(PlotFigure3.mutantTypeAdded-1,n1)) # 0 : profitable prey, nbGoodPrey = length(which(preyCommunity==0)) nbBadPrey = length(which(preyCommunity==1)) nbVeryBadPrey = length(which(preyCommunity==2)) } # 1 : toxic prey (c1), 2 : toxic prey (c2) # forward iteration print(floor(PlotFigure3.propCommunityExperienced*length(preyCommunity))) PlotFigure3.bigN = length(preyCommunity) A = SamplingStrategyMatrix(PlotFigure3.alpha, floor(PlotFigure3.propCommunityExperienced*length(preyCommunity)), PlotFigure3.f1, -PlotFigure3.f2, -PlotFigure3.f3) for(rep in 1:PlotFigure3.nbRepForward){ preyCommunityOnePred <- sample(preyCommunity,size=floor(PlotFigure3.propCommunityExperienced*length(preyCommunity))) # encounter in random order predKnowledge.ri1 = 0 # profitable prey predKnowledge.ri2 = 0 # toxic prey (c1) predKnowledge.ri3 = 0 # toxic prey (c2) nbPreyEaten = 0 nbPreyEaten1 = 0 nbPreyEaten2 = 0 nbPreyEaten3 = 0 for(prey in preyCommunityOnePred){ optBehavior = A[predKnowledge.ri1+1, predKnowledge.ri2+1, predKnowledge.ri3+1] if(optBehavior==1){ # attack nbPreyEaten = nbPreyEaten + 1 if(prey==0){ predKnowledge.ri1 = predKnowledge.ri1 + 1 nbPreyEaten1 = nbPreyEaten1 + 1 }else if(prey==1){ predKnowledge.ri2 = predKnowledge.ri2 + 1 nbPreyEaten2 = nbPreyEaten2 + 1 }else if(prey==2){ predKnowledge.ri3 = predKnowledge.ri3 + 1 nbPreyEaten3 = nbPreyEaten3 + 1 } } } popIndexVec <- c(popIndexVec, pop) n1Vec <- c(n1Vec, (PlotFigure3.n1-n1)/(PlotFigure3.n1)) # n1Vec <- c(n1Vec, PlotFigure3.n1-n1) attackRateVec <- c(attackRateVec, nbPreyEaten/PlotFigure3.bigN) if(PlotFigure3.mutantTypeAdded==1){ attackRateVec2 <- c(attackRateVec2, nbPreyEaten1/nbGoodPrey) }else if(PlotFigure3.mutantTypeAdded==2){ attackRateVec2 <- c(attackRateVec2, nbPreyEaten2/nbBadPrey) } attackRateVec3 <- c(attackRateVec3, nbPreyEaten3/nbVeryBadPrey) NbPreyEatenVec <- c(NbPreyEatenVec,nbPreyEaten) NbPreyEatenVec2 <- c(NbPreyEatenVec2,nbPreyEaten2) NbPreyEatenVec3 <- c(NbPreyEatenVec3,nbPreyEaten3) } } } tableAttackRate = data.frame(n1=n1Vec,attackRate=attackRateVec,attackRate2=attackRateVec2,attackRate3=attackRateVec3,NbPreyEaten=NbPreyEatenVec,NbPreyEaten2=NbPreyEatenVec2,NbPreyEaten3=NbPreyEatenVec3,popIndex=popIndexVec) # tableAttackRate <- tableAttackRate[tableAttackRate$popIndex==1,] tableAttackRateMean = ddply(tableAttackRate, .(n1,popIndex), summarize, mean = mean(attackRate,na.rm=T), sd = sd(attackRate,na.rm=T), mean2 = mean(attackRate2,na.rm=T), sd2 = sd(attackRate2,na.rm=T), mean3 = mean(attackRate3,na.rm=T), sd3 = sd(attackRate3,na.rm=T), meanNb = mean(NbPreyEaten,na.rm=T), sdNb = sd(NbPreyEaten,na.rm=T), meanNb2 = mean(NbPreyEaten2,na.rm=T), sdNb2 = sd(NbPreyEaten2,na.rm=T), meanNb3 = mean(NbPreyEaten3,na.rm=T), sdNb3 = sd(NbPreyEaten3,na.rm=T)) tableAttackRateMean$sampleNb = NA for(i in 1:dim(tableAttackRateMean)[1]){ tableAttackRateMean$sampleNb[i] = length(which(tableAttackRate$n1==tableAttackRateMean$n1[i])) } tableAttackRateMean$confidenceInterval <- qt(0.975,df=tableAttackRateMean$sampleNb-1)*tableAttackRateMean$sd/sqrt(tableAttackRateMean$sampleNb) tableAttackRateMean$confidenceInterval2 <- qt(0.975,df=tableAttackRateMean$sampleNb-1)*tableAttackRateMean$sd2/sqrt(tableAttackRateMean$sampleNb) tableAttackRateMean$confidenceInterval3 <- qt(0.975,df=tableAttackRateMean$sampleNb-1)*tableAttackRateMean$sd3/sqrt(tableAttackRateMean$sampleNb) tableAttackRateMean$confidenceIntervalNb <- qt(0.975,df=tableAttackRateMean$sampleNb-1)*tableAttackRateMean$sdNb/sqrt(tableAttackRateMean$sampleNb) tableAttackRateMean$confidenceIntervalNb2 <- qt(0.975,df=tableAttackRateMean$sampleNb-1)*tableAttackRateMean$sdNb2/sqrt(tableAttackRateMean$sampleNb) tableAttackRateMean$confidenceIntervalNb3 <- qt(0.975,df=tableAttackRateMean$sampleNb-1)*tableAttackRateMean$sdNb3/sqrt(tableAttackRateMean$sampleNb) # tableAttackRateMean$n1 = tableAttackRateMean$n1/PlotFigure3.n1 # attack rate with confidence interval labelAppearance = c("model A and\nmimic A", "model B and\nmimic B") xLabel <- as.numeric(c(PlotFigure3.qMutant[1],PlotFigure3.qMutant[length(PlotFigure3.qMutant)],0.5)) # xLabel = c(0,0.5,1) # HERE2 attPlot1 <- ggplot(tableAttackRateMean)+ geom_point(aes(n1,mean,color=as.factor(popIndex),shape=as.factor(popIndex)),size=3)+ # geom_errorbar(aes(n1,ymin=mean-confidenceInterval/2,ymax=mean+confidenceInterval/2,color=as.factor(popIndex)),width=0.07)+ geom_line(aes(n1,mean,color=as.factor(popIndex),lty=as.factor(popIndex)),size=0.5)+ xlab("Proportion of mimic species\nthat mimic the model B ( )")+ ylab("Per capita predation rate")+ ylim(PlotFigure3.yLim)+ scale_colour_manual(name="Appearance :", breaks=c(2, 1), labels=labelAppearance, values=c("black","#A1A1A1"))+ scale_shape_manual(name="Appearance :", breaks=c(2, 1), labels=labelAppearance, values=c(16,16))+ scale_linetype_manual(name="Appearance :", breaks=c(2, 1), labels=labelAppearance, values=c(1,1))+ theme(panel.background = element_rect(fill='white', colour='black'), axis.title.x = element_text(size=20), axis.title.y = element_text(size=20), legend.text=element_text(size=15), axis.text = element_text(size=13), panel.grid.major = element_line(colour = "#DFDFDF",size=0.003), panel.grid.minor = element_blank(),legend.key = element_blank(), legend.title=element_text(size=15), legend.text.align=0, legend.key.width=unit(2,"line"), legend.key.height=unit(3,"line"))+ scale_x_continuous(breaks=xLabel) nbSampledPlot1 <- ggplot(tableAttackRateMean)+ geom_point(aes(n1,meanNb,color=as.factor(popIndex),shape=as.factor(popIndex)),size=3)+ geom_line(aes(n1,meanNb,color=as.factor(popIndex),lty=as.factor(popIndex)),size=0.5)+ xlab("Proportion of mimic species\nthat mimic the model B ( )")+ ylab("Number of prey sampled")+ ylim(PlotFigure3.yLimNbSampled)+ scale_colour_manual(name="Appearance :", breaks=c(1, 2), labels=labelAppearance, values=c("black","#A1A1A1"))+ scale_shape_manual(name="Appearance :", breaks=c(1, 2), labels=labelAppearance, values=c(16,16))+ scale_linetype_manual(name="Appearance :", breaks=c(1, 2), labels=labelAppearance, values=c(1,1))+ theme(panel.background = element_rect(fill='white', colour='black'), axis.title.x = element_text(size=20), axis.title.y = element_text(size=20), legend.text=element_text(size=15), axis.text = element_text(size=13), panel.grid.major = element_line(colour = "#DFDFDF",size=0.003), panel.grid.minor = element_blank(),legend.key = element_blank(), legend.title=element_text(size=15), legend.text.align=0, legend.key.width=unit(2,"line"))+ scale_x_continuous(breaks=xLabel) # # nbSampledModelPlot1 <- # ggplot(tableAttackRateMean[tableAttackRateMean$popIndex==1,])+ # geom_point(aes(n1,meanNb3,color=as.factor(popIndex),shape=as.factor(popIndex)),size=3)+ # geom_line(aes(n1,meanNb3,color=as.factor(popIndex),lty=as.factor(popIndex)),size=0.5)+ # xlab("Proportion of mimic species\nthat mimic the model 2 ( )")+ # ylab("Number of models sampled")+ # ylim(PlotFigure2.yLimNbSampledModel)+ # scale_colour_manual(name="Appearance :", # breaks=c(3, popNb), # labels=labelPhenotype, # values=c("black","#A1A1A1"))+ # scale_shape_manual(name="Appearance :", # breaks=c(3, popNb), # labels=labelPhenotype, # values=c(16,16))+ # scale_linetype_manual(name="Appearance :", # breaks=c(3, popNb), # labels=labelPhenotype, # values=c(1,1))+ # theme(legend.position="none", # panel.background = element_rect(fill='white', colour='black'), # axis.title.x = element_text(size=20), # axis.title.y = element_text(size=20), # legend.text=element_text(size=15), # axis.text = element_text(size=13), # panel.grid.major = element_line(colour = "#DFDFDF",size=0.003), # panel.grid.minor = element_blank(),legend.key = element_blank(), # legend.title=element_text(size=15), # legend.text.align=0, # legend.key.width=unit(2,"line"))+ # scale_x_continuous(breaks=xLabel) # if(ArrowFig3==TRUE){ dataArrow <- data.frame(yArrowFig31I,xArrowFig3F,yArrowFig31E) attPlot1 <- attPlot1 + geom_segment(data=dataArrow, aes(x = yArrowFig31I, y = xArrowFig3F , xend = yArrowFig31E, yend = xArrowFig3F), color=colorArrowFig31, arrow = arrow(length = unit(0.4, "cm")))+ geom_hline(yintercept=LineFig3) } if(PointFig3==TRUE){ dataPoint <- data.frame(XPoint,YPoint) attPlot1 <- attPlot1 + geom_point(data=dataPoint,aes(x=XPoint,y=YPoint),size=5) } if(ArrowFig3Bis==TRUE){ dataArrow2 <- data.frame(xArrowFig3FBis,yArrowFig31IBis,yArrowFig31EBis) attPlot1 <- attPlot1 + geom_segment(data=dataArrow2, aes(x = xArrowFig3FBis, y = yArrowFig31IBis, xend = xArrowFig3FBis, yend = yArrowFig31EBis), color=colorArrowFig32, arrow = arrow(length = unit(0.4, "cm"))) } grid.arrange(attPlot1 + theme(legend.position="none")) # attPlot1 # grid.arrange(attPlot1) # grid.arrange(nbSampledPlot1) # mylegend <- g_legend(attPlot1) # grid.arrange(arrangeGrob(attPlot1 + theme(legend.position="none"), # nrow=1), # mylegend, nrow=2,heights=c(3, 1)) } if(PlotFigure4==TRUE){ setwd("~/ProjetsCEFE/OptimalSamplingStrategyDifferentToxicity") PlotFigure4.f3Vec = seq(PlotFigure4.f3Min,PlotFigure4.f3Max,by=(PlotFigure4.f3Max - PlotFigure4.f3Min)/PlotFigure4.nbStepsGraph) stepf3 = (PlotFigure4.f3Max - PlotFigure4.f3Min)/PlotFigure4.nbStepsGraph # PlotFigure4.f2f3Vec = seq(PlotFigure4.f2f3Min,PlotFigure4.f2f3Max,by=(PlotFigure4.f2f3Max - PlotFigure4.f2f3Min)/PlotFigure4.nbStepsGraphY) PlotFigure4.f2f3Vec = seq(PlotFigure4.f2f3Max,PlotFigure4.f2f3Min,by=-0.02375) #HERE stepf2f3 =(PlotFigure4.f2f3Max - PlotFigure4.f2f3Min)/PlotFigure4.nbStepsGraphY if(PlotFigure4.loadTable==TRUE){ if(length(PlotFigure4.fileName)==1){ tableDataEvolutionaryOutcome <- read.csv(paste(PlotFigure4.fileName,".csv",sep=""),header=TRUE)[,-1] }else{ tableDataEvolutionaryOutcome <- read.csv(paste(PlotFigure4.fileName[1],".csv",sep=""),header=TRUE)[,-1] tableDataEvolutionaryOutcome <- tableDataEvolutionaryOutcome[!is.na(tableDataEvolutionaryOutcome$numberSwitch),] for(i in 2:length(PlotFigure4.fileName)){ tableDataEvolutionaryOutcomeOne <- read.csv(paste(PlotFigure4.fileName[i],".csv",sep=""),header=TRUE)[,-1] tableDataEvolutionaryOutcomeOne <- tableDataEvolutionaryOutcomeOne[!is.na(tableDataEvolutionaryOutcomeOne$numberSwitch),] tableDataEvolutionaryOutcome <- rbind(tableDataEvolutionaryOutcome, tableDataEvolutionaryOutcomeOne) } } }else{ if(PlotFigure4.rePlot==FALSE){ numberSwitchVec <- c() numberStableSwitchVec <- c() numberInstableSwitchVec <- c() EvolutionaryDynamicsVec <- c() equilibriumInvasibilityVec <- c() f3Vec <- c() f2f3Vec <- c() for(PlotFigure4.f3 in PlotFigure4.f3Vec){ for(f2f3 in PlotFigure4.f2f3Vec){ if(f2f3<=limf2f3Max & f2f3>limf2f3Min){ PlotFigure4.f2 = PlotFigure4.f3 * f2f3 f1=rnorm(n=PlotFigure4.nbCombinationf,mean=PlotFigure4.f1,sd=PlotFigure4.fSd) f2=rnorm(n=PlotFigure4.nbCombinationf,mean=PlotFigure4.f2,sd=PlotFigure4.fSd) f3=rnorm(n=PlotFigure4.nbCombinationf,mean=PlotFigure4.f3,sd=PlotFigure4.fSd) f1[f1>10]=10 f2[f2>10]=10 f3[f3>10]=10 f1[f1< -10]=-10 f2[f2< -10]=-10 f3[f3< -10]=-10 if(PlotFigure4.f1>0){ f1[f1<0.0]=0.0 }else{ f1[f1>0.0]=0.0 } if(PlotFigure4.f2>0){ f2[f2<0.0]=0.0 }else{ f2[f2>0.0]=0.0 } if(PlotFigure4.f3>0){ f3[f3<0.0]=0.0 }else{ f3[f3>0.0]=0.0 } n1Vec <- c() attackRateVec <- c() attackRateVec2 <- c() attackRateVec3 <- c() NbPreyEatenVec <- c() NbPreyEatenVec2 <- c() NbPreyEatenVec3 <- c() popIndexVec <- c() # if 1 : bad prey, 2: very bad prey pop PlotFigure4.nMutant = seq(0,PlotFigure4.n1,by=PlotFigure4.stepNbMutant) # number mutants bad prey's phenotype to very bad prey's phenotype for(n1 in PlotFigure4.nMutant){ for(pop in 1:2){ if(pop==1){ preyCommunity <- c(rep(2,PlotFigure4.n2),rep(PlotFigure4.mutantTypeAdded-1,PlotFigure4.n1-n1)) # 0 : profitable prey, nbGoodPrey = length(which(preyCommunity==0)) nbBadPrey = length(which(preyCommunity==1)) nbVeryBadPrey = length(which(preyCommunity==2)) }else{ preyCommunity <- c(rep(PlotFigure4.mutantTypeAdded-1,n1),rep(2,PlotFigure4.n2)) # 0 : profitable prey, nbGoodPrey = length(which(preyCommunity==0)) nbBadPrey = length(which(preyCommunity==1)) nbVeryBadPrey = length(which(preyCommunity==2)) } PlotFigure4.bigN = length(preyCommunity) listA <- list() for(predType in 1:length(f1)){ # print(paste(f3[predType],f2[predType],f1[predType],sep=" ")) A = SamplingStrategyMatrix(PlotFigure4.alpha, floor(PlotFigure4.propCommunityExperienced*length(preyCommunity)), f1[predType], -f2[predType], -f3[predType]) listA[[length(listA)+1]] <- A } for(rep in 1:PlotFigure4.nbRepForward){ predType = sample(1:length(f1),1,replace=T) preyCommunityOnePred <- sample(preyCommunity,size=floor(PlotFigure4.propCommunityExperienced*length(preyCommunity))) # encounter in random order predKnowledge.ri1 = 0 # profitable prey predKnowledge.ri2 = 0 # toxic prey (c1) predKnowledge.ri3 = 0 # toxic prey (c2) nbPreyEaten = 0 nbPreyEaten1 = 0 nbPreyEaten2 = 0 nbPreyEaten3 = 0 for(prey in preyCommunityOnePred){ optBehavior = listA[[predType]][predKnowledge.ri1+1, predKnowledge.ri2+1, predKnowledge.ri3+1] if(optBehavior==1){ # attack nbPreyEaten = nbPreyEaten + 1 if(prey==0){ predKnowledge.ri1 = predKnowledge.ri1 + 1 nbPreyEaten1 = nbPreyEaten1 + 1 }else if(prey==1){ predKnowledge.ri2 = predKnowledge.ri2 + 1 nbPreyEaten2 = nbPreyEaten2 + 1 }else if(prey==2){ predKnowledge.ri3 = predKnowledge.ri3 + 1 nbPreyEaten3 = nbPreyEaten3 + 1 } } } popIndexVec <- c(popIndexVec, pop) n1Vec <- c(n1Vec, PlotFigure4.n2-n1) attackRateVec <- c(attackRateVec, nbPreyEaten/PlotFigure4.bigN) if(PlotFigure4.mutantTypeAdded==1){ attackRateVec2 <- c(attackRateVec2, nbPreyEaten1/nbGoodPrey) }else if(PlotFigure4.mutantTypeAdded==2){ attackRateVec2 <- c(attackRateVec2, nbPreyEaten2/nbBadPrey) } attackRateVec3 <- c(attackRateVec3, nbPreyEaten3/nbVeryBadPrey) NbPreyEatenVec <- c(NbPreyEatenVec,nbPreyEaten) NbPreyEatenVec2 <- c(NbPreyEatenVec2,nbPreyEaten2) NbPreyEatenVec3 <- c(NbPreyEatenVec3,nbPreyEaten3) } } } tableAttackRate = data.frame(n1=n1Vec,attackRate=attackRateVec,attackRate2=attackRateVec2,attackRate3=attackRateVec3,NbPreyEaten=NbPreyEatenVec,NbPreyEaten2=NbPreyEatenVec2,NbPreyEaten3=NbPreyEatenVec3,popIndex=popIndexVec) tableAttackRateMean = ddply(tableAttackRate, .(n1,popIndex), summarize, mean = mean(attackRate,na.rm=T)) diffAttackRate = tableAttackRateMean$mean[tableAttackRateMean$popIndex==1] - tableAttackRateMean$mean[tableAttackRateMean$popIndex==2] diffAttackRate = diffAttackRate[diffAttackRate!=0] switchdiffAttackRate = rep(0, length(diffAttackRate)) if(length(diffAttackRate)!=0){ equilibriumFound = FALSE for(i in 1:(length(diffAttackRate)-1)){ if(diffAttackRate[i] * diffAttackRate[i+1]<0){ if(diffAttackRate[i]<0){ switchdiffAttackRate[i] = 1 if(equilibriumFound==FALSE){ equilibriumInvasibility = (PlotFigure4.stepNbMutant * i + PlotFigure4.stepNbMutant * (i+1))/2/PlotFigure4.n1 } }else{ switchdiffAttackRate[i] = -1 if(equilibriumFound==FALSE){ equilibriumInvasibility = 0 } } equilibriumFound = TRUE } } numberSwitch = sum(abs(switchdiffAttackRate)) numberStableSwitch = length(which(switchdiffAttackRate==1)) numberInstableSwitch = length(which(switchdiffAttackRate==-1)) if(numberSwitch==1){ if(diffAttackRate[1]<0){ EvolutionaryDynamics = "Polymorphism" }else{ EvolutionaryDynamics = "Monomorphism" } }else{ EvolutionaryDynamics = "Unknown" } }else{ equilibriumInvasibility = NA numberSwitch = 0 numberStableSwitch = 0 numberInstableSwitch = 0 EvolutionaryDynamics = "Drift" } print(paste(PlotFigure4.f3,PlotFigure4.f2," : ",numberSwitch, numberStableSwitch, numberInstableSwitch, EvolutionaryDynamics,equilibriumInvasibility,sep=" ")) equilibriumInvasibilityVec <- c(equilibriumInvasibilityVec, equilibriumInvasibility) # print(paste(equilibriumInvasibility,equilibriumInvasibilityVec,length(equilibriumInvasibilityVec),sep=" ; ")) numberSwitchVec <- c(numberSwitchVec, numberSwitch) numberStableSwitchVec <- c(numberStableSwitchVec, numberStableSwitch) numberInstableSwitchVec <- c(numberInstableSwitchVec,numberInstableSwitch) EvolutionaryDynamicsVec <- c(EvolutionaryDynamicsVec,EvolutionaryDynamics) f3Vec <- c(f3Vec, PlotFigure4.f3) f2f3Vec <- c(f2f3Vec, f2f3) } } #END } tableDataEvolutionaryOutcome = data.frame(f3=f3Vec, f2f3=f2f3Vec, numberSwitch=numberSwitchVec, numberStableSwitch=numberStableSwitchVec, numberInstableSwitch=numberInstableSwitchVec, EvolutionaryDynamics=EvolutionaryDynamicsVec, equilibriumInvasibility=equilibriumInvasibilityVec) } } if(PlotFigure4.saveTable==TRUE){ write.csv(tableDataEvolutionaryOutcome, file = paste(PlotFigure4.fileName,".csv",sep="")) } # plotEvolutionaryOutcome <- ggplot(tableDataEvolutionaryOutcome)+ # geom_rect(mapping=aes(xmin= f3-stepf3/2, xmax= f3+stepf3/2, ymin=f2f3-stepf2f3/2, ymax=f2f3+stepf2f3/2, # fill=EvolutionaryDynamics, color=EvolutionaryDynamics)) + # scale_fill_discrete(name="Evolutionary outcome :")+ # scale_color_discrete(name="Evolutionary outcome :")+ # xlab(expression(f[3]))+ # ylab(expression(f[2]/f[3])) # # plotStableEquilibrium<- ggplot(tableDataEvolutionaryOutcome)+ # geom_rect(mapping=aes(xmin= f3-stepf3/2, xmax= f3+stepf3/2, ymin=f2f3-stepf2f3/2, ymax=f2f3+stepf2f3/2, # fill=numberStableSwitch, color=numberStableSwitch)) + # scale_fill_gradient2(name="Number of\nstable equilibrium :", low="red", high="blue", limits=c(0, NA))+ # scale_color_gradient2(name="Number of\nstable equilibrium :", low="red", high="blue", limits=c(0, NA))+ # xlab(expression(f[3]))+ # ylab(expression(f[2]/f[3])) # # plotUnstableEquilibrium<- ggplot(tableDataEvolutionaryOutcome)+ # geom_rect(mapping=aes(xmin= f3-stepf3/2, xmax= f3+stepf3/2, ymin=f2f3-stepf2f3/2, ymax=f2f3+stepf2f3/2, # fill=numberInstableSwitch, color=numberInstableSwitch)) + # scale_fill_gradient2(name="Number of\nunstable equilibrium :", low="red", high="blue", limits=c(0, NA))+ # scale_color_gradient2(name="Number of\nunstable equilibrium :", low="red", high="blue", limits=c(0, NA))+ # xlab(expression(f[3]))+ # ylab(expression(f[2]/f[3])) xLabplotEquilibriumInvasibility = bquote(atop("Profitability of most", "unprofitable prey (" ~ f[1]~")")) # xLabplotEquilibriumInvasibility = bquote(atop("Profitability of most unprofitable prey (" ~ f[1]~")")) yLabplotEquilibriumInvasibility = bquote(atop("Relative unprofitability of", "moderately unprofitable prey (" ~ f[2]/f[1] ~")")) # oneLine=data.frame(f3=1,f2f3=1,numberSwitch=NA,numberStableSwitch=NA,numberInstableSwitch=NA,EvolutionaryDynamics=NA,equilibriumInvasibility=0) # for(PlotFigure4.f3 in PlotFigure4.f3Vec){ # for(f2f3 in PlotFigure4.f2f3Vec){ # if( # f2f3>0.5 & # dim(tableDataEvolutionaryOutcome[tableDataEvolutionaryOutcome$f3==PlotFigure4.f3 & # tableDataEvolutionaryOutcome$f2f3 ==f2f3,])[1]==0){ # oneLine$f3 = PlotFigure4.f3 # oneLine$f2f3 = f2f3 # tableDataEvolutionaryOutcome <- rbind(tableDataEvolutionaryOutcome,oneLine) # } # } # } tableDataEvolutionaryOutcome$equilibriumInvasibility[tableDataEvolutionaryOutcome$equilibriumInvasibility>0.5] <- 0.5 summary(tableDataEvolutionaryOutcome) plotEquilibriumInvasibility <- ggplot(tableDataEvolutionaryOutcome)+ geom_rect(mapping=aes(xmin= f3-stepf3/2, xmax= f3+stepf3/2, ymin=f2f3-stepf2f3/2, ymax=f2f3+stepf2f3/2, fill=equilibriumInvasibility, color=equilibriumInvasibility)) + scale_fill_gradient(name="q at equilibrium\n", low=PlotFigure4.colorLow, high=PlotFigure4.colorHigh, limits=c(0, PlotFigure4.topLimit))+ scale_color_gradient(name="q at equilibrium\n", low=PlotFigure4.colorLow, high=PlotFigure4.colorHigh, limits=c(0, PlotFigure4.topLimit))+ xlab(xLabplotEquilibriumInvasibility)+ ylab(yLabplotEquilibriumInvasibility)+ theme(panel.background = element_blank(), axis.title.x = element_text(size=15), axis.title.y = element_text(size=15), legend.position="bottom") # multiplot(plotEvolutionaryOutcome) # multiplot(plotStableEquilibrium,plotUnstableEquilibrium,cols=2) # multiplot(plotStableEquilibrium,plotUnstableEquilibrium,plotEvolutionaryOutcome,plotEquilibriumInvasibility,cols=2) # grid.arrange(plotEquilibriumInvasibility + theme(legend.position="none")) grid.arrange(plotEquilibriumInvasibility + xlab("") + ylab("")+ theme(legend.position="none")) # grid.arrange(plotEquilibriumInvasibility) } print(proc.time() - ptm)