#Generate some data
set.seed(123)
lambda=4; gamma=0.5; omega=0.8; p=0.5
M <- 100; T <- 5
y <- array(NA, c(M, 3, T))
N <- matrix(NA, M, T)
S <- G <- matrix(NA, M, T-1)
for(i in 1:M) {
N[i,1] <- rpois(1, lambda)
y[i,1,1] <- rbinom(1, N[i,1], p) # Observe some
Nleft1 <- N[i,1] - y[i,1,1] # Remove them
y[i,2,1] <- rbinom(1, Nleft1, p) # ...
Nleft2 <- Nleft1 - y[i,2,1]
y[i,3,1] <- rbinom(1, Nleft2, p)
for(t in 1:(T-1)) {
S[i,t] <- rbinom(1, N[i,t], omega)
G[i,t] <- rpois(1, gamma)
N[i,t+1] <- S[i,t] + G[i,t]
y[i,1,t+1] <- rbinom(1, N[i,t+1], p) # Observe some
Nleft1 <- N[i,t+1] - y[i,1,t+1] # Remove them
y[i,2,t+1] <- rbinom(1, Nleft1, p) # ...
Nleft2 <- Nleft1 - y[i,2,t+1]
y[i,3,t+1] <- rbinom(1, Nleft2, p)
}
}
y=matrix(y, M)
#Create some random covariate data
sc <- data.frame(x1=rnorm(100))
if (FALSE) {
#Create unmarked frame
umf <- unmarkedFrameMMO(y=y, numPrimary=5, siteCovs=sc, type="removal")
#Fit model
(fit <- multmixOpen(~x1, ~1, ~1, ~1, K=30, data=umf))
#Compare to truth
cf <- coef(fit)
data.frame(model=c(exp(cf[1]), cf[2], exp(cf[3]), plogis(cf[4]), plogis(cf[5])),
truth=c(lambda, 0, gamma, omega, p))
#Predict
head(predict(fit, type='lambda'))
#Check fit with parametric bootstrap
pb <- parboot(fit, nsims=15)
plot(pb)
# Empirical Bayes estimates of abundance for each site / year
re <- ranef(fit)
plot(re, layout=c(10,5), xlim=c(-1, 10))
}
Run the code above in your browser using DataLab