R Markdown

Calcul de l’Abondance avec N-mixture

library(readxl)
library(unmarked)

Abondance 2026

all2026 <- matrix(c(1,0,0,0,0,2,1,1,3,2,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0, 2,1,0,2,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,1,0,0,
  0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), nrow = 25,byrow = TRUE)

Terriers = c(1,0,0,1,0,1,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0)
  
Veg = c(1,0,0,1,1,1,1,1,1,1,1,1,0,0,1,0,1,0,1,0,1,1,1,0,1)
  
GA = c(0,1,0,1,0,1,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0)

BM = c(1,0,0,1,0,0,0,1,0,0,1,1,0,0,1,1,1,1,1,0,1,1,1,0,1)
## 1. Modèle de Royle avec p VARIABLE (p=proba de détection) ----

nSites <-25
nVisits <-5
visitMat <- matrix(as.character(1:nVisits), nSites, nVisits, byrow=TRUE)
umf1 <- unmarkedFramePCount(y=all2026, siteCovs=NULL, obsCovs=list(visit=visitMat))
## Warning: obsCovs contains characters. Converting them to factors.
# Fit a model
fm1 <- pcount(~visit-1 ~ 1, umf1, K=10)
fm1
## 
## Call:
## pcount(formula = ~visit - 1 ~ 1, data = umf1, K = 10)
## 
## Abundance (log-scale):
##  Estimate    SE     z P(>|z|)
##    -0.872 0.354 -2.46  0.0137
## 
## Detection (logit-scale):
##        Estimate    SE      z P(>|z|)
## visit1   0.2976 0.744  0.400  0.6892
## visit2  -1.4417 0.814 -1.770  0.0767
## visit3  -1.4417 0.814 -1.770  0.0767
## visit4  -0.0872 0.701 -0.124  0.9010
## visit5  -0.9104 0.725 -1.256  0.2092
## 
## AIC: 95.53172 
## Number of sites: 25
plogis(coef(fm1, type="det"))
## p(visit1) p(visit2) p(visit3) p(visit4) p(visit5) 
## 0.5738500 0.1912823 0.1912823 0.4782089 0.2869202
# Empirical Bayes estimation of random effects
(fm1re <- ranef(fm1))
##             Mean Mode 2.5% 97.5%
##  [1,] 1.04337166    1    1     2
##  [2,] 3.27689440    3    3     4
##  [3,] 0.04337166    0    0     1
##  [4,] 0.04337166    0    0     1
##  [5,] 0.04337166    0    0     1
##  [6,] 2.18129887    2    2     3
##  [7,] 0.04337166    0    0     1
##  [8,] 0.04337166    0    0     1
##  [9,] 1.04337166    1    1     2
## [10,] 0.04337166    0    0     1
## [11,] 1.04337166    1    1     2
## [12,] 0.04337166    0    0     1
## [13,] 0.04337166    0    0     1
## [14,] 0.04337166    0    0     1
## [15,] 1.04337166    1    1     2
## [16,] 0.04337166    0    0     1
## [17,] 0.04337166    0    0     1
## [18,] 0.04337166    0    0     1
## [19,] 0.04337166    0    0     1
## [20,] 0.04337166    0    0     1
## [21,] 0.04337166    0    0     1
## [22,] 0.04337166    0    0     1
## [23,] 0.04337166    0    0     1
## [24,] 0.04337166    0    0     1
## [25,] 0.04337166    0    0     1
plot(fm1re, subset=site %in% 1:8, xlim=c(-1,10))

sum(bup(fm1re)) # Estimated population
## [1] 10.45574
colSums(confint(fm1re)) # and 95% CI
##  2.5% 97.5% 
##     9    34
## Modèle de Royle avec p CONSTANTE ----

