A phase-type distribution is the distribution of the time until absorption in a finite state-space time-homogeneous Markov jump process, with one absorbing state and the rest being transient. These distributions are mathematically tractable and conceptually attractive to model physical phenomena due to their interpretation in terms of a hidden Markov structure. Three recent extensions of regular phase-type distributions give rise to models which allow for heavy tails: discrete- or continuous-scaling; fractional-time semi-Markov extensions; and inhomogeneous time-change of the underlying Markov process. In this paper, we present a unifying theory for heavy-tailed phase-type distributions for which all three approaches are particular cases. Our main objective is to provide useful models for heavy-tailed phase-type distributions, but any other tail behavior is also captured by our specification. We provide relevant new examples and also show how existing approaches are naturally embedded. Subsequently, two multivariate extensions are presented, inspired by the univariate construction which can be considered as a matrix version of a frailty model. We provide fully explicit EM-algorithms for all models and illustrate them using synthetic and real-life data.