Causality-Consistent PH Model

• The causality-consistent PH model is a joint lifetime framework that uses time-inhomogeneous multivariate phase-type distributions to mediate dependency via a shared latent initial state.
• It provides closed-form joint survival and density functions through cumulative sub-intensity matrices and explicit time-inhomogeneity, enhancing causal interpretations.
• A custom EM-type estimation approach integrates multinomial regression for covariate-driven initial distributions while efficiently handling right-censored data.

The causality-consistent PH model is a joint lifetime modeling framework based on time-inhomogeneous multivariate phase-type (PH) distributions, specifically the multivariate inhomogeneous PH (mIPH) class as formulated by Albrecher, Bladt, and Müller (Hansjörg et al., 2022). This construction provides a direct causal interpretation of dependence in bivariate survival data, with all association mediated exclusively through a shared latent initial state. The mIPH approach contrasts with copula-based methods by imposing time-order and explicit mediation via covariates, ensuring full causal consistency.

1. Model Specification and State-Space Structure

The causality-consistent PH model operates on a common state-space $E=\{1,2,...,p,p+1\}$, where states $1$ to $p$ are transient and $p+1$ is absorbing. The observed random vectors $(X,Y)$, interpreted as lifetimes, are constructed as absorption times of two time-inhomogeneous pure-jump Markov chains $\{J_t^{(1)}\}_{t\ge0}$ and $\{J_t^{(2)}\}_{t\ge0}$. Both chains commence from a random initial state $J_0\sim\alpha=(\alpha_1,...,\alpha_p)$, with the initial distribution potentially varying across subjects.

For each chain $i=1,2$, the instantaneous evolution is governed by a sub-intensity matrix $T_i(t)\in\mathbb{R}^{p\times p}$, encoding transition rates between phases, and an exit vector $t_i(t) = -T_i(t)e$ denotes the rate of absorption to state $p+1$. The complete generator takes the form:

$$\Lambda_i(t) = \begin{pmatrix} T_i(t) & t_i(t)\ 0 & 0 \end{pmatrix}$$

with lifetimes defined as first hitting times of the absorbing state:

$$X = \inf\{t > 0 : J_t^{(1)} = p+1\},\quad Y = \inf\{t > 0 : J_t^{(2)} = p+1\}$$

This construction yields the joint law $$(X,Y)\sim \mathrm{mIPH}\bigl(\alpha, \{T_1(\cdot),T_2(\cdot)\}\bigr)$$, where both marginals are matrix-exponential distributions.

2. Joint Survival Functions and Density Formulation

The model admits closed-form expressions for the joint survival and density functions, conditional on the latent starting state. For $S_{X,Y}(x,y)$, the joint survival probability over times $x,y$:

$$S_{X,Y}(x, y) = \mathbb{P}(X > x, Y > y) = \sum_{j=1}^p \alpha_j \left[ e_j^{\mathsf{T}} M_1(x) e \right] \left[ e_j^{\mathsf{T}} M_2(y) e \right]$$

where $M_i(t) = \exp\left( \int_0^t T_i(u) du \right)$ is the cumulative sub-intensity matrix, and $e_j$ is the $j$th unit vector.

The joint density for observed lifetimes is:

$$f_{X,Y}(x,y) = \sum_{j=1}^p \alpha_j \left[ e_j^{\mathsf{T}} M_1(x) t_1(x) \right] \left[ e_j^{\mathsf{T}} M_2(y) t_2(y) \right]$$

with $t_i(t) = -T_i(t)e$ encoding exit rates from each phase.

3. Causal Mediation and Independence Structure

Causal consistency in the mIPH model is enforced by mediating all dependence via the single latent starting state $J_0$. The model defines a sharp causal structure:

• Covariates $A$ (such as ages at issue, health indicators) determine the initial distribution $\alpha_j^{(m)} = \mathbb{P}(J_0 = j | A^{(m)})$ through multinomial logistic regression.
• Conditional on $J_0$, the two chains and their lifetimes $(X, Y)$ evolve independently, and hazards at future times depend only on elapsed time and $J_0$, not on the partner's observed outcome.

A formal statement:

$$J_0 \rightarrow X,\qquad J_0 \rightarrow Y, \qquad X|J_0\text{ and }Y|J_0\text{ are independent}$$

There is no edge $Y \to X$ nor $X \to Y$, and updating survivor information simply refines the posterior over $J_0$, adhering to a pure causal mediation paradigm. For conditional hazard estimation:

