Multistate models provide a useful framework for modelling complex event history data in clinical settings and have recently been extended to the joint modelling framework to appropriately handle endogenous longitudinal covariates, such as repeatedly measured biomarkers, which are informative about health status and disease progression. However, the practical application of such joint models faces considerable computational challenges. Motivated by a longitudinal multimorbidity analysis of large-scale UK health records, we introduce novel Bayesian inference approaches for these models that are capable of handling complex multistate processes and large datasets with straightforward implementation. These approaches decompose the original estimation task into smaller inference blocks, leveraging parallel computing and facilitating flexible model specification and comparison. Using extensive simulation studies, we show that the proposed approaches achieve satisfactory estimation accuracy, with notable gains in computational efficiency compared to the standard Bayesian estimation strategy. We illustrate our approaches by analysing the coevolution of routinely measured systolic blood pressure and the progression of three important chronic conditions, using a large dataset from the Clinical Practice Research Datalink Aurum database. Our analysis reveals distinct and previously lesser-known association structures between systolic blood pressure and different disease transitions.
