# 1 Introduction

## 1.1 What is ubms?

ubms is an R package for fitting models of wildlife occurrence and abundance with datasets where animals are not individually marked. It provides a nearly identical interface to the popular unmarked package [1]. Instead of using maximum likelihood to fit models (as with unmarked), models are fit in a Bayesian framework using Stan [2]. It is generally expected that you are already familiar with unmarked when using ubms. You can download ubms, report issues, or help with development on Github.

## 1.2 Why should you use it?

There are several advantages to using ubms over unmarked. First, it is possible to include random effects in ubms models, which is not currently possible in unmarked. These are specified using the familiar syntax of lme4 [3]. Second, ubms generates posterior distributions for model parameters, including latent occupancy and abundance parameters. These can be useful for post-hoc analyses and diagnostics. Finally, fitting models with Stan gives you access to the large ecosystem of Stan-related tools, such as LOO (leave-one-out-cross validation; [4]).

Another alternative to ubms would be to fit models in an existing modeling language such as BUGS, JAGS, or directly in Stan. ubms abstracts away the complex process of defining and writing custom models in these languages which need to be updated each time you make changes to your model. It also provides many useful helper functions (e.g. predict) which would otherwise require custom code. Finally, because of Stan’s efficient sampler [5] and because the underlying likelihoods in ubms are marginalized, ubms will often fit models much faster than equivalent models in BUGS or JAGS [6].

## 1.3 What are the drawbacks?

Relative to unmarked, ubms has fewer types of models available. For example, models that incorporate temporary emigration (like gdistsamp) [7] are currently not available in ubms. Models should run much faster in unmarked as well. If you do not need one of the specific benefits of ubms described above, it makes sense to stick with unmarked. Even if you do plan to use ubms, it makes sense to test the models in unmarked first. The similar interface between the two packages makes this very easy, as you will see in the next section.

Relative to BUGS/JAGS/Stan, ubms is less flexible because you cannot customize your model structure. You are limited to the provided model types. Furthermore, you cannot currently customize prior distributions (although I plan to add this in the future in some form). Finally, writing your own BUGS/JAGS model can be valuable for gaining a deeper understanding of how a model works; ubms, like unmarked, is essentially a black box.

To summarize the advantages and disadvantages: I see ubms as an intermediate step along the continuum from unmarked to custom models in BUGS/JAGS. It is not meant to replace either approach but rather to supplement them, for situations when a Bayesian framework is needed and “off-the-shelf” model structures are adequate.

# 2 Example: Fitting a single-season occupancy model

Occupancy models estimate the probability $$\psi$$ that a species occupies a site, while accounting for detection probability $$p < 1$$ [8]. In order to estimate both $$p$$ and $$\psi$$, repeated observations (detection/non-detection data) at each site are required.

## 2.1 Set up the data

First, load the dataset we’ll be using, which comes with unmarked:

library(unmarked)
data(crossbill)

The crossbill dataset is a data.frame with many columns. It contains detection/non-detection data for the European crossbill (Loxia curvirostra) in Switzerland [9].

dim(crossbill)
## [1] 267  58
names(crossbill)
##  [1] "id"      "ele"     "forest"  "surveys" "det991"  "det992"  "det993"
##  [8] "det001"  "det002"  "det003"  "det011"  "det012"  "det013"  "det021"
## [15] "det022"  "det023"  "det031"  "det032"  "det033"  "det041"  "det042"
## [22] "det043"  "det051"  "det052"  "det053"  "det061"  "det062"  "det063"
## [29] "det071"  "det072"  "det073"  "date991" "date992" "date993" "date001"
## [36] "date002" "date003" "date011" "date012" "date013" "date021" "date022"
## [43] "date023" "date031" "date032" "date033" "date041" "date042" "date043"
## [50] "date051" "date052" "date053" "date061" "date062" "date063" "date071"
## [57] "date072" "date073"

Check ?crossbill for details about each column. The first three columns id, ele, and forest are site covariates.

site_covs <- crossbill[,c("id", "ele", "forest")]

The following 27 columns beginning with det are the binary detection/non-detection data; 9 years with 3 observations per year. For this example we want to fit a single-season occupancy model; thus we will use only the first three columns (year 1) of det as our response variable y.

y <- crossbill[,c("det991","det992","det993")]
head(y)
##   det991 det992 det993
## 1      0      0      0
## 2      0      0      0
## 3     NA     NA     NA
## 4      0      0      0
## 5      0      0      0
## 6     NA     NA     NA

Note that missing values are possible. The final 27 columns beginning with date are the Julian dates for each observation. As with y we want only the first three columns corresponding to year 1.

date <- crossbill[,c("date991","date992","date993")]

Finally, we build our unmarkedFrame object holding our detection/non-detection data, site covariates, and observation covariates. Since we will conduct a single-season occupancy analysis, we need to use unmarkedFrameOccu specifically. The resulting unmarkedFrame can be used by both unmarked and ubms.

