Site Occupancy Lapin
08/08
1 Intallation des packages
library(dplyr);library(ggplot2);library(sf);library(unmarked)
theme_set(theme_light(base_size = 16))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 84Total <- 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 317 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -4.798138 ymin: 47.35903 xmax: -1.023689 ymax: 48.86557
## Geodetic CRS: WGS 84VariablesSite <- 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(year = lubridate::year(date))3 Creation de la base de donnee Detection - Non Detection
detection_df <- 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))
detection_df[,-1] <- ifelse(detection_df[,-1]>0, 1, 0)
noms <- detection_df$Code_10km
detection_matrice <- as.matrix(detection_df[, -1])
rownames(detection_matrice) <- noms4 Covariables de site
site_info <- VariablesSite %>%
sf::st_drop_geometry() %>%
dplyr::arrange(Code_10km)
noms <- site_info$Code_10km
rownames(site_info) <- noms
site_info <- site_info %>%
dplyr::select(-Code_10km) %>%
dplyr::rename(Dnst_Cultures = `Dnst_Cultures `)%>%
dplyr::mutate(across(where(is.numeric),
scale))site_info %>%
dplyr::select(-air_km2) %>%
dplyr::select_if(is.numeric) %>%
GGally::ggcorr(label = TRUE, hjust = 0.75)
names(site_info)## [1] "Ind_Diversite" "Dnst_Cultures" "Dist_EcotoneArbore" "Dist_Littoral" "Dist_Eau"
## [6] "Famille_paysage" "Nom_paysage" "air_km2" "X_10km" "Y_10km"df <- site_info%>%
rename(X10km = X_10km,
Y10km = Y_10km) %>%
dplyr:: select(-air_km2) %>%
dplyr::select_if(is.numeric)
df <- as.data.frame(lapply(df, as.numeric))GGally::ggpairs(as.data.frame(lapply(df, as.numeric)))
On voit dans ces graphiques une forte correlation lineaire entre la densite de cultures et l’indice de diversite. On remarque aussi une correlation lineaire entre X et la distance au littoral. Graphiquement, on ne voit pas forcement de correlation particuliers de correlation dans les nuage de points, apart X et Y dans contre la distance au littoral, mais c’est pas forcement un probleme.
Comme ils le font dans l’etude avec les ours, on va se se retirer les covariables avec une correlation > 0.5 et regarder les diff combinaisons.
5 Implementation du Model
Mise sous format pour la modelisation
umf <- unmarked::unmarkedFrameOccu(y = detection_matrice,
obsCovs = list(periode = periodes_df),
siteCovs = site_info
)unmarked::plot(umf, main="Detection/Non Detection de lapins sur les sites depuis 2005")
On voit bien qu’on a tres peu de sites avant 2005, non pas parce qu’il n’y avais pas de presence de lapins, mais plus tres peu de detection. On a tres peu de donnees dans l’ensemble, sur toutes les especes.
On retirera ensuite les annees [2005 - 2009]. J’ai longtemps hesite entre garder 2009, ou 2025. Les 2 ont peu de donnees par rapport aux reste (2025 ca je n’ai pas l’annee entiere, et 2009 car c’est tôt). Les deux participerait a augementer la probabilite de detection.
Cependant, pour garder quelque chose de plus jolie avec des saisons de 5 ans, j’ai garde 2025.
Lapins <- Total %>%
sf::st_join(VariablesSite) %>%
filter(cd_nom == 61714,
technique_observation == "Vu",
date >= as.Date("2010-01-01"),
date < as.Date("2025-01-01")) %>%
# select(-c(insee_dept, lib_dept, lib_dept, FID, surf, CD_SIG)) %>%
select(-nom_valide, -nom_vernaculaire, -ordre, -technique_observation) %>%
mutate(year = lubridate::year(date))
detection_df <- 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))
detection_df[,-1] <- ifelse(detection_df[,-1]>0, 1, 0)
noms <- detection_df$Code_10km
detection_matrice <- as.matrix(detection_df[, -1])
rownames(detection_matrice) <- noms
umf <- unmarked::unmarkedFrameOccu(y = detection_matrice,
siteCovs = site_info
)
unmarked::plot(umf, main="Detection/Non Detection de lapins sur les sites de 2010 a 2024 inclus")
summary(umf)## unmarkedFrame Object
##
## 317 sites
## Maximum number of observations per site: 15
## Mean number of observations per site: 15
## Sites with at least one detection: 307
##
## Tabulation of y observations:
## 0 1
## 2577 2178
##
## Site-level covariates:
## Ind_Diversite.V1 Dnst_Cultures.V1 Dist_EcotoneArbore.V1
## Min. :-2.865708445210000 Min. :-2.78214762389e+00 Min. :-1.177548563760
## 1st Qu.:-0.679780380756000 1st Qu.:-6.39147940409e-01 1st Qu.:-0.538861888256
## Median :-0.103681984768000 Median : 6.52166412662e-02 Median :-0.234385970959
## Mean : 0.000000000000001 Mean : 2.00000000000e-16 Mean : 0.000000000000
## 3rd Qu.: 0.590728352478000 3rd Qu.: 6.59758341030e-01 3rd Qu.: 0.149328232852
## Max. : 3.387273586330000 Max. : 2.36743641082e+00 Max. : 7.826110849440
##
## Dist_Littoral.V1 Dist_Eau.V1 Famille_paysage
## Min. :-1.070027865720 Min. :-2.486983274260 Paysage cultive a ragosses :90
## 1st Qu.:-0.915859784196 1st Qu.:-0.653550370660 Paysage boise et de bosquets :57
## Median :-0.266683896501 Median :-0.230018369783 Paysage cultive avec talus :52
## Mean : 0.000000000000 Mean : 0.000000000000 Paysage de bocage a maille elargie :41
## 3rd Qu.: 0.699371185739 3rd Qu.: 0.584178367602 Paysage de bocage dense sur collines:34
## Max. : 2.957617628370 Max. : 6.943084217390 Paysage de littoral urbanise :27
## (Other) :16
## Nom_paysage air_km2.V1 X_10km.V1
## Reliefs des Landes de Lanvaux: 36 Min. :-3.199398376490 Min. :-1.91668816575e+00
## Plateau de Penthievre : 20 1st Qu.: 0.131669308234 1st Qu.:-7.91009617370e-01
## Bassin de Pontivy-Loudeac : 18 Median : 0.522356248052 Median : 1.43352452421e-03
## Armor morbihannais : 17 Mean : 0.000000000000 Mean : 1.00000000000e-16
## Bassin de Rennes : 17 3rd Qu.: 0.522356248055 3rd Qu.: 8.11993642375e-01
## Cornouaille interieure : 16 Max. : 0.522356248061 Max. : 1.82993231695e+00
## (Other) :193
## Y_10km.V1
## Min. :-2.36529227433e+00
## 1st Qu.:-7.80490736624e-01
## Median : 1.54633653399e-02
## Mean : 1.07000000000e-14
## 3rd Qu.: 8.34068493645e-01
## Max. : 1.96182815611e+00
## 6 Resultats
6.1 Modeles d’occupation a une seule saison
Un model de base sans variables
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.43 0.322 10.6 1.92e-26
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## -0.108 0.0295 -3.67 0.000242
##
## AIC: 6463.181
## Number of sites: 317backTransform(model_null, type="state")## Backtransformed linear combination(s) of Occupancy estimate(s)
##
## Estimate SE LinComb (Intercept)
## 0.969 0.00982 3.43 1
##
## Transformation: logisticbackTransform(model_null, type="det")## Backtransformed linear combination(s) of Detection estimate(s)
##
## Estimate SE LinComb (Intercept)
## 0.473 0.00736 -0.108 1
##
## Transformation: logisticconfint(backTransform(model_null, type="state"))## 0.025 0.975
## 0.9424205 0.9830019confint(backTransform(model_null, type="det"))## 0.025 0.975
## 0.4585328 0.487377coef(model_null)## psi(Int) p(Int)
## 3.4263979 -0.1083766En transformant, on a donc - La probabilite qu’un site est occupe par un ou plusieurs lapins est de 0,97. (En considerent toute la periode 2005-2025) - La probabilite qu’on detecte un lapin est de 0,37
On a egalement les intervalles de confiances qui sont assez serres autour de ces valeurs a 95%.
model_eau <- unmarked::occu(~1 ~ Dist_Eau, data = umf)
unmarked::summary(model_eau)##
## Call:
## unmarked::occu(formula = ~1 ~ Dist_Eau, data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 3.435 0.325 10.585 3.5e-26
## Dist_Eau -0.133 0.285 -0.468 6.4e-01
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## -0.108 0.0295 -3.67 0.000242
##
## AIC: 6464.98
## Number of sites: 317noms_site_info <- c("Ind_Diversite", "Dnst_Cultures", "Dist_EcotoneArbore", "Dist_Littoral",
"Dist_Eau", "air_km2", "X_10km", "Y_10km")
for (i in noms_site_info) {
formule <- as.formula(paste("~1 ~", i))
model <- unmarked::occu(formule, data = umf)
grid <- seq(min(site_info[[i]], na.rm=TRUE),
max(site_info[[i]], na.rm=TRUE),
length=50)
occ_pred <- plogis(coef(model)[1] + coef(model)[2] * grid)
# Plot
plot(grid,
occ_pred,
type = "l",
lwd = 3,
xlab = i,
ylab = "Probabilite d'occupation estimee")
}
On voit la relation entre nos variables de sites et la probabilite d’occupancy.
Comparaison de plusieurs models
model_null <- unmarked::occu(~1 ~ 1, data = umf)
model_eau <- unmarked::occu(~1 ~ Dist_Eau, data = umf)
model_litt <- unmarked::occu(~1 ~ Dist_Littoral, data = umf)
model_litt_eau <- unmarked::occu(~1 ~ Dist_Littoral + Dist_Eau,
data = umf)
fmList <- unmarked::fitList('{psi, p}' = model_null,
'{psi(dist_eau), p}' = model_eau,
'{psi(dist_littoral), p}' = model_litt,
'{psi(dist_littoral + dist_eau), p}' = model_litt_eau)
unmarked::modSel(fmList)## nPars AIC delta AICwt cumltvWt
## {psi(dist_littoral), p} 3 6459.21 0.00 0.576 0.58
## {psi(dist_littoral + dist_eau), p} 4 6460.43 1.22 0.313 0.89
## {psi, p} 2 6463.18 3.97 0.079 0.97
## {psi(dist_eau), p} 3 6464.98 5.77 0.032 1.00D’apres Occupancy Models to study wildlife USGS, il est interessant de regarder les models avec \(\Delta AIC < 2\). C’est a dire, ici on regarderais les 2 premiers models. Mais on va plutot regarder les models suivants comme c’est ce qui nous interesse.
6.2 Modeles d’occupation a plusieurs saisons
Au sein d’une saison, les lapins sont soit présent soit absent mais ne changent pas. Entre les saisons, les lapins peuvent coloniser un site (ie: il étiaent pas là avant et mtn ils le sont), ou disparaitre = extinction locale. C’est une chaine de Marcov.
Ici on suppose, qu’il n’y a pas de faux positif.
On a
psi \(\psi\) : prob occupation
p: prob detection
gamma \(\gamma\) : prob de colonisation
epsilon \(\epsilon\) : prob d’extinction
On separe les observations en 3 periodes de 4 ans :
2010 - 2014,
2015 - 2019,
2020 - 2024.
annees <- as.numeric(colnames(detection_matrice))
quali_periodes_vector <- dplyr::case_when(annees < 2015 ~ "periode1",
annees < 2020 ~ "periode2",
annees < 2025 ~ "periode3")
quali_periodes_matrice <- matrix(rep(quali_periodes_vector,
each = nrow(detection_df)),
nrow = nrow(detection_df),
byrow = FALSE)
colnames(quali_periodes_matrice) <- colnames(detection_matrice)
rownames(quali_periodes_matrice) <- rownames(detection_matrice)
quali_periodes_df <- as.data.frame(quali_periodes_matrice)quanti_periodes_vector <- dplyr::case_when(annees < 2015 ~ 1,
annees < 2020 ~ 2,
annees < 2025 ~ 3)
quanti_periodes_matrice <- matrix(rep(quanti_periodes_vector,
each = nrow(detection_df)),
nrow = nrow(detection_df),
byrow = FALSE)
colnames(quanti_periodes_matrice) <- colnames(detection_matrice)
rownames(quanti_periodes_matrice) <- rownames(detection_matrice)
quanti_periodes_df <- as.data.frame(quanti_periodes_matrice)site_covs_periodes <- list(
quanti_periodes = quanti_periodes_matrice,
quali_periodes = quali_periodes_matrice
)On a les periodes sous la forme qualitiative et quantitative, pour tester les 2 façons.
detections_per_year = colSums(detection_matrice)
print(detections_per_year)## 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023 2024
## 95 165 179 158 166 123 140 134 147 143 158 153 152 143 122On a bcp de detections sur toutes les années en vrai..
umf <- unmarked::unmarkedMultFrame(y = detection_matrice,
siteCovs = site_info,
yearlySiteCovs = site_covs_periodes,
numPrimary = 15)
summary(umf)## unmarkedFrame Object
##
## 317 sites
## Maximum number of observations per site: 15
## Mean number of observations per site: 15
## Number of primary survey periods: 15
## Number of secondary survey periods: 1
## Sites with at least one detection: 307
##
## Tabulation of y observations:
## 0 1
## 2577 2178
##
## Site-level covariates:
## Ind_Diversite.V1 Dnst_Cultures.V1 Dist_EcotoneArbore.V1
## Min. :-2.865708445210000 Min. :-2.78214762389e+00 Min. :-1.177548563760
## 1st Qu.:-0.679780380756000 1st Qu.:-6.39147940409e-01 1st Qu.:-0.538861888256
## Median :-0.103681984768000 Median : 6.52166412662e-02 Median :-0.234385970959
## Mean : 0.000000000000001 Mean : 2.00000000000e-16 Mean : 0.000000000000
## 3rd Qu.: 0.590728352478000 3rd Qu.: 6.59758341030e-01 3rd Qu.: 0.149328232852
## Max. : 3.387273586330000 Max. : 2.36743641082e+00 Max. : 7.826110849440
##
## Dist_Littoral.V1 Dist_Eau.V1 Famille_paysage
## Min. :-1.070027865720 Min. :-2.486983274260 Paysage cultive a ragosses :90
## 1st Qu.:-0.915859784196 1st Qu.:-0.653550370660 Paysage boise et de bosquets :57
## Median :-0.266683896501 Median :-0.230018369783 Paysage cultive avec talus :52
## Mean : 0.000000000000 Mean : 0.000000000000 Paysage de bocage a maille elargie :41
## 3rd Qu.: 0.699371185739 3rd Qu.: 0.584178367602 Paysage de bocage dense sur collines:34
## Max. : 2.957617628370 Max. : 6.943084217390 Paysage de littoral urbanise :27
## (Other) :16
## Nom_paysage air_km2.V1 X_10km.V1
## Reliefs des Landes de Lanvaux: 36 Min. :-3.199398376490 Min. :-1.91668816575e+00
## Plateau de Penthievre : 20 1st Qu.: 0.131669308234 1st Qu.:-7.91009617370e-01
## Bassin de Pontivy-Loudeac : 18 Median : 0.522356248052 Median : 1.43352452421e-03
## Armor morbihannais : 17 Mean : 0.000000000000 Mean : 1.00000000000e-16
## Bassin de Rennes : 17 3rd Qu.: 0.522356248055 3rd Qu.: 8.11993642375e-01
## Cornouaille interieure : 16 Max. : 0.522356248061 Max. : 1.82993231695e+00
## (Other) :193
## Y_10km.V1
## Min. :-2.36529227433e+00
## 1st Qu.:-7.80490736624e-01
## Median : 1.54633653399e-02
## Mean : 1.07000000000e-14
## 3rd Qu.: 8.34068493645e-01
## Max. : 1.96182815611e+00
##
##
## Yearly-site-level covariates:
## quanti_periodes quali_periodes
## Min. :1 periode1:1585
## 1st Qu.:1 periode2:1585
## Median :2 periode3:1585
## Mean :2
## 3rd Qu.:3
## Max. :3fm0 <- colext(
psiformula = ~ 1, # initial occupancy
gammaformula = ~ 1, # colonization
epsilonformula = ~ 1, # extinction
pformula = ~ 1, # detection
data = umf, # data
control = list(trace = 1)
) ## initial value 3760.705819
## iter 10 value 3003.431775
## iter 10 value 3003.431771
## final value 3003.431771
## convergedsummary(fm0)##
## Call:
## colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~1,
## pformula = ~1, data = umf, control = list(trace = 1))
##
## Initial (logit-scale):
## Estimate SE z P(>|z|)
## -0.251 0.176 -1.43 0.154
##
## Colonization (logit-scale):
## Estimate SE z P(>|z|)
## -0.97 0.0782 -12.4 2.55e-35
##
## Extinction (logit-scale):
## Estimate SE z P(>|z|)
## -1.56 0.136 -11.5 2.19e-30
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## 1.23 0.133 9.22 3.05e-20
##
## AIC: 6014.864
## Number of sites: 317colonisation
unmarked::backTransform(fm0, type="col")## Backtransformed linear combination(s) of Colonization estimate(s)
##
## Estimate SE LinComb (Intercept)
## 0.275 0.0156 -0.97 1
##
## Transformation: logisticIl a ici une probabilite de \(0,27\) de coloniser un site ou le lapin n’est pas present.
unmarked::backTransform(fm0, type="ext")## Backtransformed linear combination(s) of Extinction estimate(s)
##
## Estimate SE LinComb (Intercept)
## 0.174 0.0195 -1.56 1
##
## Transformation: logisticOna une probabilite d’extinction de \(0,17\), donc comme \(0.17<0.27\), les lapins devraient être sur de plus en plus du territoire.
unmarked::backTransform(fm0, type="psi")## Backtransformed linear combination(s) of Initial estimate(s)
##
## Estimate SE LinComb (Intercept)
## 0.438 0.0433 -0.251 1
##
## Transformation: logisticunmarked::backTransform(fm0, type="det")## Backtransformed linear combination(s) of Detection estimate(s)
##
## Estimate SE LinComb (Intercept)
## 0.774 0.0234 1.23 1
##
## Transformation: logisticLa probabilite d’etre present sur un territoire initialement, donc avant toute colonisation/extinction, est de \(0,44\) et celle de detection est \(0,77\).
6.3 Periodes qualitatives.
On essaye plusieurs sans la variance, puis selon l’AIC, on calcule l’AIC a la fin avec le modele qu’on a choisi.
fm0 <- unmarked::colext(
psiformula = ~ 1, # initial occupancy
gammaformula = ~ 1, # colonization
epsilonformula = ~ quali_periodes, # extinction
pformula = ~ 1, # detection
data = umf, # data
control = list(trace = 1),
se = FALSE)## initial value 3760.705819
## iter 10 value 3001.871608
## final value 3001.802875
## convergedunmarked::backTransform(fm0, type="col")## Backtransformed linear combination(s) of Colonization estimate(s)
##
## Estimate SE LinComb (Intercept)
## 0.272 NA -0.983 1
##
## Transformation: logisticcoef_ext <- coef(fm0, type = "ext")
plogis(coef_ext[1])## ext(Int)
## 0.1641366plogis(coef_ext[2])## ext(quali_periodesperiode2)
## 0.484936plogis(coef_ext[3])## ext(quali_periodesperiode3)
## 0.5618445unmarked::backTransform(fm0, type="psi")## Backtransformed linear combination(s) of Initial estimate(s)
##
## Estimate SE LinComb (Intercept)
## 0.442 NA -0.234 1
##
## Transformation: logisticunmarked::backTransform(fm0, type="det")## Backtransformed linear combination(s) of Detection estimate(s)
##
## Estimate SE LinComb (Intercept)
## 0.771 NA 1.21 1
##
## Transformation: logisticEncore a travailler la compréhension du truc…
- Probabilité de
summary(fm0)##
## Call:
## unmarked::colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~quali_periodes,
## pformula = ~1, data = umf, se = FALSE, control = list(trace = 1))
##
## Initial (logit-scale):
## Estimate SE z P(>|z|)
## -0.234 NA NA NA
##
## Colonization (logit-scale):
## Estimate SE z P(>|z|)
## -0.983 NA NA NA
##
## Extinction (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -1.6278 NA NA NA
## quali_periodesperiode2 -0.0603 NA NA NA
## quali_periodesperiode3 0.2487 NA NA NA
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## 1.21 NA NA NA
##
## AIC: 6015.606
## Number of sites: 317coef(fm0)## psi(Int) col(Int) ext(Int)
## -0.23418956 -0.98327873 -1.62776644
## ext(quali_periodesperiode2) ext(quali_periodesperiode3) p(Int)
## -0.06027438 0.24865144 1.21341901Maintenant, on ajoute les variables de sites, et on regarde pour l’instant uniquement l’AIC, afin de faire un tri rapide. On fait une elimination forwards manuelle.
fm0 <- unmarked::colext(
psiformula = ~ 1, # initial occupancy
gammaformula = ~ 1, # colonization
epsilonformula = ~ quali_periodes, # extinction
pformula = ~ 1, # detection
data = umf, # data
control = list(trace = 1),
se = FALSE)## initial value 3760.705819
## iter 10 value 3001.871608
## final value 3001.802875
## convergedfm1 <- unmarked::colext(
psiformula = ~ Y_10km, # initial occupancy
gammaformula = ~ 1, # colonization
epsilonformula = ~ quali_periodes, # extinction
pformula = ~ 1, # detection
data = umf, # data
control = list(trace = 1),
se = FALSE)## initial value 3760.705819
## iter 10 value 3002.865691
## final value 3001.800739
## convergedfm2 <- unmarked::colext(
psiformula = ~ X_10km, # initial occupancy
gammaformula = ~ 1, # colonization
epsilonformula = ~ quali_periodes, # extinction
pformula = ~ 1, # detection
data = umf, # data
control = list(trace = 1),
se = FALSE)## initial value 3760.705819
## iter 10 value 2997.882784
## final value 2997.313157
## convergedfm3 <- unmarked::colext(
psiformula = ~ Dist_Eau, # initial occupancy
gammaformula = ~ 1, # colonization
epsilonformula = ~ quali_periodes, # extinction
pformula = ~ 1, # detection
data = umf, # data
control = list(trace = 1),
se = FALSE)## initial value 3760.705819
## iter 10 value 3002.411072
## final value 3001.717944
## convergedfm4 <- unmarked::colext(
psiformula = ~ Dist_Littoral, # initial occupancy
gammaformula = ~ 1, # colonization
epsilonformula = ~ quali_periodes, # extinction
pformula = ~ 1, # detection
data = umf, # data
control = list(trace = 1),
se = FALSE)## initial value 3760.705819
## iter 10 value 2996.766703
## final value 2996.389556
## convergedfm5 <- unmarked::colext(
psiformula = ~ Dist_EcotoneArbore, # initial occupancy
gammaformula = ~ 1, # colonization
epsilonformula = ~ quali_periodes, # extinction
pformula = ~ 1, # detection
data = umf, # data
control = list(trace = 1),
se = FALSE)## initial value 3760.705819
## iter 10 value 3001.515376
## final value 3001.130867
## convergedfm6 <- unmarked::colext(
psiformula = ~ Dnst_Cultures, # initial occupancy
gammaformula = ~ 1, # colonization
epsilonformula = ~ quali_periodes, # extinction
pformula = ~ 1, # detection
data = umf, # data
control = list(trace = 1),
se = FALSE)## initial value 3760.705819
## iter 10 value 2995.347057
## final value 2995.043137
## convergedfm7 <- unmarked::colext(
psiformula = ~ Ind_Diversite, # initial occupancy
gammaformula = ~ 1, # colonization
epsilonformula = ~ quali_periodes, # extinction
pformula = ~ 1, # detection
data = umf, # data
control = list(trace = 1),
se = FALSE)## initial value 3760.705819
## iter 10 value 2996.033979
## final value 2995.771111
## convergedfmList <- unmarked::fitList(fm1, fm3, fm4, fm5, fm6, fm7)## Your list was unnamed, so model names were added automaticallyunmarked::modSel(fmList)## Hessian is singular.
## Hessian is singular.
## Hessian is singular.
## Hessian is singular.
## Hessian is singular.
## Hessian is singular.## nPars AIC delta AICwt cumltvWt
## fm6 7 6004.09 0.00 0.57216 0.57
## fm7 7 6005.54 1.46 0.27629 0.85
## fm4 7 6006.78 2.69 0.14886 1.00
## fm5 7 6016.26 12.18 0.00130 1.00
## fm3 7 6017.44 13.35 0.00072 1.00
## fm1 7 6017.60 13.52 0.00066 1.00Ici on garde la variable de fm6 ($ \text{Dnst_Cultures} $), et on ajoute les autres.
fm0 <- fm6
fm1 <- unmarked::colext(
psiformula = ~ Dnst_Cultures + Y_10km, # initial occupancy
gammaformula = ~ 1, # colonization
epsilonformula = ~ quali_periodes, # extinction
pformula = ~ 1, # detection
data = umf, # data
control = list(trace = 1),
se = FALSE)## initial value 3760.705819
## iter 10 value 2996.010385
## final value 2994.353904
## convergedfm2 <- unmarked::colext(
psiformula = ~ Dnst_Cultures + X_10km, # initial occupancy
gammaformula = ~ 1, # colonization
epsilonformula = ~ quali_periodes, # extinction
pformula = ~ 1, # detection
data = umf, # data
control = list(trace = 1),
se = FALSE)## initial value 3760.705819
## iter 10 value 2994.082308
## final value 2993.236118
## convergedfm3 <- unmarked::colext(
psiformula = ~ Dnst_Cultures + Dist_Eau, # initial occupancy
gammaformula = ~ 1, # colonization
epsilonformula = ~ quali_periodes, # extinction
pformula = ~ 1, # detection
data = umf, # data
control = list(trace = 1),
se = FALSE)## initial value 3760.705819
## iter 10 value 2998.136949
## final value 2994.983323
## convergedfm4 <- unmarked::colext(
psiformula = ~ Dnst_Cultures + Dist_Littoral, # initial occupancy
gammaformula = ~ 1, # colonization
epsilonformula = ~ quali_periodes, # extinction
pformula = ~ 1, # detection
data = umf, # data
control = list(trace = 1),
se = FALSE)## initial value 3760.705819
## iter 10 value 2993.915555
## final value 2992.347208
## convergedfm5 <- unmarked::colext(
psiformula = ~ Dnst_Cultures + Dist_EcotoneArbore, # initial occupancy
gammaformula = ~ 1, # colonization
epsilonformula = ~ quali_periodes, # extinction
pformula = ~ 1, # detection
data = umf, # data
control = list(trace = 1),
se = FALSE)## initial value 3760.705819
## iter 10 value 2995.573273
## final value 2993.765644
## convergedfm6 <- unmarked::colext(
psiformula = ~ Dnst_Cultures, # initial occupancy
gammaformula = ~ 1, # colonization
epsilonformula = ~ quali_periodes, # extinction
pformula = ~ 1, # detection
data = umf, # data
control = list(trace = 1),
se = FALSE)## initial value 3760.705819
## iter 10 value 2995.347057
## final value 2995.043137
## convergedfmList <- unmarked::fitList(fm0, fm1, fm3, fm4, fm5, fm6)## Your list was unnamed, so model names were added automaticallyunmarked::modSel(fmList)## Hessian is singular.
## Hessian is singular.
## Hessian is singular.
## Hessian is singular.
## Hessian is singular.
## Hessian is singular.## nPars AIC delta AICwt cumltvWt
## fm4 8 6000.69 0.00 0.551 0.55
## fm5 8 6003.53 2.84 0.133 0.68
## fm0 7 6004.09 3.39 0.101 0.79
## fm6 7 6004.09 3.39 0.101 0.89
## fm1 8 6004.71 4.01 0.074 0.96
## fm3 8 6005.97 5.27 0.039 1.00On continue ensuite avec \(\text{Dnst_Cultures} + \text{Dist_Littoral}\).
fm0 <- fm4
fm1 <- unmarked::colext(
psiformula = ~ Dnst_Cultures + Dist_Littoral + Y_10km, # initial occupancy
gammaformula = ~ 1, # colonization
epsilonformula = ~ quali_periodes, # extinction
pformula = ~ 1, # detection
data = umf, # data
control = list(trace = 1),
se = FALSE)## initial value 3760.705819
## iter 10 value 2994.338042
## final value 2992.345760
## convergedfm2 <- unmarked::colext(
psiformula = ~ Dnst_Cultures + Dist_Littoral + X_10km, # initial occupancy
gammaformula = ~ 1, # colonization
epsilonformula = ~ quali_periodes, # extinction
pformula = ~ 1, # detection
data = umf, # data
control = list(trace = 1),
se = FALSE)## initial value 3760.705819
## iter 10 value 2994.758898
## final value 2992.119507
## convergedfm3 <- unmarked::colext(
psiformula = ~ Dnst_Cultures + Dist_Littoral + Dist_Eau, # initial occupancy
gammaformula = ~ 1, # colonization
epsilonformula = ~ quali_periodes, # extinction
pformula = ~ 1, # detection
data = umf, # data
control = list(trace = 1),
se = FALSE)## initial value 3760.705819
## iter 10 value 2993.494454
## final value 2992.217760
## convergedfm5 <- unmarked::colext(
psiformula = ~ Dnst_Cultures + Dist_Littoral + Dist_EcotoneArbore, # initial occupancy
gammaformula = ~ 1, # colonization
epsilonformula = ~ quali_periodes, # extinction
pformula = ~ 1, # detection
data = umf, # data
control = list(trace = 1),
se = FALSE)## initial value 3760.705819
## iter 10 value 2992.206353
## final value 2990.443357
## convergedfmList <- unmarked::fitList(fm0, fm1, fm3, fm4, fm5)## Your list was unnamed, so model names were added automaticallyunmarked::modSel(fmList)## Hessian is singular.
## Hessian is singular.
## Hessian is singular.
## Hessian is singular.
## Hessian is singular.## nPars AIC delta AICwt cumltvWt
## fm5 9 5998.89 0.00 0.47 0.47
## fm0 8 6000.69 1.81 0.19 0.66
## fm4 8 6000.69 1.81 0.19 0.85
## fm3 9 6002.44 3.55 0.08 0.93
## fm1 9 6002.69 3.80 0.07 1.006.3.1 Periodes quantitatif.
On fait de meme avec les variables en quantitatif.
fm0 <- unmarked::colext(
psiformula = ~ 1, # initial occupancy
gammaformula = ~ 1, # colonization
epsilonformula = ~ quanti_periodes, # extinction
pformula = ~ 1, # detection
data = umf, # data
control = list(trace = 1),
se = FALSE)## initial value 3760.705819
## iter 10 value 3151.902215
## iter 20 value 3070.685628
## iter 30 value 3002.492023
## final value 3002.474185
## convergedsummary(fm0)##
## Call:
## unmarked::colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~quanti_periodes,
## pformula = ~1, data = umf, se = FALSE, control = list(trace = 1))
##
## Initial (logit-scale):
## Estimate SE z P(>|z|)
## -0.232 NA NA NA
##
## Colonization (logit-scale):
## Estimate SE z P(>|z|)
## -0.98 NA NA NA
##
## Extinction (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -1.816 NA NA NA
## quanti_periodes 0.125 NA NA NA
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## 1.21 NA NA NA
##
## AIC: 6014.948
## Number of sites: 317unmarked::backTransform(fm0, type="col")## Backtransformed linear combination(s) of Colonization estimate(s)
##
## Estimate SE LinComb (Intercept)
## 0.273 NA -0.98 1
##
## Transformation: logisticcoef_ext <- coef(fm0, type = "ext")
plogis(coef_ext[1])## ext(Int)
## 0.1399136plogis(coef_ext[2])## ext(quanti_periodes)
## 0.5311488plogis(coef_ext[3])## <NA>
## NAunmarked::backTransform(fm0, type="psi")## Backtransformed linear combination(s) of Initial estimate(s)
##
## Estimate SE LinComb (Intercept)
## 0.442 NA -0.232 1
##
## Transformation: logisticunmarked::backTransform(fm0, type="det")## Backtransformed linear combination(s) of Detection estimate(s)
##
## Estimate SE LinComb (Intercept)
## 0.77 NA 1.21 1
##
## Transformation: logistic6.4 Modeles d’occupation avec des fausses positives
Pas super pertinant je penses pour nous. Chaque detection, souvent est lie a plusieurs detections distinctes de multiples personnes. De pus, les donnees inserees sont des donnees de naturalites et sont verifies. Etant donnee que nous avons garde les donnes bien, ca nous concerne pas forcement, mais aurait pu etre utile avec les donnees diro.
7 Tests d’hypotheses
Hypothèses:
- occupancy status at each unit does not change over the survey season; i.e., units are ???closed??? to changes in occupancy -> probablement pas mais bon
- the probability of occupancy is constant across units, or differences in occu pancy probability are modeled using covariates ->
- the probability of detection is constant across all units and surveys or is a function of unit-survey covariates; there is no unmodeled heterogeneity in detection probabilities ->
- detection of species and detection histories at each location are independent -> non
- species are correctly identified; i.e., no false positives -> okay