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)
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
