In recent years many studies have reported large differences in the use of medical treatments and procedures across geographic regions, hospitals, and individual health care providers. Beyond reporting on the extent of observed variations, these studies examine the role of contributing factors including patient, regional, and provider characteristics. In addition, they may assess the relation between health care processes and outcomes, such as patient mortality, morbidity, and functioning. Studies of variations in health care utilization and outcomes involve the analysis of multilevel clustered data; for example, data on patients clustered by hospital and/or geographic region. The goals of the analysis include the estimation of cluster-specific adjusted responses, covariate effects, and components of variance. The analytic strategy needs to account for correlations induced by clustering and to handle the presence of large variations in cluster size. In this article we formulate a broad class of hierarchical generalized linear models (HGLMs) and discuss their applications to the analysis of health care utilization data. The models can incorporate covariates at each level of the hierarchical data structure, can account for greater variation than what is allowed by the variance in a one-parameter exponential family, and permit the use of heavy-tailed distributions for the random effects. We develop a Bayesian approach to fitting HGLMs using Markov chain Monte Carlo methods and discuss several methods for model checking. The HGLM analysis is presented in the context of two examples of applications to the study of variations in the utilization of medical procedures for elderly Medicare beneficiaries who sustained a heart attack. The first example involves the analysis of clustered longitudinal data with binomial responses and examines geographic and temporal trends in the utilization of coronary angiography across the United States during the 4-year period 1987–1990. The second example involves the analysis of multilevel, clustered data with Poisson responses and examines hospital variations in the utilization of coronary artery bypass graft surgery in 1990. The HGLM analysis incorporates state-level and hospital-level covariates and makes it possible to estimate covariate effects and cluster-specific rates of utilization for both hospitals and states.