1 Intallation des packages

library(dplyr);library(ggplot2);library(sf);library(unmarked)

2 Importation des bases de donnees

Total <- sf::st_read(paste0(wd$data, "derived/Total.shp"),
                     options = "ENCODING=UTF8")
## options:        ENCODING=UTF8 
## Reading layer `Total' from data source 
##   `/home/onyxia/work/AnalyseDonneesOpportunisteGMB/data/derived/Total.shp' using driver `ESRI Shapefile'
## Simple feature collection with 160754 features and 14 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -4.791511 ymin: 47.44135 xmax: -1.034921 ymax: 48.88095
## Geodetic CRS:  WGS 84
Total <- Total%>%
  transform_Total()

VariablesSite <- sf::st_read(paste0(wd$data, "derived/VariablesSite.shp"),
                             options = "ENCODING=UTF8")
## options:        ENCODING=UTF8 
## Reading layer `VariablesSite' from data source 
##   `/home/onyxia/work/AnalyseDonneesOpportunisteGMB/data/derived/VariablesSite.shp' using driver `ESRI Shapefile'
## Simple feature collection with 328 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -4.798138 ymin: 47.35903 xmax: -1.013009 ymax: 48.88628
## Geodetic CRS:  WGS 84
VariablesSite <- VariablesSite%>%
  transform_VarSites()
Lapins <- Total%>%
  sf::st_join(VariablesSite)%>%
  filter(cd_nom == 61714,
         technique_observation == "Vu",
         date > as.Date("2005-01-01"))%>%
  # select(-c(insee_dept, lib_dept, lib_dept, FID, surf, CD_SIG)) %>%
  select(-nom_valide, -nom_vernaculaire, -ordre, -technique_observation)%>%
  mutate(mois = lubridate::month(date),
         mois = case_when(mois == 1 ~ "Jan", 
                          mois == 2 ~ "Fev",
                          mois == 3 ~ "Mars",
                          mois == 4 ~ "Avr",
                          mois == 5 ~ "Mai",
                          mois == 6 ~ "Juin",
                          mois == 7 ~ "Juil",
                          mois == 8 ~ "Aout",
                          mois == 9 ~ "Sept",
                          mois == 10 ~ "Oct",
                          mois == 11 ~ "Nov",
                          mois == 12 ~ "Dec"),
         year = lubridate::year(date))

3 Creation de la base de donnee Detection - Non Detection

t_lapins <- Lapins %>%
  sf::st_drop_geometry()%>%
  dplyr::mutate(annee = lubridate::year(date))%>%
  dplyr::select(Code_10km, annee)%>%
  dplyr::count(Code_10km, annee) %>%
  tidyr::complete(Code_10km, 
                  annee = tidyr::full_seq(annee, 1), 
                  fill = list(n = 0)) %>%
  tidyr::pivot_wider(names_from = annee, 
                     values_from = n) %>%
  dplyr::arrange(Code_10km) %>%
  filter(!is.na(Code_10km))

t_lapins[,-1] <- ifelse(t_lapins[,-1]>0, 1, 0)

noms <- t_lapins$Code_10km
t_lapins <- as.matrix(t_lapins[, -1])
rownames(t_lapins) <- noms

Observation des 10 premiers sites pour verifier

head(t_lapins)
##          2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022
## E012N679    0    0    0    0    0    0    0    1    0    1    1    1    1    1    1    0    1    1
## E012N680    0    0    0    0    0    0    0    0    0    0    0    1    1    0    1    0    1    1
## E012N683    0    0    0    0    0    0    0    0    1    0    1    1    1    1    1    1    0    1
## E012N684    0    0    0    0    0    0    0    0    0    0    0    0    1    0    1    1    1    1
## E012N685    0    0    0    0    0    0    0    0    1    1    1    0    0    0    1    1    1    1
## E013N679    0    1    0    0    0    0    0    1    1    0    1    0    0    1    1    1    1    1
##          2023 2024 2025
## E012N679    1    1    0
## E012N680    0    0    0
## E012N683    1    1    0
## E012N684    1    1    1
## E012N685    1    1    1
## E013N679    1    1    1

4 Covariables de observations

On separe les observations en 2 periodes de 10 ans : Avant 2015 et Apres 2015

annees <- as.numeric(colnames(t_lapins))
periodes_vector <- ifelse(annees < 2015, "periode1", "periode2")

periodes_matrix <- matrix(rep(periodes_vector, each = nrow(t_lapins)), 
                          nrow = nrow(t_lapins), 
                          byrow = FALSE)
