In jsocolar/flocker: Flexible Occupancy Estimation with Stan

Here we show how we can use flocker to fit nonlinear occupancy models via brms. In most occupancy models, occupancy and detection probabilities are modeled as logit-linear combinations of covariates. In some models (e.g. those with splines or Gaussian processes), probabilities are modeled as the sum of more flexible functions of covariates. These are straightforward to fit in flocker using the brms functions s(), t2(), and gp(); see the flocker tutorial vignette for details.

This vignette focuses on more complicated nonlinear models that require the use of special nonlinear brms formulas. We showcase two models. The first fits a parametric nonlinear predictor. The second fits a model with a spatially varying coefficient that is given a gaussian process prior.

Parametric nonlinear predictor

In this scenario, we consider a model where the response is a specific nonlinear parametric function whose parameters are fitted and might or might not depend on covariates. Suppose for example that an expanding population of a territorial species undergoes logistic growth, and also that some unknown proportion of territories are unsuitable due to an unobserved factor, such that occupancy asymptotes at some probability less than one. Thus, occupancy probability changes through time as $\frac{L}{1 + e^{-k(t-t_0)}}$, where $L$ is the asymptote, $k$ is a growth rate, $t$ is time, and $t_0$ is the timing of the inflection point. At multiple discrete times, we randomly sample several sites to survey, and survey each of those sites over several repeat visits.

library(flocker); library(brms)
set.seed(3)

L <- 0.5
k <- .1
t0 <- -5
t <- seq(-15, 15, 1)
n_site_per_time <- 30
n_visit <- 3
det_prob <- .3

data <- data.frame(
  t = rep(t, n_site_per_time)
)

data$psi <- L/(1 + exp(-k*(t - t0)))
data$Z <- rbinom(nrow(data), 1, data$psi)
data$v1 <- data$Z * rbinom(nrow(data), 1, det_prob)
data$v2 <- data$Z * rbinom(nrow(data), 1, det_prob)
data$v3 <- data$Z * rbinom(nrow(data), 1, det_prob)

fd <- make_flocker_data(
  obs = as.matrix(data[,c("v1", "v2", "v3")]),
  unit_covs = data.frame(t = data[,c("t")]),
  event_covs <- list(dummy = matrix(rnorm(n_visit*nrow(data)), ncol = 3))
)

We wish to fit an occupancy model that recovers the unknown parameters $L$, $k$, and $t_0$. We can achieve this using the nonlinear formula syntax provided by brms via flocker.

flocker will always assume that the occupancy formula is provided on the logit scale. Thus, we need to convert our nonlinear function giving the occupancy probability to a function giving the logit occupancy probability. A bit of simplification via Wolfram Alpha and we arrive at $\log(\frac{L}{1 + e^{-k(t - t_0)} - L})$. We then write a brms formula representing occupancy via this function. To specify a formula wherein a distributional parameter (occ in this case, referring to occupancy) is nonlinear we need to use brms::set_nl() rather than merely providing the nl = TRUE argument to brms::bf().

flocker's main fitting function flock() accepts brmsformula inputs to its f_det argument. When supplying a brmsformula to f_det (rather than the typical one-sided detection formula), the following behaviors are triggered:

Several input checks are turned off. For example, flocker no longer checks to ensure that event covariates are absent from the occupancy formula. flocker also no longer explicitly checks that formulas are provided for all of the required distributional terms for a given family (detection, occupancy, colonization, extinction, and autologistic terms, depending on the family).
All inputs to f_occ, f_col, f_ex, f_auto are silently ignored. It is obligatory to pass the entire formula for all distributional parameters as a single brmsformula object. This means in turn that the user must be familiar with flocker's internal naming conventions for all of the relevant distributional parameters (det and one or more of occ, colo, ex, autologistic, Omega). If fitting a data-augmented model, it will be required to pass the Omega ~ 1 formula within the brmsformula (When passing the traditional one-sided formula to f_det, flocker includes the formula for Omega internally and automatically).
Nonlinear formulas that involve data that are required to be positive might fail! Internally, some irrelevant data positions get filled with -99, but these positions might still get evaluated by the nonlinear formula, even though they make no contribution to the likelihood.

With all of that said, we can go ahead and fit this model!

fit <- flock(f_det = brms::bf(
                 det ~ 1 + dummy,
                 occ ~ log(L/(1 + exp(-k*(t - t0)) - L)),
                 L ~ 1,
                 k ~ 1,
                 t0 ~ 1
               ) +
               brms::set_nl(dpar = "occ"),
             prior = 
               c(
                 prior(normal(0, 5), nlpar = "t0"),
                 prior(normal(0, 1), nlpar = "k"), 
                 prior(beta(1, 1), nlpar = "L", lb = 0, ub = 1)
                ),
             flocker_data = fd, 
             control = list(adapt_delta = 0.9),
             cores = 4)

