Simulare ocupanță

install.packages(unmarked)

library(unmarked)

set.seed(1) # Setăm seed-ul pentru reproducibilitate 
M <- 100 # Selectăm numărul de situri
J <- 3 # Numărul de măsurători de detecție/nondetecție
y <- matrix(NA, nrow = M, ncol = J) # Creăm matricea care să conțină datele observațiilor

vegHt <- sort(runif(M, -1, 1)) # Sortăm datele pentru a fi mai ușor de vizualizat

beta0 <- 0 # Interceptare la scară logaritmică
beta1 <- 3 # Panta pe scara logaritmică pentru vegHt
psi <- plogis(beta0 + beta1 * vegHt) # Probabilitatea ocupanței 
plot(vegHt, psi, ylim = c(0,1), type = "l", lwd = 3) # Plotăm relația psi 

z <- rbinom(M, 1, psi) # Generăm date de detecție/nondetecție reale
# Simulăm detecția/nondetecția reală pentru fiecare sit pe baza probabilității de ocupare (psi)

table(z) # Creăm un tabel de frecvență pentru a vedea câte situri sunt ocupate (1) și câte nu (0)
z # Afișăm vectorul z pentru a vedea detecția/nondetecția reală la fiecare sit

plot(vegHt, z, xlab="Vegetation height", ylab="True presence/absence (z)", frame = F, cex = 1.5) 
# Plotăm înălțimea vegetației (vegHt) în funcție de prezența/absența reală (z) pentru fiecare sit
plot(function(x) plogis(beta0 + beta1*x), -1, 1, add=T, lwd=3, col = "red")
# Adăugăm la plot curba logistică care arată relația dintre înălțimea vegetației și ocupanță

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

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

# Generăm datele de detecție/nondetecție pentru fiecare sit și vizită pe baza 
# probabilității de detecție (p) și a stării reale de ocupare (z)
for(j in 1:J) {
  y[,j] <- rbinom(M, z, p[,j]) }
    sum(apply(y, 1, max))

plot(wind, y, xlab="Wind", ylab="Observed det./nondetection data (y)", frame = F, cex = 1.5)
plot(function(x) plogis(alpha0 + alpha1*x), -1, 1, add=T, lwd=3, col = "red")

cbind(psi=round(psi,2), z=z, y1=y[,1], y2=y[,2], y3=y[,3])
# Combinăm într-un tabel valorile pentru probabilitatea de ocupare (psi), starea reală (z), și datele de detecție/nondetecție (y)

time <- matrix(rep(as.character(1:J), M), ncol = J, byrow = TRUE)
hab <- c(rep("A", 33), rep("B", 33), rep("C", 34)) # Trebuie să aveți M = 100

umf <- unmarkedFrameOccu(y = y, siteCovs = data.frame(vegHt = vegHt, hab = hab),obsCovs = list(wind = wind, time = time))
# covariate specifice sitului: vegHt = înălțimea vegetației
# covariate specifice observației: obsCovs = list(wind = wind, time = time)), wind = viteza vântului, time = data observației

summary(umf)

# Ajustăm modelul de ocupanță folosind covariatele de detecție (wind) și covariatele de ocupanță (vegHt)
summary(fm1.occ <- occu(~wind ~vegHt, data=umf)) # Covariatele de detecție urmează prima sedilă, apoi covariatele de ocupanță

# Creăm un nou set de date pentru predicții, variabila vegetație (vegHt) variind între -1 și 1, cu pasul de 0.01
newdat <- data.frame(vegHt=seq(-1, 1, 0.01))
# Obținem predicțiile pentru probabilitatea de ocupare (state) folosind noul set de date
pred.occ <- predict(fm1.occ, type="state", newdata=newdat)
# Creăm un alt set de date pentru predicții, de data aceasta pentru variabila vânt (wind)
newdat <- data.frame(wind=seq(-1, 1, 0.1))
# Obținem predicțiile pentru probabilitatea de detecție (det) folosind noul set de date
pred.det <- predict(fm1.occ, type="det", newdata=newdat)
# Creăm un set de date pentru predicții cu valori pentru wind 
newdat <- data.frame(wind=seq(-1, 1, , 5))
# Obținem predicțiile
predict(fm1.occ, type="det", newdata=newdat, append = T) 

# Ajustăm un GLM binomial folosind prezența observată și covariata vegHt
summary(fm.glm <- glm(apply(y, 1, max) ~ vegHt, family=binomial))
plot(vegHt, apply(y, 1, max), xlab="Vegetation height", ylab="Observed occurrence ('ever observed?')", frame = F, cex = 1.5)
plot(function(x) plogis(beta0 + beta1*x), -1, 1, add=T, lwd=3, col = "red") 
lines(vegHt, predict(fm.glm,,"response"), type = "l", lwd = 3)
lines(vegHt, predict(fm1.occ, type="state")[,1], col = "blue", lwd = 3)
legend(-1, 0.9, c("Truth", "'LR' with p", "LR without p"), col=c("red", "blue", "black"), lty = 1, lwd=3, cex = 1.2)