colnames(periodes_matrix) <- colnames(t_lapins)
rownames(periodes_matrix) <- rownames(t_lapins)

periodes_df <- as.data.frame(periodes_matrix)
head(periodes_df)
##              2005     2006     2007     2008     2009     2010     2011     2012     2013     2014
## E012N679 periode1 periode1 periode1 periode1 periode1 periode1 periode1 periode1 periode1 periode1
## E012N680 periode1 periode1 periode1 periode1 periode1 periode1 periode1 periode1 periode1 periode1
## E012N683 periode1 periode1 periode1 periode1 periode1 periode1 periode1 periode1 periode1 periode1
## E012N684 periode1 periode1 periode1 periode1 periode1 periode1 periode1 periode1 periode1 periode1
## E012N685 periode1 periode1 periode1 periode1 periode1 periode1 periode1 periode1 periode1 periode1
## E013N679 periode1 periode1 periode1 periode1 periode1 periode1 periode1 periode1 periode1 periode1
##              2015     2016     2017     2018     2019     2020     2021     2022     2023     2024
## E012N679 periode2 periode2 periode2 periode2 periode2 periode2 periode2 periode2 periode2 periode2
## E012N680 periode2 periode2 periode2 periode2 periode2 periode2 periode2 periode2 periode2 periode2
## E012N683 periode2 periode2 periode2 periode2 periode2 periode2 periode2 periode2 periode2 periode2
## E012N684 periode2 periode2 periode2 periode2 periode2 periode2 periode2 periode2 periode2 periode2
## E012N685 periode2 periode2 periode2 periode2 periode2 periode2 periode2 periode2 periode2 periode2
## E013N679 periode2 periode2 periode2 periode2 periode2 periode2 periode2 periode2 periode2 periode2
##              2025
## E012N679 periode2
## E012N680 periode2
## E012N683 periode2
## E012N684 periode2
## E012N685 periode2
## E013N679 periode2

5 Covariables de site

site_info <- VariablesSite%>%
  dplyr::arrange(Code_10km)

noms <- site_info$Code_10km
rownames(site_info) <- noms

site_info <- site_info%>%
  dplyr::select(-Code_10km)%>%
  sf::st_drop_geometry()
head(site_info)
##          Ind_Diversite Dnst_Cultures  Dist_EcotoneArbore Dist_Littoral Dist_Eau
## E012N679      2.398615     0.09440944          256.87302     0.7661493 358.3517
## E012N680      2.374172     0.24340895          142.15814     0.4119377 257.5322
## E012N683      2.113717     0.35821721           75.82455     1.1937126 215.3183
## E012N684      2.323193     0.33730915           70.87636     1.4593335 219.9087
## E012N685      2.232987     0.28420287           92.34334     0.7645159 300.4883
## E013N679      2.373326     0.29707381           79.00361     1.3386052 182.1931
##                     Famille_paysage     Nom_paysage   air_km2    X_10km   Y_10km
## E012N679 Paysage cultive avec talus Baie d'Audierne 13.412471 -4.692649 48.03723
## E012N680 Paysage cultive avec talus Baie d'Audierne  7.146477 -4.689630 48.05803
## E012N683 Paysage cultive avec talus Plateau leonard 44.667218 -4.743275 48.36667
## E012N684 Paysage cultive avec talus Plateau leonard 45.037581 -4.752107 48.44740
## E012N685 Paysage cultive avec talus Plateau leonard 15.914374 -4.751892 48.51838
## E013N679 Paysage cultive avec talus Baie d'Audierne 51.271751 -4.584359 48.03197

6 Implementation du Model

head(t_lapins)
##          2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022
## E012N679    0    0    0    0    0    0    0    1    0    1    1    1    1    1    1    0    1    1
## E012N680    0    0    0    0    0    0    0    0    0    0    0    1    1    0    1    0    1    1
## E012N683    0    0    0    0    0    0    0    0    1    0    1    1    1    1    1    1    0    1
## E012N684    0    0    0    0    0    0    0    0    0    0    0    0    1    0    1    1    1    1
## E012N685    0    0    0    0    0    0    0    0    1    1    1    0    0    0    1    1    1    1
## E013N679    0    1    0    0    0    0    0    1    1    0    1    0    0    1    1    1    1    1
##          2023 2024 2025
## E012N679    1    1    0
## E012N680    0    0    0
## E012N683    1    1    0
## E012N684    1    1    1
## E012N685    1    1    1
## E013N679    1    1    1

Mise sous format pour la modélisation

umf <- unmarked::unmarkedFrameOccu(y = t_lapins, 
                                   obsCovs = list(periode = periodes_df),
                                   siteCovs = site_info
                                   )
