generalised N-mixture models (state-space form)

[MilaniC]

Hi Paul

Preamble:
Further to previous conversation (and the following depends on whether brms is configured to allow for 2 separate linear predictor functions on the same response value (hence the state-space formulation)). Stan can do this but can brms be brought to bear on formulating the model for Stan. If it can then that provides a very powerful frontend to Stan for implementing a wide range of demographic models with the full Bayesian inference machinery and modelling support like various likelihoods (ziP, ZiNB, etc) and nonparametric (GAM) covariate functional form on either or both the state and observation components.

The Issue:
With the bf(,part 1, part 2, ..) formula form for brm(), would it be possible to fit a state-space model formulation like in the Stan discussion here:

https://groups.google.com/forum/#!topic/stan-users/9mMsp1oB69g

part1= the latent or hidden state component, part 2 = the observation component

These so-called N-mixtures models are implemented in R at the moment via the “unmarked” package on CRAN:

Fiske, Chandler (2011) unmarked: An R package for fitting hierarchical models of wildlife occurrence and abundance. Journal of Statistical Software, 43(10):1–23, 2011

These models are used to estimate various wildlife population metrics like abundance for unmarked animals sampled at multiple sites over multiple occasions while accounting for imperfect detection.

The demographic model most people are after is the generalised N-mixture model (also known as the Dail-Madsen model and there are recent extensions). In “unmarked” that is the “pcountOpen() function"

We already fit some forms of them using JAGS/rjags but in Stan via brms might be a way better approach.

Very simple population dynamics state-space models are already fitted using Stan as in our script below for a Australia flatback marine turtle nesting population in Western Australia.

