##### READING DATA, THIS ALL ARE IN THE VINCENZI_ET_AL_JFB13.ZIP FILE ##### ##Unzip the file in the current directory####### L_age = read.table('Lage.txt') S_age = read.table('Sage.txt') l_agethree_med_evo = read.table("l_agethree_med_evo.txt") l_agethree_med_fix = read.table("l_agethree_med_fix.txt") is.na(l_agethree_med_evo) = l_agethree_med_evo==0 is.na(l_agethree_med_fix) = l_agethree_med_fix==0 timetoext_evo = read.table("timetoext_evo.txt") timetoext_fix = read.table("timetoext_fix.txt") propage1_fix = read.table("propage1_fix.txt") propage1_evo = read.table("propage1_evo.txt") n_etauno_evo = read.table("n_etauno_evo.txt") n_etadue_evo = read.table("n_etadue_evo.txt") n_etatre_evo = read.table("n_etatre_evo.txt") n_etaquattro_evo = read.table("n_etaquattro_evo.txt") n_etauno_fix = read.table("n_etauno_fix.txt") n_etadue_fix = read.table("n_etadue_fix.txt") n_etatre_fix = read.table("n_etatre_fix.txt") n_etaquattro_fix = read.table("n_etaquattro_fix.txt") n_tot_evo = read.table("n_tot_evo.txt") n_tot_fix = read.table("n_tot_fix.txt") flood_evo = read.table("flood_evo.txt") flood_fix = read.table("flood_fix.txt") phimaxevo = read.table("phimaxevo.txt") phiminevo = read.table("phiminevo.txt") rangephievo = read.table("rangephievo.txt") phimaxfix = read.table("phimaxfix.txt") phiminfix = read.table("phiminfix.txt") rangephifix = read.table("rangephifix.txt") require(MASS) ####### INCREASE OF MEDIAN LENGTH AT AGE 3 ############ synth_ts = function(l3med,text,year_init = 150,year_final = 300) { ### This computes the slope of the linear model with # median length at age 3 as response and year as independent variable #l3med = l3med[,which(text==0)] l3med = l3med[,-which(text>0 & text0 & textthresh_up) { min_up_prov = min(which(n_tot_pers[(min_down[cont_mindown]:rowmax),aa]>thresh_up)) if (min_up_prov<10){ min_up[cont_minup] = min_up_prov cont_minup =cont_minup+1 } } cont_mindown = cont_mindown+1 } } min_up = min_up[min_up!=0] return(min_up) } res_fix1050 = synth_res(n_tot_fix, timetoext_fix,thresh_down = 10,thresh_up = 50,year_start = 150) # year from 10 fish to 50 fo FIX res_evo1050 = synth_res(n_tot_evo, timetoext_evo,thresh_down = 10,thresh_up = 50,year_start = 150) # year from 10 fish to 50 fo EVO res_fix2050 = synth_res(n_tot_fix, timetoext_fix,thresh_down = 20,thresh_up = 50,year_start = 150) # year from 10 fish to 50 fo FIX res_evo2050 = synth_res(n_tot_evo, timetoext_evo,thresh_down = 20,thresh_up = 50,year_start = 150) # year from 10 fish to 50 fo EVO res_fix550 = synth_res(n_tot_fix, timetoext_fix,thresh_down = 5,thresh_up = 50,year_start = 150) # year from 10 fish to 50 fo FIX res_evo550 = synth_res(n_tot_evo, timetoext_evo,thresh_down = 5,thresh_up = 50,year_start = 150) # year from 10 fish to 50 fo EVO wilcox.test(res_fix550,res_evo550, correct=FALSE) # some wilcoxon tests to check for differences between FIX and EVO scenarios wilcox.test(res_fix2050,res_evo2050, correct=FALSE) wilcox.test(res_fix1050,res_evo1050, correct=FALSE) #####PLOTS###### # ALL PLOTS ARE OPTIMIZED FOR EPS FORMAT, DIMENSIONS REFER TO THE WIDHT AND HEIGHT OF THE RSTUDIO PLOT WINDOW IN MY IMAC ####### #### FIG 2 #### #Length at age (a) and Survival at age (b) This is Fig. 2 dimensions are 1000 and 600 par(mar=c(5,4,2,1),mfrow=c(1,2),oma = c(1,1.3,0,0)) seq.x = seq(1,25,0.1) seq.x.s = seq(0,25,1) lwd.l = 3 lty.l = c(1,2,3,4) cex.lab.l =1 cex.axis.l = 1 xlab.l = "Age (year)" # Length for (i in 1:ncol(L_age)) { if (i == 1){ #dev.new(width=9, height=6) plot(L_age[,i] ~ seq.x ,xlab="",ylab="",type="l", lwd=lwd.l,lty = lty.l[i], ylim=c(0,400),cex.axis = cex.axis.l,cex.lab =cex.lab.l) } else { points(L_age[,i] ~ seq.x,type="l", lwd=lwd.l,lty = lty.l[i]) } } text(3,400,"a",cex=1.3) legend(15,120,c("0.5","1","1.5","2.5"),lty=c(1,2,3,4),lwd = c(3,3,3,3),y.intersp = 1.2,seg.len=3, bty ='n',cex = 1.1) ## Age for (i in 1:ncol(S_age)) { if (i == 1){ #dev.new(width=9, height=6) plot(S_age[,i] ~ seq.x.s ,xlab="",ylab=expression(paste(sigma)),type="l", lwd=lwd.l,lty = lty.l[i], ylim=c(0,1),cex.axis = cex.axis.l,cex.lab =cex.lab.l,cex.lab = 1.3) } else { points(S_age[,i] ~ seq.x.s,type="l", lwd=lwd.l,lty = lty.l[i]) } } text(3,1,"b",cex=1.3) mtext(xlab.l, outer=T,side=1,cex=1.1,padj=-1.5) mtext("Length (mm)", outer=T,side=2,cex=1.1,padj=1, adj =0.55) ##FIG. 6 ###### ylab.prop = expression(paste(beta[prop]," ",10^4)) ylab.length = expression(paste(beta[L])) ylab.ntot = "Number of fish" #all together #### #dimensions are 1000 800 labels etc optimized for eps par(mar=c(5,6.2,2,1.5),mfrow=c(2,3)) ylab.prop = expression(paste(beta['prop,200']," ",10^4)) boxplot(synthpropage_fix_200[,2]*10^4,synthpropage_evo_200[,2]*10^4,names=c("FL-FIX","FL-EVO"), ylab=ylab.prop,cex.lab=1.4,cex.axis=1) text(0.6,35,"a",cex=1.3) ylab.length = expression(paste(beta['L,200'])) boxplot(synthl3_fix_200[,2],synthl3_evo_200[,2],names=c("FL-FIX","FL-EVO"), ylab=ylab.length,cex.lab=1.4,cex.axis=1) text(0.6,0.95,"b",cex=1.3) ylab.ntot = "Number of fish at year 200" boxplot(ntot_med_fix_200,ntot_med_evo_200,names=c("FL-FIX","FL-EVO"), ylab=ylab.ntot,cex.lab=1.2,cex.axis=1,ylim = c(30,180),yaxt='n') text(0.6,177,"c",cex=1.3) axis(side = 2, at = seq(40,180,20),cex.axis = 1.) ylab.prop = expression(paste(beta['prop,300']," ",10^4)) boxplot(synthpropage_fix_300[,2]*10^4,synthpropage_evo_300[,2]*10^4,names=c("FL-FIX","FL-EVO"), ylab=ylab.prop,cex.lab=1.4,cex.axis=1.2) text(0.6,9,"d",cex=1.3) ylab.length = expression(paste(beta['L,300'])) boxplot(synthl3_fix_300[,2],synthl3_evo_300[,2],names=c("FL-FIX","FL-EVO"), ylab=ylab.length,cex.lab=1.4,cex.axis=1) text(0.6,0.25,"e",cex=1.3) ylab.ntot = "Number of fish at year 300" boxplot(ntot_med_fix_300,ntot_med_evo_300,names=c("FL-FIX","FL-EVO"), ylab=ylab.ntot,cex.lab=1.2,cex.axis=1,yaxt='n',ylim = c(30,180)) text(0.6,177,"f",cex=1.3) axis(side = 2, at = seq(40,180,20),cex.axis = 1) ######### Fig. 7 ###### ###########Median length at age 3 from year 250 to 300, Median Proportion of age-1 from year 250 to 300 ######## #dimensions of plot are 900 500 Fig 7 par(mfrow=c(1,2),mar=c(3,5,2,2),oma=c(0,0,0,0)) boxplot(apply(l_agethree_med_fix[250:300,],2,median,na.rm=T),apply(l_agethree_med_evo[250:300,],2,median,na.rm=T), names=c("FL-FIX","FL-EVO"),ylab="Median length at age 3 (mm)",cex.lab=1,cex.axis = 0.9) text(0.6,268,"a",cex=1.2) ap_propage_fix_300 = apply(propage1_fix[250:300,],2,median,na.rm=T) ap_propage_fix_300 = as.matrix(ap_propage_fix_300[ap_propage_fix_300>0.2]) ap_propage_evo_300 = apply(propage1_evo[250:300,],2,median,na.rm=T) ap_propage_evo_300 = as.matrix(ap_propage_evo_300[ap_propage_evo_300>0.2]) #ap_propage_300 = data.frame(evop = ap_propage_evo_300,fixp = ap_propage_fix_300) boxplot(as.numeric(ap_propage_fix_300),as.numeric(ap_propage_evo_300),names=c("FL-FIX","FL-EVO"), ylab="Median proportion of age-1 fish",cex.lab=1,cex.axis = 0.9) text(0.6,0.64,"b",cex=1.2) ## some statistics ##### t.test(ap_propage_fix_300,ap_propage_evo_300) t.test(apply(l_agethree_med_fix[250:300,],2,median,na.rm=T),apply(l_agethree_med_evo[250:300,],2,median,na.rm=T)) t.test(ntot_med_fix_300,ntot_med_evo_300) t.test(ntot_med_fix_200,ntot_med_evo_200) ###### FIG. 3 ########## ###4 simulations, 2 for Evo, 2 for fixed #### #####function to plot a number of simulations together as graylines#### plot_nad = function(ntot, start_sim=1, end_sim=100,cex.axis=1) for (bb in start_sim:end_sim) { if (bb == start_sim){ #dev.new(width=9, height=6) plot(n_tot_fix[,bb],xlab="",ylab="",type="l", lwd=lwed.traj.multi,col="gray90",ylim=c(0,300),cex.axis = cex.axis) } else { points(n_tot_fix[,bb],type="l", lwd=lwed.traj.multi, col="gray80",cex.axis = cex.axis) } } ######## #### moving average for the simulations number 32 and 167 for scenario FIX and 20 and 166 for scenario EVO ###### f10 =rep(1/20,20) seq_filter=c(141:300) n_tot_evo20.ts = ts(n_tot_evo[,20]) n_tot_evo20.filter=filter(n_tot_evo20.ts[141:300],f10,method="convolution") n_tot_fix32.ts = ts(n_tot_evo[,32]) n_tot_fix32.filter=filter(n_tot_fix32.ts[141:300],f10,method="convolution") n_tot_evo166.ts = ts(n_tot_evo[,166]) n_tot_evo166.filter=filter(n_tot_evo166.ts[141:300],f10,method="convolution") n_tot_fix167.ts = ts(n_tot_fix[,167]) n_tot_fix167.filter=filter(n_tot_fix167.ts[141:300],f10,method="convolution") ## start plotting ###### lwd.traj = 0.8 lwd.mov = 3 lwed.traj.multi = 0.5 par(mfrow=c(3,2),oma=c(4,5,0,0),mar=c(3,3,2,2)) ylimntot = c(0,300) plot(n_tot_evo[,20],xlab="",ylab="",type="l", lwd=lwd.traj,ylim=ylimntot,cex.axis = 0.8) points(n_tot_evo20.filter ~ seq_filter, lwd=lwd.mov,type="l") text(30,280,"a",cex=1.5) #mtext("FL-EVO",outer=F,cex=1.2,at=c(40),side=1,padj=-1) mtext("FL-EVO",outer=F,cex=1,at=c(40),side=1,padj=-2) abline(v=150,lty=2) plot(n_tot_fix[,32],xlab="",ylab="",type="l", lwd=lwd.traj,ylim=ylimntot,cex.axis = 0.8) points(n_tot_fix32.filter ~ seq_filter, lwd=lwd.mov,type="l") text(30,280,"b",cex=1.5) mtext("FL-FIX",outer=F,cex=1,at=c(40),side=1,padj=-2) abline(v=150,lty=2) plot(n_tot_evo[,166],xlab="",ylab="",type="l", lwd=lwd.traj,ylim=ylimntot,cex.axis = 0.8) points(n_tot_evo166.filter ~ seq_filter, lwd=lwd.mov,type="l") text(30,280,"c",cex=1.5) mtext("FL-EVO",outer=F,cex=1,at=c(40),side=1,padj=-2) abline(v=150,lty=2) plot(n_tot_fix[,167],xlab="",ylab="",type="l", lwd=lwd.traj,ylim=ylimntot,cex.axis = 0.8) points(n_tot_fix167.filter ~ seq_filter, lwd=lwd.mov,type="l") text(30,280,"d",cex=1.5) mtext("FL-FIX",outer=F,cex=1,at=c(40),side=1,padj=-2) abline(v=150,lty=2) plot_nad(n_tot_evo,start_sim=1, end_sim=10,cex.axis = 0.8) text(30,280,"e",cex=1.5) mtext("FL-EVO",outer=F,cex=1,at=c(40),side=1,padj=-2) abline(v=150,lty=2) plot_nad(n_tot_fix,start_sim=1, end_sim=10,cex.axis = 0.8) text(30,280,"f",cex=1.5) mtext("FL-FIX",outer=F,cex=1,at=c(40),side=1,padj=-2) abline(v=150,lty=2) ylab.not = "Number of trout" xlab.not = "Year of simulation" mtext(ylab.not, outer=T,side=2,cex=1.1,padj=-1.5) mtext(xlab.not, outer=T,side=1,cex=1.1,padj=1.5) ### FIG 5 ######## ###lenght a age 3 fix and evo ########## ### rolling average for median length at age 3, same simulations as before ##### require(zoo) rolling_width = 20 l_agethree_med_evo20.roll = rollapply(l_agethree_med_evo20.ts,rolling_width,mean, na.rm = TRUE) l_agethree_med_fix32.roll = rollapply(l_agethree_med_fix32.ts,rolling_width,mean, na.rm = TRUE) l_agethree_med_evo166.roll = rollapply(l_agethree_med_evo166.ts,rolling_width,mean, na.rm = TRUE) l_agethree_med_fix167.roll = rollapply(l_agethree_med_fix167.ts,rolling_width,mean, na.rm = TRUE) #Fig 5 dimnesions 1000 650 lwd.traj = 0.5 lwd.mov = 4 par(mfrow=c(2,2),oma=c(4,4,0,0),mar=c(3,2.5,2,2)) ylimntot = c(120,300) plot(l_agethree_med_evo[,20],xlab="",ylab="",type="l", lwd=lwd.traj,ylim=ylimntot) #points(l_agethree_med_evo20.filter ~ seq_filter, lwd=2.3,type="l") points(l_agethree_med_evo20.roll[150:length(l_agethree_med_evo20.roll)] ~ c(150:length(l_agethree_med_evo20.roll)) ,lwd=lwd.mov,type="l") text(30,280,"a",cex=1.2) abline(v=150,lty=2) mtext("FL-EVO",outer=F,cex=1,at=c(50),side=1,padj=-2) plot(l_agethree_med_fix[,32],xlab="",ylab="",type="l", lwd=lwd.traj,ylim=ylimntot) abline(v=150,lty=2) #points(l_agethree_med_fix32.filter ~ seq_filter, lwd=2.3,type="l") points(l_agethree_med_fix32.roll[150:length(l_agethree_med_fix32.roll)] ~ c(150:length(l_agethree_med_fix32.roll)) ,lwd=lwd.mov,type="l") text(30,280,"b",cex=1.2) mtext("FL-FIX",outer=F,cex=1,at=c(50),side=1,padj=-2) plot(l_agethree_med_evo[,166],xlab="",ylab="",type="l", lwd=lwd.traj,ylim=ylimntot) #points(l_agethree_med_evo166.filter ~ seq_filter, lwd=2.3,type="l") text(30,280,"c",cex=1.2) abline(v=150,lty=2) #points(l_agethree_med_fix166.filter ~ seq_filter, lwd=2.3,type="l") points(l_agethree_med_evo166.roll[150:length(l_agethree_med_evo166.roll)] ~ c(150:length(l_agethree_med_evo166.roll)) ,lwd=lwd.mov,type="l") mtext("FL-EVO",outer=F,cex=1,at=c(50),side=1,padj=-2) abline(v=150,lty=2) plot(l_agethree_med_fix[,167],xlab="",ylab="",type="l", lwd=lwd.traj,ylim=ylimntot) #points(l_agethree_med_fix167.filter ~ seq_filter, lwd=2.3,type="l") points(l_agethree_med_fix167.roll[150:length(l_agethree_med_fix167.roll)] ~ c(150:length(l_agethree_med_fix167.roll)) ,lwd=lwd.mov,type="l") text(30,280,"d",cex=1.2) abline(v=150,lty=2) mtext("FL-FIX",outer=F,cex=1,at=c(50),side=1,padj=-2) ylab.not = "Median length at age 3 (mm)" xlab.not = "Year of simulation" mtext(ylab.not, outer=T,side=2,cex=1,padj=-1.5) mtext(xlab.not, outer=T,side=1,cex=1,padj=1.5) ### FIG 5 ######## ###proportion of age 1 ########## ### rolling average for proportion of age-1 fish at age 3, same simulations as before ##### propage1_fix propage1_evo require(zoo) propage1_evo20.ts = ts(propage1_evo[,20]) propage1_fix32.ts = ts(propage1_fix[,32]) propage1_evo166.ts = ts(propage1_evo[,166]) propage1_fix167.ts = ts(propage1_fix[,167]) rolling_width = 20 propage1_evo20.roll = rollapply(propage1_evo20.ts,rolling_width,mean, na.rm = TRUE) propage1_fix32.roll = rollapply(propage1_fix32.ts,rolling_width,mean, na.rm = TRUE) propage1_evo166.roll = rollapply(propage1_evo166.ts,rolling_width,mean, na.rm = TRUE) propage1_fix167.roll = rollapply(propage1_fix167.ts,rolling_width,mean, na.rm = TRUE) #### Fig 4 dimensions 1000 650 lwd.traj = 0.5 lwd.mov = 4 par(mfrow=c(2,2),oma=c(4,4,0,0),mar=c(3,2.5,2,2)) ylimntot = c(0,1) plot(propage1_evo[,20],xlab="",ylab="",type="l", lwd=lwd.traj,ylim=ylimntot) points(propage1_evo20.roll[150:length(propage1_evo20.roll)] ~ c(150:length(propage1_evo20.roll)) ,lwd=lwd.mov,type="l") text(30,0.95,"a",cex=1.2) abline(v=150,lty=2) mtext("FL-EVO",outer=F,cex=1,at=c(40),side=1,padj=-2) plot(propage1_fix[,32],xlab="",ylab="",type="l", lwd=lwd.traj,ylim=ylimntot) points(propage1_fix32.roll[150:length(propage1_fix32.roll)] ~ c(150:length(propage1_fix32.roll)) ,lwd=lwd.mov,type="l") text(30,0.95,"b",cex=1.2) abline(v=150,lty=2) mtext("FL-FIX",outer=F,cex=1,at=c(40),side=1,padj=-2) plot(propage1_evo[,166],xlab="",ylab="",type="l", lwd=lwd.traj,ylim=ylimntot) points(propage1_evo166.roll[150:length(propage1_evo166.roll)] ~ c(150:length(propage1_evo166.roll)) ,lwd=lwd.mov,type="l") text(30,0.95,"c",cex=1.2) abline(v=150,lty=2) mtext("FL-EVO",outer=F,cex=1,at=c(40),side=1,padj=-2) abline(v=150,lty=2) plot(propage1_fix[,167],xlab="",ylab="",type="l", lwd=lwd.traj,ylim=ylimntot) points(propage1_fix167.roll[150:length(propage1_fix167.roll)] ~ c(150:length(propage1_fix167.roll)) ,lwd=lwd.mov,type="l") mtext("FL-FIX",outer=F,cex=1,at=c(40),side=1,padj=-2) text(30,0.95,"d",cex=1.2) abline(v=150,lty=2) ylab.not = "Proportion of age-1 fish" xlab.not = "Year of simulation" mtext(ylab.not, outer=T,side=2,cex=1,padj=-1.5) mtext(xlab.not, outer=T,side=1,cex=1,padj=1.5) ### FIG. 8 ######## dev.new(width=9,height=10) require(MASS) xlimphimax = c(1.5,3.5) xlimphimin = c(0,2) xlimphirange = c(0,3.2) ylimphimax = c(0,2) ylimphimin = c(0,2) ylimphirange = c(0,2) xlabphimax= expression(paste("max ",phi)) xlabphimin= expression(paste("min ",phi)) xlabphirange= expression(paste("range of ",phi)) #Fig 8 dimensions 800 600 par(mfrow=c(2,3),oma=c(4,5,2,0),mar=c(4,3,2.5,2)) truehist(phimaxevo,nbins=5,col="white", lwd=2,xlim=xlimphimax,ylim=ylimphimax,xlab=xlabphimax,cex.lab=1.2,h=0.4,ylab ='') mtext("FL-EVO",outer=F,cex=0.7,at=c(1.9),side=1,padj=-15.5) truehist(phiminevo,nbins=5,col="white",lwd=2,xlim=xlimphimin, ylim=ylimphimin,ylab="",xlab=xlabphimin,cex.lab=1.2,h=0.4) truehist(rangephievo,nbins=5,col="white", lwd=2,xlim=xlimphirange, ylim = ylimphirange,ylab="",xlab=xlabphirange,cex.lab=1.2,h=0.4) truehist(phimaxfix,nbins=5,col="white", lwd=2,xlim=xlimphimax, ylim=ylimphimax,ylab="Frequency",xlab=xlabphimax,cex.lab=1.2,h=0.4) mtext("FL-FIX",outer=F,cex=0.7,at=c(1.9),side=1,padj=-15.5) truehist(phiminfix,nbins=5,col="white", lwd=2,xlim=xlimphimin, ylim=ylimphimin,ylab="",xlab=xlabphimin,cex.lab=1.2,h=0.4) truehist(rangephifix,nbins=5,col="white", lwd=2,xlim=xlimphirange,ylim = ylimphirange,ylab="",xlab=xlabphirange,cex.lab=1.2,h=0.4) mtext("Frequency", outer=T,side=2,cex=1,padj=-1.5) #### FIG. 9 ####### ### Proportion of replicate that went quasi-extinct (1 to 20 fish remaining at a certain point in the population) in scenarios EVO and FIX ## min_ext = 2 max_ext = 21 cont_quasiext = 1 ext_var_evo = rep(0,length(min_ext:max_ext)) minyearext = 150 for (aa in min_ext:max_ext) { ext_var_evo[cont_quasiext] = sum(apply(n_tot_evo[minyearext:300,],2,min)