## R Code relative to the paper #"Food abundance, kittiwake life histories, and colony dynamics in the Northeastern Pacific: implications of climate change and regime shifts #Simone Vincenzia,b,*, Marc Mangelb,c #a Dipartimento di Elettronica, Informazione e Bioingegneria, Politecnico di Milano, Via Ponzio 34/5, I-20133 Milan, Italy #b Center for Stock Assessment Research and Department of Applied Mathematics and Statistics, University of California, Santa Cruz, CA 95064, USA #c Department of Biology, University of Bergen, P.O. Box 7803, 5020 Bergen, Norway #*Corresponding author. Email: simon.vincenz@gmail.com, Tel: +1 831-420-3936 #load some libraries require(Rlab) require(MASS) require(parallel) # Function for multiplot # Multiple plot function # # ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects) # - cols: Number of columns in layout # - layout: A matrix specifying the layout. If present, 'cols' is ignored. # # If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE), # then plot 1 will go in the upper left, 2 will go in the upper right, and # 3 will go all the way across the bottom. # multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) { require(grid) # Make a list from the ... arguments and plotlist plots <- c(list(...), plotlist) numPlots = length(plots) # If layout is NULL, then use 'cols' to determine layout if (is.null(layout)) { # Make the panel # ncol: Number of columns of plots # nrow: Number of rows needed, calculated from # of cols layout <- matrix(seq(1, cols * ceiling(numPlots/cols)), ncol = cols, nrow = ceiling(numPlots/cols)) } if (numPlots==1) { print(plots[[1]]) } else { # Set up the page grid.newpage() pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout)))) # Make each plot, in the correct location for (i in 1:numPlots) { # Get the i,j matrix positions of the regions that contain this subplot matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE)) print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row, layout.pos.col = matchidx$col)) } } } ### Here below I report the functions and data that were later used # for the model of colony dynamics #Simulating an AR(1) with exponential innovations #Let us consider simulating the AR(1) process {X(t)} defined by X(t) = a*X(t-1) +Z(t) where a is a #constant and the {Z(t)} are independently distributed Exp(1) random variables. #We will start the process off at X(1)=0. #Implementation in R arsim = function (n, alpha,rate,vec.part) { vec <- vector("numeric", n) vec[1] <- vec.part for (i in 2:n) { vec[i] <- alpha * vec[i - 1] + rexp(1, rate) } return(vec) } ### Cumulative Probability or recruting for chicks B and A for ages 3-10 prob_recruit_B = c(0.0474151,0.1243792,0.3960766,0.7257886,0.8836074,0.9785819,0.9974608,0.9999925) prob_recruit_A = c(0.03373131,0.08955652,0.29969203,0.59906623,0.78113076,0.93379513,0.98532027,0.99976028) ############Function to simulate log(CPUE) ################# #correlated = 0 --> drawing from uniforma distribution of log10CPUE between 1 and 3.5 #correlated = 1 --> ARIMA with exponential innovations #correlated = 2 --> drawing from uniforma distribution of log10CPUE between 1 and 3.5/2 log10CPUE.f = function(iter,cpue.min = 1, cpue.max = 3.5, correlated = 0, start.ar = 2, alpha.ar = 0.85, rate.ar = 3) { burn=1000 if (correlated == 0) { CPUE.gen = runif(iter,cpue.min,cpue.max) } else if (correlated ==1) { CPUE.gen = arsim(burn + iter,alpha.ar,rate.ar,start.ar)[(burn+1):(burn+iter)] } else if (correlated ==2) { CPUE.gen = runif(iter,cpue.min,cpue.max/2) } CPUE.gen[CPUE.gen>cpue.max] = cpue.max return(CPUE.gen) } #log10CPUE.f(100,1,3.5,1,0,0.85,3) ########################## Function tha return the new breeders for the year ### kitti.sexualmat.f = function(area.mat,prob_recruit_B,prob_recruit_A, breeding = 6, age = 1, social = 3) { #print(1) potential.breeders = which(area.mat[breeding,]==0 & as.numeric(area.mat[age,])>=3) #& as.numeric(area.mat[age,]!=100) & as.numeric(area.mat[age,])<=10) #extract the potential breeders (columns #print(2) if (length(potential.breeders) > 0 ){ #print(potential.breeders) potential.breeders.A = potential.breeders[which(area.mat[social,potential.breeders]=="A")] #print(potential.breeders.A ) potential.breeders.B = potential.breeders[which(area.mat[social,potential.breeders]=="B")] #print(potential.breeders.B) potential.breeders.rand.A = runif(length(potential.breeders.A)) potential.breeders.rand.B = runif(length(potential.breeders.B)) ###random uniform numbers for bernoulli trial if (length(potential.breeders.A)>0) { potential.breeder.prob.A = prob_recruit_A[as.numeric(area.mat[age,potential.breeders.A])]#predict the probability of breeding for non-breeders A.point = which(potential.breeders.rand.A < potential.breeder.prob.A) if (length(A.point)>0) { new.breeders.A = potential.breeders.A[which(potential.breeders.rand.A < potential.breeder.prob.A)] } else { new.breeders.A = 0 } } if (length(potential.breeders.B)>0) { potential.breeder.prob.B = prob_recruit_B[as.numeric(area.mat[age,potential.breeders.B])] B.point = which(potential.breeders.rand.B < potential.breeder.prob.B) if (length(B.point)>0) { new.breeders.B = potential.breeders.B[which(potential.breeders.rand.B < potential.breeder.prob.B)] } else { new.breeders.B = 0 } } ##where are the new breeders? (winners of bernoulli trial) if (length(potential.breeders.A)>0 & length(potential.breeders.B)>0) { new.breeders = c(new.breeders.A,new.breeders.B) #print(new.breeders.A) #print(new.breeders.B) } if (length(potential.breeders.A)>0 & length(potential.breeders.B)==0) { new.breeders = (new.breeders.A) #print(new.breeders.A) } if (length(potential.breeders.A)==0 & length(potential.breeders.B)>0) { new.breeders = (new.breeders.B) #print(new.breeders.B) } } else { new.breeders = 0 } if (max(new.breeders) == 0) { new.breeders = 0 } return(new.breeders) #return the new breeders (columns of matrix) #return(potential.breeders) } # kitti.sexualmat.f(area,prob_recruit_B,prob_recruit_A) ### function for random drawing of CORT value given prey_avail (that is log10CPUE) kitti.cort.f = function(prey_avail) { #cort.mean = (1 - ((prey_avail/pstar)^gamma_pa))*(cort.max - cort.min) + cort.min cort.meanlog10 = 1.5 + (-0.36) * prey_avail + rnorm(1,0,0.1253) if (cort.meanlog10<0.1) { cort.meanlog10 = 0.1 } cort.sd = runif(1,0.5,2.9) cort.mean.log = log(10^cort.meanlog10) - 0.5 * log(1+(cort.sd/(10^cort.meanlog10)^2)) cort.sd.log = sqrt(log(1+(cort.sd/(10^cort.meanlog10))^2)) cort.return = list(c.mean = cort.mean.log, c.sd = cort.sd.log,c.mean.nat=10^cort.meanlog10,c.mean.log10 = cort.meanlog10) return(cort.return) } # Function for mean colony-wide productivity given max and min productivity of colonies for max #and min log10CORT prod.scenario = function(cortdata,scenario,min_high = 0.3,max_high = 1, min_med = 0.2,max_med = 0.6,min_low = 0.1,max_low = 0.2) { beta_scen1 = (min_high - max_high)/log10(19.95262/1.58) beta_scen2 =(min_med - max_med)/log10(19.95262/1.58) beta_scen3 =(min_low - max_low)/log10(19.95262/1.58) if (scenario == 1) { if (cortdata$c.mean.nat >19.95262) { prod.y = min_high } else if (cortdata$c.mean.nat <=1.58) { prod.y = max_high } else { prod.y = max_high + beta_scen1*log10(cortdata$c.mean.nat/1.58) } } if (scenario == 2) { if (cortdata$c.mean.nat >19.95262) { prod.y = min_med } else if (cortdata$c.mean.nat <=1.58) { prod.y = max_med } else { prod.y = max_med + beta_scen2*log10(cortdata$c.mean.nat/1.58) } } if (scenario == 3) { if (cortdata$c.mean.nat >19.95262) { prod.y = min_low } else if (cortdata$c.mean.nat <=1.58) { prod.y = max_low } else { prod.y = max_low + beta_scen3*log10(cortdata$c.mean.nat/1.58) } } return(prod.y) } #####Function for getting CORT at the individual level given a mean CORT for that year (area is #the matrix of information - every column is a bird) kitti.cort.ind.f = function(area,cort_year) { col.alive = which(area[1,]!=100) area[cort,col.alive] = rlnorm(length(col.alive),cort_year$c.mean,cort_year$c.sd) return(area) } #### Function for getting a possible variance of productivity (values in the range) in a colony given # (1) a mean colony-wide productivity (prod.gen) # (2) p0 + p1 + p2 = 1 # (3) Var(X) = E[X^2] - (E[X])^2 range.var = function(prod.gen) { p0.plus = - prod.gen^2 + 3*prod.gen-2 p1.minus = 2*prod.gen - prod.gen^2 p0p1.plus = prod.gen - prod.gen^2 return(c(max(p0.plus,p0p1.plus),p1.minus)) } #Function that given # (1) mean colony-wide productivity, (2) variance in productivity, (3) number of pairs # computes probability of having 0,1,2 chicks # the function then assign 0,1,2 chicks to kittiwake pairs kitti.prod.ind.f = function(prod.gen,var.prod,couples) { chicks = rep(100,2000) #var.prod = sd.prod^2 for (pos in 1:100) { if (pos==1) { prod.gen.prov = prod.gen var.prod.prov = var.prod } #### computing probability of having 0,1,2 chicks a.mat = matrix(c(-4,-2,-3,-1),2,2) b.mat = c(prod.gen.prov^2 + var.prod.prov - 4, prod.gen.prov - 2) p01 = solve(a.mat,b.mat) p0 = p01[1] p1 = p01[2] p2 = 1-p0-p1 if (min(p0,p1,p2)>0) { break } prod.gen.prov = prod.gen + runif(1,prod.gen-0.1,prod.gen+0.1) var.prod.prov = var.prod + runif(1,var.prod-0.1,var.prod+0.1) } if (min(couples)>0) { for (i in 1:dim(couples)[2]) { n.chicks = sample(c(0,1,2),p=c(p0,p1,p2),1) chicks[i] = n.chicks } chicks = chicks[chicks!=100] } else {chicks = 0} return(chicks) } ##### ##### Function for getting probability of survival of adult birds given CORT cort.survival = function(cort.value,exp.alfa = -2.399, exp.beta = 0.0803){ cort = exp(exp.alfa + exp.beta*cort.value)/(1+exp(exp.alfa + exp.beta*cort.value)) return(cort) } ##### survival of adult and prebreeders given g_max for m = males and f = females, maximum mortality # for pre-breeder m_max #cort.survival for adults is computed within the function #the function returns the position of dead birds kitti.survival.f = function (area.mat,gmax_m,gmax_f,m_max) { gmin = 7 m_min = 0.2 beta_f = -(log(m_max/m_min)/log(gmax_f/gmin)) beta_m = -(log(m_max/m_min)/log(gmax_m/gmin)) breeders = which(area.mat[breeding,]==1) ##breeders nonbreeders = which(area.mat[breeding,]==0 & area.mat[age,]!=100) ## non-breeders breeders.rand = runif(length(breeders)) ###random uniform numbers for bernoulli trial nonbreeders.rand = runif(length(nonbreeders)) ###random uniform numbers for bernoulli trial #print(as.numeric(area.mat[cort,breeders])) breeders.mort = cort.survival(as.numeric(area.mat[cort,breeders])) #print(breeders.mort) if (length(nonbreeders)>0) { non.breeders.mort = rep(0,length(nonbreeders)) for (i in 1:length(nonbreeders)){ if(as.numeric(area.mat[body_growth,nonbreeders[i]]) < 7) { non.breeders.mort[i] = 1 } else if (area.mat[sex,nonbreeders[i]] == "F" & as.numeric(area.mat[body_growth,nonbreeders[i]])gmax_f) { non.breeders.mort[i] = m_min } else if (area.mat[sex,nonbreeders[i]] == "M" & as.numeric(area.mat[body_growth,nonbreeders[i]])gmax_m) { non.breeders.mort[i] = m_min } } } dead =c((breeders[which(breeders.rand0) { num.dead = length(dead_partner) couples[dead_partner] = 2.5 pos.dead = which(couples==2.5,arr.ind=T) pos.dead.males = pos.dead[pos.dead[,1]==1,] ## matrix of row(1st col),col (2nd col) of position of dead males pos.dead.females = pos.dead[pos.dead[,1]==2,] # matrix of row(1st col),col (2nd col) of position of dead females fem.dead = length(which(pos.dead[,1]==2)) male.dead = num.dead - fem.dead } new.breeders.pot_col = which(area.mat[1,]!=100 & area.mat[breeding,]==1 & area.mat[attempts,]==0) ###columns with potential breeder, ie kitti alive and ready to breed new.breeders.pot = new.breeders.pot_col if (length(new.breeders.pot)>1) { males = new.breeders.pot[which(area.mat[sex,c(new.breeders.pot)]=="M")] #print("b") if (length(males)>0) { males = sample(males, length(males),replace=F) #reshuffle to make random the choice of males } else { males = 0 } #print("c") females = new.breeders.pot[which(area.mat[sex,new.breeders.pot]=="F")]#extract females from potential breeders if (length(females)>0) { females = sample(females, length(females),replace=F)#reshuffle to make random the choice of females } else { females = 0 } } else { males = 0 females=0 } #print(males) #print(females) if (length(dead_partner)>0 & min(males)!=0 ) { if (length(males)>male.dead) { males.subs = males[1:male.dead] males = males[(male.dead+1):length(males)] } else { males.subs = males males = 0 } couples[pos.dead.males[1:length(males.subs)]] = males.subs } if (length(dead_partner)>0 & min(females)!=0) { if (length(females)>fem.dead) { females.subs = females[1:fem.dead] females = females[(fem.dead+1):length(females)] } else { females.subs = females females = 0 } couples[pos.dead.females[1:length(females.subs)]] = females.subs } length.minsex = min(length(males),length(females)) if (length.minsex>0 & min(length(males),length(females))>0) { new.couples = rbind(males[1:length.minsex],females[1:length.minsex]) } if (exists("new.couples")==T) { #print(new.couples) } if (exists("new.couples")==T & length(couples)>1) { couples = cbind(couples,new.couples) if(length(couples)>1) { row1.del = c(which(couples[1,]==2.5),which(is.na(couples[1,])),which(couples[1,]==0)) row2.del = c(which(couples[2,]==2.5),which(is.na(couples[2,])),which(couples[1,]==0)) if (min(length(row1.del),length(row2.del))>0) { couples = couples[,-unique(c(row1.del,row2.del))] } #area.mat[attempts,couples] = as.numeric(area.mat[attempts,couples])+1 } } else if (exists("new.couples")==T & length(couples)<=1){ couples = new.couples if(length(couples)>1) { row1.del = c(which(couples[1,]==2.5),which(is.na(couples[1,])),which(couples[1,]==0)) row2.del = c(which(couples[2,]==2.5),which(is.na(couples[2,])),which(couples[1,]==0)) if (min(length(row1.del),length(row2.del))>0) { couples = couples[,-unique(c(row1.del,row2.del))] } } #area.mat[attempts,couples] = as.numeric(area.mat[attempts,couples])+1 } else if (exists("new.couples")==F & length(couples)<=1){ couples = 0 } return(couples) } ### couples = kitti.pairs(area,0) #### Function that returns the matrix of birds with the immigrants (immigrant occupy columns # in the matrix of birds, each column represents a single bird) kitti.immigration = function(area,max.pop.size,n.immigrants,cort.min,cort.max) { available_spots = max.pop.size - length(which(area[1,]!=100)) max.immigrants = min(available_spots,n.immigrants) #print(max.immigrants) if (max.immigrants > 0) { immigrant.pos = sample(which(area[1,]==100),max.immigrants,replace=F) #print(immigrant.pos) immigrant.length = length(immigrant.pos) #print(immigrant.length) immigrants.age = round(runif(immigrant.length,min=2,max=7)) immigrants.sex = sample(c("M","F"),immigrant.length,replace=T) immigrants.social=sample(c("A","B"),immigrant.length,replace=T) immigrants.quality = runif(immigrant.length,0,1) immigrants.body.growth = rnorm(immigrant.length,15,1.5) immigrants.breeding = rep(1,immigrant.length) immigrants.cort = runif(immigrant.length,cort.min,cort.max) immigrants.attempts = rep(0,immigrant.length) immigrants.vitals = rbind(immigrants.age,immigrants.sex,immigrants.social,immigrants.quality, immigrants.body.growth,immigrants.breeding,immigrants.cort,immigrants.attempts) area[,immigrant.pos] = immigrants.vitals } return(area) } ###prova.immi = kitti.immigration(area,50,10,5,20) ################################ #### Ind.Kitti is the model of population dynamics that pieces together the functions above #The function returns the list ris.list that you find at the bottom of the function Ind.Kitti <- function (max.pop.size=500,initial.pop.size=250,iter=50,corrCPUE=0, plotorno=1,scenario.sim = 1, mort_max = 0.7, immigration = 1, Imin = 1, Imax = 20, seed.sim = 2,av.prod.gen = 0.4) { initial.pop = sample(1:max.pop.size,initial.pop.size,replace=F) cort.min = 5 cort.max = 20 pairs_year = rep(2.5,iter) extinct = 0 prob_recruit_B = c(0.0474151,0.1243792,0.3960766,0.7257886,0.8836074, 0.9785819,0.9974608,0.9999925) prob_recruit_A = c(0.03373131,0.08955652,0.29969203,0.59906623, 0.78113076,0.93379513,0.98532027,0.99976028) n.immigrants = round(runif(iter,Imin,Imax)) #####assign row number to model variables#### age = 1 sex = 2 social = 3 quality = 4 body_growth = 5 breeding = 6 cort = 7 attempts = 8 ################################build intial population################ area = data.frame(rbind(rep(100,max.pop.size),matrix(0,7,max.pop.size))) area[age,initial.pop] = as.numeric(round(runif(initial.pop.size,1,15))) area[sex,initial.pop] = as.factor(sample(c("M","F"),initial.pop.size,replace=T)) area[social,initial.pop] = as.factor(sample(c("A","B"),initial.pop.size,replace=T)) area[quality,initial.pop] = as.numeric(runif(initial.pop.size,0,1)) area[body_growth,initial.pop] = as.numeric(runif(initial.pop.size,10,20)) area[breeding,initial.pop] = as.numeric(round(runif(initial.pop.size,0,1))) area[cort,initial.pop] = as.numeric(round(runif(initial.pop.size,0,20))) colnames(area) = (1:dim(area)[2]) couples = 0 log.growth = function(c1,c2,x) { growth.mean = (c1*(1-exp(-c2*x))) return(growth.mean) } #########################mean and sd of productivity (chicks fledged per nest) prod.gen.year.cons = rep(0,iter) var.prod.year.cons = rep(0,iter) growth_conditions.cons = rep(0,iter) chicks_recruited = rep(0,iter) gmax_male = 18 gmax_female = 17 startCPUE = 0 pCPUE = 0.85 rateCPUE = 3 minCPUE = 1 maxCPUE = 3.5 log10CPUE.ott = 2.25 set.seed(seed.sim) log10CPUE.year = log10CPUE.f(iter,minCPUE,maxCPUE,corrCPUE,startCPUE,pCPUE,rateCPUE) #prod.gen.year = kitti.prod.gen.f(iter,correlated,prod.min, prod.max , #start.ar , burn = 1000, alpha.ar,rate.ar) #var.prod.year = rep(0,iter) #for (vprod.cont in 1:length(prod.gen.year)) { #var.prod.year[vprod.cont] = runif(1,range.var(prod.gen.year[vprod.cont])[1]+0.01, #range.var(prod.gen.year[vprod.cont])[2]-0.01) #} ##########################start the cycle for iter year############### for(i in 1:iter) { #print(i) if(length(area[1,area[1,]!=100])<2){ print("extinct") ##is population extinct? #open for to check for extinct print(i) #prints the year of extinction extinct = 1 break } if (immigration == 1) { # if immigration = 1, runif(Imin,Imax) arrive from other colonies area = kitti.immigration(area,max.pop.size,n.immigrants[i],cort.min,cort.max) ###add immigrants } new.breeders.area = kitti.sexualmat.f(area,prob_recruit_B,prob_recruit_A) #new breeders given the age-specific prob of recruitment for B and A chicks if (min(new.breeders.area)!=0) { # if there are new breeders area[breeding,new.breeders.area] = 1 ## new breeders have breeding value = 1 } cort_year = kitti.cort.f(log10CPUE.year[i]) # get mean cort for the year given log10CPUE.year area = kitti.cort.ind.f(area,cort_year) #assign cort at each individuals ####get the productivty given the scenario and mean cort for the year prod.gen.year = prod.scenario(cortdata = cort_year, scenario = scenario.sim, min_high = 0.3, max_high = 1,min_med = 0.2,max_med = 0.6,min_low = 0.1,max_low = 0.2) prod.gen.year.cons[i] = prod.gen.year ####get the variance of productivty for the year var.prod.year = runif(1,range.var(prod.gen.year)[1]+0.01, range.var(prod.gen.year)[2]-0.01) var.prod.year.cons[i] = var.prod.year couples = kitti.pairs(area,couples) #form the new pairs for the year #print(couples) if (length(couples)<= 1) { #save the number of pairs in the colony for the year pairs_year[i] = 0 } else { pairs_year[i] = dim(couples)[2] } #print("before chicks") chicks.tot = kitti.prod.ind.f(prod.gen.year,var.prod.year,couples) #compute the total number of chicks fledging given mean prod, variance of prod #and couples #print(paste("chicks",sum(chicks.tot))) #print the number of chicks fledging if (sum(chicks.tot)>0) { # if there are chicks fledging chicks.A = length(chicks.tot[chicks.tot==1]) + length(chicks.tot[chicks.tot==2]) #A-chicks (including singletons) chicks.B = length(chicks.tot[chicks.tot==2]) #B chicks chicks.social = sample(c(rep("A",chicks.A),rep("B",chicks.B)),sum(chicks.tot),replace=F) #create the vector with the social status of chicks if (prod.gen.year > av.prod.gen) { # if the mean productivity is greater than av.prod.gen chicks.sex = sample(c("M","F"),p=c(0.55,0.45),sum(chicks.tot),replace=T) # we have more males } else { chicks.sex = sample(c("M","F"),p=c(0.45,0.55),sum(chicks.tot),replace=T) #otherwise more females } #print("C") chicks.AM = which(chicks.social=="A")[which(chicks.social=="A")%in%(which(chicks.sex=="M"))] # A-chicks males chicks.BM = which(chicks.social=="B")[which(chicks.social=="B")%in%(which(chicks.sex=="M"))] #B-chicks males chicks.AF = which(chicks.social=="A")[which(chicks.social=="A")%in%(which(chicks.sex=="F"))] #A-chicks females chicks.BF = which(chicks.social=="B")[which(chicks.social=="B")%in%(which(chicks.sex=="F"))] #B-chicks females #sort(c(chicks.AM,chicks.BM,chicks.AF,chicks.BF)) chicks.growth = rep(0,sum(chicks.tot)) growth_conditions = log10CPUE.year[i]/log10CPUE.ott growth_conditions.cons[i] = growth_conditions chicks.growth[chicks.AM] = rnorm(length(chicks.AM),log.growth(20,2,growth_conditions),1) chicks.growth[chicks.BM] = rnorm(length(chicks.BM),log.growth(20,1.7,growth_conditions),1) chicks.growth[chicks.AF] = rnorm(length(chicks.AF),log.growth(19,2,growth_conditions),1) chicks.growth[chicks.BF] = rnorm(length(chicks.BF),log.growth(19,1.7,growth_conditions),1) chicks.growth[chicks.growth<5] = 5 chicks.quality = runif(sum(chicks.tot),0,1) chicks.age = rep(0,sum(chicks.tot)) chicks.breeding = rep(0,sum(chicks.tot)) chicks.cort = rep(0,sum(chicks.tot)) chicks.attempts = rep(0,sum(chicks.tot)) if (max(as.numeric(area[1,])) == 100) { spots_for_juv = (which(area[1,]==100)) #print(paste("spots_for_juv",spots_for_juv)) #print(paste("sum(chicks.tot)",sum(chicks.tot))) #print(paste("length(spots_for_juv)",sum(length(spots_for_juv)))) chicks.rec = min(sum(chicks.tot),length(spots_for_juv)) chicks_recruited[i] = chicks.rec #print(paste("chicks.rec",sum(chicks.rec))) area[age,spots_for_juv[1:chicks.rec]] = chicks.age[1:chicks.rec] area[sex,spots_for_juv[1:chicks.rec]] = chicks.sex[1:chicks.rec] area[social,spots_for_juv[1:chicks.rec]] = chicks.social[1:chicks.rec] area[quality,spots_for_juv[1:chicks.rec]] = chicks.quality[1:chicks.rec] area[body_growth,spots_for_juv[1:chicks.rec]] = chicks.growth[1:chicks.rec] area[breeding,spots_for_juv[1:chicks.rec]] = chicks.breeding[1:chicks.rec] area[cort,spots_for_juv[1:chicks.rec]] = chicks.cort[1:chicks.rec] area[attempts,spots_for_juv[1:chicks.rec]] = chicks.attempts[1:chicks.rec] } } kitti.dead = kitti.survival.f(area,gmax_m = gmax_male,gmax_f = gmax_female, m_max = mort_max) if(length(kitti.dead)>0) { #print(paste("dead",length(kitti.dead))) area[1,kitti.dead] = 100 area[2:dim(area)[1],kitti.dead] = 0 } alive = which(area[age,]!=100) area[age,alive] = as.numeric(area[age,alive]) + 1 } pairs_year = pairs_year[pairs_year!=2.5] if (extinct==0) { median_pairs = median(pairs_year) min_pairs = min(pairs_year) max_pairs = max(pairs_year) range_pairs = max(pairs_year) - min(pairs_year) if (min(pairs_year) < 20) { quasi_ext = 1 } else { quasi_ext = 0 } } else { median_pairs = 0 min_pairs = 0 max_pairs = 0 range_pairs = 0 quasi_ext = 1 } if (plotorno >0) { par(mfrow=c(1,2)) plot(pairs_year,type="b",pch=16,xlab="Year",ylab="Number of breeding pairs",box=NULL, ylim=c(0,max(pairs_year)+10),cex.lab=1.5,bty="n",cex.axis=1.6) plot(prod.gen.year.cons,type="b",pch=16,xlab="Year",ylab="Mean productivity",box=NULL, xlim = c(0,length(pairs_year)),ylim=c(0,max(prod.gen.year.cons)+0.2),cex.lab=1.5,bty="n",cex.axis=1.6) } ris.list = list("extinct"=extinct,"area" = area, "prod.sd" = var.prod.year.cons, "prod.mean"= prod.gen.year.cons , "pairs" = pairs_year,"growth_cond" = growth_conditions.cons, "corr" = corrCPUE,"scenario" = scenario.sim, "mort_max" = mort_max, "immigration" = immigration, "minImm" = Imin, "maxImm" = Imax, "seed" = seed.sim,"median_pairs" = median_pairs,"min_pairs" = min_pairs, "max_pairs" = max_pairs, "range_pairs" = range_pairs, "quasi_ext" = quasi_ext,"CPUE" = log10CPUE.year,"recruit" = chicks_recruited ) return(ris.list) }# close all # Example of simulation # #ris = Ind.Kitti(max.pop.size=300,initial.pop.size=100,iter=50,corrCPUE=0, #plotorno=1,scenario.sim = 1, mort_max = 0.5, immigration = 1, # Imin = 1, Imax = 5, seed.sim = 2,av.prod.gen = 0.4) #### The results of the simulations are in the file outris.all.RDS #read the data outris.all <- readRDS("outris.all.RDS") # synth.df is the data frame with the results of the simulation. Each simulations gives a list of # results as output. Here below I extract the relevant info from each list. synth.df = data.frame(ext = sapply(outris.all, with, extinct),med_pairs = sapply(outris.all, with,median_pairs), min_pairs = sapply(outris.all, with,min_pairs), max_pairs = sapply(outris.all, with,max_pairs), range_pairs = sapply(outris.all, with,range_pairs),quasi_ext = sapply(outris.all, with, quasi_ext), max_imm = sapply(outris.all, with, maxImm),sc.corr = sapply(outris.all, with, corr), scenario = sapply(outris.all, with, scenario), mort_max = sapply(outris.all, with, mort_max),seed = sapply(outris.all, with, seed)) # ext = population went extinct (1) or not (0) # med_paris = median number of pairs # min_pairs = min number of pairs # max_pairs = max number of pairs # range_pairs = max(pairs_year) - min(pairs_year) # quasi_ext = 1 if number of pairs went below 20 pairs at any point during simulation time # max_imm = max number of immigrants each year (min is always 0) # sc.corr = 0 for scenario of food availability RF, 1 for AF and 2 for PF # scenario = high sensitivity (1), med sensitivity (2), min sensitivity (3) # mort_max = max mortality for pre-breeders # seed = seed for random number generator # load the package rms for the function lrm for logistic regression require(rms) #### extinction risk with the different scenarios of food availability (logistic regression) quasi_ext2.glm= with(subset(synth.df, sc.corr == 2), lrm(quasi_ext ~ as.factor(scenario) + mort_max + max_imm)) print(quasi_ext2.glm) quasi_ext1.glm= with(subset(synth.df, sc.corr == 1), lrm(quasi_ext ~ as.factor(scenario) + mort_max + max_imm + as.factor(seed))) print(quasi_ext1.glm) quasi_ext0.glm= with(subset(synth.df, sc.corr == 0), lrm(quasi_ext ~ as.factor(scenario) + mort_max + max_imm)) print(quasi_ext0.glm) ### linear models for median number of pairs for the various scenarios of food availability median_corr1.lm = with(subset(synth.df, sc.corr ==1),lm(log(med_pairs) ~ as.factor(scenario) + mort_max + max_imm + as.factor(seed))) ### Let's try a Poisson regression -- results are the same as with lm with response variable # log-transformed after rounding median_corr1.pois = with(subset(synth.df, sc.corr ==1),glm(round(med_pairs) ~ as.factor(scenario) + mort_max + max_imm + as.factor(seed),family = 'poisson')) ##### median_corr2.lm = with(subset(synth.df, sc.corr ==2),lm(log(med_pairs) ~ as.factor(scenario) + mort_max + max_imm + as.factor(seed))) median_corr0.lm = with(subset(synth.df, sc.corr==0),lm(log(med_pairs) ~ as.factor(scenario) + mort_max + max_imm + as.factor(seed))) #####semi-partial correlation coefficient require(rockchalk) sr.sc.corr_0 = getDeltaRsquare(median_corr0.lm) sr.sc.corr_1 = getDeltaRsquare(median_corr1.lm) sr.sc.corr_2 = getDeltaRsquare(median_corr2.lm) ### PLOTS # read the results for quasi-extinction quasi_ext.df = readRDS("quasi_ext.df.RDS") # read the results for median number of pairs med.df = read.table("med.txt",header = T) #### load the data from # Kitaysky AS, Piatt JF, Hatch SA, Kitaiskaia E V., # Benowitz-Fredericks ZM, Shultz MT, Wingfield JC (2010) # Food availability and population processes: severity of nutritional stress # during reproduction predicts survival of long-lived seabirds. Funct Ecol 24:625–637 kita2010.df = read.table("kita2010.txt",header = T) with(kita2010.df,plot(cort ~ cpue,xlab="log10(CPUE)",ylab="log10(CORT)",pch=16,xlim=c(0,4), ylim=c(0,1.8),cex = 1.5,cex.lab = 1.2,cex.axis = 1.1)) abline(with(kita2010.df,lm(cort ~ cpue),lwd = 2)) ############ Plot for quasi-extinction - different panel for each immigration scenario require(ggplot2) require(grid) labelh = expression(italic("Col")["H"]) labelm = expression(italic("Col")["M"]) labell = expression(italic("Col")["L"]) labelr = expression(italic("R")["F"]) labela = expression(italic("A")["F"]) labelp = expression(italic("P")["F"]) p1 = ggplot(subset(quasi_ext.df,imm==5), aes(x=sc.corr, y=sc.mean, group=scenario)) + # geom_line(size=0.7) + geom_point(aes(shape=scenario), # Shape depends on cond size = 6) + geom_line(aes(linetype=scenario), # Shape depends on cond size = 1.2) + scale_x_discrete(breaks=c("1","0","2"), labels=c(labela, labelr, labelp), name = "Temporal scenario") + scale_y_continuous(name="Quasi-extinction risk") + scale_shape_discrete(name="Colony", breaks=c("1", "2", "3"), labels=c(labelh,labelm,labell)) + theme( legend.position="none", #axis.title.x = element_text(size=20), #axis.text.x = element_text(size=16), axis.title.x = element_blank(), axis.text.x = element_text(size=16), axis.title.y = element_text(size=20), axis.text.y = element_text(size=16), axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank()) + geom_text(aes(3, .2, label="5"),size = 7) p2 = ggplot(subset(quasi_ext.df,imm==10), aes(x=sc.corr, y=sc.mean, group=scenario)) + #geom_line(size=0.7) + geom_point(aes(shape=scenario), # Shape depends on cond size = 6)+ geom_line(aes(linetype=scenario), # Shape depends on cond size = 1.2) + scale_x_discrete(breaks=c("1","0","2"), labels=c(labela, labelr, labelp), name = "Temporal scenario") + scale_y_continuous(name="Quasi-extinction risk") + scale_shape_discrete(name="Colony", breaks=c("1", "2", "3"), labels=c(labelh,labelm,labell)) + theme(legend.position="none", #axis.title.x = element_text(size=20), # axis.text.x = element_text(size=16), axis.title.y = element_blank(), axis.text.y = element_blank(), axis.title.x = element_text(size=20), axis.text.x = element_text(size=16), axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), plot.margin = unit(c(1,1,0.05,1), "cm"))+ geom_text(aes(3, .2, label="10"),size = 7) p3 = ggplot(subset(quasi_ext.df,imm==15), aes(x=sc.corr, y=sc.mean, group=scenario)) + geom_point(aes(shape=scenario), # Shape depends on cond size = 6)+ geom_line(aes(linetype=scenario), # Shape depends on cond size = 1.2) + scale_x_discrete(breaks=c("1","0","2"), labels=c(labela, labelr, labelp), name = "Temporal scenario") + scale_y_continuous(name="Quasi-extinction risk") + scale_shape_discrete(name="Colony", breaks=c("1", "2", "3"), labels=c(labelh,labelm,labell)) + scale_linetype_discrete(name = "Colony",breaks=c("1", "2", "3"), labels=c(labelh,labelm,labell)) + theme(legend.title = element_text(size=20, face="bold"), legend.text = element_text(size = 18), legend.justification=c(0,1), legend.position=c(0,1), axis.title.x = element_blank(), axis.text.x = element_text(size=16), axis.title.y = element_blank(), axis.text.y = element_blank(), axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), legend.key.size=unit(2.5,"lines")) + geom_text(aes(3, .2, label="15"),size = 7) #plot the multiplot multiplot(p1,p2,p3,cols=3) ############plot for median ######### labelh = expression(italic("Col")["H"]) labelm = expression(italic("Col")["M"]) labell = expression(italic("Col")["L"]) pd <- position_dodge(.1) p4 = ggplot(subset(med.df,imm==5), aes(x=as.factor(sc.corr), y=sc.mean, group=as.factor(scenario))) + geom_point(aes(shape=as.factor(scenario)),size = 6,position=pd) + # Shape depends on cond geom_line(aes(linetype=as.factor(scenario)), # Shape depends on cond size = 1.2) + geom_errorbar(aes(ymin=sc.mean-sc.sd, ymax=sc.mean+sc.sd), width=.1,position = pd) + scale_x_discrete(breaks=c("1","0","2"), labels=c(labela, labelr, labelp), name = "Temporal scenario") + scale_y_continuous(name="Median number of pairs",limits = c(0,100)) + scale_shape_discrete(name="Colony", breaks=c("1", "2", "3"), labels=c(labelh,labelm,labell)) + scale_linetype_discrete(name = "Colony",breaks=c("1", "2", "3"), labels=c(labelh,labelm,labell)) + theme( legend.position="none", #axis.title.x = element_text(size=20), #axis.text.x = element_text(size=16), axis.title.x = element_blank(), axis.text.x = element_text(size=16), axis.title.y = element_text(size=20), axis.text.y = element_text(size=16), axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank()) + geom_text(aes(3, .2, label="5"),size = 7) p5 = ggplot(subset(med.df,imm==10), aes(x=as.factor(sc.corr), y=sc.mean, group=scenario)) + #geom_line(size=0.7,position=pd) + geom_point(aes(shape=as.factor(scenario)), # Shape depends on cond size = 6,position=pd)+ geom_line(aes(linetype=as.factor(scenario)), # Shape depends on cond size = 1.2) + geom_errorbar(aes(ymin=sc.mean-sc.sd, ymax=sc.mean+sc.sd), width=.1,position = pd) + scale_x_discrete(breaks=c("1","0","2"), labels=c(labela, labelr, labelp), name = "Temporal scenario") + scale_y_continuous(name="Median number of pairs") + scale_shape_discrete(name="Colony", breaks=c("1", "2", "3"), labels=c(labelh,labelm,labell)) + theme(legend.position="none", #axis.title.x = element_text(size=20), # axis.text.x = element_text(size=16), axis.title.y = element_blank(), axis.text.y = element_blank(), axis.title.x = element_text(size=20), axis.text.x = element_text(size=16), axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), plot.margin = unit(c(1,1,0.05,1), "cm"))+ geom_text(aes(3, .2, label="10"),size = 7) p6 = ggplot(subset(med.df,imm==15), aes(x=as.factor(sc.corr), y=sc.mean, group=scenario)) + #geom_line(size=0.7,position=pd) + geom_point(aes(shape=as.factor(scenario)), # Shape depends on cond size = 6,position=pd)+ geom_line(aes(linetype=as.factor(scenario)), # Shape depends on cond size = 1.2) + geom_errorbar(aes(ymin=sc.mean-sc.sd, ymax=sc.mean+sc.sd), width=.1,position = pd) + scale_x_discrete(breaks=c("1","0","2"), labels=c(labela, labelr, labelp), name = "Temporal scenario") + scale_y_continuous(name="Median number of pairs") + scale_shape_discrete(name="Colony", breaks=c("1", "2", "3"), labels=c(labelh,labelm,labell)) + theme(legend.position="none", #axis.title.x = element_text(size=20), # axis.text.x = element_text(size=16), axis.title.y = element_blank(), axis.text.y = element_blank(), axis.title.x = element_text(size=20), axis.text.x = element_text(size=16), axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), legend.key.size=unit(2.2,"lines")) + geom_text(aes(3, .2, label="15"),size = 7) multiplot(p4,p5,p6,cols=3)