model{
  
  ## indexing
  # i = microcensus enumeration zone
  # t = residential type
  # k = regression coefficient
  
  for(i in 1:ni){
    
    # N = true number of people
    N[i] ~ dpois(D[i]*A[i])
    
    # D = true density of people
    D[i] ~ dlnorm(Dbar[i], pow(sigmaD[type[i],region[i],state[i],local[i]],-2))
    
    # Dbar = expected density
    Dbar[i] <- alpha[type[i],region[i],state[i],local[i]] + beta[1]*wpglobal[i] + beta[2]*schools[i] + beta[3]*hhsize[i] + beta[4]*settled[i] + beta[5]*residential[i] + beta[6]*nonresidential[i]
    
    # posterior predictions
    Nhat[i] ~ dpois(Dhat[i]*A[i])
    Dhat[i] ~ dlnorm(Dbar[i], pow(sigmaD[type[i],region[i],state[i],local[i]],-2))
  }
  
  # alpha = average population density (log-scale) for each settlement type per LGA (hierarchical by type/region/state/lga)
  mu_alpha_national ~ dnorm(0, 1e-3)
  sigma_alpha_national ~ dunif(0, 1e3)
  
  alpha0_national ~ dnorm(mu_alpha_national, pow(sigma_alpha_national, -2))
  
  for(t in 1:ntype){
    alpha0_type[t] ~ dnorm(mu_alpha_type[t], pow(sigma_alpha_type[t],-2))
    
    mu_alpha_type[t] ~ dnorm(mu_alpha_national, pow(sigma_alpha_national,-2))
    sigma_alpha_type[t] ~ dunif(0, sigma_alpha_national)
    
    for(r in 1:nregion){
      alpha0_region[t,r] ~ dnorm(mu_alpha_region[t,r], pow(sigma_alpha_region[t,r],-2))
      
      mu_alpha_region[t,r] ~ dnorm(mu_alpha_type[t], pow(sigma_alpha_type[t],-2))
      sigma_alpha_region[t,r] ~ dunif(0, sigma_alpha_type[t])
      
      for(s in 1:nstate[r]){
        alpha0_state[t,r,s] ~ dnorm(mu_alpha_state[t,r,s], pow(sigma_alpha_state[t,r,s],-2))
        
        mu_alpha_state[t,r,s] ~ dnorm(mu_alpha_region[t,r], pow(sigma_alpha_region[t,r],-2))
        sigma_alpha_state[t,r,s] ~ dunif(0, sigma_alpha_region[t,r])
        
        for(l in 1:nlocal[r,s]){
          alpha[t,r,s,l] ~ dnorm(mu_alpha_state[t,r,s], pow(sigma_alpha_state[t,r,s],-2))
        }
      }
    }
  }
  
  # beta = fixed effects of covariates
  for(b in 1:nbeta){
    beta[b] ~ dnorm(0, pow(5,-2))
  }
  
  # sigmaD = residual variance in population density (log-scale) for an LGA (hierarchical by type/region/state/lga)
  mu_sigmaD_national ~ dnorm(0, 1e-3) I(0,)
  sigma_sigmaD_national ~ dunif(0, 3)
  
  sigmaD0_national ~ dnorm(mu_sigmaD_national, pow(sigma_sigmaD_national, -2)) T(0,)
  
  for(t in 1:ntype){
    sigmaD0_type[t] ~ dnorm(mu_sigmaD_type[t], pow(sigma_sigmaD_type[t],-2)) T(0,)
    
    mu_sigmaD_type[t] ~ dnorm(mu_sigmaD_national, pow(sigma_sigmaD_national,-2)) T(0,)
    sigma_sigmaD_type[t] ~ dunif(0, sigma_sigmaD_national)
    
    for(r in 1:nregion){
      sigmaD0_region[t,r] ~ dnorm(mu_sigmaD_region[t,r], pow(sigma_sigmaD_region[t,r],-2)) T(0,)
      
      mu_sigmaD_region[t,r] ~ dnorm(mu_sigmaD_type[t], pow(sigma_sigmaD_type[t],-2)) T(0,)
      sigma_sigmaD_region[t,r] ~ dunif(0, sigma_sigmaD_type[t])
      
      for(s in 1:nstate[r]){
        sigmaD0_state[t,r,s] ~ dnorm(mu_sigmaD_state[t,r,s], pow(sigma_sigmaD_state[t,r,s],-2)) T(0,)
        
        mu_sigmaD_state[t,r,s] ~ dnorm(mu_sigmaD_region[t,r], pow(sigma_sigmaD_region[t,r],-2)) T(0,)
        sigma_sigmaD_state[t,r,s] ~ dunif(0, sigma_sigmaD_region[t,r])
        
        for(l in 1:nlocal[r,s]){
          sigmaD[t,r,s,l] ~ dnorm(mu_sigmaD_state[t,r,s], pow(sigma_sigmaD_state[t,r,s],-2)) T(0,)
        }
      }
    }
  }
}