Bayesian occupancy -> abundance problem (WinBUGS)

[James Russell]

Hi everyone,

I am facing a problem I can't quite put my finger on. In fact I think its a BUGS problem but again not so sure, and given expertise in this forum think maybe someone could help out. Here goes:

I'm trying to scale from an occupancy (0/1) to abundance (0,1,...) model. Here is some stripped down code and data to higlight the problem.

The occupancy code works, the individual (change prior on N and z[i] distribution) does not. Increasing the prior on N (dunif(0,100) say) does work, presumably because it removes zeros from the binomial distribution. My problem is why would these zeroes be a problem in the abundance but not species case?

I would be indebted for a solution (even if I've missed something stupidly obvious as is often my case).

------< warning code below!! >-------

SPECIES MODEL & DATA

model {

p~dunif(0,1)
N~dunif(0,1)

for (i in 1:nsite) {
z[i] ~ dbern(N)
for(k in 1:K) {
x[i,k] ~ dbin(p,z[i])
} } }
}

list(nsite=4, K=4,
x = structure(.Data = c(1,0,1,1,1,1,0,1,1,1,0,1,1,1,0,1), .Dim=c(4,4)))

INDIVIDUAL MODEL & DATA

model {

p~dunif(0,1)
N~dunif(0,10)

for (i in 1:nsite) {
z[i] ~ dpois(N)
for(k in 1:K) {
x[i,k] ~ dbin(p,z[i])
} } }
}

list(nsite=4, K=4,
x = structure(.Data = c(1,0,1,5,1,1,0,4,1,2,0,1,2,1,0,1), .Dim=c(4,4)))

[sbonner]

Hi James,

What's the error message you get in WinBUGS? In OpenBUGS I get an error that the slice sampler updater of z[1] ran out of iterations.

I think that the problem may be that the model, as written, is not identifiable. You have z[i]~Pois(N) and x[i,k]|z[i] ~ Bin(z[i],p), which implies that the unconditional distribution of x[i,k] is x[i,k] ~ Pois(Np). Because the only data given are the x[i,k]'s there is no information to differentiate between N and p. Given your uniform priors, it is equally likely that x[i,k]=1 because N=9 and p=1/9 or because N=1 and p=1.

Essentially, there is no information to distinguish between the expected number of individuals present and the capture probability. Few individuals may be observed because there aren't many or because the capture probability is low.

I think that BUGS is running into difficulty because it's trying to find a better value for z[i] and there isn't one; all possible z[i]'s are equally probable. My understanding of the slice sampler is that it generates many proposals for the value and then chooses the best of these. I've encountered similar problems in the past, and my guess is that because the model is not identifiable the sampler can't find a better value and eventually quits.

Unfortunately, I'm not an expert on abundance estimation, but my guess is that you need to provide more information about either N or p in some way.

Hope that helps!

Cheers...

[sbonner]

Forgot to mention, if I set N to a fixed value (N <- 100 for example) then it all seems to work well.

Cheers...

[darryl]

Hi James,
Simon's right, basically you need to have a link between detection of the species (p) and abundance. One way of doing that is to redefine p as:

1-(1-r)^z[i] where r is now the probability of detecting an individual of the species, ie 1-p = Pr(don't detect any of the z[i] individuals at the site) = (1-r)^z[i]

Of course there's a few more assumptions in there like individuals are detected independently.

This is model presented by Royle and Nichols (2003) Ecology.

My personal position though is that trying to estimate abundance from just presence/absence data is potentially dangerous and that you have to look very carefully at all the assumptions and make sure they are being met (or at least not too badly violated). However, if you're just using this types of methods to account for heterogeneity in p to get a better estimate of occupancy (fraction of places where N>=1) rather than abundance per se, violations of assumptions might not so critical.

Cheers
Darryl

[James Russell]

Thanks Simon and Darryl,

Indeed I was having a non-identifiability issue due to a lapse in thought process. We have data on the number of (un-marked) animals at each survey, which I was hoping to exploit, but of course scaling from species to individual detectability would require a CMR framework on the individuals to determine Pr(individual detectability). The abundance distribution alone which I currently have provides no extra information since within that I can't seperate detection and density.

Perhaps we could bridge to the MARK forum and work on a supermodel for CMR for individuals and resurveys, etc (that's a field work issue though, very data hungry).

Here is the correct code:

model {

r ~dunif(0,1)
N~dunif(0,10)

for (i in 1:nsite) {
z[i] ~ dpois(N)
p[i] <- 1-pow((1-r),z[i]) # density model goes to P(detect) species
for(k in 1:K) {
x[i,k] ~ dbern(p[i])
} } }
}

list(nsite=4, K=4,
x = structure(.Data = c(1,0,1,1,1,1,0,1,1,1,0,1,1,1,0,1), .Dim=c(4,4)))