Predicting species interaction term in occumulti

159 views
Skip to first unread message

Divyajyoti Ganguly

unread,
Aug 19, 2022, 4:56:34 PM8/19/22
to unmarked
Hi everyone,

As part of my study I am trying to look at patterns of co-occurrence between species pairs of mesocarnivores. 

For this purpose, I have used the Rota et al. 2016 model on unmarked using the occumulti() function. The reason behind choosing this model is that it estimates species co-occurrence as a separate parameter and allows for checking the effects of covariates on co-occurrence. 

Now I want to plot the co-occurrence parameters for each species pair against relevant covariates to be able to visualize the change along a gradient. I tried using the predict() function to achieve this, but it is throwing up an error: "Supplied species name not found". The same function works while predicting the occupancy for individual species. 

How can I achieve this? Any help would be much appreciated.

Best,
Divyajyoti Ganguly,
National Centre for Biological Sciences, Bangalore,
India 

Ken Kellner

unread,
Aug 20, 2022, 8:40:32 AM8/20/22
to unmarked
Hi,

Can you show the summary output of your model and the code you are trying to use to predict?

Ken

Divyajyoti Ganguly

unread,
Aug 20, 2022, 9:05:15 AM8/20/22
to unma...@googlegroups.com
Thank you for your response, Dr. Ken. Please find the code and the output below:
#Dog Interactions
y6OF_3 <- c('~OS_Z+AG_Z','~OS_Z','~Dog_Z') # OS and AG: habitat types proportion within grid-cell; Dog: free-ranging dog capture rates
yDF <- c('~TrepZ:OC','~TrepZ:OC') # Trep: effort; OC: detection type

> y6_3 <- occuMulti(yDF,y6OF_3,data6, se = T, silent = T, penalty = 1)
Bootstraping covariance matrix
> #Look at output
> summary(y6_3)

Call:
occuMulti(detformulas = yDF, stateformulas = y6OF_3, data = data6,
    penalty = 1, se = T, silent = T, maxOrder = 2L)

Occupancy (logit-scale):
                      Estimate    SE     z  P(>|z|)
[sp1] (Intercept)       -0.607 0.398 -1.52 1.27e-01
[sp1] OS_Z               0.594 0.380  1.56 1.18e-01
[sp1] AG_Z              -0.412 0.245 -1.68 9.23e-02
[sp2] (Intercept)       -1.599 0.342 -4.67 2.99e-06
[sp2] OS_Z               0.524 0.182  2.88 4.00e-03
[sp1:sp2] (Intercept)    0.554 0.459  1.21 2.27e-01
[sp1:sp2] Dog_Z         -0.488 0.339 -1.44 1.49e-01

Detection (logit-scale):
                  Estimate     SE       z  P(>|z|)
[sp1] (Intercept) -3.58432 0.2259 -15.866 1.10e-56
[sp1] TrepZ:OCCT   1.14998 0.2370   4.851 1.23e-06
[sp1] TrepZ:OCSS   0.18588 0.0971   1.913 5.57e-02
[sp2] (Intercept) -2.21978 0.2682  -8.277 1.27e-16
[sp2] TrepZ:OCCT   0.66255 0.3093   2.142 3.22e-02
[sp2] TrepZ:OCSS  -0.00302 0.1890  -0.016 9.87e-01

AIC: 701.6538
Number of sites: 163
optim convergence code: 0
optim iterations: 48
Bootstrap iterations: 30

 range(SC$Dog_Z)
[1] -1.581462  1.158784
> seqD <- seq(-1.581462, 1.158784, length.out=100)
> nd <- data.frame(Dog_Z = seqD)
> predict(y6_3, type = "state", species = "sp1:sp2", newdata = nd)
Error in name_to_ind(species, names(object@data@ylist)) :
  Supplied species name not found

Best,
Divyajyoti

--
*** 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 a topic in the Google Groups "unmarked" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/unmarked/2N0G4I60dMY/unsubscribe.
To unsubscribe from this group and all its topics, send an email to unmarked+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/unmarked/8b989296-d2ed-4abb-a72d-f3d0df2ace24n%40googlegroups.com.


--
Divyajyoti Ganguly (he/him),
PG program in Wildlife Biology and Conservation,

Ken Kellner

unread,
Aug 20, 2022, 3:55:54 PM8/20/22
to unmarked
If you want the probability that both sp1 and sp2 occupy a site over a range of the covariate value, you can get that from predict with a slight change to your code:

predict(y6_3, type = "state", species = c("sp1", "sp2"), newdata = nd)

If you want to plot the value of the actual interaction term/natural parameter for a given range of covariate values (f[12] eg as in equation 1 of Rota et al. 2016), predict can't do that, and you will have to calculate it manually using the estimated [sp1:sp2] intercept and [sp1:sp2] slope. You can extract those values using the coef() function on the model output.

Ken

Divyajyoti Ganguly

unread,
Aug 21, 2022, 9:19:26 AM8/21/22
to unma...@googlegroups.com
Yes, I want to predict the natural parameter (f[12]). This makes sense, thanks a lot for your help.

Best,
Divyajyoti

Reply all
Reply to author
Forward
0 new messages