Predicting species interaction term in occumulti

Skip to first unread message

Divyajyoti Ganguly

Aug 19, 2022, 4:56:34 PMAug 19
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.

Divyajyoti Ganguly,
National Centre for Biological Sciences, Bangalore,

Ken Kellner

Aug 20, 2022, 8:40:32 AMAug 20
to unmarked

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


Divyajyoti Ganguly

Aug 20, 2022, 9:05:15 AMAug 20
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)

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

[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


*** 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
To unsubscribe from this group and all its topics, send an email to
To view this discussion on the web visit

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

Ken Kellner

Aug 20, 2022, 3:55:54 PMAug 20
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.


Divyajyoti Ganguly

Aug 21, 2022, 9:19:26 AMAug 21
Yes, I want to predict the natural parameter (f[12]). This makes sense, thanks a lot for your help.


Reply all
Reply to author
0 new messages