umf2 <- unmarkedFramePCount(y=all2026, siteCovs=NULL, obsCovs=NULL)
summary(umf2)
## unmarkedFrame Object
## 
## 25 sites
## Maximum number of observations per site: 5 
## Mean number of observations per site: 5 
## Sites with at least one detection: 6 
## 
## Tabulation of y observations:
##   0   1   2   3 
## 113   7   4   1
# Fit a model
fm2 <- pcount(~1 ~ 1, umf2, K=10)
fm2
## 
## Call:
## pcount(formula = ~1 ~ 1, data = umf2, K = 10)
## 
## Abundance (log-scale):
##  Estimate    SE     z P(>|z|)
##    -0.818 0.363 -2.25  0.0244
## 
## Detection (logit-scale):
##  Estimate    SE     z P(>|z|)
##    -0.725 0.416 -1.74  0.0815
## 
## AIC: 93.08232 
## Number of sites: 25
plogis(coef(fm2, type="det"))
##    p(Int) 
## 0.3261904
# Empirical Bayes estimation of random effects
(fm2re <- ranef(fm2))
##             Mean Mode 2.5% 97.5%
##  [1,] 1.06131328    1    1     2
##  [2,] 3.37663640    3    3     5
##  [3,] 0.06131328    0    0     1
##  [4,] 0.06131328    0    0     1
##  [5,] 0.06131328    0    0     1
##  [6,] 2.24934448    2    2     3
##  [7,] 0.06131328    0    0     1
##  [8,] 0.06131328    0    0     1
##  [9,] 1.06131328    1    1     2
## [10,] 0.06131328    0    0     1
## [11,] 1.06131328    1    1     2
## [12,] 0.06131328    0    0     1
## [13,] 0.06131328    0    0     1
## [14,] 0.06131328    0    0     1
## [15,] 1.06131328    1    1     2
## [16,] 0.06131328    0    0     1
## [17,] 0.06131328    0    0     1
## [18,] 0.06131328    0    0     1
## [19,] 0.06131328    0    0     1
## [20,] 0.06131328    0    0     1
## [21,] 0.06131328    0    0     1
## [22,] 0.06131328    0    0     1
## [23,] 0.06131328    0    0     1
## [24,] 0.06131328    0    0     1
## [25,] 0.06131328    0    0     1
plot(fm2re, subset=site %in% 1:8, xlim=c(-1,10))

sum(bup(fm2re)) # Estimated population size
## [1] 11.03619
colSums(confint(fm2re)) # and 95% CI
##  2.5% 97.5% 
##     9    35
## choix des modèles

cbind(AIC.fm2=fm2@AIC, AIC.fm1=fm1@AIC)
##       AIC.fm2  AIC.fm1
## [1,] 93.08232 95.53172
deltaAIC=fm2@AIC-fm1@AIC
deltaAIC # 2026, diff significative entre proba constante et variable, p constante
## [1] -2.449405
## 2. Modèle de Royle avec p CONSTANTE (p VARIABLE peut être choisie aussi) et covariables gîtes anthropiques et terriers ----

