Electronic health records (EHRs) systematically capture patients' health information in a digital format. To accurately characterize individual disease comorbidity and susceptibility, it is necessary to integrate information of biomarkers over time. However, due to observational data nature, the biomarkers measured over time are multivariate with mixed types. In addition, for each patient, the biomarker data are sparsely measured and the measurement times are irregular and potentially informative. In this paper, we propose multivariate generalized linear models to model different types of biomarkers. We also incorporate latent multivariate Gaussian processes to account for between-biomarker dependence over time. We allow the covariate effects and covariance matrix of the latent processes to be time-dependent. For parameters' estimations, we apply kernel-weight estimating equations based on the first and second moments of the assumed biomarker processes, where our kernel weights account for the heterogeneous intensity of the measurement times. We illustrate our method through extensive simulation studies and then apply our method to analyze an EHR for Type 2 Diabetes (T2D) patients.