model {

 # Priors for random effects precisions
 yearInterceptSD ~  dunif(0.0001, 5)
 plantInterceptSD ~ dunif(0.0001, 5)
 plantSlopeSD ~     dunif(0.0001, 5)

 ## When we use gamma priors for the precision, the lines above
 ## would be replaced by (for example):
 ## plantSlopeSD <- 1/sqrt(plantSlopePrecision)

 yearInterceptPrecision <- 1/(yearInterceptSD * yearInterceptSD)
 plantInterceptPrecision <- 1/(plantInterceptSD * plantInterceptSD)
 plantSlopePrecision <- 1/(plantSlopeSD * plantSlopeSD)

 ## ... and the above lines would be replaced by:
 ## 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])
 }
}