data <- data.frame(Terriers = Terriers, GA = GA, Veg = Veg, BM = BM)
data
##    Terriers GA Veg BM
## 1         1  0   1  1
## 2         0  1   0  0
## 3         0  0   0  0
## 4         1  1   1  1
## 5         0  0   1  0
## 6         1  1   1  0
## 7         0  0   1  0
## 8         0  1   1  1
## 9         0  1   1  0
## 10        0  0   1  0
## 11        1  0   1  1
## 12        1  0   1  1
## 13        0  0   0  0
## 14        0  0   0  0
## 15        0  0   1  1
## 16        0  0   0  1
## 17        0  0   1  1
## 18        0  0   0  1
## 19        0  0   1  1
## 20        0  0   0  0
## 21        0  0   1  1
## 22        0  0   1  1
## 23        0  1   1  1
## 24        0  0   0  0
## 25        0  0   1  1
# GITES ANTHROPIQUES :
nSites <- 25
nVisits <- 5
umf3 <- unmarkedFramePCount(y=all2026, siteCovs=data,obsCovs=NULL)
summary(umf3)
## unmarkedFrame Object
## 
## 25 sites
## Maximum number of observations per site: 5 
## Mean number of observations per site: 5 
## Sites with at least one detection: 6 
## 
## Tabulation of y observations:
##   0   1   2   3 
## 113   7   4   1 
## 
## Site-level covariates:
##     Terriers         GA            Veg             BM      
##  Min.   :0.0   Min.   :0.00   Min.   :0.00   Min.   :0.00  
##  1st Qu.:0.0   1st Qu.:0.00   1st Qu.:0.00   1st Qu.:0.00  
##  Median :0.0   Median :0.00   Median :1.00   Median :1.00  
##  Mean   :0.2   Mean   :0.24   Mean   :0.68   Mean   :0.56  
##  3rd Qu.:0.0   3rd Qu.:0.00   3rd Qu.:1.00   3rd Qu.:1.00  
##  Max.   :1.0   Max.   :1.00   Max.   :1.00   Max.   :1.00
# Fit a model
fm3 <- pcount(~1 ~ GA, umf3, K=5)
fm3
## 
## Call:
## pcount(formula = ~1 ~ GA, data = umf3, K = 5)
## 
## Abundance (log-scale):
##             Estimate    SE     z P(>|z|)
## (Intercept)    -1.63 0.597 -2.73 0.00640
## GA              2.06 0.696  2.96 0.00311
## 
## Detection (logit-scale):
##  Estimate    SE     z P(>|z|)
##    -0.955 0.446 -2.14  0.0323
## 
## AIC: 85.86947 
## Number of sites: 25
plogis(coef(fm3, type="det"))
##    p(Int) 
## 0.2778633
# Empirical Bayes estimation of random effects
(fm3re <- ranef(fm3))
##             Mean Mode 2.5% 97.5%
##  [1,] 1.03855813    1    1     2
##  [2,] 4.08888471    4    3     5
##  [3,] 0.03855813    0    0     1
##  [4,] 0.30166476    0    0     2
##  [5,] 0.03855813    0    0     1
##  [6,] 2.92538047    3    2     5
##  [7,] 0.03855813    0    0     1
##  [8,] 0.30166476    0    0     2
##  [9,] 1.30159241    1    1     3
## [10,] 0.03855813    0    0     1
## [11,] 1.03855813    1    1     2
## [12,] 0.03855813    0    0     1
## [13,] 0.03855813    0    0     1
## [14,] 0.03855813    0    0     1
## [15,] 1.03855813    1    1     2
## [16,] 0.03855813    0    0     1
## [17,] 0.03855813    0    0     1
## [18,] 0.03855813    0    0     1
## [19,] 0.03855813    0    0     1
## [20,] 0.03855813    0    0     1
## [21,] 0.03855813    0    0     1
## [22,] 0.03855813    0    0     1
## [23,] 0.30166476    0    0     2
## [24,] 0.03855813    0    0     1
## [25,] 0.03855813    0    0     1
plot(fm3re, subset=site %in% 1:8, xlim=c(-1,10))

