library(dplyr);library(ggplot2);library(sf);library(unmarked)
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))
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
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
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
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
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
#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.