Hi all,
I am trying to fit a hierarchical binomial model to some spider foraging data that I have and am running into some issues with divergent transitions and low bulk effective and tail effective sample sizes.
For some background, I am trying to fit spider feeding rates (or functions describing how spider feeding rates change other variables) to feeding surveys of spiders eating midges on walls. The way this method works is by combining the fraction of feeding spiders during a survey and the time over which feeding events are detectable to estimate the feeding rate. Skipping some details, the proportion of spiders feeding in a survey p should be equal to the product of the feeding rate of the spiders during that survey f and the detection time of spiders feeding on their prey d, or
p = fd.
Now, we want to be able to use a whole series of surveys to be able to estimate the parameters of a function describing how f changes with other variables such as midge densities. The classic function for this is
f = \frac{aR}{1 + ahR}
where a is the so-called space clearance rate of the predator and h is the handling time which for our convenience here we will assume is equal to the detection time, d, and R is the resource density. For each survey, we have estimated the detection times using data from another experiment and we also have estimates of the resource densities in each survey.
Finally, to the actual model. Originally, we simply fit the following the model
y_{i} \sim Binomial(n = n_{t,i} , p_{i} = \frac{aR_{i}}{1 + ad_{i}R_{i}}d_{i})
where y_{i} is the number of spiders feeding in survey i and n_{t,i} is the total number of spiders in survey i. Fitting this model everything worked fine.
However, an astute reviewer noticed that we had surveyed certain walls multiple times and suggested that we account for this non-independence by using a hierarchical model. I figured this wouldn’t be a big deal and tried to fit the following model
y_{ji} \sim Binomial(n = n_{t,ji}, p_{i} = \frac{a_{j}R_{ji}}{1 + a_{j}d_{ji}R_{ji}}d_{ji})
a_{j} \sim Normal(\mu_{a}, \sigma_{a})
where y_{ji} is now the number of feeding spiders on wall j in survey i, a_{j} is wall-specific space clearance rate for wall j, and so on. \mu_{a} and \sigma_{a} are hyperpriors describing desitrbution of wall-specific space clearance rates, a.
Now, however, I get a bunch of warnings about divergent transitions and low bulk effective and tail effective sample sizes. I’m looking for any tips on how to get rid of these issues. This is likely the simplest of the models I will fit to the data since I will also be fitting models in which, for example, the a parameter is temperature-dependent. My worry is that if I’m running into issues with the hierarchical models already, fitting those more complex models may not be possible.
Below, I give my R code to fit the model to the data, the Stan code to define the model, and a link to the .csv file containing the data on GitHub,
R Code:
### load required packages
library(dplyr); library(rstan);
### set some options for Stan
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
### load data
spider <- read.csv('FeedingData_Manipulated_f.csv', stringsAsFactors = TRUE)
### drop surveys with no midges
spider <- spider %>% filter(NumberMidges > 0)
### set up data for Stan
notemp_noint_data <- list(N = nrow(spider),
y = spider$NumberFeed,
n_t = spider$NumberPred,
K = length(unique(spider$BuildingWall)),
d = spider$DetectionEst,
R = spider$NumberMidges,
j = as.numeric(spider$BuildingWall))
notemp_noint_fit <- stan(file = 'Female_notemp_BDnone_hierarchical.stan', data = notemp_noint_data)
The contents of the file ‘Female_notemp_BDnone_hierarchical.stan’
// Stan model to fit type II FR
//
// Input data:
data {
int<lower=0> N; // Number of surveys
int<lower=0> y[N]; // Number of spiders feeding
int<lower=0> n_t[N]; // Number of spiders observed
int<lower=1> K; // Number of different building/wall combinations
vector[N] d; // Detection times
vector[N] R; // Prey densities
int<lower=1> j[N]; // building/wall index for each observation
}
// The parameters accepted by the model.
parameters {
real<lower = 0> aj[K];
real<lower = 0> a_mean;
real<lower = 0> a_sd;
}
transformed parameters {
vector[N] f;
vector<lower = 0, upper = 1>[N] p;
for (i in 1:N) {
f[i] = (aj[j[i]] * R[i])/(1 + (aj[j[i]] * R[i] * d[i]));
p[i] = f[i] * d[i];
}
}
// The model to be estimated. y is binomial with the probability
// equal to the feeding rate times the detection time
model {
// Likelihood
for (i in 1:N) {
y[i] ~ binomial(n_t[i], p[i]);
}
// Need priors and hyper priors
for (k in 1:K) {
aj[k] ~ normal(a_mean, a_sd);
}
a_mean ~ normal(30, 15);
a_sd ~ normal(0, 10);
}
generated quantities {
vector[N] log_lik;
for (n in 1:N){
log_lik[n] = binomial_lpmf(y[n] | n_t[n], (aj[j[n]]*R[n])/(1 + aj[j[n]]*d[n]*R[n]) * d[n]);
}
}
Link to the .csv file containing the data on GitHub: ZebraSpiderFR/FeedingData_Manipulated_f.csv at main · KyleCoblentz/ZebraSpiderFR · GitHub
I really appreciate any tips or help that anyone can offer. Thank you in advance and please let me know if there are any questions you have.