sum(bup(fm3re)) # Estimated population size
## [1] 12.95346
colSums(confint(fm3re)) # and 95% CI
##  2.5% 97.5% 
##     9    41
# TERRIERS
nSites <- 25
nVisits <- 5
umf4 <- unmarkedFramePCount(y=all2026, siteCovs=data,obsCovs=NULL)
summary(umf4)
## unmarkedFrame Object
## 
## 25 sites
## Maximum number of observations per site: 5 
## Mean number of observations per site: 5 
## Sites with at least one detection: 6 
## 
## Tabulation of y observations:
##   0   1   2   3 
## 113   7   4   1 
## 
## Site-level covariates:
##     Terriers         GA            Veg             BM      
##  Min.   :0.0   Min.   :0.00   Min.   :0.00   Min.   :0.00  
##  1st Qu.:0.0   1st Qu.:0.00   1st Qu.:0.00   1st Qu.:0.00  
##  Median :0.0   Median :0.00   Median :1.00   Median :1.00  
##  Mean   :0.2   Mean   :0.24   Mean   :0.68   Mean   :0.56  
##  3rd Qu.:0.0   3rd Qu.:0.00   3rd Qu.:1.00   3rd Qu.:1.00  
##  Max.   :1.0   Max.   :1.00   Max.   :1.00   Max.   :1.00
# Fit a model
fm4 <- pcount(~1 ~ Terriers, umf4, K=20)
fm4
## 
## Call:
## pcount(formula = ~1 ~ Terriers, data = umf4, K = 20)
## 
## Abundance (log-scale):
##             Estimate    SE     z P(>|z|)
## (Intercept)    -1.18 0.471 -2.51  0.0122
## Terriers        1.21 0.666  1.82  0.0692
## 
## Detection (logit-scale):
##  Estimate    SE     z P(>|z|)
##     -0.76 0.429 -1.77  0.0767
## 
## AIC: 92.09703 
## Number of sites: 25
plogis(coef(fm4, type="det"))
##    p(Int) 
## 0.3187396
# Empirical Bayes estimation of random effects
(fm4re <- ranef(fm4))
##             Mean Mode 2.5% 97.5%
##  [1,] 1.15118825    1    1     2
##  [2,] 3.28667321    3    3     4
##  [3,] 0.04507392    0    0     1
##  [4,] 0.15118825    0    0     1
##  [5,] 0.04507392    0    0     1
##  [6,] 2.54667410    2    2     4
##  [7,] 0.04507392    0    0     1
##  [8,] 0.04507392    0    0     1
##  [9,] 1.04507392    1    1     2
## [10,] 0.04507392    0    0     1
## [11,] 1.15118825    1    1     2
## [12,] 0.15118825    0    0     1
## [13,] 0.04507392    0    0     1
## [14,] 0.04507392    0    0     1
## [15,] 1.04507392    1    1     2
## [16,] 0.04507392    0    0     1
## [17,] 0.04507392    0    0     1
## [18,] 0.04507392    0    0     1
## [19,] 0.04507392    0    0     1
## [20,] 0.04507392    0    0     1
## [21,] 0.04507392    0    0     1
## [22,] 0.04507392    0    0     1
## [23,] 0.04507392    0    0     1
## [24,] 0.04507392    0    0     1
## [25,] 0.04507392    0    0     1
plot(fm4re, subset=site %in% 1:8, xlim=c(-1,10))

sum(bup(fm4re)) # Estimated population size
## [1] 11.2945
colSums(confint(fm4re)) # and 95% CI
##  2.5% 97.5% 
##     9    35
# TERRIERS + GITES ANTHROPIQUES
nSites <- 25
nVisits <- 5
umf5 <- unmarkedFramePCount(y=all2026, siteCovs=data,obsCovs=NULL)
summary(umf5)
## unmarkedFrame Object
## 
## 25 sites
## Maximum number of observations per site: 5 
## Mean number of observations per site: 5 
## Sites with at least one detection: 6 
## 
## Tabulation of y observations:
##   0   1   2   3 
## 113   7   4   1 
## 
## Site-level covariates:
##     Terriers         GA            Veg             BM      
##  Min.   :0.0   Min.   :0.00   Min.   :0.00   Min.   :0.00  
##  1st Qu.:0.0   1st Qu.:0.00   1st Qu.:0.00   1st Qu.:0.00  
##  Median :0.0   Median :0.00   Median :1.00   Median :1.00  
##  Mean   :0.2   Mean   :0.24   Mean   :0.68   Mean   :0.56  
##  3rd Qu.:0.0   3rd Qu.:0.00   3rd Qu.:1.00   3rd Qu.:1.00  
##  Max.   :1.0   Max.   :1.00   Max.   :1.00   Max.   :1.00
# Fit a model
fm5 <- pcount(~1 ~ GA+Terriers, umf5, K=20)
fm5
## 
## Call:
## pcount(formula = ~1 ~ GA + Terriers, data = umf5, K = 20)
## 
## Abundance (log-scale):
##             Estimate    SE     z P(>|z|)
## (Intercept)   -1.732 0.654 -2.65 0.00806
## GA             1.984 0.729  2.72 0.00652
## Terriers       0.756 0.656  1.15 0.24940
## 
## Detection (logit-scale):
##  Estimate    SE    z P(>|z|)
##     -1.12 0.588 -1.9   0.058
## 
## AIC: 86.25589 
## Number of sites: 25
plogis(coef(fm5, type="det"))
##    p(Int) 
## 0.2469176
# Empirical Bayes estimation of random effects
(fm5re <- ranef(fm5))
##             Mean Mode 2.5% 97.5%
##  [1,] 1.09120465    1    1     2
##  [2,] 4.35385811    4    3     6
##  [3,] 0.04283634    0    0     1
##  [4,] 0.66325976    0    0     3
##  [5,] 0.04283634    0    0     1
##  [6,] 3.66927624    3    2     6
##  [7,] 0.04283634    0    0     1
##  [8,] 0.31151505    0    0     2
##  [9,] 1.31151505    1    1     3
## [10,] 0.04283634    0    0     1
## [11,] 1.09120465    1    1     2
## [12,] 0.09120465    0    0     1
## [13,] 0.04283634    0    0     1
## [14,] 0.04283634    0    0     1
## [15,] 1.04283634    1    1     2
## [16,] 0.04283634    0    0     1
## [17,] 0.04283634    0    0     1
## [18,] 0.04283634    0    0     1
## [19,] 0.04283634    0    0     1
## [20,] 0.04283634    0    0     1
## [21,] 0.04283634    0    0     1
## [22,] 0.04283634    0    0     1
## [23,] 0.31151505    0    0     2
## [24,] 0.04283634    0    0     1
## [25,] 0.04283634    0    0     1
plot(fm5re, subset=site %in% 1:8, xlim=c(-1,10))

