This could be done relatively straightforwardly for the logitnormal (and cloglognormal...) case by using standard random-effects style additive normal random effects for each level.
For two levels there would be just one additional unknown (with the $\sigma$ replaced with $\sigma_1$ and $\sigma_2$ for the two levels of clustering$). There would be only one additional set of integrals (compared to the existing two) however, they would be 2-d integrals.
Could use something like cubature to do the multidimensional integration efficiently and accurately without resorting to monte carlo.
This could be done relatively straightforwardly for the logitnormal (and cloglognormal...) case by using standard random-effects style additive normal random effects for each level.
For two levels there would be just one additional unknown (with the $\sigma$ replaced with $\sigma_1$ and $\sigma_2$ for the two levels of clustering$). There would be only one additional set of integrals (compared to the existing two) however, they would be 2-d integrals.
Could use something like cubature to do the multidimensional integration efficiently and accurately without resorting to monte carlo.