# -------------------------------------------------------------------------------------------- # # - FILE NAME: IPMFunctions.R # - DATE: 30/07/2017 # - TITLE: Warming impacts on early life stages increase the vulnerability and delay the # population recovery of a long-lived habitat-forming macroalga # - AUTHORS: Pol Capdevila, Bernat Hereu, Roberto Salguero-Gómez, # Gracielˇla Rovira, Alba Medrano, Emma Cebrián, Joaquim Garrabou, # Diego K. Kersting, Cristina Linares # - SCRIPT: P. Capdevila (pcapdevila.pc@zoo.ox.ac.uk) # - JOURNAL: Journal of Ecology # -------------------------------------------------------------------------------------------- # # Script Content------------------------------------------------ # Description of the growth, survival and fecundity functions used # to build stochastich and density-dependent IPMs. # Additionally, it is provided a function to project the density- # dependent IPM. # Note that this code must be loaded in the CzIPM.R script ## Growth function, given you are size z now returns the pdf of size z1 next time g_z1z <- function(z, z1, m.par) { p.den.grow <- dnorm(z1, mean = m.par["grow.int"] + m.par["grow.z"] * z , sd = m.par["grow.sd"]) return(p.den.grow) } ## Survival function, logistic regression s_z <- function(z, m.par) { u <- exp(m.par["surv.int"] + m.par["surv.z"] * z ) return((u/(u+1))) } ## Reproduction function, logistic regression pb_z <- function(z, m.par) { u <- exp(m.par["flow.int"] + m.par["flow.z"] * z) return((u/(1+u))) } ## Recruits/adult pr_z <- function(Nt, mod=mod.sum){ sum.ad <- Nt mu.fec <- ifelse(sum.ad <= 0.1, 0,{ preds.sum <- data.frame(Adult=sum.ad, q=1, AdOffset=1) predict(mod, preds.sum, "response") }) return(mu.fec)} #Settlement probability load(paste(DataPath,"/sett.RData", sep = "")) pr_s <- function(Temp){ if(Temp==24|Temp=="24"){ settlement <- rtnorm(1, mean=mean(sett$sett), sd=sd(sett$sett), upper=max(sett$sett), lower=min(sett$sett)) } else{ settlement <- 1 } return(settlement) } #Recruit survival load(paste(DataPath,"/survrec.RData", sep = "")) r_s <- function(Temp){ if(Temp==24|Temp=="24"){ s <-rtnorm(1, mean=mean(survrec[[2]]), sd=sd(survrec[[2]]), upper=1, lower=min(survrec[[2]])) } else if(Temp==20|Temp=="20"){ s <- rtnorm(1, mean=mean(survrec[[1]]), sd=sd(survrec[[1]]), upper=1, lower=min(survrec[[1]])) } else s <- 1 return(s) } ## Recruit size function c_z1z <- function(z1, z, m.par) { recr.sizes <- dnorm(z1, m.par["c.mean"],m.par["c.sd"])/(1-pnorm(m.par["c.mean"],m.par["c.sd"])) return(recr.sizes) } #Kernel functions--------------------------------------------------------------------- ## Define the survival kernel P_z1z <- function (z1, z, m.par) { return( s_z(z, m.par) * g_z1z(z1, z, m.par) ) } # Define the reproduction kernel F_z1z <- function (z1, z, Nt, m.par, mod, Temp) { return( pb_z(z, m.par) * pr_z(Nt, mod) *pr_s(Temp)* r_s(Temp)*c_z1z(z1, z, m.par)) } ## Build the iteration matrix for the discretized kernel mk_K <- function(Nt, m=10, m.par, mod=mod.sum, Temp, L, U) { b=L+c(0:m)*(U-L)/m # boundary points y=0.5*(b[1:m]+b[2:(m+1)]) # mesh points h=y[2]-y[1] # step size #Create P and F matrices P <- h * (t(outer(y, y, P_z1z, m.par = m.par))) Fec <- h * (outer(y, y, F_z1z, Nt = Nt, m.par = m.par, mod=mod.sum, Temp=Temp)) K <- P+Fec return(list(K = K, meshpts = y, P = P, F = Fec, h=h)) } #Function to project the IPM--------------------------------------------------- dd.projection <- function(N, m, Temp, Time, mort=0, prob=0){ n <- matrix(ncol=Time, nrow=length(N)) #Empty matrix: colums are year rows the size of the kernel lambd <-numeric() #Store lambda IPM.list <- list() #Store IPM n[,1] <- N y<-mk_K(Nt=sum(N), m=m, m.par=m.par, Temp=Temp, L=L, U=U)$meshpts for(t in 2:Time){ u <- rbinom(1, 1, prob)# 1: disturbance; 0: no disturbance n[which(y>1), t-1] <- n[which(y>1), t-1]*(1 - u*mort)# if u = 0 or mort = 0, N[t-1] doesn't change Nt<- sum(n[which(y>1),t-1]) params.t=m.par + qnorm(runif(length(m.par),0.001,0.999))*m.par.sd.est #11 is the total number of parameters IPM<-mk_K(Nt=Nt, m=m, m.par=params.t, Temp=Temp, L=L, U=U)# Create dens-dep kernel IPM.list[[t-1]] <- IPM n[,t] <- IPM$K%*%n[, t-1] lambd[t-1] <- sum(n[,t])/sum(n[,t-1]) } return(list(n=n, lambda=lambd)) }