Interpreting OccuMulti results - logit scale conversion

[Emily Collins]

Hi all, I am wondering if anyone can advise. These are the results of my OccuMulti model:

Occupancy:
                                       Estimate    SE              z          P(>|z|)
[sp1] (Intercept)           -0.355        0.297     -1.196    2.32e-01
[sp2] (Intercept)           -2.101        0.517     -4.068    4.74e-05
[sp1:sp2] (Intercept)    0.185         0.817     0.227     8.21e-01

In reporting these results, I have been using 'plogis' to convert the intercept estimates out of the logit scale. 

For example, intercept estimate for sp1 

> plogis (-0.355)
[1] 0.4121705

which I am interpreting as probability of 0.4 that sp 1 will be present at a site.

Should I do this for the sp1:sp2 intercept as well? 

And am I correct in reporting that there was non significant positive interaction between sp1 and sp2 based on these results?

Thanks!

Emily

[Ken Kellner]

Hi Emily,

For the multispecies model you don't want to use plogis on the intercepts, because they are not on the logit scale but rather the multinomial logit. The intercepts represent estimates of the natural [f] parameters. See Rota et al. 2016 equations 1-2 for the math here. If you want marginal occupancy probability for a given species from the model use the predict() function, e.g.

predict(model, type='state', species='sp1')

See ?occuMulti help for many examples of the use of predict.

And yes you are correct that the interaction term here is not significantly different from 0.

Ken

[Emily Collins]

Thanks for the response. I am confused because the output of my multi-species model reads: Occupancy (logit-scale).

Call:
occuMulti(detformulas = detFormulas, stateformulas = occFormulas,
    data = umf.multi2)

Occupancy (logit-scale):

                      Estimate    SE      z  P(>|z|)
[sp1] (Intercept)       -0.355 0.297 -1.196 2.32e-01
[sp2] (Intercept)       -2.101 0.517 -4.068 4.74e-05
[sp1:sp2] (Intercept)    0.185 0.817  0.227 8.21e-01

In comparison, when I use a single species model I get:

Occupancy (logit-scale):

 Estimate    SE     z P(>|z|)

   -0.332 0.281 -1.18   0.237

# to convert:

> predict (mod.ch1, type="state")
    Predicted         SE     lower    upper
1   0.4176463 0.06843316 0.2923682 0.554538

> plogis(-0.332)
[1] 0.4177541

versus the marginal occupancy for sp1 only in the multispecies model:

Occupancy (logit-scale):

                      Estimate    SE      z  P(>|z|)
[sp1] (Intercept)       -0.355 0.297 -1.196 2.32e-01

# to convert:

> predict(fit, type="state", species="sp1")
Bootstrapping confidence intervals with 100 samples
    Predicted         SE     lower     upper
1   0.4174412 0.06973618 0.2911559 0.5446485

> plogis(-0.355)
[1] 0.4121705

**this was similar enough to the predict() value in both single and multispecies models that I thought plogis was working here, and could be used to convert the [sp1:sp2]  intercept estimate as well.

Basically I just want to make sure that I am reporting the intercept for the [sp1:sp2] interaction term on the right scale (even though it is not significantly different from 0). 

Should the intercept term for occupancy of sp1 and sp2 just be reported as 0.185, or should this be converted somehow, as the intercepts for single species were converted using predict(). 

[sp1:sp2] (Intercept)    0.185 0.817  0.227 8.21e-01

I tried using predict(fit, type="state", species="sp1:sp2") just to see, but it didn't work.

Thanks so much,

Emily

--
*** Three hierarchical modeling email lists ***
(1) unmarked (this list): for questions specific to the R package unmarked
(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology: for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2021) and Schaub & Kéry (2022)
---
You received this message because you are subscribed to the Google Groups "unmarked" group.
To unsubscribe from this group and stop receiving emails from it, send an email to unmarked+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/unmarked/92398e7a-902f-4d71-a63d-8496a1d3dfd6n%40googlegroups.com.

[Ken Kellner]

Short answers:

The summary really should say multinomial logit, not logit.

The calculation of occupancy probability for a single species model is related to, but not identical to the calculation for the multispecies model. In your case, species 2 is very rare, so the species 1 intercepts end up looking similar between the single- and multispecies models. However despite this you cannot isolate and transform any individual parameter from the multispecies model, you have to consider them all together when calculating probabilities of occupancy (see below).

Thus the interaction term should be reported on the original scale (0.185). It doesn't make sense to transform it.

Ken

-----

Suppose the intercept for the single species model is b, the intercepts for the multi species models are a1 and a2 for species 1 and 2 respectively, and the interaction term is a12.

The probability of occupancy psi[1] for based on the single species model (i.e., plogis) is

psi = exp[b] / (1 + exp[b])

Based on the two-species model, the probability both species occupy a site psi[11] is

psi[11] = exp(a1 + a2 + a12) / (1 + exp(a1) + exp(a2) + exp(a1 + a2 + a12))

Notice if you remove the species 2 term and the interaction this simplifies back to the single-species formula.

The probability only species 1 occupies the site is

psi[10] = exp(a1) / (1 + exp(a1) + exp(a2)+ exp(a1 + a2 + a12))

The marginal probability that species 1 occupies a site is

psi[11] + psi[10]

So you can't isolate any one of the parameters in the multispecies model, since each occupancy state probability depends on all model parameters.

[Emily Collins]

Hi Ken, thanks so much for all of your help. This makes sense now.

Much appreciated,

Emily

To view this discussion on the web visit https://groups.google.com/d/msgid/unmarked/41ecbc3f-31e5-48f2-801e-ca6f764d3babn%40googlegroups.com.