if (FALSE) {
# Simulate data
M = 1000 # Number of sites
nrep <- 3 # Number of visits per site
Tmax = 5 # Max duration of a visit
alpha1 = -1 # Covariate on rate
beta1 = 1 # Covariate on density
mu.lambda = 1 # Rate at alpha1 = 0
mu.dens = 1 # Density at beta1 = 0
covDet <- matrix(rnorm(M*nrep),nrow = M,ncol = nrep) #Detection covariate
covDens <- rnorm(M) #Abundance/density covariate
dens <- exp(log(mu.dens) + beta1 * covDens)
sum(N <- rpois(M, dens)) # Realized density per site
lambda <- exp(log(mu.lambda) + alpha1 * covDet) # per-individual detection rate
ttd <- NULL
for(i in 1:nrep) {
ttd <- cbind(ttd,rexp(M, N*lambda[,i])) # Simulate time to first detection per visit
}
ttd[N == 0,] <- 5 # Not observed where N = 0; ttd set to Tmax
ttd[ttd >= Tmax] <- 5 # Crop at Tmax
#Build unmarked frame
umf <- unmarkedFrameOccuTTD(y = ttd, surveyLength=5,
siteCovs = data.frame(covDens=covDens),
obsCovs = data.frame(covDet=as.vector(t(covDet))))
#Fit model
fit <- nmixTTD(~covDens, ~covDet, data=umf, K=max(N)+10)
#Compare to truth
cbind(coef(fit), c(log(mu.dens), beta1, log(mu.lambda), alpha1))
#Predict abundance/density values
head(predict(fit, type='state'))
}
Run the code above in your browser using DataLab