sum(bup(fm5re)) # Estimated population size
## [1] 14.57993
colSums(confint(fm5re)) # and 95% CI
##  2.5% 97.5% 
##     9    44
# VEGETATION : 
nSites <- 25
nVisits <- 5
umf6 <- unmarkedFramePCount(y=all2026, siteCovs=data,obsCovs=NULL)
summary(umf6)
## unmarkedFrame Object
## 
## 25 sites
## Maximum number of observations per site: 5 
## Mean number of observations per site: 5 
## Sites with at least one detection: 6 
## 
## Tabulation of y observations:
##   0   1   2   3 
## 113   7   4   1 
## 
## Site-level covariates:
##     Terriers         GA            Veg             BM      
##  Min.   :0.0   Min.   :0.00   Min.   :0.00   Min.   :0.00  
##  1st Qu.:0.0   1st Qu.:0.00   1st Qu.:0.00   1st Qu.:0.00  
##  Median :0.0   Median :0.00   Median :1.00   Median :1.00  
##  Mean   :0.2   Mean   :0.24   Mean   :0.68   Mean   :0.56  
##  3rd Qu.:0.0   3rd Qu.:0.00   3rd Qu.:1.00   3rd Qu.:1.00  
##  Max.   :1.0   Max.   :1.00   Max.   :1.00   Max.   :1.00
# Fit a model
fm6 <- pcount(~1 ~ Veg, umf6, K=20)
fm6
## 
## Call:
## pcount(formula = ~1 ~ Veg, data = umf6, K = 20)
## 
## Abundance (log-scale):
##             Estimate    SE      z P(>|z|)
## (Intercept)   -0.721 0.602 -1.198   0.231
## Veg           -0.140 0.704 -0.198   0.843
## 
## Detection (logit-scale):
##  Estimate   SE     z P(>|z|)
##    -0.732 0.42 -1.74  0.0816
## 
## AIC: 95.04373 
## Number of sites: 25
plogis(coef(fm6, type="det"))
##    p(Int) 
## 0.3248282
# Empirical Bayes estimation of random effects
(fm6re <- ranef(fm6))
##             Mean Mode 2.5% 97.5%
##  [1,] 1.05935512    1    1     2
##  [2,] 3.41332176    3    3     5
##  [3,] 0.06824072    0    0     1
##  [4,] 0.05935512    0    0     1
##  [5,] 0.05935512    0    0     1
##  [6,] 2.24208834    2    2     3
##  [7,] 0.05935512    0    0     1
##  [8,] 0.05935512    0    0     1
##  [9,] 1.05935512    1    1     2
## [10,] 0.05935512    0    0     1
## [11,] 1.05935512    1    1     2
## [12,] 0.05935512    0    0     1
## [13,] 0.06824072    0    0     1
## [14,] 0.06824072    0    0     1
## [15,] 1.05935512    1    1     2
## [16,] 0.06824072    0    0     1
## [17,] 0.05935512    0    0     1
## [18,] 0.06824072    0    0     1
## [19,] 0.05935512    0    0     1
## [20,] 0.06824072    0    0     1
## [21,] 0.05935512    0    0     1
## [22,] 0.05935512    0    0     1
## [23,] 0.05935512    0    0     1
## [24,] 0.06824072    0    0     1
## [25,] 0.05935512    0    0     1
plot(fm6re, subset=site %in% 1:8, xlim=c(-1,10))