## Warning: obsCovs contains characters. Converting them to factors.
head(umf); summary(umf)
## Data frame representation of unmarkedFrame object.
##          y.1 y.2 y.3 y.4 y.5 y.6 y.7 y.8 y.9 y.10 y.11 y.12 y.13 y.14 y.15 y.16 y.17 y.18 y.19 y.20
## E012N679   0   0   0   0   0   0   0   1   0    1    1    1    1    1    1    0    1    1    1    1
## E012N680   0   0   0   0   0   0   0   0   0    0    0    1    1    0    1    0    1    1    0    0
## E012N683   0   0   0   0   0   0   0   0   1    0    1    1    1    1    1    1    0    1    1    1
## E012N684   0   0   0   0   0   0   0   0   0    0    0    0    1    0    1    1    1    1    1    1
## E012N685   0   0   0   0   0   0   0   0   1    1    1    0    0    0    1    1    1    1    1    1
## E013N679   0   1   0   0   0   0   0   1   1    0    1    0    0    1    1    1    1    1    1    1
## E013N680   0   0   0   0   0   0   0   0   0    0    0    0    1    1    1    0    1    1    1    1
## E013N681   0   0   0   0   1   1   1   1   1    1    1    1    1    1    0    0    1    1    0    0
## E013N682   0   0   0   0   1   1   1   0   0    0    0    0    1    1    1    1    1    1    1    1
## E013N683   0   0   0   0   0   0   0   1   0    1    1    1    1    0    1    1    1    1    1    1
##          y.21 Ind_Diversite Dnst_Cultures. Dist_EcotoneArbore Dist_Littoral  Dist_Eau
## E012N679    0      2.398615     0.09440944          256.87302     0.7661493  358.3517
## E012N680    0      2.374172     0.24340895          142.15814     0.4119377  257.5322
## E012N683    0      2.113717     0.35821721           75.82455     1.1937126  215.3183
## E012N684    1      2.323193     0.33730915           70.87636     1.4593335  219.9087
## E012N685    1      2.232987     0.28420287           92.34334     0.7645159  300.4883
## E013N679    1      2.373326     0.29707381           79.00361     1.3386052  182.1931
## E013N680    1      2.023110     0.36020261          115.81535     0.9278021  217.1809
## E013N681    0      2.371093     0.00770415          766.33655     0.2196287 1012.0054
## E013N682    0      2.245646     0.02477821          231.63570     0.4566141  397.3353
## E013N683    1      1.907151     0.45765102           48.66792     3.2717547  218.9354
##                     Famille_paysage           Nom_paysage   air_km2    X_10km   Y_10km periode.1
## E012N679 Paysage cultive avec talus       Baie d'Audierne 13.412471 -4.692649 48.03723  periode1
## E012N680 Paysage cultive avec talus       Baie d'Audierne  7.146477 -4.689630 48.05803  periode1
## E012N683 Paysage cultive avec talus       Plateau leonard 44.667218 -4.743275 48.36667  periode1
## E012N684 Paysage cultive avec talus       Plateau leonard 45.037581 -4.752107 48.44740  periode1
## E012N685 Paysage cultive avec talus       Plateau leonard 15.914374 -4.751892 48.51838  periode1
## E013N679 Paysage cultive avec talus       Baie d'Audierne 51.271751 -4.584359 48.03197  periode1
## E013N680 Paysage cultive avec talus       Baie d'Audierne 25.225340 -4.592079 48.06757  periode1
## E013N681 Paysage cultive avec talus Crozon-Aulne maritime  6.678999 -4.556092 48.20397  periode1
## E013N682 Paysage cultive avec talus Crozon-Aulne maritime 21.927737 -4.587397 48.27703  periode1
## E013N683 Paysage cultive avec talus Crozon-Aulne maritime 78.753249 -4.639600 48.37835  periode1
##          periode.2 periode.3 periode.4 periode.5 periode.6 periode.7 periode.8 periode.9 periode.10
## E012N679  periode1  periode1  periode1  periode1  periode1  periode1  periode1  periode1   periode1
## E012N680  periode1  periode1  periode1  periode1  periode1  periode1  periode1  periode1   periode1
## E012N683  periode1  periode1  periode1  periode1  periode1  periode1  periode1  periode1   periode1
## E012N684  periode1  periode1  periode1  periode1  periode1  periode1  periode1  periode1   periode1
## E012N685  periode1  periode1  periode1  periode1  periode1  periode1  periode1  periode1   periode1
## E013N679  periode1  periode1  periode1  periode1  periode1  periode1  periode1  periode1   periode1
## E013N680  periode1  periode1  periode1  periode1  periode1  periode1  periode1  periode1   periode1
## E013N681  periode1  periode1  periode1  periode1  periode1  periode1  periode1  periode1   periode1
## E013N682  periode1  periode1  periode1  periode1  periode1  periode1  periode1  periode1   periode1
## E013N683  periode1  periode1  periode1  periode1  periode1  periode1  periode1  periode1   periode1
##          periode.11 periode.12 periode.13 periode.14 periode.15 periode.16 periode.17 periode.18
## E012N679   periode2   periode2   periode2   periode2   periode2   periode2   periode2   periode2
## E012N680   periode2   periode2   periode2   periode2   periode2   periode2   periode2   periode2
## E012N683   periode2   periode2   periode2   periode2   periode2   periode2   periode2   periode2
## E012N684   periode2   periode2   periode2   periode2   periode2   periode2   periode2   periode2
## E012N685   periode2   periode2   periode2   periode2   periode2   periode2   periode2   periode2
## E013N679   periode2   periode2   periode2   periode2   periode2   periode2   periode2   periode2
## E013N680   periode2   periode2   periode2   periode2   periode2   periode2   periode2   periode2
## E013N681   periode2   periode2   periode2   periode2   periode2   periode2   periode2   periode2
## E013N682   periode2   periode2   periode2   periode2   periode2   periode2   periode2   periode2
## E013N683   periode2   periode2   periode2   periode2   periode2   periode2   periode2   periode2
##          periode.19 periode.20 periode.21
## E012N679   periode2   periode2   periode2
## E012N680   periode2   periode2   periode2
## E012N683   periode2   periode2   periode2
## E012N684   periode2   periode2   periode2
## E012N685   periode2   periode2   periode2
## E013N679   periode2   periode2   periode2
## E013N680   periode2   periode2   periode2
## E013N681   periode2   periode2   periode2
## E013N682   periode2   periode2   periode2
## E013N683   periode2   periode2   periode2
## unmarkedFrame Object
## 
## 328 sites
## Maximum number of observations per site: 21 
## Mean number of observations per site: 21 
## Sites with at least one detection: 316 
## 
## Tabulation of y observations:
##    0    1 
## 4446 2442 
## 
## Site-level covariates:
##  Ind_Diversite   Dnst_Cultures      Dist_EcotoneArbore Dist_Littoral        Dist_Eau      
##  Min.   :1.445   Min.   :0.007704   Min.   : 23.80     Min.   : 0.2196   Min.   :  77.86  
##  1st Qu.:1.923   1st Qu.:0.288156   1st Qu.: 40.38     1st Qu.: 2.7277   1st Qu.: 204.78  
##  Median :2.051   Median :0.388762   Median : 48.53     Median :15.1991   Median : 237.46  
##  Mean   :2.076   Mean   :0.376411   Mean   : 58.12     Mean   :20.9589   Mean   : 257.46  
##  3rd Qu.:2.218   3rd Qu.:0.466718   3rd Qu.: 59.67     3rd Qu.:34.7019   3rd Qu.: 293.17  
##  Max.   :3.061   Max.   :0.694237   Max.   :766.34     Max.   :79.8731   Max.   :1399.09  
##                                                                                           
##                              Famille_paysage                        Nom_paysage     air_km2       
##  Paysage cultive a ragosses          :93     Reliefs des Landes de Lanvaux: 36   Min.   :  2.294  
##  Paysage boise et de bosquets        :57     Plateau de Penthievre        : 22   1st Qu.: 78.606  
##  Paysage cultive avec talus          :56     Armor morbihannais           : 19   Median :100.000  
##  Paysage de bocage a maille elargie  :41     Bassin de Pontivy-Loudeac    : 18   Mean   : 84.662  
##  Paysage de bocage dense sur collines:34     Bassin de Rennes             : 17   3rd Qu.:100.000  
##  Paysage de littoral urbanise        :29     Cornouaille interieure       : 16   Max.   :100.000  
##  (Other)                             :18     (Other)                      :200                    
##      X_10km           Y_10km     
##  Min.   :-4.752   Min.   :47.41  
##  1st Qu.:-3.668   1st Qu.:47.92  
##  Median :-2.907   Median :48.18  
##  Mean   :-2.896   Mean   :48.18  
##  3rd Qu.:-2.140   3rd Qu.:48.45  
##  Max.   :-1.021   Max.   :48.87  
##                                  
## 
## Observation-level covariates:
##      periode    
##  periode1:3280  
##  periode2:3608

