A few approaches have been proposed to learn healthcare delivery network using electronic health records (EHR) data, that is to infer the relationship among different types of medical encounters (e.g. disease diagnosis, treatments, and procedures), which might in turn improve the healthcare system. However, these approaches do not fully account for the properties of the EHR data, including the longitudinal nature and patient heterogeneity. In this article, we propose a flexible model, built upon multivariate Hawkes process, for HDN construction. Our model allows for patient-specific time-varying background intensity functions via random effects, which can also adjust for effects of important covariates. We adopt a penalized approach to select fixed effects, yielding a sparse network structure, and to remove unnecessary random effects from the model. Through extensive simulation studies and a dataset of type 2 diabetes, we show that our proposed method performs well in recovering the network structure and that it is essential to account for patient heterogeneities.