data {
  int<lower = 1> N_obs;                 // Number of observations
  int<lower = 1> N_var;                 // Number of variables including intercept
  int<lower = 1> N_S;                   // Number of sites
  int<lower = 1, upper = N_S> S[N_obs]; // Site
  int<lower = 1> N_B;                   // Maximum number of blocks
  int<lower = 1, upper = N_B> B[N_obs]; // Index of blocks
  int<lower = 1> N_Q;                   // Maximum number of quadrats
  int<lower = 1, upper = N_Q> Q[N_obs]; // Index of quadrats
  int<lower = 0> Y[N_obs];              // Observations
  matrix[N_obs, N_var] X;               // explanatory variables
}
parameters {
  vector[N_var] beta;         // coefficients
  vector<lower = 0>[3] sigma; // scale parameters
  real epsilon_S[N_S];
  real epsilon_B[N_S, N_B];
  real epsilon_Q[N_Q, N_B, N_Q];
}
model {
  vector[N_obs] lambda;
  for (i in 1:N_obs)
    lambda[i] = exp(X[i, ] * beta + epsilon_S[S[i]] + epsilon_B[S[i], B[i]]
                      + epsilon_Q[S[i], B[i], Q[i]]);
  Y ~ poisson(lambda);

  // Random effects
  for (s in 1:N_S) {
    epsilon_S[s] ~ normal(0, sigma[1]);
    for (b in 1:N_B) {
      epsilon_B[s, b] ~ normal(0, sigma[2]);
        for (q in 1:N_Q)
          epsilon_Q[s, b, q] ~ normal(0, sigma[3]);
    }
  }

  // Priors
  beta ~ normal(0, 100);
  sigma ~ normal(0, 10);
}