sum(bup(fm6re)) # Estimated population size
## [1] 11.08278
colSums(confint(fm6re)) # and 95% CI
##  2.5% 97.5% 
##     9    35
#Bois mort
nSites <- 25
nVisits <- 5
umf7 <- unmarkedFramePCount(y=all2026, siteCovs=data,obsCovs=NULL)
summary(umf7)
## unmarkedFrame Object
## 
## 25 sites
## Maximum number of observations per site: 5 
## Mean number of observations per site: 5 
## Sites with at least one detection: 6 
## 
## Tabulation of y observations:
##   0   1   2   3 
## 113   7   4   1 
## 
## Site-level covariates:
##     Terriers         GA            Veg             BM      
##  Min.   :0.0   Min.   :0.00   Min.   :0.00   Min.   :0.00  
##  1st Qu.:0.0   1st Qu.:0.00   1st Qu.:0.00   1st Qu.:0.00  
##  Median :0.0   Median :0.00   Median :1.00   Median :1.00  
##  Mean   :0.2   Mean   :0.24   Mean   :0.68   Mean   :0.56  
##  3rd Qu.:0.0   3rd Qu.:0.00   3rd Qu.:1.00   3rd Qu.:1.00  
##  Max.   :1.0   Max.   :1.00   Max.   :1.00   Max.   :1.00
# Fit a model
fm7 <- pcount(~1 ~ BM, umf7, K=20)
fm7
## 
## Call:
## pcount(formula = ~1 ~ BM, data = umf7, K = 20)
## 
## Abundance (log-scale):
##             Estimate    SE      z P(>|z|)
## (Intercept)   -0.295 0.461 -0.639   0.523
## BM            -1.070 0.709 -1.509   0.131
## 
## Detection (logit-scale):
##  Estimate   SE     z P(>|z|)
##    -0.819 0.46 -1.78   0.075
## 
## AIC: 92.63187 
## Number of sites: 25
plogis(coef(fm7, type="det"))
##    p(Int) 
## 0.3059415
# Empirical Bayes estimation of random effects
(fm7re <- ranef(fm7))
##             Mean Mode 2.5% 97.5%
##  [1,] 1.04113705    1    1     2
##  [2,] 3.66147554    3    3     5
##  [3,] 0.11994850    0    0     1
##  [4,] 0.04113705    0    0     1
##  [5,] 0.11994850    0    0     1
##  [6,] 2.45028385    2    2     4
##  [7,] 0.11994850    0    0     1
##  [8,] 0.04113705    0    0     1
##  [9,] 1.11994850    1    1     2
## [10,] 0.11994850    0    0     1
## [11,] 1.04113705    1    1     2
## [12,] 0.04113705    0    0     1
## [13,] 0.11994850    0    0     1
## [14,] 0.11994850    0    0     1
## [15,] 1.04113705    1    1     2
## [16,] 0.04113705    0    0     1
## [17,] 0.04113705    0    0     1
## [18,] 0.04113705    0    0     1
## [19,] 0.04113705    0    0     1
## [20,] 0.11994850    0    0     1
## [21,] 0.04113705    0    0     1
## [22,] 0.04113705    0    0     1
## [23,] 0.04113705    0    0     1
## [24,] 0.11994850    0    0     1
## [25,] 0.04113705    0    0     1
plot(fm7re, subset=site %in% 1:8, xlim=c(-1,10))

