class: center, middle, title-slide # Population Projection Models (PPMs) ### Sarah Cubaynes for the Team ### last updated: 2022-03-18 --- # What we've learned so far ### Obj. 1: Assess current and past trends in population abundance - see class 1 -- ### Obj. 2: Estimate demographic parameters and identify causes of variation - see class 2 -- ### Obj. 3: Evaluate population viability to inform decision about management actions - some hints about this now! --- # Population Viability Analysis (PVA) Morris et al. (1999). *A practical handbook for population viability analysis*. The Nature Conservancy. + Use of **quantitative methods** to predict the **likely future status** of a population or collection of populations of conservation concern + **Tentative assessments** based upon what we now know **rather than as iron-clad predictions **of population fate --- # Why is PVA useful ? + Quantify rate of population change over time -- + Estimate extinction risks (used by IUCN) -- + Identify key parameters for population management -- + Evaluate and compare relative impact of population management actions --- ## Can we assess population viability from counts ? <img src="img/pop_size.png" width="60%" style="display: block; margin: auto;" /> --- class: center, middle # 1- Count-based PVA --- # Finite population growth rate `\(\lambda\)` .pull-left[ <img src="matrixmodels_files/figure-html/unnamed-chunk-2-1.svg" style="display: block; margin: auto;" /> ] .pull-right[ - Exponential growth or decay at a constant rate of change - `\(\lambda = \displaystyle{\frac{N_{t+1}}{N_{t}}}\)` gives the **proportional change** in population size - After `\(t\)` time steps, population size will be `\(N_{t} = N_{0} \cdot \lambda^{t}\)` - `\(\lambda\)` is log-normally distributed ] --- # Growth rate versus intrinsic rate of increase - **Population growth rate** `\(\lambda \sim \mbox{Lognormal}(\mu,\sigma^{2})\)` - Easier to work with the **intrinsic rate of increase** `$$r = \log(\lambda) = \log\left(\frac{N_{t+1}}{N_{t}}\right) = \log(N_{t+1}) - \log(N_{t})$$` -- - Which is normally distributed `\(r \sim N(\mu,\sigma^{2})\)` -- - `\(\mu\)` is the **mean rate of increase** - `\(\sigma^{2}\)` is the **environmental variance** --- # Environmental variance - If the mean rate of increase `\(\mu < 0\)`, extinction will certainly occur. -- - A population can still decline or go extinct even if the mean rate of increase `\(\mu>0\)`, because of environmental variance `\(\sigma^{2}\)` -- - Variable environments increase extinction risks --- ## Example of the Yellowstone grizzly bear population From Morris & Doak (2002). *Quantitative conservation biology: Theory and practice of population viability analysis*. Massachusetts, USA. .pull-left[ <img src="matrixmodels_files/figure-html/unnamed-chunk-3-1.svg" width="85%" /> ] .pull-right[ <img src="img/bear.jpg" width="60%" style="display: block; margin: auto;" /> ] --- # Step 1: Calculate `\(\mu\)` and `\(\sigma^{2}\)` from the data ```r # rate of increase over years logN <- log(N[-1]/N[-length(N)]) # log(Nt+1) - log(Nt) #mean rate of increase *mu <- mean(logN) #environmental variance *sigma2 <- var(logN) ``` ``` mu sigma2 1 0.02134027 0.01305092 ``` -- - `\(\mu >0\)` so on average the population is growing -- - `\(\sigma^{2} = 0.013\)` reflects low inter-annual variance in the rate of increase --- ## or using linear regression for unequal time intervals ```r x <- sqrt(grizzly$year[-1] - grizzly$year[-length(grizzly$year)]) # sqrt time intervals y <- logN / x mod <- lm(y ~ 0 + x) # forcing a intercept of zero mod ``` ``` Call: lm(formula = y ~ 0 + x) Coefficients: x 0.02134 ``` ```r mu <- coef(mod) # slope = mean intrinsic rate of increase ``` --- ## or using linear regression for unequal time intervals ```r # get an estimate for sigma2 anova(mod) ``` ``` Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) x 1 0.01731 0.017306 1.326 0.2569 Residuals 37 0.48288 0.013051 ``` ```r sigma2 <- anova(mod)[["Mean Sq"]][2] # environmental variance ``` --- ## Get confidence intervals for `\(\mu\)` and `\(\sigma_{2}\)` ```r ## Confidence interval for mu confint(mod, 1) ``` ``` 2.5 % 97.5 % x -0.01620969 0.05889023 ``` - Confidence interval of mean rate of increase encompasses 0, therefore we cannot rule out a potential risk of decline! -- ```r ## Confidence interval for sigma 2 df1 <- length(logN) - 1 df1 * sigma2 / qchisq(c(.975, .025), df = df1) ``` ``` [1] 0.008674359 0.021844393 ``` --- ## Back-transform to get finite population growth rate `\(\bar{\lambda}\)` ```r lambda <- exp(mu) lambda # average growth rate ``` ``` x 1.02157 ``` -- - Here `\(\bar{\lambda} > 1\)`, so the grizzly population is growing on average - It does not rule out the possibility of a decline owing to the chance occurrence of a sequence of bad years (remember confidence interval) --- # Step 2: Project the population - Expected population size using mean rate of increase and ignoring environmental variance (not recommended) `$$N_{t} = N_{0} \cdot \mu^{t}$$` `$$\ln(N_{t}) = \ln(N_{0}) + \mu \cdot t$$` -- * Initial population `\(N_{0}=44\)` * Time steps to project over: `\(t = 38\)` years (nb. years - 1) * Mean rate of increase `\(\mu = 0.021\)` --- # Step 2: Project the population <img src="matrixmodels_files/figure-html/unnamed-chunk-12-1.svg" style="display: block; margin: auto;" /> --- # Step 2: Project the population ## Let's account for observed variation in growth rate - First, set the initial population, number of time steps, and number of repetitions: ```r n0 <- grizzly$N[1] # initial pop. T <- 50 # time iterations to project over runs <- 500 # number of repetitions (pop. trajectories) stoch.pop <- matrix(NA,T,runs) # to store resuts stoch.pop[1,] <- n0 # initiate ``` --- # Step 2: Project the population ## Let's account for observed variation in growth rate - Then set a quasi-extinction threshold ```r Ne <- 30 # threshold for minimum viable pop. ``` -- - 1 female or a minimum viable population (genetic drift, demographic stochasticity) - Can also be the lowest level of abundance at which it remains feasible to attempt intervention to prevent further decline. --- ## Now run the projections ```r # let's project the population *for (i in 1:runs){ * for (t in 2:T){ # Draw r from normal using estimates of mu and sigma2 r <- rnorm(n = 1, mean = mu, sd = sqrt(sigma2)) # back-transform to get lambda and get pop. size lambda <- exp(r) #project one time step from the current pop size stoch.pop[t,i] <- stoch.pop[(t-1),i] * lambda # leave the loop if pop <= threshold if(stoch.pop[t,i] <= Ne){ stoch.pop[t,i] <- 0 i < i+1} } } ``` --- ## Now run the projections ```r # let's project the population for (i in 1:runs){ # loop over repetitions for (t in 2:T){ # loop over years # Draw r from normal using estimates of mu and sigma2 * r <- rnorm(n = 1, mean = mu, sd = sqrt(sigma2)) # back-transform to get lambda and get pop. size lambda <- exp(r) #project one time step from the current pop size stoch.pop[t,i] <- stoch.pop[(t-1),i] * lambda # leave the loop if pop <= threshold if(stoch.pop[t,i] <= Ne){ stoch.pop[t,i] <- 0 i < i+1} } } ``` --- ## Now run the projections ```r # let's project the population for (i in 1:runs){ # loop over repetitions for (t in 2:T){ # loop over years # Draw r from normal using estimates of mu and sigma2 r <- rnorm(n = 1, mean = mu, sd = sqrt(sigma2)) # back-transform to get lambda and get pop. size * lambda <- exp(r) #project one time step from the current pop size * stoch.pop[t,i] <- stoch.pop[(t-1),i] * lambda # leave the loop if pop <= threshold if(stoch.pop[t,i] <= Ne){ stoch.pop[t,i] <- 0 i < i+1} } } ``` --- ## Now run the projections ```r # let's project the population for (i in 1:runs){ # loop over repetitions for (t in 2:T){ # loop over years # Draw r from normal using estimates of mu and sigma2 r <- rnorm(n = 1, mean = mu, sd = sqrt(sigma2)) # back-transform to get lambda and get pop. size lambda <- exp(r) #project one time step from the current pop size stoch.pop[t,i] <- stoch.pop[(t-1),i] * lambda # leave the loop if pop <= threshold * if(stoch.pop[t,i] <= Ne){ * stoch.pop[t,i] <- 0 * i < i+1} } } ``` --- ## Step 3: Examine the results <img src="matrixmodels_files/figure-html/unnamed-chunk-19-1.svg" style="display: block; margin: auto;" /> --- ## Step 3: Examine the results Plot population size at the last time step: <img src="matrixmodels_files/figure-html/unnamed-chunk-20-1.svg" style="display: block; margin: auto;" /> --- ## Step 3: Examine the results Plot prediction together with observed counts: <img src="matrixmodels_files/figure-html/unnamed-chunk-21-1.svg" style="display: block; margin: auto;" /> --- # Step 4: Quantify extinction risks - The average population growth rate doesn’t do a good job at predicting what most population realizations will do -- - What are the chances that the population will go extinct at various times? + Extinction risk + Time to extinction --- # Step 4: Quantify extinction risks - Ultimate extinction probability = percentage of trajectories (over the 500 runs) reaching the extinction threshold at some point (over T years) ```r Pr.ext <- sum(lastN <= Ne) / runs # prob. to reach the extinction threshold Pr.ext ``` ``` [1] 0.22 ``` --- # Step 4: Quantify extinction risks - Cumulative extinction risk over the years <img src="matrixmodels_files/figure-html/unnamed-chunk-23-1.svg" style="display: block; margin: auto;" /> --- # Step 5: Quantify time to extinction - Mean time to extinction is an overestimation (because of few pop. growing fast) - Median time to extinction is a better measure ```r # Time to reach extinction for extinct pop. maxt <- NULL # empty vector to store results for (i in 1:runs){ # loop over repetitions N <- stoch.pop[,i] # max time N > threshold * maxt[i] <- max(which(N>0)) } # time at extinction for pseudo-extinct populations *time.ext <- maxt[maxt < T] median(time.ext) ``` ``` [1] 11 ``` --- # Step 5: Quantify time to extinction <img src="matrixmodels_files/figure-html/unnamed-chunk-25-1.svg" style="display: block; margin: auto;" /> --- # Step 6: Perturb and run the model ### Interesting to evaluate the sensitivity of the results to changes in: - Initial population size - Extinction threshold - Amount of environmental variance - Number of time steps --- # Live demo on grizzly bears From Morris & Doak (2002). *Quantitative conservation biology: Theory and practice of population viability analysis*. Massachusetts, USA. <img src="img/bear.jpg" width="40%" style="display: block; margin: auto;" /> --- # Count-based extinction analyses are based on strong hypotheses - Exhaustive counts (no sampling error) -- - No density-dependence (exponential growth) -- - Only source of variation is environmental stochasticity (no demographic stochasticity, no trends in mean or variance over time, uncorrelated environment among successive years) -- - Moderate environmental variability (no catastrophe, no bonanzas) --- # Advantages - Simplicity (data requirement at least 10 censuses and calculation) - Work relatively well when assumptions are met - Assess model quality by hindcasting --- # Limitations - `\(\lambda\)` is only a summary of the population dynamics -- - No info about the mechanisms governing the dynamics -- - No hints about which management action might be most efficient? - Is it better to act on survival? Fecundity? That of adults? Of juveniles? --- # Can we assess population viability from demographic parameters ? <img src="img/PPM.png" width="60%" style="display: block; margin: auto;" /> --- class: center, middle # 2- PVA using Matrix Projection Models (MPMs) --- # Let's assume no migrations for now - Survival and fecundity rates are enough to fully describe the population dynamics `$$N_{t+1} = N_{t} * F + N_{t} * S$$` + F = fecundity + S = survival (= 1- mortality) --- # Demographic parameters (see class 2) - Demographic parameters are often **heterogeneous across the life cycle**: survival and fecundity vary with age and/or stage in many species. <img src="img/oak-tree.jpg" width="60%" style="display: block; margin: auto;" /> --- # Matrix Population Models (MPMs) - Incorporate vital rates that are heterogeneous across the life cycle - Project the population based on the matrix summarizing the age- or stage- dependent demographic parameters --- # The Barn swallow (*Hirundo rustica*) example <img src="img/barnswallowlifecycle.png" width="60%" style="display: block; margin: auto;" /> --- # Demographic parameters Survival rates: - `\(S_{0}\)` chick survival - `\(S_{1}\)` juvenile survival (1 yo) - `\(S_{2}\)` adult survival (2+ yo) Fecundity: - `\(F_{1}\)` number of females produced by a juvenile female - `\(F_{2}\)` number of females produced by an adult female --- # Timing of data collection ? <img src="img/PreBreeding.png" width="100%" style="display: block; margin: auto;" /> + Pre-breeding census --- # Timing of data collection ? <img src="img/PostBreeding.png" width="100%" style="display: block; margin: auto;" /> + Post-breeding census --- # Step 1: write the agenda of events + Let's first consider a pre-breeding census <img src="img/AgendaPreB1.png" width="100%" style="display: block; margin: auto;" /> --- # Step 1: write the agenda of events + Let's first consider a pre-breeding census <img src="img/AgendaPreB2.png" width="100%" style="display: block; margin: auto;" /> --- # Step 1: write the agenda of events + Let's first consider a pre-breeding census <img src="img/AgendaPreB3.png" width="100%" style="display: block; margin: auto;" /> -- + Note that newborn are not observed directly! --- # Step 2: Translate into life cycle graph or 'Caswell representation' - A trick is to go 'up the arrows' <img src="img/preBLifeCycle.png" width="35%" style="display: block; margin: auto;" /> --- # Step 3: Translate into equations - Link `\(N_{(t+1)}\)` to `\(N_{(t)}\)` via survival and fertility rates - A trick is to read the parameters going up the arrows `\(N_{(1,t+1)} = F_{1} \cdot S_{0} \cdot N_{(1,t)}+ F_{2} \cdot S_{0} \cdot N_{(2,t)}\)` `\(N_{(2,t+1)} = S_{1} . N_{(1,t)} + S_{2} . N_{(2,t)}\)` --- # Step 4: Arrange in a matrix format + Called the **transition matrix**, or the **projection matrix** `$$N_{(t+1)} = \begin{bmatrix}F_{1} \cdot S_{0} & F_{2} \cdot S_{0}\\ S_{1} & S_{2} \end{bmatrix} \cdot N_{(t)}$$` --- # Step 4: Arrange in a matrix format + With `\(S_{0} = 0.2\)`, `\(S_{1} = 0.5\)` and `\(S_{2} = 0.65\)` + `\(F_{1} = 3/2\)` and `\(F_{2} = 6/2\)` `$$N_{(t+1)} = \begin{bmatrix}0.3 & 0.6\\ 0.5 & 0.65 \end{bmatrix} \cdot N_{(t)}$$` --- ## What is the difference with a post-breeeding census ? -- <img src="img/PostBreeding.png" width="100%" style="display: block; margin: auto;" /> + Post-breeding census --- # Step 1: write the agenda of events <img src="img/AgendaPostB1.png" width="100%" style="display: block; margin: auto;" /> --- # Step 1: write the agenda of events <img src="img/AgendaPostB2.png" width="100%" style="display: block; margin: auto;" /> --- # Step 1: write the agenda of events <img src="img/AgendaPostB3.png" width="100%" style="display: block; margin: auto;" /> -- + Note that newborn are now observed! --- # Step 2: Translate into life cycle graph or 'Caswell representation' <img src="img/postBLifeCycle.png" width="70%" style="display: block; margin: auto;" /> --- # Step 3: Translate into equations `\(N_{(0,t+1)} = S_{0} \cdot F_{1} \cdot N_{(0,t)} + S_{1} \cdot F_{2} \cdot N_{(1,t)}+ S_{2} \cdot F_{2} \cdot N_{(2,t)}\)` `\(N_{(1,t+1)} = S_{0} \cdot N_{(0,t)}\)` `\(N_{(2,t+1)} = S_{1} \cdot N_{(1,t)} + S_{2} \cdot N_{(2,t)}\)` -- # Step 4: Arrange in matrix format `$$N_{(t+1)} = \begin{bmatrix} S_{0} \cdot F_{1} & S_{1} \cdot F_{2} & S_{2} \cdot F _{2} \\ S_{0} & 0 & 0 \\ 0 & S_{1} & S_{2}\\ \end{bmatrix} \cdot N_{(t)}$$` --- # Why is the matrix format interesting ? - Easier to read than multiple equations - Intrinsic numeric features (back to it later) - Work the same way for complex life cycles --- # Examples <img src="img/Compadre.png" width="80%" style="display: block; margin: auto;" /> --- ## Several sites: Barn swallow example <img src="img/twosites.png" width="80%" style="display: block; margin: auto;" /> --- ## Variable age at first reproduction: Slender-billed gull example <img src="img/recruitment.png" width="80%" style="display: block; margin: auto;" /> --- ## Transition among stages: Peony example <img src="img/stages.png" width="80%" style="display: block; margin: auto;" /> --- # The projection matrix - How many age classes / stages in the life cycle? - Age at first reproduction? - Maximum age at death fixed or not? - Pre- or post-breeding census? <img src="img/Leslie.png" width="60%" style="display: block; margin: auto;" /> --- class: center, middle # 1- Deterministic MPMs --- # Step 5: Project the population + Write the transition matrix in R: ```r A.swallow <- matrix(c(0.3, 0.6, 0.5, 0.65), # pre-breeding leslie matrix nrow = 2, byrow = TRUE) A.swallow ``` ``` [,1] [,2] [1,] 0.3 0.60 [2,] 0.5 0.65 ``` --- # Step 5: Project the population + Start from an initial population `\(n_{0}\)` at time `\(t=0\)`: ```r n0 <- c(50,30) # vector with initial population n0 ``` ``` [1] 50 30 ``` --- # Step 5: Project the population + Project to the next time step: ```r n1 <- A.swallow %*% n0 # matrix product n1 ``` ``` [,1] [1,] 33.0 [2,] 44.5 ``` --- # Step 5: Project the population + Project the population over 10 years: ```r require(matrixcalc) t <- 10 # matrix product *n10 <- matrix.power(A.swallow,t) %*% n0 n10 ``` ``` [,1] [1,] 53.82434 [2,] 67.28043 ``` --- # Step 5: Project the population + Or using built-in functions from package `popbio`: ```r library(popbio) # load package t <- 11 # project the population *results <- pop.projection(A.swallow,n0,iterations = t) results ``` --- ## Step 6: Examine the results Let's plot the projection <img src="matrixmodels_files/figure-html/unnamed-chunk-51-1.svg" style="display: block; margin: auto;" /> --- ## Convergence to a stable distribution Transient dynamics: + Depends on the initial population + Damping ratio measures how fast the population converges toward equilibrium Stationary phase: + Independent of initial conditions + Depends on the transition matrix only + Constant growth rate = stable or **asymptotic population growth rate** ( `\(\lambda\)` ) + Constant proportion of individuals per age/stage = **stable age/stage structure** --- ## Step 6: Examine the results ```r results <- pop.projection(A.swallow, n0, iterations = t) names(results) ``` ``` [1] "lambda" "stable.stage" "stage.vectors" "pop.sizes" [5] "pop.changes" ``` --- ## Step 6: Examine the results First element contains the average growth rate in stationary phase = **asymptotic growth rate** ```r results$lambda ``` ``` [1] 1.05 ``` --- ## Step 6: Examine the results Second element contains the **stable age/stage structure** = % of each age/stage in the population in stationary phase ```r results$stable.stage ``` ``` [1] 0.4444444 0.5555556 ``` --- ## Step 6: Examine the results Third element contains pop. sizes per age/stage class at each time step ```r results$stage.vectors ``` ``` 0 1 2 3 4 5 6 7 8 [1,] 50 33.0 36.600 38.23500 40.16625 42.17261 44.28144 46.49549 48.82027 [2,] 30 44.5 45.425 47.82625 50.20456 52.71609 55.35177 58.11937 61.02533 9 10 [1,] 51.26128 53.82434 [2,] 64.07660 67.28043 ``` --- ## Step 6: Examine the results Fourth element contains total pop. size at each time step ```r results$pop.sizes ``` ``` [1] 80.00000 77.50000 82.02500 86.06125 90.37081 94.88870 99.63320 [8] 104.61486 109.84560 115.33788 121.10477 ``` --- ## Step 6: Examine the results Fifth element is `\(\lambda_{t}\)` the rate of change at each time step ```r results$pop.changes ``` ``` [1] 0.968750 1.058387 1.049208 1.050076 1.049993 1.050001 1.050000 1.050000 [9] 1.050000 1.050000 ``` --- ## Quantities in stationary phase Calculated directly from the transition matrix: ```r lambda(A.swallow) #stable population growth rate ``` ``` [1] 1.05 ``` --- ## Quantities in stationary phase Calculated directly from the transition matrix: ```r stable.stage(A.swallow) # stable age/stage structure ``` ``` [1] 0.4444444 0.5555556 ``` --- ## Quantities in stationary phase Calculated directly from the transition matrix: ```r # relative contribution of each age/stage to the next generation reproductive.value(A.swallow) ``` ``` [1] 1.0 1.5 ``` --- ## Quantities in stationary phase Calculated directly from the transition matrix: ```r generation.time(A.swallow) # average time between generations ``` ``` [1] 4.150924 ``` ```r # average age of mothers at birth of their daughters ``` --- # Step 7: Sensitivity analysis Let's perturb the model: What happens if female adult survival is reduced by 50%? ```r A.swallow.modified <- matrix(c(0.3, 0.6, 0.5, 0.65/2), # pre-breeding leslie matrix nrow = 2, byrow = TRUE) A.swallow.modified ``` ``` [,1] [,2] [1,] 0.3 0.600 [2,] 0.5 0.325 ``` --- # Step 7: Sensitivity analysis ```r n0 <- c(100,100) t <- 10 results.modified <- pop.projection(A.swallow.modified, n0, iterations = t) ``` --- # Step 7: Sensitivity analysis <img src="matrixmodels_files/figure-html/unnamed-chunk-64-1.svg" style="display: block; margin: auto;" /> --- # Step 7: Sensitivity analysis ```r lambda(A.swallow.modified) ``` ``` [1] 0.8603652 ``` + The population is now declining + Adult survival is important to population growth --- # Step 7: Sensitivity analysis + What happens if juvenile survival is reduced ? + What happens if fecundity or chick survival is reduced ? -- + Which demographic parameter contributes most to population dynamics ? --- # Step 7: Sensitivity analysis + Measuring the **impact of a change in a specific demographic parameter on population dynamics** + **Sensitivity** measures **absolute change** (e.g. -0.1 in parameter) `$$\frac{\delta \lambda}{\delta \theta}$$` -- + **Elasticity** measures **relative change** (e.g. -0.1% change in parameter) `$$\frac{\delta \lambda}{\delta \theta} \cdot \frac{\theta}{\lambda}$$` -- + Both are useful, elasticity is better to compare parameters that are on different scales (change in survival versus fertility) --- # Step 7: Sensitivity analysis ```r swallow.param <- list(s0 = 0.20, s1 = 0.5, s2 = 0.65, f1 = 3/2, f2 = 6/2) swallow.equation <- expression( s0 * f1, s0 * f2, s1, s2) VS <- vitalsens(swallow.equation, swallow.param) ``` --- # Step 7: Sensitivity analysis ```r VS <- vitalsens(swallow.equation, swallow.param) VS ``` ``` estimate sensitivity elasticity s0 0.20 1.82608696 0.34782609 s1 0.50 0.52173913 0.24844720 s2 0.65 0.65217391 0.40372671 f1 1.50 0.06956522 0.09937888 f2 3.00 0.08695652 0.24844720 ``` --- # Step 7: Sensitivity analysis <img src="matrixmodels_files/figure-html/unnamed-chunk-68-1.svg" style="display: block; margin: auto;" /> --- # Implications for management + Identify key parameters for population management - Best strategy here is to reduce adult mortality + Evaluate the impact of relative management actions - Actions focused on fecundity (e.g. nest protection) will have little impact --- # Live demo on the barn swallow <img src="img/swallow.jpg" width="328" style="display: block; margin: auto;" /> --- ## Assumptions and limitations of deterministic MPMs + One sex model: Females drive the demography, males are not limiting + Synchronous breeding: discrete time + No density-dependence (exponential growth or decay) -- + Demographic parameters are constant in time + No environmental stochasticity + No demographic stochasticity --- # Interest of deterministic models + Species living in stable environments (such as protected areas) + Exponential growth (recolonization, recovering from over-exploitation) -- + Sensitivity analyses provide useful information to identify key parameters for population management --- class: center, middle # 2- Stochastic MPMs --- # Environmental stochasticity - Unpredictable fluctuations in environmental conditions -- - Demographic rates vary over the years -- - Stochastic fluctuations in survival and reproductive rates reduce long-run population growth rate -- - The chance of occurrence of consecutive periods of unfavourable environmental conditions can drive small populations to extinctions -- - **MPMs in stochastic environments** project the population using a series of annual stochastic transition matrices (instead of one transition matrix) --- # The crested newt example High inter-annual variation in demographic parameters: - Mean adult survival = 0.52 varied between 0.22 et 0.74 in 8 years of monitoring - Mean adult fecundity = 3.07 juvenile females / adult female varied form 0.31 to 5.40 in 10 years - Bad years with dry pond in spring (about 1 out of 3 years) induces quasi-complete failure of reproduction --- # Crested newt (*Triturus cristatus*) life cycle and transtion matrix <img src="img/tritionlifecycle.png" width="2560" style="display: block; margin: auto;" /> --- ## Several options exist to include annual variability on demographic parameters - Random annual variation around mean values - Catastrophic events --- ## Several options exist to include annual variability on demographic parameters - **Random annual variation around mean values** - Catastrophic events --- ## Step 1: Start from the deterministic transition matrix ``` [,1] [,2] [,3] [1,] 0.00 1.03152 1.2894 [2,] 0.48 0.00000 0.0000 [3,] 0.00 0.52000 0.5200 ``` --- ## Step 2: Draw annual values for each demographic parameters from a normal distribution ```r s0stoch <- rnorm(n=1500,mean=s0,sd=0.15) s1stoch <- rnorm(n=1500,mean=s1,sd=0.15) s2stoch <- rnorm(n=1500,mean=s2,sd=0.15) s3stoch <- rnorm(n=1500,mean=s3,sd=0.15) fstoch <- rnorm(n=1500,mean=f,sd=0.5) ``` --- ## Step 2: Draw annual values for each demographic parameters from a normal distribution <img src="matrixmodels_files/figure-html/unnamed-chunk-73-1.svg" width="40%" style="display: block; margin: auto;" /> --- ## Step 3: Pile up the stochastic transition matrices ```r # Create a list of stochastic transition matrices A.newtSE <- list() # fill in by sampling from distrib. of demo. param. for(i in 1:ns){ A.newtSE[[i]] <- matrix( c( 0, * sample(alpha2,1) * sample(fstoch,1) * sample(s0stoch,1), * sample(alpha3,1) * sample(fstoch,1) * sample(s0stoch,1), * sample(s1stoch,1), 0, 0, * 0, sample(s2stoch,1), sample(s3stoch,1) ) , nrow = 3, ncol = 3) } ``` --- ## Step 4: Project the population using the stochastic transition matrices ```r T <- 30 *runSE <- stoch.projection(A.newtSE, n0 = c(50, 50, 50), tmax = T, nreps = 1000, verbose = FALSE) runSE ``` --- ## Step 4: Project the population using the stochastic transition matrices ```r runSE <- stoch.projection(A.newtSE, # initial population * n0 = c(50, 50, 50), tmax = T, nreps = 1000, verbose = FALSE) runSE ``` --- ## Step 4: Project the population using the stochastic transition matrices ```r runSE <- stoch.projection(A.newtSE, n0 = c(50, 50, 50), # number of time steps to project over * tmax = T, nreps = 1000, verbose = FALSE) runSE ``` --- ## Step 4: Project the population using the stochastic transition matrices ```r runSE <- stoch.projection(A.newtSE, n0 = c(50, 50, 50), tmax = T, # number of repetitions * nreps = 1000, verbose = FALSE) head(runSE) ``` ``` [,1] [,2] [,3] [1,] 59.83276 88.17996 91.30614 [2,] 92.12902 196.82366 274.23069 [3,] 616.31127 1282.68469 1301.33291 [4,] 71.46646 97.31048 223.07366 [5,] 999.37276 686.58327 978.99009 [6,] 294.61040 1126.87119 1255.92738 ``` --- ## Step 5: Examine the results Distribution of population sizes after `\(T=20\)` years? <img src="matrixmodels_files/figure-html/unnamed-chunk-79-1.svg" width="40%" style="display: block; margin: auto;" /> --- ## Step 5: Examine the results Long-run stochastic growth rate `\(\lambda_{s}\)`: ```r lambdastoch <- stoch.growth.rate(A.newtSE, maxt = 5000, verbose = FALSE) names(lambdastoch) ``` ``` [1] "approx" "sim" "sim.CI" ``` --- ## Step 5: Examine the results Long-run stochastic growth rate `\(\lambda_{s}\)` on a log-scale: ```r lambdastoch$approx # by Tuljapukar's approximation ``` ``` [1] 0.0365452 ``` ```r lambdastoch$sim # by simulation ``` ``` [1] 0.0335274 ``` ```r lambdastoch$sim.CI # with confidence interval ``` ``` [1] 0.02813619 0.03891860 ``` --- ## Step 5: Examine the results Long-run stochastic growth rate `\(\lambda_{s}\)`: ```r exp(lambdastoch$approx) # exponentiate to get stochastic growth rate ``` ``` [1] 1.037221 ``` --- ## Step 5: Examine the results Probability of extinction: ```r proba.ext <- stoch.quasi.ext(A.newtSE, # initial population size * n0 = c(50,50, 50), Nx = 30, nreps = 1000, tmax=50, maxruns = 10, verbose = FALSE) ``` --- ## Step 5: Examine the results Probability of extinction: ```r proba.ext <- stoch.quasi.ext(A.newtSE, # initial population size * n0 = c(50, 50, 50), Nx = 30, nreps = 1000, tmax = 50, maxruns = 10, verbose = FALSE) ``` --- ## Step 5: Examine the results Probability of extinction: ```r proba.ext <- stoch.quasi.ext(A.newtSE, n = c(50, 50, 50), # quasi-extinction threshold * Nx = 30, nreps = 1000, tmax = 50, maxruns = 10, verbose = FALSE) ``` --- ## Step 5: Examine the results Probability of extinction: ```r proba.ext <- stoch.quasi.ext(A.newtSE, n = c(50, 50, 50), Nx = 30, # number of runs * nreps = 1000, tmax = 50, maxruns = 10, verbose = FALSE) ``` --- ## Step 5: Examine the results Probability of extinction: ```r proba.ext <- stoch.quasi.ext(A.newtSE, n = c(50, 50, 50), Nx = 30, nreps = 1000, # number of time steps * tmax = 50, maxruns = 10, verbose = FALSE) ``` --- ## Step 5: Examine the results Probability of extinction: ```r proba.ext <- stoch.quasi.ext(A.newtSE, n = c(50, 50, 50), Nx = 30, nreps = 1000, tmax = 50, # repeat to get robust estimates * maxruns = 10, verbose = FALSE) ``` --- ## Step 5: Examine the results Probability of extinction: <img src="matrixmodels_files/figure-html/unnamed-chunk-89-1.svg" width="40%" style="display: block; margin: auto;" /> --- ## Step 5: Examine the results Probability of extinction: ```r proba.ext.mean <- apply(proba.ext,1,mean) proba.ext.mean[20] # in 20 years ``` ``` [1] 0.0016 ``` --- ## Several options to generate stochastic transition matrices Include annual variability on demographic parameters: - Random annual variation around mean values - **Catastrophic events** --- ## Step 1: Define transition matrix in good vs bad years ```r fgood <- 3.07 # fecundity in normal years fbad <- 0 # fecundity in years with Spring dryness of the pond ``` --- ## Step 1: Define transition matrix in good vs bad years ```r # Transition matrix in normal years A.newtGood <- matrix( c(0, s0 * alpha2 * fgood, s0 * alpha3 * fgood, s1, 0, 0, 0, s2, s3) , nrow = 3, ncol = 3, byrow = TRUE) # Transition matrix in bad years A.newtBad <- matrix( c(0, s0 * alpha2 * fbad, s0 * alpha3 * fbad, s1, 0, 0, 0, s2, s3) , nrow = 3, ncol = 3, byrow = TRUE) ``` --- ## Step 1: Define transition matrix in good vs bad years ``` $good [,1] [,2] [,3] [1,] 0.00 1.03152 1.2894 [2,] 0.48 0.00000 0.0000 [3,] 0.00 0.52000 0.5200 $bad [,1] [,2] [,3] [1,] 0.00 0.00 0.00 [2,] 0.48 0.00 0.00 [3,] 0.00 0.52 0.52 ``` --- ## Step 2: Define frequency of catastrophic events ```r # Spring dyness of the pond occurs every 3 years in average freqbad <- 1/3 ``` --- ## Step 3: Project the population Project using the 'good' or 'bad' transition matrix with probabillity defined above: ```r *stochCATA <- stoch.projection(A.newtCATA, prob = c( (1-freqbad), (freqbad)), n0 = c(50, 50, 50), tmax = 100, nreps = 1000, verbose = FALSE) head(stochCATA) ``` ``` [,1] [,2] [,3] [1,] 0.0006161299 0.0000000000 0.0002779838 [2,] 0.0020843707 0.0000000000 0.0009356701 [3,] 0.0029163790 0.0000000000 0.0011761417 [4,] 0.0000000000 0.0000000000 0.0218283887 [5,] 0.0008811366 0.0004034784 0.0003951277 [6,] 0.0089805546 0.0041055611 0.0040284151 ``` --- # Step 3: Project the population Project using the 'good' or 'bad' transition matrix with probabillity defined above: ```r stochCATA <- stoch.projection(A.newtCATA, # frequency of castastrophic events * prob = c( (1-freqbad), (freqbad)), n0 = c(50, 50, 50), tmax = 100, nreps = 1000, verbose = FALSE) stochCATA ``` --- # Step 3: Project the population Project using the 'good' or 'bad' transition matrix with probabillity defined above: ```r stochCATA <- stoch.projection(A.newtCATA, prob = c( (1-freqbad), (freqbad)), # initial population * n0 = c(50, 50, 50), tmax = 100, nreps = 1000, verbose = FALSE) stochCATA ``` --- # Step 3: Project the population Project using the 'good' or 'bad' transition matrix with probabillity defined above: ```r stochCATA <- stoch.projection(A.newtCATA, prob = c( (1-freqbad), (freqbad)), n0 = c(50, 50, 50), # number of time steps * tmax = 100, nreps = 1000, verbose = FALSE) stochCATA ``` --- # Step 3: Project the population Project using the 'good' or 'bad' transition matrix with probabillity defined above: ```r stochCATA <- stoch.projection(A.newtCATA, prob = c( (1-freqbad), (freqbad)), n0 = c(50, 50, 50), tmax = 100, # number of replicates * nreps = 1000, verbose = FALSE) stochCATA ``` --- # Step 4: Examine the results Long-run stochastic growth rate: ```r lambdaCATA <- stoch.growth.rate(A.newtCATA, prob = c( (1-freqbad), (freqbad)), maxt = 5000, verbose = FALSE) exp(lambdaCATA$approx) ``` ``` [1] 0.8871908 ``` --- # Step 4: Examine the results Probability of extinction: ```r proba.extCATA <- stoch.quasi.ext(A.newtCATA, prob = c( (1-freqbad), (freqbad)), n = c(50, 15, 50), # initial population Nx = 30, # Quasi-extinction threshold nreps = 1000, # nb. of replicates tmax = 50, # nb. of time steps maxruns = 10, # nb of repetitions verbose = TRUE) ``` ``` Calculating extinction probability for run 1 ``` ``` Calculating extinction probability for run 2 ``` ``` Calculating extinction probability for run 3 ``` ``` Calculating extinction probability for run 4 ``` ``` Calculating extinction probability for run 5 ``` ``` Calculating extinction probability for run 6 ``` ``` Calculating extinction probability for run 7 ``` ``` Calculating extinction probability for run 8 ``` ``` Calculating extinction probability for run 9 ``` ``` Calculating extinction probability for run 10 ``` ```r proba.extCATA ``` ``` [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 [2,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 [3,] 0.041 0.036 0.043 0.036 0.035 0.045 0.039 0.030 0.042 0.025 [4,] 0.117 0.118 0.121 0.101 0.109 0.104 0.112 0.100 0.116 0.105 [5,] 0.174 0.178 0.192 0.165 0.164 0.175 0.179 0.165 0.187 0.185 [6,] 0.225 0.239 0.261 0.222 0.223 0.243 0.251 0.236 0.252 0.256 [7,] 0.291 0.289 0.328 0.271 0.305 0.305 0.316 0.310 0.318 0.318 [8,] 0.341 0.342 0.379 0.328 0.363 0.356 0.358 0.374 0.373 0.382 [9,] 0.394 0.408 0.432 0.390 0.422 0.405 0.417 0.424 0.413 0.428 [10,] 0.435 0.471 0.467 0.440 0.470 0.458 0.464 0.467 0.475 0.486 [11,] 0.492 0.543 0.535 0.507 0.528 0.517 0.531 0.515 0.525 0.533 [12,] 0.553 0.607 0.584 0.563 0.582 0.565 0.577 0.562 0.586 0.593 [13,] 0.596 0.654 0.631 0.611 0.632 0.619 0.613 0.607 0.644 0.639 [14,] 0.638 0.692 0.666 0.649 0.674 0.671 0.650 0.652 0.692 0.678 [15,] 0.691 0.726 0.699 0.677 0.710 0.706 0.695 0.694 0.731 0.715 [16,] 0.709 0.755 0.728 0.707 0.736 0.735 0.730 0.721 0.762 0.745 [17,] 0.735 0.785 0.759 0.740 0.770 0.755 0.762 0.749 0.786 0.770 [18,] 0.766 0.812 0.786 0.771 0.794 0.782 0.780 0.774 0.805 0.803 [19,] 0.797 0.831 0.813 0.801 0.819 0.810 0.813 0.808 0.827 0.831 [20,] 0.823 0.852 0.834 0.831 0.840 0.834 0.835 0.822 0.843 0.857 [21,] 0.844 0.870 0.855 0.852 0.864 0.854 0.856 0.854 0.862 0.881 [22,] 0.865 0.890 0.880 0.874 0.879 0.872 0.879 0.875 0.878 0.894 [23,] 0.881 0.901 0.890 0.890 0.893 0.886 0.894 0.896 0.890 0.907 [24,] 0.891 0.912 0.900 0.908 0.904 0.901 0.908 0.909 0.908 0.920 [25,] 0.907 0.928 0.915 0.925 0.919 0.913 0.921 0.919 0.916 0.932 [26,] 0.921 0.937 0.921 0.940 0.925 0.923 0.930 0.929 0.929 0.940 [27,] 0.930 0.942 0.935 0.947 0.938 0.931 0.936 0.943 0.936 0.949 [28,] 0.938 0.950 0.942 0.950 0.946 0.942 0.943 0.949 0.941 0.956 [29,] 0.947 0.952 0.950 0.955 0.955 0.948 0.953 0.958 0.951 0.962 [30,] 0.953 0.963 0.956 0.960 0.960 0.955 0.956 0.964 0.954 0.968 [31,] 0.962 0.968 0.961 0.966 0.965 0.957 0.962 0.966 0.959 0.972 [32,] 0.966 0.970 0.965 0.968 0.973 0.961 0.968 0.971 0.965 0.974 [33,] 0.972 0.977 0.968 0.973 0.978 0.965 0.974 0.977 0.967 0.979 [34,] 0.972 0.980 0.976 0.976 0.978 0.968 0.977 0.978 0.972 0.980 [35,] 0.975 0.983 0.977 0.980 0.982 0.971 0.979 0.980 0.976 0.983 [36,] 0.976 0.985 0.981 0.985 0.984 0.975 0.979 0.984 0.980 0.985 [37,] 0.978 0.986 0.981 0.987 0.987 0.979 0.981 0.985 0.983 0.986 [38,] 0.983 0.988 0.982 0.987 0.991 0.984 0.985 0.987 0.984 0.986 [39,] 0.985 0.990 0.985 0.989 0.991 0.984 0.988 0.990 0.987 0.986 [40,] 0.987 0.992 0.985 0.991 0.992 0.984 0.988 0.990 0.988 0.988 [41,] 0.990 0.993 0.986 0.994 0.994 0.986 0.991 0.992 0.990 0.991 [42,] 0.991 0.993 0.988 0.995 0.994 0.987 0.993 0.994 0.992 0.991 [43,] 0.991 0.995 0.993 0.996 0.996 0.989 0.993 0.994 0.996 0.991 [44,] 0.993 0.996 0.993 0.996 0.996 0.990 0.993 0.995 0.997 0.993 [45,] 0.993 0.996 0.993 0.996 0.996 0.993 0.995 0.995 0.998 0.995 [46,] 0.994 0.997 0.994 0.997 0.996 0.993 0.996 0.995 0.998 0.996 [47,] 0.995 0.998 0.994 0.998 0.997 0.993 0.996 0.995 0.998 0.997 [48,] 0.995 0.998 0.996 0.998 0.998 0.994 0.996 0.996 0.998 0.997 [49,] 0.995 0.998 0.996 0.998 0.998 0.995 0.996 0.996 0.998 0.998 [50,] 0.995 0.998 0.997 0.998 0.998 0.995 0.996 0.996 0.998 0.999 ``` --- # Step 4: Examine the results Probability of extinction: <img src="matrixmodels_files/figure-html/unnamed-chunk-102-1.svg" width="40%" style="display: block; margin: auto;" /> --- # Step 4: Examine the results Probability of extinction: ```r proba.ext.mean <- apply(proba.extCATA,1,mean) proba.ext.mean[20] # in 20 years ``` ``` [1] 0.8371 ``` --- # Step 5: Interesting to evaluate the sensitivity of the results to changes in: - Amount of environmental variance - Frequency of catastrophic events -- - Initial population size - Extinction threshold - Number of time steps --- # Live demo on crested newts <img src="img/triton.jpg" width="50%" style="display: block; margin: auto;" /> --- # Assumptions and limitations - One sex model: Females drive the demography, males are not limiting - Synchronous breeding: discrete time - No density-dependence (exponential growth or decay) - No demographic stochasticity - No trend or change of mean in environmental variance --- # Assumptions and limitations - Keep a critical eyes on results: efficiency of PVAs debated in the litterature - We are making assumptions about demographic parameters in the future - Short-term studies = underestimation of variance in demographic rates while extinction risks increase with increased temporal variance in pop size --- # In which cases constant MPMs are useful despite its limitations? - Good knowledge of the species biology, life cycle and estimate of demographic parameters are critical - Sensitivity analyses remain a powerful tool to identify key demographic parameters and evalute management actions. - Especially for big populations (little impact of demographic stochasticity) - When survival and fertility are best structured by age or stage --- # Other (more complicated) models exist - MPMs including demographic stochasticity for small populations - Density-dependent MPMs - Two-sex models for species with skewed reproductive success - Integral projection models for size-structured demography - Continuous-time models - Multi-species models... --- # Useful references Morris, William, et al. (1999). *A practical handbook for population viability analysis.* The Nature Conservancy. Caswell, H. (2000). *Matrix population models (Vol. 1)*. Sunderland, MA: Sinauer. Brook, Barry W., et al. (2000). Predictive accuracy of population viability analysis in conservation biology. *Nature* 404.6776: 385-387. Beissinger, Steven R., and Dale R. McCullough, eds. (2002). *Population viability analysis*. University of Chicago Press. Reed, J. Michael, et al. (2002). Emerging issues in population viability analysis. *Conservation Biology* 16.1: 7-19.