It uses some Stan code from here (https://github.com/stan-dev/example-models/tree/master/BPA)

But it is far too simple EGSS model (useful nonetheless) but importantly is limited to Gaussian likelihood - it really should be Poisson on NB likelihood (and that is the problem discussed at:

https://groups.google.com/forum/#!topic/stan-users/9mMsp1oB69g

Anyway, the N-mixture state-space population model discussed above (https://groups.google.com/forum/#!topic/stan-users/9mMsp1oB69g) could be a new frontier for brms!!

Not sure if this is possible with brms as the frontend to Stan though you are in a far far better position to know.

cheers

Milani

## Mundabullangana nester time series Bayesian state-space model using STAN in R

## based on code from Kery & Schaub 2012 BPA WinBUGS book but translated into STAN …
## https://github.com/stan-dev/example-models/tree/master/BPA

setwd("/ MC Files/Flatback RLA")
require(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
set.seed(123)
require(coda)
require(ggmcmc)
require(gridExtra)
require(ggthemes)

## specify model as character vector using model_code instead of calling the separate 'file = "ssm2.stan" '
mda.model.code<-"
data {
  int<lower=0> T;
  vector[T] y;
  int<lower=0> pyears;                // Number of years for prediction
}

parameters {
  real logN_est1;                     // Initial log population size
  real mean_r;                        // Mean growth rate
  real<lower=0,upper=1> sigma_proc;   // SD of state process
  real<lower=0,upper=1> sigma_obs;    // SD of observation process
  vector[T-1] r;
}

transformed parameters {
  vector[T] logN_est;                 // Log population size

  // State process
  logN_est[1] <- logN_est1;
  for (t in 1:(T - 1))
    logN_est[t + 1] <- logN_est[t] + r[t];
}

model {
  // Priors and constraints
  logN_est1 ~ normal(5.6, 10);
  mean_r ~ normal(1, sqrt(1000));
  sigma_proc ~ uniform(0, 1);
  sigma_obs ~ uniform(0, 1);

  // Likelihood
  r ~ normal(mean_r, sigma_proc);

  // Observation process
  y ~ normal(logN_est, sigma_obs);
}

generated quantities {
  // Population sizes on real scale
  real sigma2_proc;
  real sigma2_obs;
  vector[pyears] pr;
  vector[pyears] plogN_est;              // Predicted log population size
  vector<lower=0>[T + pyears] N_est;     // Population size

  sigma2_proc <- square(sigma_proc);
  sigma2_obs <- square(sigma_obs);

  pr[1] <- normal_rng(mean_r, sigma_proc);
  plogN_est[1] <- logN_est[T] + pr[1];  
  for (t in 2:pyears) {
    pr[t] <- normal_rng(mean_r, sigma_proc);
    plogN_est[t] <- plogN_est[t - 1] + pr[t];
  }
  for (t in 1:T)
    N_est[t] <- exp(logN_est[t]);
  for (t in 1:pyears)
    N_est[T + t] <- exp(plogN_est[t]);
}
"

## Munda data
pyears <- 3 # Number of future years for predictions
## mda <- c(863,1361,1712,2101,NA,1336,3150,1826,2166,1939,1619,2314,1609,2176,2227,2298,2366,2192,2073,1872)
## replaced NA with 1750 because STAN cannot handle NAs in the y var!!!!!
## JAGS estimated missing 1999 nester abundance == mean=1747 sd=297 2.5%=1253 median=1712 97.5%=2434
mda <- c(863,1361,1712,2101,1750,1336,3150,1826,2166,1939,1619,2314,1609,2176,2227,2298,2366,2192,2073,1872)
year <- seq(from=1995,to=2014)

## JAGS estimated missing 1999 nester abundance == mean=1747 sd=297 2.5%=1253 median=1712 97.5%=2434

## Bundle data
mda_data <- list(y=log(mda),T=length(year),pyears=pyears)

## Parameters monitored
params <- c("r", "mean_r", "sigma2_obs", "sigma2_proc",
            "N_est")

# MCMC settings
ni <- 200000
nt <- 6
nb <- 100000
nc <- 3

## Initial values
inits <- lapply(1:nc, function(i) {
    list(sigma_proc = runif(1, 0, 1),
         mean_r = rnorm(1),
         sigma_obs = runif(1, 0, 1))})

## Call Stan from R
mda.ssm <- stan(model_code=mda.model.code,
               data = mda_data, init = inits, pars = params,
               chains = nc, thin = nt, iter = ni, warmup = nb,
               seed = 1,
               control = list(adapt_delta = 0.999),
               open_progress = FALSE)


## Summarize posteriors
print(mda.ssm, probs=c(0.025,0.50,0.975), digits = 3)

## parameter      mean    sd       2.5%     50%     97.5%
## mean_r         0.031   0.042   -0.051    0.029    0.125
## sigma2_obs     0.046   0.028    0.004    0.043    0.114
## sigma2_proc    0.032   0.037    0.000    0.019    0.134

## extract estimated params
out<-extract(mda.ssm)
Nt<-out$N_est
nsims = length(out$sigma)

## subset of output …
trend<-summary(mda.ssm)$summary[23:45,c(1,4,8)]  ## extract estimated mean and CRIs for each year

[paul-buerkner]

Thanks for the detailed explanations! I will dive deeper into this model within the next week or so. As soon as I understand it thoroughly, I can tell you more if / how we can make this possible with brms.

[paul-buerkner]

Nearly one year ago, I promised to respond in a week or so, which apparently I did not keep. Sorry for that.

I have taken a closer look at such models, and while I really like them, I see the following problems when trying to implement them in brms:

(1) The models require sampling of discrete parameters such as the latent abundance, which Stan cannot naturally handle. Their are at least two options, which are both discussed in the forum post you linked to above: Marginalizing over the discrete parameters and approximation through continuous distribution (e.g. via a normal distribution on the log-scale as in your example above). I currently don't see marginalizing to be a feasible option for these models. Using continuous distributions as approximations would certainly be an option, but we need to figure out, which approximations are most generally applicable to such models. Possibly, there are already papers out there discussing this topic.

(2) The models are very specific. We need to write new family functions, for instance binomial_poisson for a N-mixture model, with binomial observered data and a poisson distribution of the latent abundance (or some continuous approximation of it). In addition, we need a new correlation structure (for argument autocor) to specify time and site variables, which will only work with these new families. This would be admittedly be a lot of work, but in general something that would be possible from the brms side. However, my concern is that such an implemention would not really satisfy anyone. Because these models are very specific, it would require much special case coding, which made brms much harder to maintain in the future. Given its current flexibility this is already a non-trivial task to date. Also the rather basic models, which could make it into brms, will likely not satisfy most researchers: They will expect and require some of the more complex forms, for instance those with colonization and extinction parameters. The number of options in unmarked::pcountOpen provides a clear picture of how complex such models can get.

In summary, I do not think implementing N-mixture models in brms would be a good idea. Rather, one should think of a separate R package (or an extention to unmarked) focusing just on such models in a Bayesian framework fitted via JAGS or Stan behind the scences. With regard to Stan, there are already quite a few R packages on CRAN, which do a great job focusing on narrow model classs suited only for a specific research setting, but providing everything users expect for such setting.

[MilaniC]

Thanks Paul We had come to a similar conclusion sometime ago. But appreciate your thoughtful response. Wishing you all the best for the festive season. Cheers Milani

…

Sent from my iPhone

On 24 Dec 2017, at 10:06 am, Paul-Christian Bürkner ***@***.***> wrote: Nearly one year ago, I promised to respond in a week or so, which apparently I did not keep. Sorry for that. I have taken a closer look at such models, and while I really like them, I see the following problems when trying to implement them in brms: (1) The models require sampling of discrete parameters such as the latent abundance, which Stan cannot naturally handle. Their are at least two options, which are both discussed in the forum post you linked to above: Marginalizing over the discrete parameters and approximation through continuous distribution (e.g. via a normal distribution on the log-scale as in your example above). I currently don't see marginalizing to be a feasible option for these models. Using continuous distributions as approximations would certainly be an option, but we need to figure out, which approximations are most generally applicable to such models. Possibly, there are already papers out there discussing this topic. (2) The models are very specific. We need to write new family functions, for instance binomial_poisson for a N-mixture model, with binomial observered data and a poisson distribution of the latent abundance (or some continuous approximation of it). In addition, we need a new correlation structure (for argument autocor) to specify time and site variables, which will only work with these new families. This would be admittedly be a lot of work, but in general something that would be possible from the brms side. However, my concern is that such an implemention would not really satisfy anyone. Because these models are very specific, it would require much special case coding, which made brms much harder to maintain in the future. Given its current flexibility this is already a non-trivial task to date. Also the rather basic models, which could make it into brms, will likely not satisfy most researchers: They will expect and require some of the more complex forms, for instance those with colonization and extinction parameters. The number of options in unmarked::pcountOpen provides a clear picture of how complex such models can get. In summary, I do not think implementing N-mixture models in brms would be a good idea. Rather, one should think of a separate R package (or an extention to unmarked) focusing just on such models in a Bayesian framework fitted via JAGS or Stan behind the scences. With regard to Stan, there are already quite a few R packages on CRAN, which do a great job focusing on narrow model classs suited only for a specific research setting, but providing everything users expect for such setting. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub, or mute the thread.