# Simulate data using the multinomial-Poisson model with a
# repeated constant-interval removal design.
n <- 100 # number of sites
T <- 4 # number of primary periods
J <- 3 # number of secondary periods
lam <- 3
phi <- 0.5
p <- 0.3
#set.seed(26)
y <- array(NA, c(n, T, J))
M <- rpois(n, lam) # Local population size
N <- matrix(NA, n, T) # Individuals available for detection
for(i in 1:n) {
N[i,] <- rbinom(T, M[i], phi)
y[i,,1] <- rbinom(T, N[i,], p) # Observe some
Nleft1 <- N[i,] - y[i,,1] # Remove them
y[i,,2] <- rbinom(T, Nleft1, p) # ...
Nleft2 <- Nleft1 - y[i,,2]
y[i,,3] <- rbinom(T, Nleft2, p)
}
y.ijt <- cbind(y[,1,], y[,2,], y[,3,], y[,4,])
umf1 <- unmarkedFrameGMM(y=y.ijt, numPrimary=T, type="removal")
(m1 <- gmultmix(~1, ~1, ~1, data=umf1, K=30))
backTransform(m1, type="lambda") # Individuals per plot
backTransform(m1, type="phi") # Probability of being avilable
(p <- backTransform(m1, type="det")) # Probability of detection
p <- coef(p)
# Multinomial cell probabilities under removal design
c(p, (1-p) * p, (1-p)^2 * p)
# Or more generally:
head(getP(m1))
# Empirical Bayes estimates of super-population size
re <- ranef(m1)
plot(re, layout=c(5,5), xlim=c(-1,20), subset=site%in%1:25)
Run the code above in your browser using DataLab