#################################################################################### # # # # # Code for the Integral Projection Model # # # # # #################################################################################### #Load Reproduction Data Rep=read.table("Rep.txt",h=T) #Load Survival Data Sur=read.table("Sur.txt",h=T) #CH include the individual capture histories CH=as.matrix(Sur[,c(6:19)]) CH[CH==0]=2 n.occasions=dim(CH)[2] #f is the the time of first capture f=Sur$f #SVL include the individual SVL at each capture SVL=as.matrix(Sur[,c(21:34)]) #To get the growth between seasons 1 and 2, we averaged observed SVL in each season SVL1=apply(SVL[,1:6],1,mean,na.rm=T) SVL2=apply(SVL[,7:14],1,mean,na.rm=T) datSVL=data.frame(SVL1,SVL2) datSVL[is.na(datSVL)]=NA ###create variables with indexes indicating the different seasons and sites #laps_season include the number of months between two sessions of captures for the three sites (in row). lapsseason=array(0,dim=c(3,13)) lapsseason[1,]=c(1,1.5,2,1,1,7,1,2,1,1,1,0.5,1) lapsseason[2,]=c(1,2.5,1,1,1,7,1,1,1,1,1,1,1) lapsseason[3,]=c(1,2.5,0.5,1,1,6,0.5,1,2,1,1,1,1) ####Index of the season for each capture session. seas=c(rep(0,6),rep(1,7)) ####there are only 6 sessions of capture in the sites 1 and 2 while 7 in site 3. seap=array(1,dim=c(3,13)) seap[1:2,13]=0 ##Indexes for SITE effect for survival and reproduction SITE=site23=as.numeric(Sur$plot)-1 SITE[SITE==0]=1 site23[site23==2]=1 SITER=siter23=Rep$plot-1 SITER[SITER==0]=1 siter23[siter23==2]=1 #function to create an initial matrix for known true states: 1: ALIVE, NA: unknwon. known.state.ms <- function(ms, notseen){ # notseen: label for 'not seen' state <- ms state[state==notseen] <- NA for (j in 1:dim(ms)[1]){ deb=min(which(ms[j,]<2)) fin=max(which(ms[j,]<2)) for (k in which(is.na(state[j,deb:fin]))){ state[j,deb+k-1]=1} } return(state) } zdata=known.state.ms(CH,2) library(jagsUI) #model in BUGS language sink("Surgrorep.jags") cat(" model { #### Priors ##survival alpha~ dnorm(0, 0.01) #intercept svl~ dnorm(0, 0.01) #effect of SVL alphasea~ dnorm(0, 0.01) #effect of season on intercept svlsea~ dnorm(0, 0.01) #interaction between SVL and season ##growth #mean K~ dnorm(0, 0.01) #growth rate between season 1 and 2 L~ dnorm(0, 0.01) #asymptotic size #variance salpha~ dnorm(0, 0.01) #intercept sb~ dnorm(0, 0.01) #effect of SVL ##reproduction alpha_rep~ dnorm(0, 0.01) #intercept b_rep~ dnorm(0, 0.01) #effect of SVL alpha_repsea~ dnorm(0, 0.01) #effect of season on intercept b_repsea~ dnorm(0, 0.01) #interaction between SVL and season tau_rep <- pow(sigma_rep, -2)#precision of the normal distribution sigma_rep ~ dunif(0, 5) #standard deviation of the normal distribution ##recapture palpha~ dnorm(0, 0.01) #intercept psvl~ dnorm(0, 0.01) #effect of SVL palphasea~ dnorm(0, 0.01) #effect of season on intercept psvlsea~ dnorm(0, 0.01) #interaction between SVL and season #site effect for the different functions for(u in 1:2){ alphasite[u]~ dnorm(0, 0.01) #on the intercept of survival Ksite[u]~ dnorm(0, 0.01) #on growth rate Lsite[u]~ dnorm(0, 0.01) #on asymptotic size salpha_site[u]~ dnorm(0, 0.01) #on the interept of the variance of growth sb_site[u]~ dnorm(0, 0.01) alpha_repsite[u]~ dnorm(0, 0.01) #on the intercept of reproduction palphasite[u]~ dnorm(0, 0.01) #on the intercept of recapture } ###Models for each vital rate for (i in 1:nind){ for (t in f[i]:(n.occasions-1)){ #survival log(phi[i,t])<- -exp(alpha+alphasea*seas[t]+svl*SVL[i,(seas[t]+1)]+svlsea*seas[t]*SVL[i,(seas[t]+1)]+alphasite[SITE[i]]*site23[i]) #recapture logit(p[i,t]) <- palpha+palphasea*seas[t]+psvl*SVL[i,(seas[t]+1)]+psvlsea*seas[t]*SVL[i,(seas[t]+1)]+palphasite[SITE[i]]*site23[i] } #t } #i #Growth for (i in 1:650){ #mean mu[i] <- exp((-K)-Ksite[SITE[i]]*site23[i])*(SVL[i,1]+(L+Lsite[SITE[i]]*site23[i])*(1/exp((-K)-Ksite[SITE[i]]*site23[i])-1)) #variance sigma[i] <-exp(salpha+salpha_site[SITE[i]]*site23[i]+sb*SVL[i,1]+sb_site[SITE[i]]*SVL[i,1]*site23[i]) tau[i] <- pow(sigma[i], -2) } #reproduction for (i in 1:ndataR){ mu_rep[i] <- alpha_rep+alpha_repsea*sear[i]+b_rep*svlr[i]+b_repsea*sear[i]*svlr[i]+alpha_repsite[SITER[i]]*siter23[i] } ###Likelihood # survival and recapture # Define state-transition and observation matrices for (i in 1:nind){ # Define probabilities of survival for (t in f[i]:(n.occasions-1)){ ps[1,i,t,1] <- phi[i,t]^sea[Site[i],t] ps[1,i,t,2] <- 1-phi[i,t]^sea[Site[i],t] ps[2,i,t,1] <- 0 ps[2,i,t,2] <- 1 # Define probabilities of recapture po[1,i,t,1] <- p[i,t]*seap[Site[i],t] po[1,i,t,2] <- 1-p[i,t]*seap[Site[i],t] po[2,i,t,1] <- 0 po[2,i,t,2] <- 1 } #t } #i for (i in 1:nind){ for (t in (f[i]+1):n.occasions){ # State process: survival z[i,t] ~ dcat(ps[z[i,t-1], i, t-1,]) # Observation process: recapture y[i,t] ~ dcat(po[z[i,t], i, t-1,]) } #t } #i #Growth for (i in 1:650){ SVL[i,2] ~ dnorm(mu[i],tau[i]) } #Reproduction for (i in 1:ndataR){ egg[i] ~ dnorm(mu_rep[i],tau_rep) } } ",fill = TRUE) sink() # Bundle data jags.data <- list(y = CH,z=zdata,Site=as.numeric(Sur$plot), f = f, nind = dim(CH)[1], n.occasions = dim(CH)[2],seas=seas, SVL=datSVL,SITE=SITE,site23=site23,sea=lapsseason,seap=seap, egg=Rep$egg,svlr=Rep$svl,sear=Rep$season-1,SITER=SITER,siter23=siter23,ndataR=length(Rep$egg)) # Initial values inits <- function(){list(alpha=0, svl=0,palpha=0,psvl=0, palphasea=0,alphasea=0, svlsea=0,psvlsea=0, K=0,L=0,salpha=0,sb=0, alpha_rep=0,b_rep=0,alpha_repsea=0,b_repsea=0, palphasite=rep(0,2),alphasite=rep(0,2),alpha_repsite=rep(0,2), Ksite=rep(0,2),Lsite=rep(0,2),salpha_site=rep(0,2),sb_site=rep(0,2) )} # Parameters monitored parameters <- c('alpha','svl','palpha','psvl',' palphasea','alphasea','svlsea','psvlsea', 'K','L','salpha','sb', 'alpha_rep','b_rep','alpha_repsea','b_repsea', 'palphasite','alphasite','alpha_repsite', 'Ksite','Lsite','salpha_site','sb_site') #run the model par<- jags(jags.data, inits, parameters, "Surgrorep.jags", n.chains = 3, n.thin = 50, n.iter = 100000, n.burnin =50000, parallel=T) print(par, digits = 3) ##### IPM part that runs with the parameters estimated in the previous model niter=(100000-50000)/50 nn=100 #Initiate population growth rate (lambda) lambda=array(0,dim=c(2,niter)) #Initiate stable age distribution (sad) and reproductive value (rv) sad=rv=array(0,dim=c(2,niter,nn)) #Initiate mean and SD of svl svl.stats=array(0,dim=c(2,niter,2)) #SVL parameters to build the meshpoints: maxsize<-65; L<-0; U<-1.1*maxsize # boundary points b and mesh points y b<-L+c(0:nn)*(U-L)/nn meshpoints<-0.5*(b[1:nn]+b[2:(nn+1)]) #parameters from the litterature inhparam=c(25,2)#Inheritance mean and variance. eggsur=0.9 #egg survival juvsur=0.574 #first-year survival #To get the uncertainty of estimated population growth rate and stable age distribution, # we used the posterior distributions of the parameters for survival , growth and reproduction estimated from the previous model. # Thus, we build an IPM for each iteration of the previous model. for (j in 1:niter){ #Store the survival, growth and reproduction parameters for each season (1 and 2) and each iteration i sparam1=c(par$sims.list$alpha0[j],par$sims.list$bm0[j]) sparam2=c(par$sims.list$alpha0[j]+par$sims.list$alphasea[j],par$sims.list$bm0[j]+par$sims.list$bmsea[j]) gparam1=c(par$sims.list$K0[j],par$sims.list$L0[j],par$sims.list$salpha_bm0[j],par$sims.list$sb_bm0[j]) rparam1=c(par$sims.list$alpha_rep0[j],par$sims.list$b_rep0[j]) rparam2=c(par$sims.list$alpha_rep0[j]+par$sims.list$alpha_repsea[j],par$sims.list$b_rep0[j]+par$sims.list$b_repsea[j]) #The following line source the functions and matrices of the IPM. source('IPM_functions.R') #For season 1 #Build the big matrix of the IPM: mat <- I%*%(R1/2) + G1%*%S1 #Estimate population growth rate, stable age distribution and reproductive value goat=get.eigen.stuff(mat) lambda[1,j] <-goat[[1]] sad[1,j,] <- as.vector(goat[[2]])/sum(as.vector(goat[[2]])) # work out mean SVL svl.stats[1,j,] <- means.and.vars(meshpoints,sad[1,j,]) #For season 2 #Build the big matrix of the IPM: mat <- I%*%(R2/2) +G1%*%S2 goat=get.eigen.stuff(mat) #Estimate population growth rate, stable age distribution and reproductive value lambda[2,j] <-goat[[1]] sad[2,j,] <- as.vector(goat[[2]])/sum(as.vector(goat[[2]])) # work out mean SVL for season one and two svl.stats[2,j,] <- means.and.vars(meshpoints,sad[2,j,]) } Res=list(lambda,sad,svl.stats)