Joint models for longitudinal and survival data have become an established tool for optimally handling scenarios when both types of data co-exist. Multivariate extensions to the classic univariate joint model have started to emerge but are typically restricted to the Gaussian case, deployed in a Bayesian framework or focused on dimension reduction. An approximate EM algorithm is utilised which circumvents the oft-lamented curse of dimensionality and offers a likelihood-based implementation which ought to appeal to clinicians and practitioners alike. The proposed method is validated in a pair of simulation studies, which demonstrate both its accuracy in parameter estimation and efficiency in terms of computational cost. Its clinical use is demonstrated via an application to primary billiary cirrhosis data. The proposed methodology for estimation of these joint models is available in R package gmvjoint.