1 Motivation

Suite à la question d’Alain du 31 mai “Sur la viabilité des petites populatiosn, type de celle sur laquelle nous travaillons pour le lynx, avons nous des chiffres, des exemples quantifiés sur l’impact de 10 à 20 ind. tués chaque année par collision…Quel impact sur une population en terme démographique à 10, 20 et 50 ans…?”, j’ai fait l’exercice. Je me suis basé sur le modèle démographique de Jean-Michel Vandel et d’Eric Marboutin qu’ils ont utilisé dans leurs papiers :

2 Modèle démographique

On considère trois stades: chaton (kitten), subadulte et adulte. Le modèle est femelle-centrée, en post-breeding. Les primipares font 1 chaton, les multipares 2 chatons. Les paramètres de survie sont ceux utilisés par Vandel et al. (2006). Je prends une survie adulte à 0.85, au milieu de l’intervalle considéré par Vandel et al. (2006).

F1 <- 1 # primiparous
F2 <- 2 # multiparous
Sk <- 0.45 # kitten survival
Ssa <- 0.53 # subadult survival
Sa <- 0.85 # adult survival
A <- matrix(c(0,  0.5 * F1, 0.5 * F2,
             Sk, 0 ,       0,
             0,  Ssa,      Sa),
           byrow = T,nrow = 3)

On calcule le taux de croissance asymptotique.

out <- eigen(A)
Re(out$values[1])
## [1] 1.09491

Avec le modèle considéré, on a une croissance annuelle de 10% en gros, à l’équilibre.

La structure en classe stable est la suivante.

Re(out$vectors[,1])
## [1] 0.7142909 0.2935684 0.6353000
round(Re(out$vectors[,1])/sum(Re(out$vectors[,1]))*100,1) # on normalise
## [1] 43.5 17.9 38.7

Les effectifs à l’équilibre se répartissent en proportions approx 2/5 pour chatons et adultes, et 1/5 pour subadultes.

Les élasticités sont les suivantes.

vr <- list(F1=1,F2=2,Sk=0.45,Ssa=0.53,Sa=0.85)
elm <- expression( 
0,  0.5*F1, 0.5*F2,
Sk, 0,     0,
0,  Ssa,    Sa)
vitalsens(elm, vr)

Sans surprises, c’est à la survie adulte que le taux de croissance est le plus sensible.

3 Taux de croissance asymptotique et effet des collisions

Ici on étudie quel est l’effet des collisions sur le taux de croissance asymptotique. On part d’une population de 100 individus qui se répartissent en fonction de la structure stable. On applique une baisse des survies subadultes et adultes de 0% à 20% et on regarde l’effet sur le taux de croissance.

n <- 100
percent <- seq(0, 0.20, length = n)
growth_rate <- rep(NA, n*n)
adult_survival <- rep(NA, n*n)
subadult_survival <- rep(NA, n*n)
index <- 1
for (i in 1:n){
  for (j in 1:n){
    A <- matrix(c(0,  0.5 * F1, 0.5 * F2,
                  Sk, 0 ,       0,
                  0,  Ssa * (1 - percent[i]),      Sa * (1 - percent[j])),
                byrow = T,nrow = 3)
    growth_rate[index] <- Re(eigen(A)$values[1])
    adult_survival[index] <- Sa * (1 - percent[j])
    subadult_survival[index] <- Ssa * (1 - percent[i])
    index <- index + 1
  }
}
df <- data.frame(growth_rate = growth_rate,
                 adult_survival = adult_survival,
                 subadult_survival = subadult_survival)
df %>%
  ggplot() + 
  aes(x = adult_survival, y = subadult_survival, z = growth_rate) +
  geom_contour_fill() +
  scale_fill_divergent(midpoint = 1, guide = ) + 
  labs(x = "survie adulte", 
       y = "survie subadulte",
       fill = "taux de croissance asymptotique")

Clairement, c’est la baisse de survie sur les adultes plus que celle sur les subadultes qui entraine une baisse marquée sur le taux de croissance.

On a regardé dans cette section un effet sur les paramètres démographiques. Si on peut en théorie calculer le nombre d’individus tués par collisions correspondant à une diminution en pourcentage des paramètres de survie, ça ne répond pas exactement à la question d’Alain. Plutôt qu’un pourcentage de baisse dans les paramètres démographiques, on tue des individus directement par collision, et on regarde le taux de croissance observé, plus l’asymptotique.

4 Taux de croissance observé, modèle déterministe

D’abord un modèle déterministe.

F1 <- 1 # primiparous
F2 <- 2 # multiparous
Sk <- 0.45 # kitten survival
Ssa <- 0.53 # subadult survival
Sa <- 0.85 # adult survival
A <- matrix(c(0,  0.5 * F1, 0.5 * F2,
             Sk, 0 ,       0,
             0,  Ssa,      Sa),
           byrow = T,nrow = 3)
