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.