nonlinear vignette

vignettes/nonlinear_models.Rmd

@@ -0,0 +1,87 @@
+---
+title: "Nonlinear models in flocker"
+author: "Jacob Socolar"
+date: "2023-10-04"
+output: html_document
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+```
+
+## R Markdown
+
+Here we show how we can use `flocker` to fit nonlinear occupancy models via `brms`. Suppose an expanding population of a territorial species undergoes logistic growth, and also that some 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.
+
+```{r true-occupancy}
+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`. 
+
+Our first challenge is that `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}{-e^{-k(t - t_0)} + L - 1})$.
+
+Our second challenge is in remembering how to write down a nonlinear brms formula where a distributional parameter is nonlinear. `brms::set_nl()` is our friend here.
+
+Our final challenge is to pass a `brmsformula` as an input function to `flocker`'s main fitting function `flock()`, which accepts `brmsformula` inputs to its `f_det` argument. When supplying a `brmsformula` to `f_det`, the following behaviors are triggered:
+
+* Several input checks are turned off.  For example, we no longer check to ensure that no event covariates are used in the occupancy formula.  We also no longer explicitly check that we have formulas for all of the required distributional terms for a given family.
+
+* All inputs to `f_occ`, `f_det`, `f_col`, `f_ex`, `f_auto` are silently ignored. It is obligatory to pass the entire formula as a `brmsformula` object. In turn this means that the user must be familiar with `flocker`'s internal naming conventions for all of the relevant distributional parameters (`det` and at least one of `occ`, `colo`, `ex`, `autologistic`, `Omega`). If fitting a data-augmented model, it will be requried to pass the `Omega ~ 1` formula within the `brmsformula` (otherwise `flocker` includes this formula 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!
+
+```{r fitting}
+fit <- flock(f_det = brms::bf(
+                 det ~ 1 + dummy,
+                 occ ~ log(-L/(-exp(-k*(t - t0)) + L - 1)),
+                 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, 
+             adapt_delta = .9,
+             backend = "cmdstanr", 
+             cores = 4)
+summary(fit)
+
+```
+
+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](https://cran.r-project.org/web/packages/brms/vignettes/brms_nonlinear.html).
+
+