Hi. It’s my first time posting here, and I’m just getting started using Stan. If someone can help me figure this out it would be a big help. I’m working on converting NIMBLE code to Stan.
I have a fairly simple twostep hurdle model of bird abundance. Model 1 models the probability that the species occurred during that observation, so it predicts an occurrence probability for each observation which I (stochastically) convert to a binary (0,1) estimate using occurred ~ dbern(p1). Then, Model 2 models abundance given occurrence, where count ~ dpois(lambda * occurred). So the models are connected in that way. Sitelevel variation is the biggest source of variation, so I include sitelevel random effects.
I can figure out converting this to Stan except for two things:

The
occurred ~ dbern(p1)
piece that I mention. Occurred ends up being a generated quantity (I think? Or something like that?) and so I don’t know how to code it in a way for Stan to understand it’s generated but also intermediate between the two submodels. 
My primary inferential objective is the derived quantity, which is at the bottom, and is the annual estimate of typical abundance, which is basically
avg_abd[t] = lambda2[t] * p1[t]
. I think this would go undertransformed parameters
? Is that correct? Or is this also a generated quantity?
I also am unsure how sumtozero random effects are handled in Stan.
Many thanks for any and all advice and guidance. Again, this would be a really big help.
#  Data 
dataList < list(
year = df1$year,
occurrence = df1$occurrence,
count = df1$count
)
#  Constants 
constList < list(
N1 = 8436,
S = 30,
nyears = 33,
years = sort(unique(df1$year)),
site = df1$site
)
#  Model code 
code < nimbleCode({
#  Priors 
# Model 1: Occurrence
int_m1 ~ dnorm(0, 0.01)
beta_m1_year ~ dnorm(0, 0.01)
# random site intercepts (sumtozero)
sigma_m1 ~ dunif(0, 10)
for (s in 1:(S1)) {
eps_site_m1[s] ~ dnorm(0, sd = sigma_m1)
}
eps_site_m1[S] < sum(eps_site_m1[1:(S1)])
# Model 2: Count
int_m2 ~ dnorm(0, 0.001)
beta_m2_year ~ dnorm(0, 0.001)
# random site intercepts (sumtozero)
sigma_m2 ~ dunif(0, 10)
for (s in 1:(S1)) {
eps_site_m2[s] ~ dnorm(0, sd = sigma_m2)
}
eps_site_m2[S] < sum(eps_site_m2[1:(S1)])
#  Likelihoods 
for(i in 1:N1) {
#  MODEL 1 
# Likelihood
occurrence[i] ~ dbern(p1[i])
logit(p1[i]) < int_m1 + (beta_m1_year * year[i]) + eps_site_m1[site[i]]
# Prediction from Model 1 to use in count model  this is a "generated quantity"?
occurred[i] ~ dbern(p1[i])
#  MODEL 2 
# Likelihood
count[i] ~ dpois(lambda2[i] * occurred[i])
log(lambda2[i]) < int_m2 + (beta_m2_year * year[i]) + eps_site_m2[site[i]]
}
#  Derived quantities 
# Typical annual abundance
for(t in 1:nyears){
logit(p[t]) < int_m1 + (beta_m1_year * years[t])
log(lambda[t]) < int_m2 + (beta_m2_year * years[t])
avg_abd[t] < p[t] * lambda[t]
}
})