How to obtain psi and p coefficients in community and species level

[Natasha Grosch Loureiro]

Hello, everyone. I'm using spOccupancy package (function stMsPGOcc), and I have acquired the values at the community and species level for each covariate, but I'm having trouble finding how to obtain the occupation and detection coefficients for the community and species level. Could anyone help me? Thank you.

[Jeffrey Doser]

Hi Natasha, 

There is no direct estimate of "psi" and "p" at the community level, it is the coefficients (the "beta.comm.samples" and the "alpha.comm.samples") that are community-level parameters. If you want to get a value of average occupancy probability across all species in the community (so an average psi across all species), you could the following, assuming that your saved model output from stMsPGOcc is called "out"

psi.avg.samples <- apply(out$psi.samples, c(1, 3, 4)

That gives you the average occupancy probability across all species at each site in each season (for each MCMC sample). 

For detection, there is one more step as you would need to use the fitted() function to get the species specific detection probability values at each site/season/replicate: 

fitted.out <- fitted(out)

# Average detection probability for each MCMC sample, site, season, and replicate

p.avg.samples <- apply(fitted.out$p.samples, c(1, 3, 4, 5))

Alternatively, if you just want a single overall value of average occupancy or detection at the community level, you can back-transform the community-level intercept that shows in the summary output. So, if your community-level occupancy intercept was 0, the value plogis(0) would be the average occupancy probability across all sites and seasons (when all covariates in your model are 0, which if you standardized them is the mean). 

Hope that helps,

Jeff

[Natasha Grosch Loureiro]

Thank you for your help! It is possible do to some procedure to estimate psi and p at the species level?

Best wishes,

Natasha

[Jeffrey Doser]

Hi Natasha,

Yes, you can get estimates of occurrence probability (psi) by species in the "out$psi.samples" component of the model object, where "out" is the name of the model object. For detection probability (p), you can get those estimates by using the fitted() function (e.g., p.samples <- fitted(out)$p.samples will give you detection probabilities by species at each of the sites/replicates in your data set). 

Jeff