7 Resultats

Ici on regarde 3 modeles seulement (car pb de colinearite?) pour l’instant. D’autres se rajouteront evidamment.

Un modele nul, on l’ont ne regarde que …

model_null <- unmarked::occu(~1 ~1, data = umf)
unmarked::summary(model_null)
## 
## Call:
## unmarked::occu(formula = ~1 ~ 1, data = umf)
## 
## Occupancy (logit-scale):
##  Estimate    SE    z  P(>|z|)
##      3.27 0.294 11.1 1.02e-28
## 
## Detection (logit-scale):
##  Estimate     SE     z   P(>|z|)
##    -0.541 0.0255 -21.2 3.93e-100
## 
## AIC: 8838.276 
## Number of sites: 328

Un modele avec la periode en covariance

model_periodes <- unmarked::occu(~periode ~1, data = umf)
unmarked::summary(model_periodes)
## 
## Call:
## unmarked::occu(formula = ~periode ~ 1, data = umf)
## 
## Occupancy (logit-scale):
##  Estimate    SE    z  P(>|z|)
##      3.27 0.295 11.1 1.12e-28
## 
## Detection (logit-scale):
##                 Estimate     SE     z   P(>|z|)
## (Intercept)       -0.924 0.0394 -23.4 2.54e-121
## periodeperiode2    0.700 0.0522  13.4  5.10e-41
## 
## AIC: 8656.4 
## Number of sites: 328
# model_eau <- unmarked::occu(~1 ~ Dist_Eau, data = umf)
# unmarked::summary(model_eau)

