rm(list=ls(all=TRUE)) library(ggplot2) library(dplyr) # set working directory (example of my directory): setwd("/Users/maria/Dropbox/Marmots/dryad") # LOAD MODEL PARAMETERS load("Bay_mod_param.rda") ### DEFINE VITAL RATE FUNCTIONS mcmc=out1J$sims.list # FUNCTIONS TO PREDICT DEMOGRAPHIC PROCESSES # Inputs are: # mass (z) # mass next (zz) # life-cycle stage # year (here repeating the chronological sequences, but in projections can be sampled from a Normal distribution with mean = 0 and sd = 'SDunabs' parameters from respective submodels) # Q.sim (latent variable) # post (posterior parameter sample) # WINTER SURVIVAL: AUGUST TO JUNE NEXT S1.fun <- function(z,stage,year,Q.sim,post) { mu.surv=exp(mcmc[["a0.surv1"]][post]+ mcmc[["bc.surv1"]][post]*z + mcmc[["a1.surv1"]][post,stage]+ mcmc[["aY.surv1"]][post,year]+ mcmc[["bq.surv1"]][post]*Q.sim) return(mu.surv/(1+mu.surv)) } # SUMMER SURVIVAL: JUNE TO AUGUST S2.fun <- function(z,stage,year,Q.sim,post) { mu.surv=exp(mcmc[["a0.surv2"]][post]+ mcmc[["bc.surv2"]][post]*z + mcmc[["a1.surv2"]][post,stage]+ mcmc[["aY.surv2"]][post,year]+ mcmc[["bq.surv2"]][post]*Q.sim) return(mu.surv/(1+mu.surv)) } # WINTER TRAIT TRANSITION AUGUST TO JUNE (we assume a contant variance) GR1.fun <- function(z,zz,stage,year,Q.sim,post){ growth.mu=(mcmc[["a0.gr"]][post]+ (mcmc[["bc.gr"]][post]+mcmc[["bcstage.gr"]][post,stage])*z + mcmc[["a1.gr"]][post,stage]+ mcmc[["aY.gr"]][post,year]+ mcmc[["bq.gr"]][post]*Q.sim) # Get residual variance var.res= (mcmc[["sigma.gr"]][post])^2 # Density distribution function of the normal distribution gr1 = sqrt(2*pi*var.res) gr2 = ((zz-growth.mu)^2)/(2*var.res) return(exp(-gr2)/gr1) } # SUMMER TRAIT TRANSITION JUNE TO AUGUST (we assume a contant variance) GR2.fun <- function(z,zz,stage,year,Q.sim,post){ growth.mu=(mcmc[["a0.gr2"]][post]+ (mcmc[["bc.gr2"]][post]+mcmc[["bcstage.gr2"]][post,stage])*z + mcmc[["a1.gr2"]][post,stage]+ mcmc[["aY.gr2"]][post,year]+ mcmc[["bq.gr2"]][post]*Q.sim) # Get residual variance var.res= (mcmc[["sigma.gr2"]][post])^2 # Density distribution function of the normal distribution gr1 = sqrt(2*pi*var.res) gr2 = ((zz-growth.mu)^2)/(2*var.res) return(exp(-gr2)/gr1) } # PROBABILITY OF REPRODUCING: PR.fun <- function(z,stage,year,Q.sim,post) { mu.rep=exp(mcmc[["a0.pr"]][post]+ mcmc[["bc.pr"]][post]*z+ mcmc[["a1.pr"]][post,stage]+mcmc[["aY.pr"]][post,year]+ mcmc[["bq.pr"]][post]*Q.sim) return(mu.rep/(1+mu.rep)) } # NUMBER OF RECRUITS: R.fun <- function(z,year,Q.sim,post) { mu.rec=exp(mcmc[["a0.rec"]][post]+ mcmc[["bc.rec"]][post]*z+ mcmc[["aY.rec"]][post,year]+ mcmc[["bq.rec"]][post]*Q.sim) return(mu.rec) } ## OFFSPRING MASS OffMass.fun <- function(z,zz,year,Q.sim,post){ growth.mu=(mcmc[["a0.off"]][post]+ mcmc[["bc.off"]][post]*z+ mcmc[["aY.off"]][post,year]+ mcmc[["bq.off"]][post]*Q.sim) # Get residual variance var.res= (mcmc[["sigma.off"]][post])^2 # Density distribution function of the normal distribution gr1 = sqrt(2*pi*var.res) gr2 = ((zz-growth.mu)^2)/(2*var.res) return(exp(-gr2)/gr1) } ############################################################################################# ############# !!!!! SIMULATIONS CAN BE LOADED (line 378) !!!!!!!! # IPMS ### Necassary IPM input parameters # Define the lower and upper integration limit L=7.791458 # minimum size U=17.07776 n=100 # bins b <- L+c(0:n)*(U-L)/n # interval that each cell of the discretized IPM z <- 0.5*(b[1:n]+b[2:(n+1)]) # midpoint h=(U-L)/n # bin width # Simulation parameters year.sim=40 # 40 years of simulations n.stage=4 # number of stages par.sub=sample(1:3000,100, replace=F)# sample posterior parameters (note that here we only sample 100 to decrease the simulation time) # Arrays to hold results in pop.growth.simJ=array(NA,c(length(par.sub),year.sim)) pop.growth.simY=array(NA,c(length(par.sub),year.sim)) pop.growth.simR=array(NA,c(length(par.sub),year.sim)) pop.growth.simN=array(NA,c(length(par.sub),year.sim)) ### INITIAL POPULATION VECTOR (Number of individuals in a given stage and mass class) load("N1.rda") N1.orig=N1 for(x in 1:length(par.sub)){ # loop through posterior samples N1=N1.orig for(i in 1:year.sim){ # loop through years IPMja=array(0,c(n*n.stage,n*n.stage)) # summer IPM IPMaj=array(0,c(n*n.stage,n*n.stage)) # winter IPM ### AUGUST - JUNE # NOTE: in projections, both the random year effect and the latent Q can be sampled from distributions as described in the main text year = i # here we just replicate the chronoligal order of years Q.sim=mcmc$Q[par.sub[x],i] # we also pick Q from derived year-specific values Sj <- diag(S1.fun(z,stage=1,year,Q.sim,par.sub[x])) # Survival juveniles Sy <- diag(S1.fun(z,stage=2,year,Q.sim,par.sub[x])) # Survival yearlings Snra <- diag(S1.fun(z,stage=3,year,Q.sim,par.sub[x])) # Survival non-reproductive adults Sra <- diag(S1.fun(z,stage=4,year,Q.sim,par.sub[x])) # Survival reproductive adults # Transition To RA or NRA Ty <- diag(PR.fun(z,stage=1,year,Q.sim,par.sub[x])) Tnra <- diag(PR.fun(z,stage=2,year,Q.sim,par.sub[x])) Tra <- diag(PR.fun(z,stage=3,year,Q.sim,par.sub[x])) # Trait transition (growth) - stage specific like for survival G <- h*t(outer(z,z,GR1.fun,stage=1,year,Q.sim,par.sub[x])) # Control for eviction: # this is equivalent to redistributing evicted sizes evenly among existing size classes G=G/matrix(as.vector(apply(G,2,sum)),nrow=n,ncol=n,byrow=TRUE) G[is.na(G)]=0 # FILL IPM IPMaj[(n+1):(2*n),1:n]=G%*%Sj # Juvenile to Yearling G <- h*t(outer(z,z,GR1.fun,stage=2,year,Q.sim,par.sub[x])) # Control for eviction: G=G/matrix(as.vector(apply(G,2,sum)),nrow=n,ncol=n,byrow=TRUE) G[is.na(G)]=0 # FILL IPM IPMaj[(2*n+1):(3*n),(n+1):(2*n)]=G%*%(Sy*diag(1-PR.fun(z,stage=1,year,Q.sim,par.sub[x]))) # Yearling to Non-Reproductive Adult IPMaj[(3*n+1):(4*n),(n+1):(2*n)]=G%*%(Sy*Ty) # Yearling to Reproductive Adult G <- h*t(outer(z,z,GR1.fun,stage=3,year,Q.sim,par.sub[x])) G=G/matrix(as.vector(apply(G,2,sum)),nrow=n,ncol=n,byrow=TRUE) G[is.na(G)]=0 # FILL IPM IPMaj[(2*n+1):(3*n),(2*n+1):(3*n)]=G%*%(Snra*diag(1-PR.fun(z,stage=2,year,Q.sim,par.sub[x]))) # Non-Reproductive Adult to Non-Reproductive Adult IPMaj[(3*n+1):(4*n),(2*n+1):(3*n)]=G%*%(Snra*Tnra) # Non-Reproductive Adult to Reproductive Adult G <- h*t(outer(z,z,GR1.fun,stage=4,year,Q.sim,par.sub[x])) G=G/matrix(as.vector(apply(G,2,sum)),nrow=n,ncol=n,byrow=TRUE) G[is.na(G)]=0 # FILL IPM IPMaj[(2*n+1):(3*n),(3*n+1):(4*n)]=G%*%(Sra*diag(1-PR.fun(z,stage=3,year,Q.sim,par.sub[x]))) # Reproductive Adult to non-Reproductive Adult IPMaj[(3*n+1):(4*n),(3*n+1):(4*n)]=G%*%(Sra*Tra) # Reproductive Adult to Reproductive Adult ######################################################## ### JUNE - AUGUST # SURVIVAL S2y <- diag(S2.fun(z,stage=1,year,Q.sim,par.sub[x]),nrow=n,ncol=n) # Survival yearlings S2nra <- diag(S2.fun(z,stage=2,year,Q.sim,par.sub[x]),nrow=n,ncol=n) # Survival non-reproductive adults S2ra <- diag(S2.fun(z,stage=3,year,Q.sim,par.sub[x]),nrow=n,ncol=n) # Survival reproductive adults # TRAIT TRANSITIONS # Yearlings Gy <- h*t(outer(z,z,GR2.fun,stage=1,year,Q.sim,par.sub[x])) # Reproductive Adults Gnra <- h*t(outer(z,z,GR2.fun,stage=2,year,Q.sim,par.sub[x])) # Non-Reproductive Adults Gra <- h*t(outer(z,z,GR2.fun,stage=3,year,Q.sim,par.sub[x])) # Control for eviction: # this is equivalent to redistributing evictd sizes evenly among existing size classes Gy=Gy/matrix(as.vector(apply(Gy,2,sum)),nrow=n,ncol=n,byrow=TRUE) Gy[is.na(Gy)]=0 Gra=Gra/matrix(as.vector(apply(Gra,2,sum)),nrow=n,ncol=n,byrow=TRUE) Gra[is.na(Gra)]=0 Gnra=Gnra/matrix(as.vector(apply(Gnra,2,sum)),nrow=n,ncol=n,byrow=TRUE) Gnra[is.na(Gnra)]=0 # Offspring mass OffMass <- h*t(outer(z,z,OffMass.fun,year,Q.sim,par.sub[x])) # Control for eviction: # this is equivalent to redistributing evictd sizes evenly among existing size classes OffMass=OffMass/matrix(as.vector(apply(OffMass,2,sum)),nrow=n,ncol=n,byrow=TRUE) OffMass[is.na(OffMass)]=0 # recruitment (as individuals that die in August still produce offspring, recruitment does not need to be multiplied by S2) Rec=(diag(R.fun(z,year,Q.sim,par.sub[x])))/2 Fkernel <- as.matrix(OffMass%*%(Rec*S2ra)) # FILL IPM IPMja[(n+1):(2*n),(n+1):(2*n)]=Gy%*%S2y # Yearling stay Yearling IPMja[(2*n+1):(3*n),(2*n+1):(3*n)]=Gnra%*%S2nra # Non-Reproductive Adults stay Non-Reproductive Adults IPMja[(3*n+1):(4*n),(3*n+1):(4*n)]=Gra%*%S2ra # Reproductive Adults stay Reproductive Adults IPMja[1:n,(3*n+1):(4*n)]=Fkernel # Adults producing juveniles #### Observed population changes N0=N1 N05=IPMaj%*%N0 N1=IPMja%*%N05 pop.growth.simJ[x,i]=sum(N1[1:100])/sum(N0[1:100]) pop.growth.simY[x,i]=sum(N1[101:200])/sum(N0[101:200]) pop.growth.simN[x,i]=sum(N1[201:300])/sum(N0[201:300]) pop.growth.simR[x,i]=sum(N1[301:400])/sum(N0[301:400]) } } # Put results in one data frame to plot pop.growth.sim.ALL=NULL pop.growth.simJ.mu=data.frame(mean=apply(pop.growth.simJ,2,mean), LB=apply(pop.growth.simJ,2,quantile,.025), UB=apply(pop.growth.simJ,2,quantile,.975), year=1976:2015, cat="J") pop.growth.sim.ALL=rbind(pop.growth.sim.ALL,pop.growth.simJ.mu) pop.growth.simY.mu=data.frame(mean=apply(pop.growth.simY,2,mean), LB=apply(pop.growth.simY,2,quantile,.025), UB=apply(pop.growth.simY,2,quantile,.975), year=1976:2015, cat="Y") pop.growth.sim.ALL=rbind(pop.growth.sim.ALL,pop.growth.simY.mu) pop.growth.simN.mu=data.frame(mean=apply(pop.growth.simN,2,mean), LB=apply(pop.growth.simN,2,quantile,.025), UB=apply(pop.growth.simN,2,quantile,.975), year=1976:2015, cat="N") pop.growth.sim.ALL=rbind(pop.growth.sim.ALL,pop.growth.simN.mu) pop.growth.simR.mu=data.frame(mean=apply(pop.growth.simR,2,mean), LB=apply(pop.growth.simR,2,quantile,.025), UB=apply(pop.growth.simR,2,quantile,.975), year=1976:2015, cat="R") pop.growth.sim.ALL=rbind(pop.growth.sim.ALL,pop.growth.simR.mu) ### Get observed changes in numbers load("pop.growth.obs.rda") pop.growth.sim.ALL$sim="Simulated" pop.growth.obs$sim="Observed" pop.growth=rbind(pop.growth.obs,pop.growth.sim.ALL[,c(4,1,5,2,3,6)]) #### PLOT pop.growth[pop.growth$sim=="Simulated"&pop.growth$cat=="R",c(2,4,5)][1,]=NA ggplot(data=pop.growth,aes(x=year,y=mean,col=sim,fill=sim))+ facet_wrap( ~ cat, ncol=2,scales="free")+ geom_ribbon(aes(ymin=LB,ymax=UB),alpha=0.2,col=NA)+ geom_line(size=1)+ scale_color_manual(name="",values=c("black","darkgreen"))+ scale_fill_manual(name="",values=c("black","darkgreen"))+ theme_bw()+ theme(legend.title = element_text(size=16, face="bold"), legend.text = element_text(size=14), legend.key.size = unit(1, "lines"), legend.background =element_rect(fill = "transparent"), legend.position = c(.1,0.9))+ guides(color=F) + xlab("") + ylab("Change in number")+ theme(panel.grid = element_blank())+ theme(axis.text = element_text(size=14,colour="black"))+ theme(axis.title = element_text(size=18))+ theme(axis.title.x = element_text(vjust=-0.4), axis.title.y = element_text(vjust=0.9), axis.text.x = element_text(angle = 65, hjust = 1), plot.title = element_text(size=12,lineheight=.8, face="bold"))+ theme(plot.margin = unit(c(0.2,0.1,0.5,0.1), "cm"))+ scale_x_continuous(breaks = seq(1975,2015,by=5))+ theme(strip.text = element_text(size=12), strip.background =element_blank()) #### Plot points colnames(pop.growth.sim.ALL)[1]="mean.sim" pop.growth.obs[,c("mean.sim","LB","UB")]=left_join(pop.growth.obs,pop.growth.sim.ALL,by=c("year","cat"))[,c("mean.sim","LB.y","UB.y")] pd <- position_dodge(0.1) pop.growth.obs=pop.growth.obs[is.finite(pop.growth.obs$mean.sim),] ggplot(data=pop.growth.obs,aes(x=mean,y=mean.sim))+ facet_wrap( ~ cat, ncol=2,scales="free")+ geom_abline(intercept = 0, slope = 1, color="grey", size=1.2)+ geom_errorbar(aes(ymin=LB, ymax=UB), width=.1, position=pd) + geom_point(size=3,position=pd)+ theme_bw()+ theme(legend.title = element_text(size=16, face="bold"), legend.text = element_text(size=14), legend.key.size = unit(1, "lines"), legend.background =element_rect(fill = "transparent"), legend.position = c(.1,0.9))+ xlab("Observed change in number") + ylab("Simulated change in number")+ theme(panel.grid = element_blank())+ theme(axis.text = element_text(size=14,colour="black"))+ theme(axis.title = element_text(size=18))+ theme(axis.title.x = element_text(vjust=-0.4), axis.title.y = element_text(vjust=0.9), plot.title = element_text(size=12,lineheight=.8, face="bold"))+ theme(plot.margin = unit(c(0.2,0.1,0.5,0.1), "cm"))+ theme(strip.text = element_text(size=12), strip.background =element_blank())