data {
  int<lower = 1> N_obs;                 // Number of observations
  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
  int<lower = 0, upper = 1> Treat[N_obs]; // explanatory variable
  real<lower = 0> Height[N_obs]; // explanatory variable
  real<lower = 0> Dist[N_obs]; // explanatory variable
}
parameters {
  vector[6] 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 {
  for (i in 1:N_obs) {
    real lambda = exp(beta[1] + beta[2] * Treat[i] + beta[3] * Height[i]
                       + beta[4] * Dist[i]
                       + epsilon_S[S[i]] + epsilon_B[S[i], B[i]]
                       + epsilon_Q[S[i], B[i], Q[i]]);
    real omega = inv_logit(beta[5] + beta[6] * Dist[i]);
    if (Y[i] == 0) {
      vector[2] lp;
      lp[1] = bernoulli_lpmf(0 | omega);
      lp[2] = bernoulli_lpmf(1 | omega)
              + poisson_lpmf(0 | lambda);
      target += log_sum_exp(lp);
    } else {
      target += bernoulli_lpmf(1 | omega)
                + poisson_lpmf(Y[i] | 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);
}