Joint models are an increasingly popular way to characterise the relationship between one or more longitudinal responses and an event of interest. However, for multivariate joint models the increased dimensionality and complexity of random effects present in the model specification are commensurate with increased computing time, hampering the implementation of many classic approaches. An approximate EM algorithm which ameliorates the so-called ‘curse of dimensionality’ is developed. The scaleability and accuracy of the proposed method are demonstrated via two simulation studies and applied to data arising from two clinical trials in the disease areas of cirrhosis and Alzheimer's disease, each with three biomarkers.