boxplot(siteCovs(umf)$Dist_Eau, main = "Dist_Eau") 

On a des gros outliers… Ce qui créer des soucis pour le model. On peux centrer reduire.

siteCovs(umf)$Dist_Eau_scaled <- scale(siteCovs(umf)$Dist_Eau)
boxplot(siteCovs(umf)$Dist_Eau_scaled, main = "Dist_Eau_scaled") 

model_eau <- occu(~1 ~ Dist_Eau_scaled, data = umf)
unmarked::summary(model_eau)
## 
## Call:
## occu(formula = ~1 ~ Dist_Eau_scaled, data = umf)
## 
## Occupancy (logit-scale):
##                 Estimate    SE      z  P(>|z|)
## (Intercept)       3.2755 0.295 11.091 1.38e-28
## Dist_Eau_scaled  -0.0747 0.239 -0.312 7.55e-01
## 
## Detection (logit-scale):
##  Estimate     SE     z   P(>|z|)
##    -0.541 0.0255 -21.2 3.81e-100
## 
## AIC: 8840.19 
## Number of sites: 328

Un avec une covariance de site et une d’observation.

model_eau_periode <- unmarked::occu(~periode ~Dist_Eau_scaled, data = umf)
unmarked::summary(model_eau_periode)
## 
## Call:
## unmarked::occu(formula = ~periode ~ Dist_Eau_scaled, data = umf)
## 
## Occupancy (logit-scale):
##                 Estimate    SE      z  P(>|z|)
## (Intercept)       3.2777 0.296 11.091 1.38e-28
## Dist_Eau_scaled  -0.0652 0.246 -0.265 7.91e-01
## 
## Detection (logit-scale):
##                 Estimate     SE     z   P(>|z|)
## (Intercept)       -0.924 0.0394 -23.4 2.87e-121
## periodeperiode2    0.699 0.0522  13.4  5.82e-41
## 
## AIC: 8658.317 
## Number of sites: 328

Comparaison des modeles uniquement base sur l’AIC (On peut pas comparer les deux parce que centrer reduire l’eau change les données de la base de donnée sur lequel on a fait les 2 modeles d’avant. A changer.)

unmarked::modSel(unmarked::fitList(aucune_variable = model_null, 
                                   deux_periodes = model_era))
##                 nPars     AIC  delta   AICwt cumltvWt
## deux_periodes       3 8656.40   0.00 1.0e+00     1.00
## aucune_variable     2 8838.28 181.88 3.2e-40     1.00
unmarked::modSel(unmarked::fitList(dp_et_eau = model_eau_periode,
                                   eau = model_eau))
##           nPars     AIC  delta   AICwt cumltvWt
## dp_et_eau     4 8658.32   0.00 1.0e+00     1.00
## eau           3 8840.19 181.87 3.2e-40     1.00

8 Tests d’hypotheses

#car::vif(lm(Dist_Eau ~ periode, data = toutesdonnees))

Sur 2 périodes, il y a une augementation des lapins. Peut être car en 2005 a 2010, on a peu de données. D’ou l’idée de créer 3 périodes.