umf <- unmarkedFrameOccu(y=y, siteCovs=site_covs, obsCovs=list(date=date))
head(umf)
## Data frame representation of unmarkedFrame object.
##    y.1 y.2 y.3 id  ele forest date.1 date.2 date.3
## 1    0   0   0  1  450      3     34     59     65
## 2    0   0   0  2  450     21     17     33     65
## 3   NA  NA  NA  3 1050     32     NA     NA     NA
## 4    0   0   0  4  950      9     29     59     65
## 5    0   0   0  5 1150     35     24     45     65
## 6   NA  NA  NA  6  550      2     NA     NA     NA
## 7    0   0   0  7  750      6     26     54     74
## 8    0   0   0  8  650     60     23     43     71
## 9    0   0   0  9  550      5     21     36     56
## 10   0   0   0 10  550     13     37     62     75

## 2.2 Fit a model

### 2.2.1 Fit the model in unmarked

First, we fit a null model (no covariates) in unmarked using the occu function. The occu function requires as input a double formula (for detection and occupancy, respectively) along with our unmarkedFrame.

(fit_unm <- occu(~1~1, data=umf))
##
## Call:
## occu(formula = ~1 ~ 1, data = umf)
##
## Occupancy:
##  Estimate    SE     z P(>|z|)
##    -0.546 0.218 -2.51  0.0121
##
## Detection:
##  Estimate    SE     z P(>|z|)
##    -0.594 0.208 -2.86 0.00426
##
## AIC: 511.2538

### 2.2.2 Fit the model in ubms

Next, we fit the same model in ubms. The equivalent to occu in ubms is stan_occu. Functions in ubms generally use this stan_ prefix, based on the approach used in package rstanarm for GLMs. We need to provide the same arguments to stan_occu. In addition, we will specify that we want 3 MCMC chains (chains=3), with 500 iterations per chain (iter=500) of which the first half will be warmup iterations. It is beyond the scope of this vignette to discuss the appropriate number or length of chains; see the Stan user’s guide for more details. Generally 4 chains of 2000 iterations each is recommended (of which 1000 per chain are warmups). Thus, 500 iterations per chain is probably not enough, but to keep things running quickly it is sufficient for this vignette. Note that if you are more familiar with BUGS or JAGS, Stan generally requires a smaller number of iterations to reach convergence thanks to its default NUTS sampler [5]. If you have a good multi-core CPU, you can run chains in parallel. Tell Stan how many cores you want to use by setting the mc.cores option.

library(ubms)
options(mc.cores=3)

(fit_stan <- stan_occu(~1~1, data=umf, chains=3, iter=500))
##
## Call:
## stan_occu(formula = ~1 ~ 1, data = umf, chains = 3, iter = 1000,
##     refresh = 0)
##
## Occupancy (logit-scale):
##  Estimate    SD   2.5%   97.5% n_eff  Rhat
##    -0.506 0.231 -0.922 -0.0334  1237 0.999
##
## Detection (logit-scale):
##  Estimate    SD  2.5%  97.5% n_eff Rhat
##    -0.625 0.208 -1.04 -0.215  1422    1
##
## LOOIC: 511.771

### 2.2.3 Compare results

The structure of the output from unmarked and ubms is intentionally similar. Estimates of the occupancy and detection parameters are also similar, but not identical. For a more direct comparison, call the coef function on both model objects:

cbind(unmarked=coef(fit_unm), stan=coef(fit_stan))
##            unmarked       stan
## psi(Int) -0.5461203 -0.5059667
## p(Int)   -0.5939612 -0.6245178

### 2.2.4 Understanding the ubms output summary

Let’s look at the output from our fit_stan model again:

fit_stan
##
## Call:
## stan_occu(formula = ~1 ~ 1, data = umf, chains = 3, iter = 1000,
##     refresh = 0)
##
## Occupancy (logit-scale):
##  Estimate    SD   2.5%   97.5% n_eff  Rhat
##    -0.506 0.231 -0.922 -0.0334  1237 0.999
##
## Detection (logit-scale):
##  Estimate    SD  2.5%  97.5% n_eff Rhat
##    -0.625 0.208 -1.04 -0.215  1422    1
##
## LOOIC: 511.771

The first part (under Call:) is the command we used to get this model output. Underneath are two tables, one per submodel, corresponding to the occupancy and detection parts of the model. Within each table there is one row per parameter in the submodel. Since fit_stan had no covariates, there is only an intercept term for each submodel. Model parameters in this summary table are always shown on the appropriate transformed scale, in this case logit. To get the corresponding probabilities, you can use the predict function, which we will demonstrate later.

For each parameter, the mean and standard deviation of the posterior distribution are given. Unlike unmarked output, there is no $$Z$$ or $$p$$-value. Instead, there is a 95% uncertainty interval provided.