out <- eigen(A)
prop <- round(Re(out$vectors[,1])/sum(Re(out$vectors[,1])),1)

Scénario pas de mortalité par collision.

initial.n <- 100 * prop # vecteur d'effectifs initiaux
num.times <- 20 # le nombre de pas de temps de projection
N <- matrix(0,3,num.times + 1)
N[,1] <- initial.n
for(year in 2:(num.times+1)) N[,year] <- A %*% N[,year-1]

Scénario 4 morts par collision par an, avec 2 subadultes et 2 adultes.

N2 <- matrix(0,3,num.times + 1)
N2[,1] <- initial.n
for(year in 2:(num.times+1)) N2[,year] <- A %*% (N2[,year-1] - c(0, 2, 2))
N2[N2<0] <- 0

Scénario 10 morts par collision par an, avec 5 subadultes et 5 adultes.

N5 <- matrix(0,3,num.times + 1)
N5[,1] <- initial.n
for(year in 2:(num.times+1)) N5[,year] <- A %*% (N5[,year-1] - c(0, 5, 5))
N5[N5<0] <- 0

On visualise les tendances d’effectifs en fonction de ces 3 scénarios.

df <- data.frame(years = c(1:(num.times+1), 
                           1:(num.times+1), 
                           1:(num.times+1)),
                 counts = c(round(apply(N,2,sum)),
                            round(apply(N2,2,sum)),
                            round(apply(N5,2,sum))),
                 scenario = c(rep("aucune_collision", num.times+1),
                              rep("4_collisions", num.times+1),
                              rep("10_collisions", num.times+1)))
df %>%
  ggplot() + 
  aes(x = years, y = counts, color = scenario) + 
  geom_line() +
  scale_color_manual(labels = c("5 subadultes et 5 adultes tués par collision",
                               "2 subadultes et 2 adultes tués par collision",
                               "pas de collision"),
                     values = c("blue", "darkgreen", "orange")) +
  labs(x = "temps (années)", y = "effectif total de lynx", color = "")

En partant de 100 lynx, la population s’éteint avec 5 lynx subadultes et 5 adultes tués par an entre 10 et 15 ans.

5 Taux de croissance observé, modèle avec stochasticité démographique

Dans ce qui précède, on a négligé la stochasticité démographique. On corrige ça ici.

On redéfinit le modèle, puis on sépare fécondité et survie puisque la stochasticité est appliquée de façon différente sur ces deux types de paramètres (binomiale pour survie, log-normale pour fécondité).

F1 <- 1 # primiparous
F2 <- 2 # multiparous
Sk <- 0.45 # kitten survival
Ssa <- 0.53 # subadult survival
Sa <- 0.85 # adult survival
A <- matrix(c(0,  0.5 * F1, 0.5 * F2,
             Sk, 0 ,       0,
             0,  Ssa,      Sa),
           byrow = T,nrow = 3)
# separe fec de survie
sep_mat <- splitA(A = A)
sep_mat
## $T
##      [,1] [,2] [,3]
## [1,] 0.00 0.00 0.00
## [2,] 0.45 0.00 0.00
## [3,] 0.00 0.53 0.85
## 
## $F
##      [,1] [,2] [,3]
## [1,]    0  0.5    1
## [2,]    0  0.0    0
## [3,]    0  0.0    0
survie <- sep_mat$T 
fertil <- sep_mat$F 
#Une trajectoire. 
#multiresultm(n = initial.n, T = survie, F = fertil)

On écrit une fonction qui va simuler plusieurs trajectoires de la population en fonction d’un horizon et un nombre d’individus tués par collision.

stochdemo <- function(reps, tmax, n0, matS, matF, ncoll){ 
  # nreps = nb trajectoires
  # tmax = horizon
  # n0 = effectifs initiaux
  # matS = survie
  # matF = fertil
  # ncoll = nb collisions sur subadultes et adultes
  totalpop <- matrix(0, tmax, reps)
  for (j in 1:reps) { 
  n <- n0
  for (i in 1:tmax) { 
    n_priorcollision <- multiresultm(n, T = matS, F = matF) # forecast à t+1
    n <- n_priorcollision - matrix(c(0, ncoll, ncoll), ncol = 1) # retranche collisions
    n[n<0] <- 0 # si les effectifs sont négatifs, on les ramène à zéro
    totalpop[i,j] <- sum(n) # effectif total 
  }
}
return(totalpop)
}

On génère plusieurs trajectoires de la population (ntraj), à un horizon de plusieurs années (horizon). Puis on visualise.

ntraj <- 100 # 100 trajectoires
horizon <- 20 # 20 années
popsizeSD <- stochdemo(reps = ntraj,
                       tmax = horizon,
                       n0 = initial.n,
                       matS = survie,
                       matF = fertil,
                       ncoll = 0)
matplot(popsizeSD, 
        las=TRUE,
        log='y',
        type = 'l', 
        xlab = 'Temps (années)', 
        ylab = 'Effectif total de lynx')

Pas d’extinction. Qu’en est-il si 2 subadultes et 2 adultes sont tués chaque année par collision.

ntraj <- 100
horizon <- 20
popsizeSD <- stochdemo(reps = ntraj,
                       tmax = horizon,
                       n0 = initial.n,
                       matS = survie,
                       matF = fertil,
                       ncoll = 2)
matplot(popsizeSD, 
        las=TRUE,
        log='y',
        type = 'l', 
        xlab = 'Temps (années)', 
        ylab = 'Effectif total de lynx')

Pas vraiment de signe d’extinction là non plus. Et si ce sont 5 subadultes et 5 adultes qui sont tués par collision chaque année.

set.seed(1)
ntraj <- 100
horizon <- 20
popsizeSD <- stochdemo(reps = ntraj,
                       tmax = horizon,
                       n0 = initial.n,
                       matS = survie,
                       matF = fertil,
                       ncoll = 5)
matplot(popsizeSD, 
        las=TRUE,
        log='y',
        type = 'l', 
        xlab = 'Temps (années)', 
        ylab = 'Effectif total de lynx')

Plusieurs trajectoires conduisent à l’extinction. Calculons le taux d’extinction. On considère qu’il y a extinction si \(N < N_\text{seuil}\) de quasi-extinction. On fixe le seuil à 10 lynx.

Nseuil <- 10
extinction <- ifelse(popsizeSD < Nseuil, 1, 0)
popextinct <- ifelse(apply(extinction, 2, sum) > 1, 1, 0) # quelles trajectoires s'éteignent ?
proba.extinction <- sum(popextinct) / ntraj
proba.extinction
## [1] 0.84

La probabilité d’extinction est de 84%. On calcule le temps moyen à l’extinction.

tempsextinction <- NULL
for(n in (which(popextinct==1))){
  tempsextinction <- c(tempsextinction, min(which(extinction[,n]==1)))
}
# temps moyen à l'extinction pour les pop éteintes
mean(tempsextinction,na.rm=TRUE)
## [1] 14.08333

Le temps moyen à l’extinction est de 14 ans.

On reprend cette analyse en faisant varier les effectifs initiaux et le nombre de lynx tués par collision. On commence par créer une fonction qui va bien.

calcul_extinction <- function(ntraj = 100, 
                              horizon = 20, 
                              ncoll = 0, 
                              Nseuil = 10,
                              n = 100){
  F1 <- 1 # primiparous
  F2 <- 2 # multiparous
  Sk <- 0.45 # kitten survival
  Ssa <- 0.53 # subadult survival
  Sa <- 0.85 # adult survival
  A <- matrix(c(0,  0.5 * F1, 0.5 * F2,
                Sk, 0 ,       0,
                0,  Ssa,      Sa),
              byrow = T,nrow = 3)
  sep_mat <- splitA(A = A)
  sep_mat
  survie <- sep_mat$T 
  fertil <- sep_mat$F 
  out <- eigen(A)
  prop <- round(Re(out$vectors[,1])/sum(Re(out$vectors[,1])),1)
  initial.n <- n * prop
  popsizeSD <- stochdemo(reps = ntraj,
                         tmax = horizon,
                         n0 = initial.n,
                         matS = survie,
                         matF = fertil,
                         ncoll = ncoll)
  extinction <- ifelse(popsizeSD < Nseuil, 1, 0)
  popextinct <- ifelse(apply(extinction, 2, sum) > 1, 1, 0) 
  proba.extinction <- sum(popextinct) / ntraj
  tempsextinction <- 0
  if (proba.extinction > 0){
    tempsextinction <- NULL
    for(n in (which(popextinct==1))){
      tempsextinction <- c(tempsextinction, min(which(extinction[,n]==1)))
    }
  }
    res <- list(pr.ext = proba.extinction, 
                time.ext = mean(tempsextinction,na.rm = TRUE))
    
  return(res)
}

Puis on fait les calculs.

