The rise in availability of electronic health record data enables us to answer more detailed clinical questions; however, the associated increased complexity raises substantial statistical and computational challenges. Recent work in the area of joint models has introduced an extended mixed effects framework, encompassing multiple outcomes of any type, each of which could be repeatedly measured (longitudinal), with any number of levels, and with any number of random effects at each level (Crowther, 2020). This allows for sharing and linking between outcome models in an extremely flexibly way, either by linking random effects directly, or the expected value of one outcome (or function of it) within the linear predictor of another. Non-linear and time-dependent effects are also seamlessly incorporated to the linear predictor through the use of splines or fractional polynomials. In this talk, I’ll present an extension to the framework to further allow modelling of variance components directly, allowing the residual level 1 variance to have its own complex linear predictor, allowing for heteroskedasticity, which in turn provides new tools for joint modelling. Throughout my talk I will illustrate an accompanying user-friendly implementation in Stata, showing how to build and estimate a joint longitudinal-survival model with complex variance components, quantifying how between-subject variation in the level 1 variance structure of a continuous biomarker (e.g., blood pressure), can be associated with survival. Dynamic predictions from such a joint model will also be derived and presented. Due to the generality of the implementation, multiple outcomes, such as multiple biomarkers or competing risks, are also immediately available.