summary(fit)
#>  Family: occupancy_single 
#>   Links: mu = identity; occ = identity 
#> Formula: ff_y | vint(ff_n_unit, ff_n_rep, ff_Q, ff_rep_index1, ff_rep_index2, ff_rep_index3) ~ 1 + dummy 
#>          occ ~ log(L/(1 + exp(-k * (t - t0)) - L))
#>          L ~ 1
#>          k ~ 1
#>          t0 ~ 1
#>    Data: data (Number of observations: 2790) 
#>   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup draws = 4000
#> 
#> Population-Level Effects: 
#>              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept       -0.91      0.13    -1.18    -0.66 1.00     2535     2318
#> L_Intercept      0.50      0.10     0.38     0.77 1.00     1241      864
#> k_Intercept      0.19      0.08     0.07     0.36 1.00     1283     1433
#> t0_Intercept    -5.90      3.38   -10.34     3.17 1.00     1248      921
#> dummy            0.00      0.08    -0.15     0.15 1.00     2620     2236
#> 
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

It works!

Note that if desired, we could fit more complicated formulas than ~ 1 for any of the nonlinear parameters. For more see the brms nonlinear model vignette.

Spatially varying coefficients via a Gaussian process

The gp() function in brms includes a Gaussian process of arbitrary dimension in the linear predictor. We can use the nonlinear formula syntax to tell brms to include a Gaussian process prior on a coefficient as well.

First we simulate some data wherein the logit of the occupancy probability depends on a covariate, and the slope of the dependency is modeled via a two-dimensional spatial Gaussian process. It turns out that we will need quite a few of data points to constrain the standard deviation of the Gaussian process, so we simulate with 2000 sites:

set.seed(1)
n <- 2000 # sample size
lscale <- 0.3 # square root of l of the gaussian kernel
sigma_gp <- 1 # sigma of the gaussian kernel
intercept <- 0 # occupancy logit-intercept
det_intercept <- -1 # detection logit-intercept
n_visit <- 4

# covariate data for the model
gp_data <- data.frame(
  x = rnorm(n), 
  y = rnorm(n),
  covariate = rnorm(n)
  )

# get distance matrix
dist.mat <- as.matrix(
  stats::dist(gp_data[,c("x", "y")])
  )

# get covariance matrix
cov.mat <- sigma_gp^2 * exp(- (dist.mat^2)/(2*lscale^2))

# simulate occupancy data
gp_data$coef <- mgcv::rmvn(1, rep(0, n), cov.mat)
gp_data$lp <- intercept + gp_data$coef * gp_data$covariate
gp_data$psi <- boot::inv.logit(gp_data$lp)
gp_data$Z <- rbinom(n, 1, gp_data$psi)

# simulate visit data
obs <- matrix(nrow = n, ncol = n_visit)
for(j in 1:n_visit){
  obs[,j] <- gp_data$Z * rbinom(n, 1, boot::inv.logit(det_intercept))
}

And here's how we can fit this model in flocker! Because we have a large number of sites, we use a Hilbert space approximate Gaussian process for computational efficiency.

fd2 <- make_flocker_data(obs = obs, unit_covs = gp_data[, c("x", "y", "covariate")])
svc_mod <- flock(
  f_det = brms::bf(
                 det ~ 1,
                 occ ~ occint + g * covariate,
                 occint ~ 1,
                 g ~ 0 + gp(x, y, scale = FALSE, k = 20, c = 1.25)
               ) +
               brms::set_nl(dpar = "occ"),
  flocker_data = fd2,
  cores = 4
)

summary(svc_mod)
#>  Family: occupancy_single_C 
#>   Links: mu = identity; occ = identity 
#> Formula: ff_n_suc | vint(ff_n_trial) ~ 1 
#>          occ ~ occint + g * covariate
#>          occint ~ 1
#>          g ~ 0 + gp(x, y, scale = FALSE, k = 20, c = 1.25)
#>    Data: data (Number of observations: 2000) 
#>   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup draws = 4000
#> 
#> Gaussian Process Terms: 
#>                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sdgp(g_gpxy)       1.66      0.75     0.73     3.66 1.00     2849     3317
#> lscale(g_gpxy)     0.23      0.11     0.08     0.50 1.00     3850     3141
#> 
#> Population-Level Effects: 
#>                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept           -1.04      0.06    -1.15    -0.94 1.00     5532     2810
#> occint_Intercept     0.09      0.09    -0.08     0.27 1.00     5736     3182
#> 
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).