$$\mathbb{P}(X \in dx | Y=y, A) = \sum_{j=1}^p \underbrace{\frac{\alpha_j(A) e_j^{\mathsf{T}} M_2(y) t_2(y)}{\sum_{k=1}^p \alpha_k(A) e_k^{\mathsf{T}} M_2(y) t_2(y)}}_{\,\text{state law at time }y}\, e_j^{\mathsf{T}} M_1(x) t_1(x) dx$$

Thus, knowing $Y$ does not directly influence $X$ apart from its effect on the latent $J_0$ distribution.

4. Parameter Estimation: EM-Type and Covariate Modeling

Maximum likelihood estimation in the causality-consistent PH model proceeds via a custom EM algorithm (ERMI), which accommodates right-censored observations and covariate-driven heterogeneity:

• E-step: Compute posterior weights $$w_j^{(m)} = \mathbb{P}(J_0 = j | y_1^{(m)}, y_2^{(m)}, A^{(m)})$$ and sufficient statistics for transitions and sojourn times.
• M-step: Update phase transition matrix $T$ via estimated transition counts and total occupation times.
• R-step: Update multinomial regression coefficients $\{\gamma_j\}$ governing covariate influence on starting probabilities.
• I-step: Optimize the time-inhomogeneity functions $\{\beta_i\}$ controlling phase progressions.

Observed-data likelihood for individual $m$ is a mixture over latent states:

$$L_m = \sum_{j=1}^p \alpha_j^{(m)} \prod_{i=1}^2 \Bigl[ e_j^{\mathsf{T}} M_i(y_i^{(m)}) t_i(y_i^{(m)}) \Bigr]^{\delta_i^{(m)}} \Bigl[ e_j^{\mathsf{T}} M_i(y_i^{(m)}) e \Bigr]^{1-\delta_i^{(m)}}$$

where $y_i^{(m)}$ is observed or censored time, $\delta_i^{(m)}$ is the censoring indicator.

5. Time-Inhomogeneity: Flexibility and Parsimony

A central feature of the mIPH model is scalar time-inhomogeneity, enabling significant flexibility and parameter parsimony. Commonly, the sub-intensity matrices factor as $T_i(t) = \lambda_i(t) T$, with $\lambda_i(t) > 0$ scalar functions and $T$ a common constant matrix, reducing the dimensionality of parameter estimation. Notably, the matrix-Gompertz specification $\lambda_i(t) = \exp(\beta_i t)$ yields a survival matrix:

$$M_i(t) = \exp\left(T \frac{e^{\beta_i t}-1}{\beta_i} \right)$$

This structure allows for high-quality fits with few phases, as illustrated with ten phases for joint-lifetime data of Frees et al. (Hansjörg et al., 2022):

$$\hat\beta_1=43.101,\quad \hat\beta_2=47.474, \quad p=10$$

6. Practical Implications and Interpretability

In practical applications, the causality-consistent PH model provides a mechanism for inferring associations between paired lifetimes (e.g., spouses, insured lives) through a latent health state determined by explicit covariates. Conditioning on one partner's survival time updates knowledge of the shared latent state, thereby shifting the hazard of the other—a direct encoding of causal mediation. The model eschews copula construction, sidestepping issues of non-identifiability and arbitrary dependence structures. All model properties correspond precisely to the directed-graph mechanism:

• Joint survival arises entirely from the shared initial condition.
• Time-inhomogeneity admits model parsimony and robustness in fitting.
• Covariates enter exclusively via the multinomial latent-state initialization.

7. Summary Table of Key Components

Component	Mathematical Notation	Description
State-space	$E = \{1,\dots, p, p+1\}$	Phases + absorbing state
Initial latent state	$J_0 \sim \alpha$	Covariate-driven start
Sub-intensity matrix	$T_i(t)$	Time-dependent rates
Joint survival function	$S_{X,Y}(x,y)$	Closed-form formula
Inhomogeneity specification	$\lambda_i(t)$	Time scaling factor
Multinomial regression	$\alpha_j^{(m)}$ via $\gamma_j$	Covariate mapping
EM estimation steps	E, M, R, I	Posterior + parameter fit

This architecture provides a fully causally interpretable ageing mechanism for joint lifetimes, with directed associations and closed-form formulas throughout (Hansjörg et al., 2022). The mIPH construction enables exact mediation modeling, right-censoring handling, and interpretable dependence consistent with real-world causal processes.