These days, the Bayesian hierarchical modeling framework is more-or-less standard for analyzing BBS data. I wrote a collection of posts about gathering BBS and remote sensing data, summarizing and importing into R, and thereby getting the data frame in shape for an analysis. From there, I would clean up the work space such that only the data frame for the analysis is left.
rm(list= ls()[!(ls() %in% c('vars'))])
Here’s an example of a formula (“test”), and how to fit it using the package R-INLA. In this model, in keeping with the example that I linked to that got us here, “MN” (mean % tree cover) is our variable of interest, and the rest of the variables are nuisance effects. In this model, year and rteobs are random effects (hyper-parameters), while yearfix, firstyear and MN are fixed effects (parameters).
test = abundance ~ f(Year,model="iid") + yearfix + f(rteobs, model="iid") + firstyear + MN library(INLA) model = inla(formula=test, data=vars, family="gpoisson", control.compute = list(waic=TRUE), verbose=TRUE)
The generalized Poisson distribution is used to incorporate over-dispersion, and I chose to calculate Watanabe-Akaike Information Criterion (WAIC) in case I was to compare this to another model. To examine the output once the model has stopped running…
summary(model) names(model)
The latter shows all that was generated in this model run. One I find useful is summary.fitted.values, which are the modeled values of the response variable. If you want to predict values, you can add or replace response variable values in your data frame with NA.