Projected occupancy and colonisation/extinction values in stan colext

[Vishnu Menon]

Hi all,

This might be a simple question but wanted to get some opinions on interpreting projected occupancy and col/ext values. So I'm running dynamic occupancy models with two seasons (prefire and postfire) across two different treatments (baited and unbaited) for brushtail possums. I'm using ubms so as to include a random effect 'region' on my multiseason model. The top model is:

model <- stan_colext (~ baiting + (1|region), ~baiting, ~1, ~ time + baiting, umf)

I use projected occupancy to estimate changes in occupancy values for time vs treatment and the resulting graph looks like this:

I then use predict to estimate effects of baiting on colonisation probability using predict for the top model:

newdat= data.frame(baiting=c("unbaited","baited"))
col_value <-predict(model, submodel="col", newdata=newdat, appendData=TRUE)

And my colonisation estimates are 0.21 (CI: 0.11 - 0.32) in baited compared to 0.13 (CI: 0.09 - 0.18) in unbaited.

Why do I see a slight increase in my projected occupancy for brushtail possums in unbaited areas post-fire even though the colonisation probability is lower in unbaited? Is it because of using projected vs predict? Is there a way to try and explain this or am I doing this incorrectly? Any advice would be appreciated. Thank you so much.

Regards,

Vishnu

[Ken Kellner]

I'm not sure exactly how you calculated the projections here (is it for a single site or across-sites mean?).

One possible explanation: projected occupancy as calculated by ubms and unmarked is a function of not just colonization but also initial occupancy and extinction : 

psi[t+1]  = psi[t]*(1 - eps) + (1-psi[t])*gamma [equation 7 in MacKenzie et al. 2003]

So all three of these together are interacting and driving the projected values, meaning that it's not as simple as just comparing colonization by itself.

If occupancy at psi[t] is relatively large, then many sites are occupied, and persistence (1-eps) will be more important driver of how many sites are occupied in time t+1 than colonization (because there are few sites that it's even possible to colonize). If  occupancy at time t is relatively small, then few sites are occupied, and colonization becomes the more important driver of sites occupied in time t+1 (because many sites could be colonized). 

As a result it turns out that the projected occupancy will trend towards a stable value, which is (if I remember correctly) calculated as gamma / (gamma + eps). This is assuming gamma and eps are constant over time. If your initial occupancy is close to this stable value, then the projected occupancy will shift a small amount towards that value in each time step. If the initial value is far away, it'll move more quickly towards the stable value. So in your example I am guessing that the initial occupancy for the baited sites is closer to the value gamma_baited / (gamma_baited + eps), relative to the equivalent for the non-baited sites. As a result projected occupancy for the non-baited sites is changing more quickly (towards the stable value) even though non-baited colonization probability is lower. It depends on your estimated extinction probability, though, which you didn't include.

Here's a quick example in R:

# initial occupancy
psi0_bait <- 0.4
psi0_nobait <- 0.2

# extinction
eps <- 0.28

# colonization
gam_bait <- 0.2
gam_nobait <- 0.13 # lower

T <- 100

# Calculate projected
psi_bait <- rep(NA, T)
psi_bait[1] <- psi0_bait

psi_nobait <- rep(NA, T)
psi_nobait[1] <- psi0_nobait

for (t in 2:T){
  psi_bait[t] <-psi_bait[t-1] * (1-eps) + (1-psi_bait[t-1])*gam_bait
  psi_nobait[t] <-psi_nobait[t-1] * (1-eps) + (1-psi_nobait[t-1])*gam_nobait
}

# differences from times 1 to 2 (greater for no bait despite lower colonization)
psi_bait[2] - psi_bait[1] # 0.008
psi_nobait[2] - psi_nobait[1] # 0.048

# asymptotes: gam / (gam + eps)
asym_bait <- gam_bait / (gam_bait + eps)
asym_nobait <- gam_nobait / (gam_nobait + eps)

# steeper slope for nobait
plot(1:2, psi_bait[1:2], type='o', ylim=c(0.15, 0.45))
lines(1:2, psi_nobait[1:2], type='o', col='red')
legend("bottomright", col=c("black", "red"), lty=1, legend=c("bait", "nobait"))
abline(h=asym_bait, col='black', lty=2)
abline(h=asym_nobait, col='red', lty=2)

# long term: stable value reached
plot(1:T, psi_bait, type='o', ylim=c(0.15, 0.45), ylab="occupancy")
lines(1:T, psi_nobait, type='o', col='red')
legend("bottomright", col=c("black", "red"), lty=1, legend=c("bait", "nobait"))
abline(h=asym_bait, col='black', lty=2)
abline(h=asym_nobait, col='red', lty=2)

[Vishnu Menon]

Hi Ken,

That makes perfect sense, thank you. And the projections I've used here is mean across sites.

Thanks,

Vishnu