ubms is an R package for fitting Bayesian hierarchical models of animal occurrence and abundance. The package has a formula-based interface compatible with unmarked, but the model is fit using MCMC with Stan instead of using maximum likelihood. Currently there are Stan versions of unmarked functions occu, occuRN, colext, occuTTD, pcount, distsamp, and multinomPois. These functions follow the stan_ prefix naming format established by rstanarm. For example, the Stan version of the unmarked function occu is stan_occu.

Advantages compared to unmarked:

  1. Obtain posterior distributions of parameters and derived parameters
  2. Include random effects in parameter formulas (same syntax as lme4)
  3. Assess model fit using WAIC and LOO via the loo package

Disadvantages compared to unmarked:

  1. MCMC is slower than maximum likelihood
  2. Not all model types are supported
  3. Potential for convergence issues


ubms is on CRAN:

Alternatively, the latest development version can be installed from Github:

# install.packages("devtools")

If you are on Windows, you can download and install a pre-compiled binary package of the latest release:

# Install dependencies
install.packages(c("unmarked", "ggplot2", "gridExtra", "lme4", "loo",
                   "Matrix", "Rcpp", "rstan", "rstantools"))

# Download and install ubms
install.packages("ubms_1.0.2.zip", repos=NULL)


Simulate occupancy data including a random effect on occupancy:


dat_occ <- data.frame(x1=rnorm(500))
dat_p <- data.frame(x2=rnorm(500*5))

y <- matrix(NA, 500, 5)
z <- rep(NA, 500)

b <- c(0.4, -0.5, 0.3, 0.5)

re_fac <- factor(sample(letters[1:26], 500, replace=T))
dat_occ$group <- re_fac
re <- rnorm(26, 0, 1.2)
re_idx <- as.numeric(re_fac)

idx <- 1
for (i in 1:500){
  z[i] <- rbinom(1,1, plogis(b[1] + b[2]*dat_occ$x1[i] + re[re_idx[i]]))
  for (j in 1:5){
    y[i,j] <- z[i]*rbinom(1,1, 
                    plogis(b[3] + b[4]*dat_p$x2[idx]))
    idx <- idx + 1

Create unmarked frame:

umf <- unmarkedFrameOccu(y=y, siteCovs=dat_occ, obsCovs=dat_p)

Fit a model with a random intercept:

options(mc.cores=3) #number of parallel cores to use
(fm <- stan_occu(~x2 ~x1 + (1|group), umf, chains=3))

Examine residuals for occupancy and detection submodels (following Wright et al. 2019). Each panel represents a draw from the posterior distribution.


Assess model goodness-of-fit with a posterior predictive check, using the MacKenzie-Bailey chi-square test:

fm_fit <- gof(fm, draws=500)

Look at the marginal effect of x2 on detection:

plot_marginal(fm, "det")

A more detailed vignette for the package is available here.