Procedures for analyzing and comparing healthcare providers' effects on health services delivery and outcomes have been referred to as provider profiling. In a typical profiling procedure, patient-level responses are measured for clusters of patients treated by providers that in turn, can be regarded as statistically exchangeable. Thus, a hierarchical model naturally represents the structure of the data. When provider effects on multiple responses are profiled, a multivariate model rather than a series of univariate models, can capture associations among responses at both the provider and patient levels. When responses are in the form of charges for healthcare services and sampled patients include non-users of services, charge variables are a mix of zeros and highly-skewed positive values that present a modeling challenge. For analysis of regressor effects on charges for a single service, a frequently used approach is a two-part model (Duan, Manning, Morris, and Newhouse 1983) that combines logistic or probit regression on any use of the service and linear regression on the log of positive charges given use of the service. Here, we extend the two-part model to the case of charges for multiple services, using a log-linear model and a general multivariate log-normal model, and employ the resultant multivariate two-part model as the within-provider component of a hierarchical model. The log-linear likelihood is reparameterized as proposed by Fitzmaurice and Laird (1993), so that regressor effects on any use of each service are marginal with respect to any use of other services. The general multivariate log-normal likelihood is constructed in such a way that variances of log of positive charges for each service are provider-specific but correlations between log of positive charges for different services are uniform across providers. A data augmentation step is included in the Gibbs sampler used to fit the hierarchical model, in order to accommodate the fact that values of log of positive charges are undefined for unused service. We apply this hierarchical, multivariate, two-part model to analyze the effects of primary care physicians on their patients' annual charges for two services, primary care and specialty care. Along the way, we also demonstrate an approach for incorporating prior information about the effects of patient morbidity on response variables, to improve the accuracy of provider profiles that are based on patient samples of limited size.
