Hey!

I’m trying to reproduce the model described in https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/1365-2664.12710 (appendix contains more complete formulations, see https://besjournals.onlinelibrary.wiley.com/action/downloadSupplement?doi=10.1111%2F1365-2664.12710&file=jpe12710-sup-0001-AppendixS1.pdf).

In doing so, I’m constructing a couple of covariance matrices (one for a spatial random effect, and another for spatiotemporal random effect).

When these are introduced, I consistently get:

Rejecting initial value:

Gradient evaluated at the initial value is not finite.

Stan can’t start sampling from this initial value.

I’ve tried adding some small noise to the diagonals and I’ve tried using small non-zero initial values.

Could some underflow occur that then creates zeros where they shouldn’t be, leading to exploding gradients?

In the below model I’ve removed the spatiotemporal effect as it demonstrates the issue sufficiently as it is.

```
data {
int<lower=0> N; // total samples
matrix[N, 2] S; // site coordinates
vector[N] m; // number of tows needed for a specific reef at a certain year
vector[N] x_inner; // indicator whether the reef is inner-shelf
vector[N] x_outer; // indicator whether the reef is outer-shelf
vector[N] tau_pre04; // time (years) since pre-2004 closure
vector[N] tau_04; // time (years) since post-2004 closure
int<lower=0> y[N]; // COTS count
}
transformed data {
real beta_sig = 10;
}
parameters {
real<lower=0> r; // overdispersion
vector<lower=0>[2] l; // length-scale parameters for spatial correlation
real<lower=0> sigma_phi; // magnitude parameter for spatial random effects
real<lower=0> sigma_g; // magnitude parameter for temporal random effects
real<lower=0> a_pre04; // half-saturation constant for pre-2004 closure
real<lower=0> a_04; // half-saturation constant for post-2004 closure
vector[N] f;
}
model {
// priors
a_pre04 ~ student_t(4, 0, 100);
a_04 ~ student_t(4, 0, 100);
l ~ student_t(4, 0, 10^3);
r ~ gamma(4, 0.1);
// covariance matrix
matrix[N, N] Sigma_phi;
for (i in 1:N) {
for (j in 1:N) {
Sigma_phi[i, j] = sigma_phi*exp(-sqrt((((S[i, 1]-S[j, 1])/l[1])^2) + ((S[i, 2]-S[j, 2])/l[2])^2));
}
Sigma_phi[i, i] += 0.0001;
}
vector[N] H_pre04 = tau_pre04 ./ (a_pre04 + tau_pre04);
vector[N] H_04 = tau_04 ./ (a_04 + tau_04);
f ~ multi_normal(rep_vector(0, N), beta_sig*identity_matrix(N) +
beta_sig*x_inner*x_inner' +
beta_sig*x_outer*x_outer' +
beta_sig*H_pre04*H_pre04' +
beta_sig*H_04*H_04' +
Sigma_phi);
y ~ neg_binomial_2(m .* exp(f), r);
}
```