The final two columns in each summary table n_eff and Rhat are MCMC diagnostics. We will discuss their meaning later.

### 2.2.5 Extracting individual parameters

To extract summary values into an R table for easy manipulation, use the summary method. Note that you have to specify which submodel you want ("state" for occupancy or "det" for detection).

sum_tab <- summary(fit_stan, "state")
sum_tab\$mean[1]
## [1] -0.5059667

To extract the entire posterior for a parameter, use the extract method. To avoid name collisions you need to use the full name of the parameter (which contains both the submodel and the parameter name) when extracting. To see a list of the available full parameter names, use the names method.

names(fit_stan)
## [1] "beta_state[(Intercept)]" "beta_det[(Intercept)]"
occ_intercept <- extract(fit_stan, "beta_state[(Intercept)]")[[1]]
hist(occ_intercept, freq=FALSE)
lines(density(occ_intercept), col='red', lwd=2)

## 2.3 Compare candidate models

Now we’ll fit a series of candidate models to the crossbill data in ubms and compare them.

### 2.3.1 Fit the models

Along with our previous null model, we’ll fit a model with site covariates (elevation and forest), a model with an observation covariate (date), and a “global” model with both site and observation covariates. This is just an example; perhaps other models should also be considered if we were preparing this analysis for publication. In our model formulas, we have normalized all covariates with scale so they have a mean of 0 and a standard deviation of 1. This can help improve model convergence and is generally a good idea.

fit_null <- fit_stan

fit_sc <- stan_occu(~1~scale(forest)+scale(ele), data=umf,
chains=3, iter=500)

fit_oc <- stan_occu(~scale(date)~1, data=umf, chains=3, iter=500)

fit_global <- stan_occu(~scale(date)~scale(forest)+scale(ele), data=umf,
chains=3, iter=500)
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#tail-ess

The fit_global model gave us some warnings about the effective sample size (n_eff) along with a suggested solution. We will ignore this warning for now but normally it is a good idea to pay close attention to these warnings.

### 2.3.2 Compare the models

First we combine the models into a fitList:

mods <- fitList(fit_null, fit_sc, fit_oc, fit_global)

Then we generate a model selection table:

round(modSel(mods), 3)
##                elpd nparam elpd_diff se_diff weight
## fit_global -237.431  7.012     0.000   0.000  0.812
## fit_sc     -241.782  5.599    -4.351   3.566  0.142
## fit_oc     -249.894  3.418   -12.463   5.158  0.046
## fit_null   -255.886  2.418   -18.454   6.622  0.000

Instead of AIC, models are compared using leave-one-out cross-validation (LOO) [4] via the loo package. Based on this cross-validation, the expected predictive accuracy (elpd) for each model is calculated. The model with the largest elpd (fit_global) performed best. The elpd_diff column shows the difference in elpd between a model and the top model; if this difference is several times larger than the standard error of the difference (se_diff), we are confident the model with the larger elpd performed better. LOO model weights, analogous to AIC weights, are also calculated. We can see that the fit_global model is clearly the best performing model.

You can obtain LOO information for a single model using the loo method:

loo(fit_global)
##
## Computed from 1500 by 245 log-likelihood matrix
##
##          Estimate   SE
## elpd_loo   -237.4 19.3
## p_loo         7.0  0.7
## looic       474.9 38.6
## ------
## Monte Carlo SE of elpd_loo is 0.2.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

The looic value is analogous to AIC.

You can also obtain the WAIC (Widely Applicable Information Criterion) if you prefer [4]:

waic(fit_global)
##
## Computed from 1500 by 245 log-likelihood matrix
##
##           Estimate   SE
## elpd_waic   -237.4 19.3
## p_waic         7.0  0.7
## waic         474.8 38.6

## 2.4 Diagnostics and model fit

We’ll define the global model as our top model:

fit_top <- fit_global
fit_top
##
## Call:
## stan_occu(formula = ~scale(date) ~ scale(forest) + scale(ele),
##     data = umf, chains = 3, iter = 1000, refresh = 0)
##
## Occupancy (logit-scale):
##               Estimate    SD   2.5% 97.5% n_eff Rhat
## (Intercept)     -0.494 0.665 -1.150  1.72  59.9 1.03
## scale(forest)    1.231 0.658  0.631  3.45  60.7 1.03
## scale(ele)       0.566 0.252  0.049  1.03 453.9 1.00
##
## Detection (logit-scale):
##             Estimate    SD   2.5%  97.5% n_eff Rhat
## (Intercept)   -0.791 0.282 -1.470 -0.339   143 1.01
## scale(date)    0.563 0.172  0.237  0.906   978 1.00
##
## LOOIC: 474.862

### 2.4.1 MCMC diagnostics

Again looking at the summary of fit_top, it appears all chains have adequately converged based on $$\hat{R}$$ values. Ideally, you want all $$\hat{R} > 1.05$$. To visualize convergence, look at the traceplots:

traceplot(fit_top)