model {
  # priors
  # SD's for simple random effects (detection: p)
  p.spp.sd ~ dt(0, 1, 1)T(0,) # SD for spp effects
  p.spp.tau <- pow(p.spp.sd, -2) # convert to precision
  # sites are uniquely nested within lakes, no context as simple RE
  
  # detection probability per sample run, given DNA is present
  p.rep ~ dunif(0, 1)
  
  # SDs and precision for complex random effects (all presumed spp-specific)
  p.spp.lake.sd ~ dt(0, 1, 1)T(0,) # spp-lake effect
  p.spp.lake.tau <- pow(p.spp.lake.sd, -2) 
  p.spp.time.sd ~ dt(0, 1, 1)T(0,) # spp-time effect
  p.spp.time.tau <- pow(p.spp.time.sd, -2) 
  p.spp.lake.site.sd ~ dt(0, 1, 1)T(0,) # SD for site effects within lakes within spp
  p.spp.lake.site.tau <- pow(p.spp.lake.site.sd, -2) 


  # revised so that species share common distributions of coefficients
  a.jd.mu ~ dnorm(0, 1/3)
  a.jd2.mu ~ dnorm(0, 1/3)
  a.depth.mu ~ dnorm(0, 1/3)
  a.strat.mu ~ dnorm(0, 1/3)
  a.clarity.mu ~ dnorm(0, 1/3)
  a.jd.sd ~ dunif(0, 1); a.jd.tau <- pow(a.jd.sd, -2)
  a.jd2.sd ~ dunif(0, 1); a.jd2.tau <- pow(a.jd2.sd, -2)
  a.depth.sd ~ dunif(0, 1); a.depth.tau <- pow(a.depth.sd, -2)
  a.strat.sd ~ dunif(0, 1); a.strat.tau <- pow(a.strat.sd, -2)
  a.clarity.sd ~ dunif(0, 1); a.clarity.tau <- pow(a.clarity.sd, -2)
  
  # species specific priors 
  for (i in 1:n.spp){
    p.spp[i] ~ dnorm(0, p.spp.tau) # spp mean detection (logit scale)
    
    # covariate effects
    a.jd[i] ~ dnorm(a.jd.mu, a.jd.tau)
    a.jd2[i] ~ dnorm(a.jd2.mu, a.jd2.tau)
    a.depth[i] ~ dnorm(a.depth.mu, a.depth.tau)
    a.strat[i] ~ dnorm(a.strat.mu, a.strat.tau)
    a.clarity[i] ~ dnorm(a.clarity.mu, a.clarity.tau)
    # species-lake effect
    for (j in 1:n.lakes){
      p.spp.lake[i,j] ~ dnorm(0, p.spp.lake.tau) # lake-by-spp adjustments
      for (s in 1:n.sites){
        p.spp.lake.site[i,j,s] ~ dnorm(p.spp.lake[i,j], p.spp.lake.site.tau) # spp-lake-site
      }
    }
  }
  
  # Ecological model for true detection
  # Observation model for replicated detection/nondetection observations
  for (i in 1:n.obs) {
    # sample.p[i] can be envisioned as sample level occupancy
    logit(sample.p[i]) <- p.spp[spp[i]] + 
      p.spp.lake[spp[i], lake[i]] + 
      p.spp.lake.site[spp[i], lake[i], site[i]] + 
      a.jd[spp[i]] * jd[i] + a.jd2[spp[i]] * jd2[i] + 
      a.clarity[spp[i]] * clarity[i] +
      a.strat[spp[i]] * mixed_perc[i] +
      a.depth[spp[i]] * depth[lake[i], site[i]]
      rep.p[i] <- 1 - pow(1 - sample.p[i], 1/3) # convert to per replicate                              
    for (k in 1:n.reps){
      # note: 0 trials if not there
      # rep.p[i] will include modeled variation due to spp, lake, site
      y[i,k] ~ dbinom(rep.p[i], occupancy[lake[i],spp[i]]) 
      
      # GOF
      ynew[i, k] ~ dbern(rep.p[i])
      presi[i,k] <- (y[i, k] - rep.p[i]) / sqrt(rep.p[i] * (1 - rep.p[i]))
      presi.new[i,k] <- (ynew[i, k] - rep.p[i]) / sqrt(rep.p[i] * (1 - rep.p[i]))
      D[i,k] <- pow(presi[i,k], 2)
      D.new[i,k] <- pow(presi.new[i,k], 2)
    }
  } # close n
  # Discrepnacy measure
  fit <- sum(D[1:n.obs, 1:n.reps])
  fit.new <- sum(D.new[1:n.obs, 1:n.reps])
} # end jags model