In this paper we propose a four-part random effects model, with application to correlated medical cost data. Four joint equations are used to model respectively: (1) the probability of seeking medical treatment, (2) the probability of being hospitalized (conditional on seeking medical treatment), and the actual amount of (3) outpatient and (4) inpatient costs. Our model simultaneously takes account of the inter-temporal (or within-cluster) correlation of each patient and the cross-equation correlation of the four equations, by means of joint linear mixed models and generalized linear mixed models. The estimation is accomplished by the high-order Laplace approximation technique in Raudenbush et al. [Raudenbush, S.W., Yang, M., Yosef, M., 2000. Maximum likelihood for generalized linear models with nested random effects via high-order, multivariate Laplace approximation. Journal of Computational and Graphical Statistics 9, 141-157] and Olsen and Schafer [Olsen, M.K., Schafer, J.L., 2001. A two-part random effects model for semicontinuous longitudinal data. Journal of the American Statistical Association 96, 730-745]. Our model is used to analyze monthly medical costs of 1397 chronic heart failure patients from the clinical data repository (CDR) at the University of Virginia.