Determining break point/threshold in species richness in relation to land cover from multi-species occupancy models in R I am conducting research in which I would like to see if I can detect a threshold response in mammal species richness in relation to the surrounding level of forest  across study sites. I plan to use multi-species occupancy modelling methods using a bayesian approach in R to determine species richness in relation to the forest cover covariate, but I am a bit stuck on how I would then proceed to determine break points/thresholds in the species richness across sites in relation to forest cover. I have done a bit of a research into methods/tools for determining breakpoints (such as Muggeo's segmented package in R); however, as I am admittedly very naive still when it comes to these methods (and R in general) I am wondering if anyone has suggestions for how I could both use these multi-species occupancy modelling methods to determine species richness in relation to the forest cover covariate, and also determine if there is a threshold relationship between species richness and forest cover in any sites using tools in R? Thanks very much!
 A: This is a great question. I'm going to address it in the context of the model as formulated by this paper: https://onlinelibrary.wiley.com/doi/full/10.1002/ece3.4821
That is to say you have some occupation probability vector $\psi_{ik}$ where k represents the species and i the site (using the notation from the linked paper). You now want to model those probabilities in terms of some covariates (e.g., forest cover).
The linked paper puts this in terms of a logistic regression: $$logit(\psi_{ik}) = \beta_{0k} + \beta_{1k}*X_{1i} + \beta_{2k}*X_{2i} + ...$$
So in this case since you only have forest cover you would not include a $X_{2i}$.
Now the thing about your particular model is that you believe that there is a breakpoint in this relationship so let's assume that you believe below some forest cover there is no effect and above it there is. Let's call the forest cover at which you think there is a change as FC. Here $X_{1i}$ is the Forest Cover at site i. Then you have effectively two models:
1) If $X_{1i}$ $\le$ FC then $logit(\psi_{ik}) = \beta_{0k}$
2) If $X_{1i}$ > FC then $logit(\psi_{ik}) = \beta_{0k} + \beta_{1k}*X_{1i}$
What you want to do is estimate $\beta_{0k}$, $\beta_{1k}$, and FC
So the way that you do that is you fit a piecewise linear regression (because our model is a linear regression but piecewise because we're essentially sticking different lines together)
Sorry for that preamble I just wanted to layout the model in math so that we're on the same page. The way to fit this is to evaluate a model that has the flexibility to move the breakpoint and then basically use a fitting algorithm to infer the optimal parameters. You could then also evaluate a model without the breakpoint and perform a comparison to see if the breakpoint is necessary in the first place. 
To implement this in R I would look at the following reference:
https://janhove.github.io/analysis/2018/07/04/bayesian-breakpoint-model <- this implements a bayesian piecewise regression in Stan and R. I particularly like this approach because you can see how the model works under the hood and because Stan is a powerful probabilistic programming language.
You could also look at the changepoint literature as what you're asking is if there is a change in the slope at some value. This means you could use the bcp package (https://www.jstatsoft.org/article/view/v023i03/v23i03.pdf) to address your problem.
The nice thing about using a changepoint analysis for this is that it would give you a measure of the confidence in the existence of that changepoint and not just where it is. Indeed the bcp package is capable of deciding whether or not the changepoint should be included in the model based on the BIC criteria.