sum(bup(fm7re)) # Estimated population size
## [1] 11.76721
colSums(confint(fm7re)) # and 95% CI
##  2.5% 97.5% 
##     9    36
#Comparaison des modèles d'abondance
cbind(AIC.fm1=fm1@AIC, AIC.fm2=fm2@AIC, AIC.fm3=fm3@AIC, AIC.fm4=fm4@AIC, AIC.fm5=fm5@AIC, AIC.fm6=fm6@AIC, AIC.fm7=fm7@AIC)
##       AIC.fm1  AIC.fm2  AIC.fm3  AIC.fm4  AIC.fm5  AIC.fm6  AIC.fm7
## [1,] 95.53172 93.08232 85.86947 92.09703 86.25589 95.04373 92.63187

Fm3 et fm5 se dissocient des autres en étant meilleurs.

Ces mêmes modèles et comparaisons ont été mis en place pour chaque année

Graphiques pour comparer l’évolution de l’abondance entre les campagnes annuelles de suivi

df_estimation <- data.frame(
  Annee = c(2015, 2018, 2020, 2022, 2024, 2026),
  Estimation = c(37.2, 29.4, 21.87, 15.15, 18.75, 12.95),
  IC_inf = c(21, 14, 14, 11, 12, 9),
  IC_sup = c(69, 61, 33, 27, 49, 41))
#graphique en histogramme
library(ggplot2)
ggplot(df_estimation, aes(x = Annee, y = Estimation)) +
  
geom_col(fill = "lightgreen",colour = "black",width = 0.7) +
  geom_errorbar(aes(ymin = IC_inf,ymax = IC_sup),width = 0.15,linewidth = 0.2,colour = "darkgreen") +
  geom_text(aes(label = round(Estimation,1)),vjust = -0.6,size = 4, colour="black") +
  scale_x_continuous(breaks = 2015:2026) +
  labs(x = "Année",y = "Abondance estimée") +
  ylim(0,75) +
  theme_classic(base_size = 12) +
  theme(plot.title = element_text(face = "bold", hjust = 0.5),
    axis.title = element_text(face = "bold"),
    axis.text = element_text(colour = "black"))

Régression linéaire

#Régression linéaire pour montrer si il existe une tendance à l'augmentation ou à la diminution des abondances estimées

mod <- lm(Estimation ~ Annee, data = df_estimation)

summary(mod)
## 
## Call:
## lm(formula = Estimation ~ Annee, data = df_estimation)
## 
## Residuals:
##       1       2       3       4       5       6 
##  2.0491  0.7278 -2.4830 -4.8838  3.0354  1.5545 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 4386.7200   765.3312   5.732  0.00459 **
## Annee         -2.1596     0.3787  -5.702  0.00467 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.405 on 4 degrees of freedom
## Multiple R-squared:  0.8905, Adjusted R-squared:  0.8631 
## F-statistic: 32.52 on 1 and 4 DF,  p-value: 0.004675
ggplot(df_estimation, aes(Annee, Estimation))+
  geom_point(size=1)+
  geom_smooth(method="lm", se=TRUE)
## `geom_smooth()` using formula = 'y ~ x'

Cela signifie que l’abondance estimée diminue en moyenne de 2,16 individus tous les deux ans (dans ton jeu de données, les suivis sont réalisés tous les deux ans), soit environ 1,08 individu par an. Cette diminution est statistiquement significative (t = -5,70 ; p = 0,0047). Le coefficient de détermination est très élevé : R2=0,89 Cela signifie que 89 % de la variabilité des estimations est expliquée par l’année de suivi, ce qui traduit une forte tendance temporelle.

