I am using a binomial mixture model to estimate abundance (N) and
detection probability (p) using simulated count data:
-Each site has a simulated abundance that follow a Poisson
distribution with lambda = 5
-There are 200 simulated sampled sites
-3 repeated counts at each site
- only 50 percent of the animals are counted during each count (i.e,
detection probability p =0.5, see codes)
We removed sites in which animals were never counted (see matrix y, in
the script)
I would like to use a zero truncated version of the Poisson
distribution (I am aware of zero-inflated binomial mixture models, but
still want to solve the problem described above).
The codes below:
(1) generate a matrix of counts (y), rows correspond to sites and
column to repeat visits at each sites. The script also removes sites
when no animals were counted during the 3 consecutive visits
(2) The second part of the script calls WinBUGS and run the binomial
mixture models on the count data. In this case the count matrix y was
converted to a vector C1 before being passed over to BUGS
Any idea how to create a zero truncated Poisson for parameter lam1
(i.e., parameter lambda of the Poisson distribution)
Thank you for your help.
#R script
#Simulated abundance data
n.site <- 200# 200 sites visited
lam <- 5 #mean number of animals per site
R <- n.site # nubmer of sites
T <- 3 # Number of replicate counts at each site
N=rpois(n = n.site, lambda = lam) #True abundance
# Simulate count data; only a fraction of N is counted which results in y
y <- array(dim = c(R, T))
for(i in 1:T){
y[,i] <- rbinom(n = n.site, size = N, prob = 0.5)
}
#truncate y matrix
y # R-by-T matrix of counts
sumy=apply(y,1,sum)
cbindysumy=cbind(y,sumy)
subsetcbindysumy=subset(cbindysumy,sumy!=0)
y=subsetcbindysumy[,1:3]# sites where no animals ever counted are removed
C1<-c(y) #vectorized matrix y
R=dim(y)[1]
site = 1:R
site.p <- rep(site, T)
#################################################################################################################################
#WinBUGS codes
#################################################################################################################################
library("R2WinBUGS") # Load R2WinBUGS package
sink("Model.txt")
cat("
model {
# Priors: new uniform priors
p0~dunif(0,1)
lam1~dgamma(.01,.01)
# Likelihood
# Biological model for true abundance
for (i in 1:R) { # Loops over R sites
N1[i] ~ dpois(lambda1[i])
lambda1[i] <- lam1
}
# Observation model for replicated counts
for (i in 1:n) { # Loops over all n observations
C1[i] ~ dbin(p1[i], N1[site.p[i]])
p1[i] <-p0
}
# Derived quantities
totalN1 <- sum(N1[]) # Estimate total population size across all sites
}
",fill=TRUE)
sink()
# Package data for WinBUGS
R = dim(y)[1] # number of sites: 200
n = dim(y)[1] * dim(y)[2]#number of observations (sites*surveys)
win.data <- list(R = R, n = n, C1 = C1, site.p = site.p)
# Inits
Nst <- apply(y, 1, max) + 1
inits <- function(){list(N1 = Nst,lam1=runif(1,1,8), p0=runif(1))}
parameters <- c( "totalN1", "p0", "lam1")
# MCMC settings
nc <- 3
nb <- 400 #need to push to 400 for convergence
ni <- 1400 #need to push to 14000 for convergence
nt <- 1
out <- bugs (win.data, inits, parameters, "Model.txt", n.chains=nc,
n.iter=ni, n.burn = nb, n.thin=nt, debug = T)