Multistage models have been effective at describing length of stay (LOS) distributions for diverse patient groups. Our study objective was to determine whether such models could be used for patient groups restricted by diagnosis, severity of illness, or hospital in order to facilitate comparisons conditioned on these factors. We performed a retrospective cohort study using data from 317,876 hospitalizations occurring over 2 years in 17 hospitals in a large, integrated health care delivery system. We estimated model parameters using data from the first year and validated them by comparing the predicted LOS distribution to the second year of data. We found that 3- and 4-stage models fit LOS data for either the entire hospital cohort or for subsets of patients with specific conditions (e.g. community-acquired pneumonia). Probability distributions were strongly influenced by the degree of physiologic derangement on admission, pre-existing comorbidities, or a summary mortality risk combining these with age, sex, and diagnosis. The distributions for groups with greater severity of illness were shifted slightly to the right, but even more notable was the increase in the dispersion, indicating the LOS is harder to predict with greater severity of illness. Multistage models facilitate computation of the hazard function, which shows the probability of imminent discharge given the elapsed LOS, and provide a unified method of fitting, summarizing, and studying the effects of factors affecting LOS distributions. Future work should not be restricted to expected LOS comparisons, but should incorporate examination of LOS probability distributions.