Simulare abundență

install.packages(unmarked)

library(unmarked)
set.seed(24) # Setăm seed-ul pentru reproducibilitate
# Selectăm dimensiunile eșantionului și pregătim matricea de date
M <- 150 # 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 conține datele observațiilor

vegHt <- sort(runif(M, -1, 1)) # Sortați datele pentru a fi mai ușor de vizualizat

beta0 <- 0 # Interceptare la scară logaritmică
beta1 <- 2 # Panta pe scara logaritmică pentru vegHt
# Calculăm abundența așteptată (lambda) pentru fiecare sit pe baza modelului log-linear
lambda <- exp(beta0 + beta1 * vegHt) # Abundență așteptată
# Plotăm relația dintre înălțimea vegetației și abundența așteptată (lambda)
plot(vegHt, lambda, type = "l", lwd = 3) # Abundență așteptată

# Generăm abundența locală (N) pentru fiecare sit utilizând o distribuție Poisson cu parametru lambda
N <- rpois(M, lambda) 
points(vegHt, N) # Adăugăm abundența realizată în grafic table(N)

plot(vegHt, N, xlab="Înălțimea vegetației", ylab="Abundență reală (N)", frame = F, cex = 1.5) 
lines(seq(-1,1,,100), exp(beta0 + beta1* seq(-1,1,,100)), lwd=3, col = "red")

# Generăm o matrice de dimensiuni M x J care conține valori aleatoare uniform distribuite între -1 și 1
wind <- array(runif(M * J, -1, 1), dim = c(M, J))

alpha0 <- -2 # Interceptul pe scara logit
alpha1 <- -3 # Panta pe scara logit pentru vânt
# Calculăm probabilitatea de detecție pentru fiecare sit și vizită utilizând funcția logistică
p <- plogis(alpha0 + alpha1 * wind)
# Plotăm relația dintre probabilitatea de detecție (p) și vânt (wind)
plot(p ~ wind, ylim = c(0,1))

# Generăm datele de detecție pentru fiecare măsurătoare (vizită) la fiecare sit
for(j in 1:J) {
  C[,j] <- rbinom(M, N, p[,j]) }

# Plotăm valorile vântului (wind) în funcție de datele observate de detecție scalate 
plot(wind, C/max(C), xlab="Wind", ylab="Scaled counts: C/max(C)", frame = F, cex = 1.5)
# Adăugăm la plot 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")

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)
# Creăm o variabilă categorică hab care indică habitatul pentru fiecare sit
hab <- c(rep("A", 50), rep("B", 50), rep("C", 50))
# Creăm un nou cadru de date unmarkedFramePCount cu covariate specifice sitului și observației
umf <- unmarkedFramePCount( y = C, siteCovs = data.frame(vegHt = vegHt, hab = hab), obsCovs = list(time = time, wind = wind))
# Rezumăm noul set de date pentru a vedea structura și un sumar statistic al acestora 
summary(umf)

summary(fm.Nmix1 <- pcount(~wind ~vegHt, data=umf, control=list(trace=T, REPORT=1))) 
fm.Nmix2 <- pcount(~wind ~vegHt, data=umf, mixture="NB", control=list(trace=TRUE, REPORT=5)) 
fm.Nmix3 <- pcount(~wind ~vegHt, data=umf, mixture="ZIP", control=list(trace=TRUE, REPORT=5))
cbind(AIC.P=fm.Nmix1@AIC, AIC.NB=fm.Nmix2@AIC, AIC.ZIP=fm.Nmix3@AIC)

# Obținem predicțiile pentru abundența (lambda) pentru fiecare sit
(lambda.hat <- predict(fm.Nmix1, type="state")) # lambda pentru fiecare sit 
# Obținem predicțiile pentru probabilitatea de detecție (p) pentru fiecare sit și vizită
(p.hat <- predict(fm.Nmix1, type="det")) # p în timpul fiecărui sondaj

# Utilizăm datele de detecție (C) combinate pe toate cele 3 vizite (replicate)
summary(fm.glm <- glm(c(C) ~ rep(vegHt, 3), family=poisson)) # p-model naiv
# Plotăm datele observate de detecție în funcție de înălțimea vegetației
matplot(vegHt, C, xlab="Vegetation height", ylab="Counts", frame = F, cex = 1.5, pch = 1, col = "black")
# Adăugăm la plot curba care arată relația reală dintre înălțimea vegetației și abundența (lambda)
lines(seq(-1,1,,100), exp(beta0 + beta1* seq(-1,1,,100)), lwd=3, col = "red")
# Adăugăm la plot curba estimată de modelul GLM Poisson naiv
curve(exp(coef(fm.glm)[1]+coef(fm.glm)[2]*x), -1, 1, type ="l", lwd=3, add=TRUE)
# Adăugăm la plot predicțiile modelului N-mixture ajustat 
lines(vegHt, predict(fm.Nmix1, type="state")[,1], col = "blue", lwd = 3)
legend(-1, 6, c("Truth", "'Poisson GLM' with p", "Poisson GLM without p"), col=c("red", "blue", "black"), lty = 1, lwd=3, cex = 1.2)