ncoll.grid <- seq(0, 10)
n.grid <- seq(10, 150, 10)
pr.ext <- rep(NA, length(ncoll.grid) * length(n.grid))
tps.ext <- rep(NA, length(ncoll.grid) * length(n.grid))
dead.lynx <- rep(NA, length(ncoll.grid) * length(n.grid))
total.lynx <- rep(NA, length(ncoll.grid) * length(n.grid))
index <- 1
for (i in 1:length(ncoll.grid)){
  for (j in 1:length(n.grid)){
    res <- calcul_extinction(ncoll = ncoll.grid[i], n = n.grid[j])
    pr.ext[index] <- res$pr.ext
    tps.ext[index] <- res$time.ext
    dead.lynx[index] <- ncoll.grid[i]
    total.lynx[index] <- n.grid[j]
    index <- index + 1
  }
}
df <- data.frame(pr.ext = pr.ext,
                 tps.ext = tps.ext,
                 dead.lynx = dead.lynx,
                 total.lynx = total.lynx)

df %>%
  ggplot() + 
  aes(x = dead.lynx, y = total.lynx, fill = pr.ext) +
  geom_raster() +
  scale_x_continuous(labels = dead.lynx,
                     breaks = dead.lynx) + 
  scale_y_continuous(labels = total.lynx,
                     breaks = total.lynx) + 
  scale_fill_viridis_c() + 
  labs(x = "nombre de lynx tués par collision (chaque année, subadulte et adulte)", 
       y = "effectif initial de la population",
       fill = "Pr(extinction à 20 ans)")

6 Avec une survie de 0.95

Jusqu’à présent, on a considéré une survie adulte médiane. On refait l’exercice avec une survie de 0.95 la borne max considérée par Vandel et al. (2006).

calcul_extinction <- function(ntraj = 100, 
                              horizon = 20, 
                              ncoll = 0, 
                              Nseuil = 10,
                              n = 100){
  F1 <- 1 # primiparous
  F2 <- 2 # multiparous
  Sk <- 0.45 # kitten survival
  Ssa <- 0.53 # subadult survival
  Sa <- 0.95 # adult survival
  A <- matrix(c(0,  0.5 * F1, 0.5 * F2,
                Sk, 0 ,       0,
                0,  Ssa,      Sa),
              byrow = T,nrow = 3)
  sep_mat <- splitA(A = A)
  sep_mat
  survie <- sep_mat$T 
  fertil <- sep_mat$F 
  out <- eigen(A)
  prop <- round(Re(out$vectors[,1])/sum(Re(out$vectors[,1])),1)
  initial.n <- n * prop
  popsizeSD <- stochdemo(reps = ntraj,
                         tmax = horizon,
                         n0 = initial.n,
                         matS = survie,
                         matF = fertil,
                         ncoll = ncoll)
  extinction <- ifelse(popsizeSD < Nseuil, 1, 0)
  popextinct <- ifelse(apply(extinction, 2, sum) > 1, 1, 0) 
  proba.extinction <- sum(popextinct) / ntraj
  tempsextinction <- 0
  if (proba.extinction > 0){
    tempsextinction <- NULL
    for(n in (which(popextinct==1))){
      tempsextinction <- c(tempsextinction, min(which(extinction[,n]==1)))
    }
  }
    res <- list(pr.ext = proba.extinction, 
                time.ext = mean(tempsextinction,na.rm = TRUE))
    
  return(res)
}

Puis on fait les calculs.

ncoll.grid <- seq(0, 10)
n.grid <- seq(10, 150, 10)
pr.ext <- rep(NA, length(ncoll.grid) * length(n.grid))
tps.ext <- rep(NA, length(ncoll.grid) * length(n.grid))
dead.lynx <- rep(NA, length(ncoll.grid) * length(n.grid))
total.lynx <- rep(NA, length(ncoll.grid) * length(n.grid))
index <- 1
for (i in 1:length(ncoll.grid)){
  for (j in 1:length(n.grid)){
    res <- calcul_extinction(ncoll = ncoll.grid[i], n = n.grid[j])
    pr.ext[index] <- res$pr.ext
    tps.ext[index] <- res$time.ext
    dead.lynx[index] <- ncoll.grid[i]
    total.lynx[index] <- n.grid[j]
    index <- index + 1
  }
}
df <- data.frame(pr.ext = pr.ext,
                 tps.ext = tps.ext,
                 dead.lynx = dead.lynx,
                 total.lynx = total.lynx)

df %>%
  ggplot() + 
  aes(x = dead.lynx, y = total.lynx, fill = pr.ext) +
  geom_raster() +
  scale_x_continuous(labels = dead.lynx,
                     breaks = dead.lynx) + 
  scale_y_continuous(labels = total.lynx,
                     breaks = total.lynx) + 
  scale_fill_viridis_c() + 
  labs(x = "nombre de lynx tués par collision (chaque année, subadulte et adulte)", 
       y = "effectif initial de la population",
       fill = "Pr(extinction à 20 ans)")

7 Discussion

Il s’agit d’un exercice fait à la va-vite, mais les éléments sont là. Ca fera un chouette complément au papier collisions. Dans les points à revoir, il y a certainement les suivants :