#corrélation de spearman
cor.test(df_estimation$Annee,df_estimation$Estimation,method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  df_estimation$Annee and df_estimation$Estimation
## S = 68, p-value = 0.01667
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.9428571

ρ=−0,943. Cette valeur est très proche de -1, ce qui indique une très forte relation monotone décroissante entre l’année et les effectifs estimés. La relation est significative : p=0,0167

#Test de Mann-Kendall
library(Kendall)
## Warning: le package 'Kendall' a été compilé avec la version R 4.5.3
MannKendall(df_estimation$Estimation)
## tau = -0.867, 2-sided pvalue =0.024171

τ=−0,867 avec p=0,024 Le Tau négatif indique une tendance décroissante. Comme la p-value est inférieure à 0,05, cette tendance est statistiquement significative.

Ces tests ne prennent pas en compte les intervalles de confiance très large Hors ici les estimations d’abondance présentent des intervalles très large donc une régression pondérée a été mis en place afin de donner plus de poid aux estimations ayant un intervalle de confiance plus étroit

#Régression pondérée

#calcul de SE
SE <- (df_estimation$IC_sup - df_estimation$IC_inf)/(2*1.96)

Cette conversion est approximative, car tes intervalles ne sont pas symétriques autour de l’estimation. Par exemple, pour 2015 : estimation : 37,2 ; distance à la borne inférieure : 16,2 ; distance à la borne supérieure : 31,8. Cela suggère que la distribution de l’abondance estimée n’est probablement pas parfaitement normale. L’utilisation d’une seule erreur standard simplifie donc l’incertitude réelle

#Régression pondéré
mod <- lm(
  Estimation ~ Annee,
  data = df_estimation,
  weights = 1/SE^2)

Une année dont l’estimation est précise, donc avec un petit SE, reçoit un poids important. À l’inverse, une estimation dont l’intervalle de confiance est large influence moins la droite. Par exemple, 2021 reçoit davantage de poids que 2015, car son intervalle de confiance est beaucoup plus étroit.

summary(mod)
## 
## Call:
## lm(formula = Estimation ~ Annee, data = df_estimation, weights = 1/SE^2)
## 
## Weighted Residuals:
##        1        2        3        4        5        6 
##  0.40320  0.27746 -0.01587 -0.65428  0.53564  0.41432 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 4189.5812  1162.4847   3.604   0.0227 *
## Annee         -2.0632     0.5751  -3.588   0.0230 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5307 on 4 degrees of freedom
## Multiple R-squared:  0.7629, Adjusted R-squared:  0.7036 
## F-statistic: 12.87 on 1 and 4 DF,  p-value: 0.02301

L’intercept de 4189,58 n’a pas d’interprétation biologique utile : il correspond à l’abondance théorique en l’an 0. Ce nombre élevé vient simplement du fait que les années sont codées 2015, 2017, etc. La valeur importante est la pente : β=−2,063 Cela signifie que l’abondance estimée diminue en moyenne de : 2,06 individus par année, soit environ 4,13 individus entre deux campagnes espacées de deux ans. La pente est négative et la p-value est inférieure à 0,05. D’après cette régression pondérée, il existe donc une diminution statistiquement significative de l’abondance estimée entre 2015 et 2025. R2 : Cela indique qu’environ 76 % de la variation des estimations pondérées est associée à l’année. Il vaut mieux écrire « associée » ou « expliquée statistiquement » plutôt que conclure que le temps provoque directement le déclin. Dans une régression avec un seul prédicteur, le test F conduit à la même conclusion que le test de la pente : le modèle contenant l’année est significativement plus informatif qu’un modèle ne contenant qu’une constante. La régression linéaire pondérée par l’inverse de la variance met en évidence une diminution significative de l’abondance estimée du Lézard ocellé entre 2015 et 2025. La pente estimée est de −2,06 individus par an, soit une diminution moyenne d’environ 4,13 individus entre deux campagnes successives espacées de deux ans (R² = 0,76 ; p = 0,023).

#Graphique de la régression linéaire pondérée
ggplot(df_estimation, aes(x = Annee, y = Estimation)) +
  geom_point(size = 1) +
  geom_smooth(
    aes(weight = 1/SE^2),
    method = "lm",
    se = TRUE)
## `geom_smooth()` using formula = 'y ~ x'