In multilevel models for binary responses, estimation is computationally challenging due to the need to evaluate intractable integrals. In this paper, we investigate the performance of a recently proposed Bayesian method for deterministic fast approximate inference, called Integrated Nested Laplace Approximation (INLA). In particular, we conducted a simulation study, comparing the results obtained via INLA with the results obtained via MCMC, i.e. the traditional estimation method in the Bayesian context, and via maximum likelihood with adaptive quadrature. Particular attention is devoted to the case of small sample size and to the specification of the prior distribution for the variance.