Obiectivele acestui exercițiu sunt următoarele:
1) să vă introducă în simularea datelor utilizate în modelarea abundenței
2) să configureze datele de abundență pentru utilizare în pachetul R “unmarked”
3) să ruleze modele N-mixture simple (single species single–season) de abundență și să evalueze rezultatele
1. Simulăm cel mai simplu model N-Mixture binomial
Cel mai bun mod de a înțelege modul în care funcționează modelarea, în general, este să simulați propriile date, apoi să rulați modele predictive și să evaluați rezultatele. Deoarece știm la ce să ne așteptăm de la seturile de date simulate acest exercițiu reprezintă o modalitate excelentă atât pentru înțelegerea structurii datelor, cât și pentru evaluarea rezultatelor modelelor complexe.
Pentru început, instalați pechetul R pe care îl vom folosi pentru implementarea modelelor de ocupanță
install.packages(unmarked)Odată ce ați instalat pachetele, nu este nevoie să le mai instalați dacă reluați exercițiul
Selectăm dimensiunile eșantionului și pregătiți matricea de date
library(unmarked)
set.seed(24) # Setăm seed-ul pentru reproducibilitate
M <- 100 # Selectăm numărul de situri
J <- 3 # Numărul de măsurători de abundență pe site (rep.counts)
C <- matrix(NA, nrow = M, ncol = J) # Matrice care va conține datele observațiilorGenerăm date locale de abundență (reale)
lambda <- 4 # Abundență așteptată
p <- 0.6 # Probabilitatea de detecțieSelectăm valorile parametrilor
N <- rpois(n = M, lambda = lambda)
# Generăm numărul real de indivizi (N) la fiecare sit folosind o distribuție Poisson cu parametru lambdaEfectuăm măsurători repetate (replicate)
# Generăm datele de detecție pentru fiecare măsurătoare la fiecare sit
for(j in 1:J) {
C[, j] <- rbinom(n = M, size = N, prob = p)
}Explorăm datele
table(N) # Distribuția reală a abundenței
N # Verificăm datele
sum(N) # Mărimea totală reală a populației la M site-uri
sum(N>0) # Numărul real de situri ocupate
mean(N) # Abundența medie reală (estimarea lambda)
table(apply(C, 1, max)) # Distribuția abundenței observată (număr maxim)
sum(apply(C, 1, max)) [1] # Dimensiunea totală a populației observată în M situri,
# Număr observat de situri ocupate
sum(apply(C, 1, max)>0) [1] # Media observată „abundență relativă” abundance"
mean(apply(C, 1, max)) [1] # Media observată „abundență relativă” abundance"
head(cbind(N=N, count1=C[,1], count2=C[,2])) # Selectăm și verificăm primele 6 situri
cor(C)[1,2]Observăm acordul cu valoarea p
# Creăm un cadru de date unmarkedFramePCount utilizând matricea de date de abundență (C)
umf <- unmarkedFramePCount(y = C) # Creăm unui cadru de date
summary(umf) # Rezumăm datele nou create
(fm1 <- pcount(~1 ~1, data = umf)) # # Definim un model de abundență folosind funcția pcount
# Modelul nu include covariate (modelul null), doar interceptul (~1 ~1)# Obținem estimările parametrilor la scară naturală (abundența și detectabilitatea)
backTransform(fm1, "state") # Estimarea abundenței
backTransform(fm1, "det") # Estimarea probabilității de detecțieExaminăm relația
# Generăm de detecție pentru fiecare vizită la fiecare sit
for(j in 1:J) {
C[,j] <- rbinom(M, N, p[,j]) }Plotăm datele observate și efectul vântului
# Plotăm valorile vântului (wind) în funcție de datele observate de detecție
plot(wind, C/max(C), xlab="Wind", ylab="Scaled counts: C/max(C)", frame = F, cex = 1.5)
# Adăugăm curba logistică care arată relația dintre vânt și probabilitatea de detecție
lines(seq(-1,1,,100), plogis(alpha0 + alpha1*seq(-1,1,,100)), lwd=3, col="red")Abundența așteptată (lambda) și cea realizată (N) și măsurătorile (C)
# Examinăm abundența așteptată (lambda) și cea realizată (N) și măsurătorile repetate (C)
cbind(lambda=round(lambda,2), N=N, C1=C[,1], C2=C[,2], C3=C[,3])
time <- matrix(rep(as.character(1:J), M), ncol = J, byrow = TRUE)
hab <- c(rep("A", 50), rep("B", 50), rep("C", 50))
umf <- unmarkedFramePCount( y = C, siteCovs = data.frame(vegHt = vegHt, hab = hab), obsCovs = list(time = time, wind = wind))
summary(umf)Codul R este disponibil pentru descărcare: