A comparison of Bayesian and likelihood-based methods for fitting multilevel models

We use simulation studies, whose design is realistic for educational and medical research (as well as other fields of inquiry), to compare Bayesian and likelihood-based methods for fitting variance-components (VC) and random-effects logistic regression (RELR) models. The likelihood (and approximate likelihood) approaches we examine are based on the methods most widely used in current applied ultilevel (hierarchical) analyses: maximum likelihood (ML) and restricted ML (REML) for Gaussian outcomes, and marginal and penalized quasi-likelihood MQL and PQL) for Bernoulli outcomes. Our Bayesian methods use Markov chain Monte Carlo (MCMC) estimation, with adaptive hybrid Metropolis-Gibbs sampling for RELR models, and several diffuse prior distributions (Gamma(-1)(epsilon,epsilon) and U(0,1/epsilon) priors for variance components). For evaluation criteria we consider bias of point estimates and nominal versus actual coverage of interval estimates in repeated sampling. In two-level VC models we find that (a) both likelihood-based and Bayesian approaches can be made to produce approximately unbiased estimates, although the automatic manner in which REML accomplishes this is an dvantage, but (b) both approaches had difficulty achieving nominal coverage in small samples and with small values of the intraclass correlation. With the three level RELR models we examine we find that (c) quasi-likelihood methods for estimating random-effects variances perform badly with respect to bias and coverage in the example we simulated, and (d) Bayesian diffuse-prior methods lead to well-calibrated point and interval RELR estimates. While it is true that the likelihood-based methods we study are considerably faster computationally than MCMC, (i) steady improvements in recent years in both hardware speed and efficiency of Monte Carlo algorithms and (ii) the lack of calibration of likelihood-based methods in some common hierarchical settings combine to make MCMC-based Bayesian fitting of multilevel models an attractive approach, even with rather large data sets. Other analytic strategies based on less approximate likelihood methods are also possible but would benefit from further study of the type summarized here.