model {

 # Priors for random effects precisions
 yearInterceptSD ~ dunif(0.01, 5)
 plantInterceptSD ~ dunif(0.01, 5)
 plantSlopeSD ~ dunif(0.01, 5)
 ## plantSlopeSD <- 1/sqrt(plantSlopePrecision)

 yearInterceptPrecision <- 1/(yearInterceptSD * yearInterceptSD)
 plantInterceptPrecision <- 1/(plantInterceptSD * plantInterceptSD)
 plantSlopePrecision <- 1/(plantSlopeSD * plantSlopeSD)
 ## plantSlopePrecision ~ dgamma(0.0001, 0.0001) 

 # Priors for fixed effects (intercepts) of previous stage
 ## WE COULD CONSIDER OTHER OPTIONS HERE...
 for(i in 1:Nstage) { 
  intercept[i] ~ dnorm(0, 1.0E-6)
 }

 # Prior for slope
 ## ... AND HERE
 slope ~ dnorm(0, 1.0E-6)

 ## random year effects:
 for (i in 1:Nyear) {
   yearInterceptEffect[i] ~ dnorm(0, yearInterceptPrecision)
 }

 ## random plant effects:
 for (i in 1:Nplant) {
   plantInterceptEffect[i] ~ dnorm(0, plantInterceptPrecision)
   plantSlopeEffect[i] ~ dnorm(0, plantSlopePrecision)
 }

 ## logit model and data:
 for(i in 1:Ndata) {
   logit(probF[i]) <- intercept[ stage[i] ] + yearInterceptEffect[ year[i] ] + plantInterceptEffect[ plant[i] ] + slope * Pods[i] + plantSlopeEffect[ plant[i] ] * Pods[i]
   toF[i] ~ dbern(probF[i])
 }
 
}