Hi Cheng,

If you have 5 habitat types, I am guessing the proportion of each is closer to 0.2 and not very tiny.

Those large beta's mean the probability of occupancy at each site is either 0 (e.g., when shrubs are present) or 1 (e.g., when conifers are present). Typically, this results when you have a small data set and too many covariates. For example, if you only have 5 detections and they all occur in areas with conifers and without shrubs, then you will get a huge beta for conifers and a very negative value for shrubs. How large is your data set? How many sites, how many surveys, how many detections? With a small data set, it is more common to get a large SE and you have a very small SE, so that is a little peculiar.

Also, if there are only 5 land use types (i.e., sub+con+shrub+plan+const=1), you should not use more than 4 in the model because the intercept accounts for the other land use type. If you have an intercept, plus every possible land use type, you will get bad estimates.

Also, if plantation and construction are also percent of land use, then your interaction terms would be very hard to interpret and another possible source of the problem.

I would say the most likely issue is fitting a complex model to a small data set, but if you want a more definitive answer, you will need to share more information about your data and models. Good luck!

Matthew J. Clement,

When I have an almost perfect detection probability (as I built below), with occupancy following clearly a gradient effect (as I built below), why the confidence interval (in occupancy prediction under the gradient) is huge (0-1) when it should be tiny. I would love, if you can share some thoughts about that.

library(unmarked)

m=matrix(rep(c(0,0,1,1),c(12,12,12,12)),4,12,byrow=T)

m[3:4,12]=0

gradient=data.frame(gradient=c(100,400,700,1000))

frame=unmarkedFrameOccu(m,gradient)

null=occu(~1 ~1,frame)

grad=occu(~1 ~gradient,frame)

f=fitList(null=null,gradient=grad)

modSel(f)

predict(grad, type="state")

grad Predicted SE lower upper

100 1.930998e-10 1.633219e-08 1.958904e-82 1

400 7.034032e-04 2.735281e-02 5.297396e-37 1

700 9.996104e-01 1.660390e-02 1.307839e-33 1

1000 1.000000e+00 9.601140e-09 3.395093e-67 1

Luiz,

Your main problem is a small sample size. In your gradient model, you are estimating two occupancy parameters from four sites. You cannot expect precise estimates from such a tiny data set.

A secondary problem is that you have some significant data separation. When the gradient is below 700, it appears the sites are never occupied. When the gradient is above 400, it appears the sites are always occupied. That creates an ambiguous result when the gradient is between 400 and 700 because those sites could be always occupied, never occupied, or something in between. That is called "data separation" and it makes it difficult to estimate the relationship between the gradient and occupancy.

But like I said, your main problem is small sample size. Try your modeling again with more data.

Thank you so much!

A last question... perceive above that SE is tiny but 95%CI is 0-1.

For example

Predicted: 0.99

SE: 0.02

CI: 0 to 1

Supposedly, for a Gaussian distribution, mean+- 1.96*SE should results in CI bounds. It does not happen here. Why?

Hi Luiz,

The distribution is Gaussian on the logit scale, not the probability scale.

In occupancy modeling, the estimation is done on the logit scale, which means when you tell R "~gradient", it does NOT model psi~gradient (where psi is probability of occupancy). Instead, it models log(psi/1-psi)~gradient, which is called the logit scale.

It is easier to look at your null model.

> predict(null, type="state")

Predicted SE lower upper

1 0.5 0.25 0.1234709 0.8765291

2 0.5 0.25 0.1234709 0.8765291

3 0.5 0.25 0.1234709 0.8765291

4 0.5 0.25 0.1234709 0.8765291

You cannot get the CI from 0.5 +/- 1.96*.25 = 0.01 and 0.99.

Let's switch to the logit scale.

> summary(null)

Call: occu(formula = ~1 ~ 1, data = frame)

Occupancy (logit-scale): Estimate SE z P(>|z|)

1.83e-14 1 1.83e-14 1

Detection (logit-scale): Estimate SE z P(>|z|)

2.4 0.739 3.25 0.00117

AIC: 23.3133

First we use Gaussian assumptions to get the CI on the logit scale,

1.83e-14 +/- 1.96*1.0 = -1.96 and 1.96

Then we transform to probability scale

plogis(-1.96) = 0.123479

plogis(1.96) = 0.876533

That gives the correct result.

For your gradient model, the math gets more complicated, but you can see on the logit scale that the intercept estimate is very imprecise (SE=101).

> summary(grad)

Call: occu(formula = ~1 ~ gradient, data = frame)

Occupancy (logit-scale): Estimate SE z P(>|z|)

(Intercept) -27.4041 101.755 -0.269 0.788

gradient 0.0504 0.182 0.277 0.782

Detection (logit-scale): Estimate SE z P(>|z|)

2.4 0.739 3.25 0.00117

AIC: 19.77031

As I said, the math gets more complicated, but for the gradient model, you should use the CI from 0 to 1, and not the tiny SE estimate it gives you.