Trying to calculate the SIF with OccuMulti

110 views
Skip to first unread message

Gabriel Patricio Andrade Ponce

unread,
Jul 16, 2021, 8:21:59 PM7/16/21
to unma...@googlegroups.com
Hi all,
I am working with the multi-species models developed by Rota et al (2016) to investigate co-occurrence patterns of three species: the bobcat, gray fox and cottontail rabbit. I fitted several models with different hypotheses of species co-occurrence and chose the best models with the AICc. My best model indicates that there is a "spatial interaction" between two species pairs Grayfox vs Cottontail and Bobcat vs Cottontail.
My model looks like this
m8c <- occuMulti(c("~Vertcover_50", "~Effort", "~Cam"), c("~Neob", "~1", "~Neob",0 , "~Neob", "~1" ,0), data=multi)
Where Neob is a two-level factor: 1 for vegetation dominated by cacti and 0 for other vegetation.
I am interested in calculating the species interaction factor (SIF) or odds ratio v (MacKenzie et al 2018, p.530).  SIF=(ψ_11/(1-ψ_11))/(ψ_10/(1-ψ_10))
Since it is a derived parameter and to get the confidence intervals, it must be estimated by the delta method or by parboot. I am trying to do it by parboot using the following function:

SIF_par <- function(fm) {

  Psi11f <- plogis(coef(fm)[["psi([gray_fox:rabbit] (Intercept))"]])

  minPsi11f <- 1-Psi11f

  Psi10f <- plogis(coef(fm)[["psi([gray_fox] (Intercept))"]])

  minPsi10f <- 1-Psi10f

  Psi11b <- plogis(coef(fm)[["psi([bob_cat:rabbit] (Intercept))"]])

  minPsi11b <- 1-Psi11b

  Psi10b <- plogis(coef(fm)[["psi([bob_cat] (Intercept))"]])

  minPsi10b <- 1-Psi10b

    SIFf <- (Psi11f/minPsi11f)/(Psi10f/minPsi10f)

  SIFfb <- (Psi11b/minPsi11b)/(Psi10b/ minPsi10b)

    total <- c(siff=SIFf, sifb=SIFfb)

  return(total)

}

trySIF <- parboot(m8c, SIF_par, nsim= 1000)

SIFCI <- cbind(trySIF@t0,  t(apply(try...@t.star, 2, quantile, probs= c(0.025, 0.975))))

I manage to get the estimates and intervals, but the values are huge. For example, for the fox and rabbit the SIF= 462 with 95%CI = 9.56 --1170482927, and for the bobcat and rabbit it is similar but sadder SIF= 10.50108; 95%CI= 0.22 -- 12071948.
Is something wrong with my function? Should I try the delta method? Is it possible to calculate the SIF in relation to my categorical covariate? Or is the function OK and my data is definitely not very good?

Thanks in advance for the help

--
Gabriel P. Andrade-Ponce
Biólogo - Universidad Nacional de Colombia
MSc y Estudiante de Doctorado -Instituto de Ecología A.C.
Xalapa, Veracruz, México

Ken Kellner

unread,
Jul 17, 2021, 8:55:00 AM7/17/21
to unmarked
Hi Gabriel,

I don't have the MacKenzie book but I assume they are fitting a co-occupancy model like Richmond 2010 model? I don't know if the SIF for that type of model has the same meaning for a Rota-style multispecies model where there isn't a "dominant" species. So I'm not certain it's appropriate to use for inference here, particularly when you have > 2 species.

If you are sure you want to calculate something like this, one issue with your function is that you are assuming logit links for the occupancy probabilities when the Rota model uses a multinomial logit. The intercept and slope terms returned by occuMulti are used to calculate what Rota calls the 'natural parameters' or "f"s. All natural parameters are then used in the multinomial logit link function to calculate the probability of each possible occupancy state. See Rota 2016, eqn 2. To get the occupancy probabilities you can do the math yourself following eqn 2 or use the predict() function. For example, something like this should work:

nd <- data.frame(Neob=0) # or whatever other covariate levels you want
pr <- predict(fm, type="state", newdata=nd, se.fit=FALSE)$Predicted

This will give you a vector of probabilities summing to 1 for each possible state in a 3-species system e.g. [111], [110], [101] and so on. You would have to add the appropriate probabilities together to simplify this down to two-species combinations, e.g. if you want [11] representing probability of sites being occupied by species 1 and 3 regardless of the occupancy state of species 2, you would add together the probabilities for states [111] and [101]. 

In addition to the problem with the link function, the gigantic values you were getting make me think you might have parameter estimate(s) in your model that have a very large absolute magnitude and/or SE which would indicate problems with the model fit.

Ken

Reply all
Reply to author
Forward
0 new messages