Bayesian tree-based heterogeneous mediation analysis with a time-to-event outcome

Abstract

Mediation analysis aims at quantifying and explaining the underlying causal mechanism between an exposure and an outcome of interest. In the context of survival analysis, mediation models have been widely used to achieve causal interpretation for the direct and indirect effects on the survival of interest. Although heterogeneity in treatment effect is drawing increasing attention in biomedical studies, none of the existing methods have accommodated the presence of heterogeneous causal pathways pointing to a time-to-event outcome. In this study, we consider a heterogeneous mediation analysis for survival data based on a Bayesian tree-based Cox proportional hazards model with shared topologies. Under the potential outcomes framework, individual-specific conditional direct and indirect effects are derived on the scale of the logarithm of hazards, survival probability, and restricted mean survival time. A Bayesian approach with efficient sampling strategies is developed to estimate the conditional causal effects through the Monte Carlo implementation of the mediation formula. Simulation studies show the satisfactory performance of the proposed method. The proposed model is then applied to an HIV dataset extracted from the ACTG175 study to demonstrate its usage in detecting heterogeneous causal pathways.

Appendices

Appendix A: Derivation of Eqs. (3) and (4)

Under the identification assumptions stated in Sect. 2.3, the targeted function of potential survival time for subject i, $\omega \{T_i(a, M_i(a'))\}$, is identified as

$$\begin{aligned}{} & {} \omega \{T_i(a,M_i(a'))\} {=} \int \omega \{T_i(a, m) | M_i(a') {=} m, \textbf{x}_i = \varvec{x}\}dF_{M_i(a')|\textbf{x}_i=\varvec{x}} \nonumber \\{} & {} \quad {\left[ \text {By Assumption } 4\right] }\nonumber \\{} & {} \quad = \int \omega \{T_i(a, m) |\textbf{x}_i = \varvec{x}\}dF_{M_i(a')|\textbf{x}_i=\varvec{x}}\nonumber \\{} & {} \quad {\left[ \text {By Assumption } 1\right] }\nonumber \\{} & {} \quad = \int \omega \{T_i(a, m) |A_i = a, \textbf{x}_i = \varvec{x}\}dF_{M_i(a')|\textbf{x}_i=\varvec{x}} \nonumber \\{} & {} \quad {\left[ \text {By Assumption } 2\right] }\nonumber \\{} & {} \quad = \int \omega \{T_i(a, m) |A_i = a, M_i = m, \textbf{x}_i = \varvec{x}\}dF_{M_i(a')|\textbf{x}_i=\varvec{x}} \nonumber \\{} & {} \quad {\left[ \text {By Assumption } 3\right] }\nonumber \\{} & {} \quad = \int \omega \{T_i(a, m) |A_i = a, M_i = m, \textbf{x}_i = \varvec{x}\}dF_{M_i(a')|A_i = a', \textbf{x}_i=\varvec{x}} \nonumber \\{} & {} \quad {\left[ \text {By Consistency}\right] }\nonumber \\{} & {} \quad = \int \omega \{T_i |A_i = a, M_i = m, \textbf{x}_i = \varvec{x}\}dF_{M_i | A_i = a', \textbf{x}_i=\varvec{x}}. \nonumber \\ \end{aligned}$$

(A.1)

The last equation follows from the basic consistency assumption of the counterfactual framework, which states that an individual’s potential outcome under their observed exposure history is exactly their observed outcome (Imai et al. 2010; Ten Have and Joffe 2012), i.e., $M_i(A_i) = M_i$ and $T_i(A_i, M_i(A_i)) = T_i$. We used the Monte Carlo method to integrate over $F_{M_i|A_i, \textbf{x}_i}$, and also adopted the comonotonic sampling strategy (Linero and Zhang 2022) to enhance the efficiency of the numerical implementation of the mediation formula. Major steps are outlined in Algorithm A1.

Algorithm 2

The Monte Carlo implementation of mediation formula

By using inverse transform sampling with shared $u_i$, $M_i^{(itr)}(1)$ and $M_i^{(itr)}(0)$ are sampled as comonotonic pairs. As explained by Linero and Zhang (2022), this approach is effective in reducing Monte Carlo errors by making the $M_i^{(itr)}(a)$ pairs as correlated as possible, thereby allowing for accurate estimation of the conditional PSEs with only one sample per iteration.

Fig. 5

Trace plots of a $\sigma$, b $\lambda _1$, and c–d two randomly selected $\tau (\textbf{x}_i)$ generated under setup (i) with $\lambda _0(t) = 1$

Table 10 Average bias, relative bias, RMSE, coverage rate and length of the 95% credible interval for the average PSEs on the scale of log hazard, survival probability, and RMST at the maximal observed time under setup (i) with $n = 1000$ and a true baseline hazard of $\lambda _0(t)=1/2t^{1/2}$

Table 11 Average MSE of mediation effects across individuals and average (in-sample) PRMSE of the mediator and the time-to-event outcome on three scales under setup(i) with $\lambda _0(t)=1/2t^{1/2}$

The group-average or sample-average PSEs are identified through an additional integration over the covariate distribution $F_{\textbf{x}}$. We followed the common practice in the mediation literature to approximate the integration using empirical distribution of the covariates based on the collected sample (Albert and Nelson 2011; Bonetti and Gelber 2004), under which circumstance Eq. 4 simplifies to

$$\begin{aligned}{} & {} {} {\bar{\omega }}\{T(a,M(a'))\} \approx \frac{1}{n}\sum _{i=1}^n \omega \{T_i(a,M_i(a'))\} \nonumber \\{}{} & {} {} = \frac{1}{n}\sum _{i=1}^n \int \omega \{T_i |A_i = a, M_i = m, {\textbf {x}}_i = \varvec{x}\} dF_{M_i | A_i = a', {\textbf {x}}_i=\varvec{x}}.\nonumber \\ \end{aligned}$$

(A.2)

Appendix B: Full conditional distributions of the tree ensembles

Denote ${\textbf{Y}} = (Y_1, Y_2, \dots , Y_n)^T$, ${\textbf{M}} = (M_1, M_2, \dots , M_n)^T$, ${\textbf{A}} = (A_1, A_2, \dots , A_n)^T$, ${\textbf{X}} = (\textbf{x}_1^T, \textbf{x}_2^T, \dots , \textbf{x}_n^T)^T$, and $\varvec{\delta } = (\delta _1, \dots , \delta _n)^T$. Let ${\mathcal {D}} = \{{\textbf{Y}},{\textbf{X}},{\textbf{A}},\varvec{\delta }, {\textbf{M}}\}$ denote the observed data. The complete data likelihood of the proposed model is expressed as

$$\begin{aligned} p({\mathcal {D}}|\cdot )= & {} \prod _{i=1}^{n}\exp \Bigg [-\Lambda _0(Y_i)\exp \{\upsilon (\textbf{x}_i, M_i)+ \tau (\textbf{x}_i)A_i\}\big ]\nonumber \\{} & {} \quad \big [\lambda _0(Y_i)\exp \{\upsilon (\textbf{x}_i,M_i)+ \tau (\textbf{x}_i)A_i\}\Bigg ]^{\delta _i}\nonumber \\{} & {} \quad \times \frac{1}{\sqrt{2\pi }\sigma }\exp \{-\frac{(M_i - \upsilon _M(\textbf{x}_i) - \tau _M(\textbf{x}_i)A_i )^2}{2\sigma ^2}\},\nonumber \\ \end{aligned}$$

(A.3)

where $\Lambda _0(t) = \int _0^t \lambda _0(u)du$ is the cumulative baseline hazard function. After assuming prior independence among the unknown parameters $\{\mathcal {T}_j, \mathcal {M}_j\}_{j=1}^J$, $\{\widetilde{\mathcal {T}}_h, \widetilde{\mathcal {M}}_h^{(M)}, \widetilde{\mathcal {M}}_h\}_{h=1}^H$, $\{\breve{\mathcal {T}}_k, \breve{\mathcal {M}}_j\}_{k=1}^K$, $\varvec{\lambda }$, and $\sigma$, and prior independence among the individual binary trees inside each tree ensembles as well, the full conditional distribution of the parameters are derived as

follows.

(I) Full conditional distribution of $\{\mathcal {T}_j, \mathcal {M}_j\}_{j=1}^J$

Let $R_{ij} = M_i - \sum _{l \ne j}g({\textbf {x}}_i;\mathcal {T}_l, \mathcal {M}_l) - \tau _M({\textbf {x}}_i)A_i$ and $\varvec{R}_j = (R_{1j}, \dots , R_{nj})^T$ denote the partial residuals in the mediator regression equation w.r.t the jth binary tree $\{\mathcal {T}_j,\mathcal {M}_j\}$. By prior independence and the i.i.d. assumption on the sample data,

$$\begin{aligned}{} & {} {} p\left( \{\mathcal {T}_j,\mathcal {M}_j\}_{j = 1}^{J}|{\mathcal {D}}, \sigma , \{\widetilde{\mathcal {T}}_h,\widetilde{\mathcal {M}}_h^{(M)}\}_{h=1}^H \right) \\{}{} & {} {} =\prod _{j=1}^{J}p(\mathcal {M}_j|\varvec{R}_j, \mathcal {T}_j, \sigma )p(\mathcal {T}_j|\varvec{R}_j,\sigma ), \end{aligned}$$

where

$$\begin{aligned}{} & {} {} p( \mathcal {M}_j|{\mathcal {D}},\mathcal {T}_j,\sigma ) \propto p(\varvec{R}_j|\mathcal {M}_j,\mathcal {T}_j, \sigma )p(\mathcal {M}_j|\mathcal {T}_j)\nonumber \\{}{} & {} {} \propto \prod _{l=1}^{b_j}\prod _{i:\mathbb {1}_{il}=1}\exp \Big \{-\frac{1}{2\sigma ^2}(R_{ij} -\mu _{jl})^2\Big \}\exp \Big \{-\frac{1}{2\sigma _\mu ^2}\mu _{jl}^2\Big \}\nonumber \\{}{} & {} {} \propto \prod _{l=1}^{b_j}\exp \Big \{-\frac{1}{2\sigma _\mu ^2}\mu _{jl}^2-\frac{1}{2\sigma ^2}\sum _{i:\mathbb {1}_{il} = 1}(R_{ij} - \mu _{jl})^2\Big \}\nonumber \\{}{} & {} {} \overset{d}{=}\ \prod _{l=1}^{b_j}N(a_l,\sigma ^2_l), \end{aligned}$$

(A.4)

with $\sigma _l^2 = \frac{\sigma ^2\sigma _\mu ^2}{\sigma ^2 + \sigma _\mu ^2\sum _i\mathbb {1}_{il}}$, $a_l = \frac{\sigma _\mu ^2\sum _{i:\mathbb {1}_{il} = 1}R_{ij}}{\sigma ^2 + \sigma _\mu ^2\sum _i\mathbb {1}_{il}}$, $\sigma ^2_\mu = \frac{c^2_1}{4Jk_1^2}$, $\mathbb {1}_{il}$ denoting that subject i is allocated to the lth terminal node of $\mathcal {T}_j$ with mean parameter $\mu _{jl}$, and ‘$\overset{d}{=}$’ representing equivalence in distribution. The conjugacy makes it feasible to draw $\mathcal {T}_j$ from its marginal posterior distribution where $\mathcal {M}_j$ is integrated out to avoid the need for reversible jumps. That is,

$$\begin{aligned}{} & {} {} p(\mathcal {T}_j| \varvec{R}_j,\sigma ) \propto p(\mathcal {T}_j)\left( \int p(\varvec{R}_j|\mathcal {M}_j, \mathcal {T}_j, \sigma )p(\mathcal {M}_j|\mathcal {T}_j)d\mathcal {M}_j\right) \nonumber \\{}{} & {} {} \propto p(\mathcal {T}_j)\prod _{l=1}^{b_j} \int \prod _{i:\mathbb {1}_{il}=1}\frac{1}{\sqrt{2\pi \sigma ^2}}\exp \left\{ -\frac{1}{2\sigma ^2}(R_{ij} - \mu _{jl})^2\right\} \nonumber \\{}{} & {} {} \quad \times \frac{1}{\sqrt{2\pi \sigma ^2_\mu }}\exp \left\{ -\frac{\mu _{jl}^2}{2\sigma ^2_\mu }\right\} d\mu _{jl}\nonumber \\{}{} & {} {} \propto p(\mathcal {T}_j)\prod _{l=1}^{b_j} \left( \frac{1}{\sqrt{2\pi \sigma ^2}}\right) ^{n}\exp \left\{ -\frac{\sum _{i=1}^n(R_{ij})^2}{2\sigma ^2}\right\} \nonumber \\{}{} & {} {} \quad \times \sqrt{\frac{\sigma _l^2}{\sigma ^2_\mu }}\exp \left\{ \frac{\left( \sum _{i:\mathbb {1}_{il} = 1}R_{ij}\right) ^2\sigma _l^2}{2\sigma ^4}\right\} . \end{aligned}$$

(A.5)

To acquire posterior samples of $\mathcal {T}_j$, we employed the Metropolis-Hastings algorithm (Metropolis et al. 1953; Hastings 1970), with a classical proposal distribution suggested by (Chipman et al. 2010) that includes four possible moves with probabilities as follows: splitting at a leaf node (0.25), pruning a pair of leaf nodes (0.25), changing the splitting rule of an internal node (0.4), and swapping the splitting rule between an internal node and its child (0.1). Given the current iteration $\mathcal {T}_j^{(itr)}$, a candidate tree $\mathcal {T}_j^*$ is proposed according to the proposal distribution $q(\cdot |\cdot )$, and is accepted as $\mathcal {T}_j^{(itr+1)}$ with a probability of $\min \left\{ 1, \frac{q\left( \mathcal {T}_j^{(itr)}\big |\mathcal {T}_j^*\right) p\left( \mathcal {T}_j^*\big |\varvec{R}^{(itr)}_j, \sigma ^{(itr)}\right) }{q\left( \mathcal {T}_j^*\big |\mathcal {T}_j^{(itr)}\right) p\left( \mathcal {T}_j^{(itr)}\big |\varvec{R}^{(itr)}_j, \sigma ^{(itr)}\right) }\right\}$.

(II) Full conditional distribution of $\{\widetilde{\mathcal {T}}_h, \widetilde{\mathcal {M}}_h, \widetilde{\mathcal {M}}_h^{(M)}\}_{h=1}^H$

Let $\widetilde{R}_{ih}^{(M)} = M_i - \upsilon _M({\textbf {x}}_i) - \sum _{l \ne h}g({\textbf {x}}_i;\widetilde{\mathcal {T}}_l, \widetilde{\mathcal {M}}_l^{(M)})A_i$ and $\varvec{\widetilde{R}}_h^{(M)} = (\widetilde{R}_{1\,h}^{(M)}, \dots , \widetilde{R}_{nh}^{(M)})^T$ denote the partial residuals in the mediator regression equation w.r.t the hth binary tree $\{\widetilde{\mathcal {T}}_h,\widetilde{\mathcal {M}}_h^{(M)}\}$. Let $\widetilde{R}_{ih} = \upsilon ({\textbf {x}}_i, M_i) + \sum _{l \ne h}g({\textbf {x}}_i;\widetilde{\mathcal {T}}_l, \widetilde{\mathcal {M}}_l)A_i$ and $\varvec{\widetilde{R}}_h = (\widetilde{R}_{1h}, \dots , \widetilde{R}_{nh})^T$ denote the partial residuals in the PH model w.r.t $\{\widetilde{\mathcal {T}}_h,\widetilde{\mathcal {M}}_h\}$. The full conditionals of $\{\widetilde{\mathcal {T}}_h, \widetilde{\mathcal {M}}_h, \widetilde{\mathcal {M}}_h^{(M)}\}_{h=1}^H$ are derived similarly, with the only difference being that the leaf node parameters are bivariate due to the shared tree topologies. Assuming independence among the leaf node parameters, we have

$$\begin{aligned}{} & {} {} p\left( \{\widetilde{\mathcal {T}}_h,\widetilde{\mathcal {M}}_h^{(M)}, \widetilde{\mathcal {M}}_h\}_{h = 1}^{H}|{\mathcal {D}}, \sigma , \{(\mathcal {T}_j,\mathcal {M}_j)\}_{j=1}^J, \{\breve{\mathcal {T}}_k,\breve{\mathcal {M}}_k\}_{k=1}^K, \varvec{\lambda } \right) \\{}{} & {} {} =\,\prod _{h=1}^{H}p\left( \widetilde{\mathcal {M}}_h^{(M)},\widetilde{\mathcal {M}}_h|\varvec{\widetilde{R}}_h^{(M)}, \varvec{\widetilde{R}}_h,\widetilde{\mathcal {T}}_h, \varvec{\lambda },\sigma \right) \\{}{} & {} {} p\left( \widetilde{\mathcal {T}}_h|\varvec{\widetilde{R}}_h^{(M)}, \varvec{\widetilde{R}}_h,\varvec{\lambda },\sigma \right) =\,\prod _{h=1}^{H}p\left( \widetilde{\mathcal {M}}_h| \varvec{\widetilde{R}}_h,\widetilde{\mathcal {T}}_h, \varvec{\lambda }\right) p\left( \widetilde{\mathcal {M}}_h^{(M)}|\varvec{\widetilde{R}}_h^{(M)}, \widetilde{\mathcal {T}}_h, \sigma \right) \\{}{} & {} {} \quad \times p\left( \widetilde{\mathcal {T}}_h|\varvec{\widetilde{R}}_h^{(M)}, \varvec{\widetilde{R}}_h,\varvec{\lambda },\sigma \right) \end{aligned}$$

where

$$\begin{aligned}{} & {} {} p\left( \widetilde{\mathcal {M}}_h| \varvec{\widetilde{R}}_h,\widetilde{\mathcal {T}}_h, \varvec{\lambda }\right) \propto p\left( \varvec{\widetilde{R}}_h|\widetilde{\mathcal {M}}_h,\widetilde{\mathcal {T}}_h,\varvec{\lambda }\right) p\left( \widetilde{\mathcal {M}}_h|\widetilde{\mathcal {T}}_h\right) \nonumber \\{}{} & {} {} \propto \prod _{l=1}^{b_h}\prod _{i:\mathbb {1}_{il}=1}\exp \left[ \delta _i\left\{ \widetilde{R}_{ih} + \widetilde{\mu }_{hl}A_i\right\} - \exp \left\{ \widetilde{R}_{ih} + \widetilde{\mu }_{hl}A_i\right\} \right. \nonumber \\{}{} & {} {} \quad \left. \times \sum _{s:i\in {\mathcal {R}}(s)}\lambda _s\delta _s\right] \exp \left[ \widetilde{\mu }_{hl}{\widetilde{\zeta }} - \exp (\widetilde{\mu }_{hl}){\widetilde{\eta }}\right] \nonumber \\{}{} & {} {} \propto \prod _{l=1}^{b_h} \exp \left[ \widetilde{\mu }_{hl}\left( \sum _{i:\mathbb {1}_{il}= 1}\delta _iA_i + {\widetilde{\zeta }}\right) \right. \nonumber \\{}{} & {} {} \quad \left. - \exp (\widetilde{\mu }_{hl})\left\{ {\widetilde{\eta }} +\sum _{i:\mathbb {1}_{il}= 1,A_i=1}\exp (\widetilde{R}_{ih})\sum _{s:i\in {\mathcal {R}}(s)}\lambda _s\delta _s\right\} \right] \nonumber \\{}{} & {} {} \overset{d}{=}\prod _{l=1}^{b_h}logGam({\widetilde{\zeta }}^*, {\widetilde{\eta }}^*), \end{aligned}$$

(A.6)

with ${\widetilde{\zeta }}^* = \sum _{i:\mathbb {1}_{il}= 1}\delta _iA_i + {\widetilde{\zeta }}$ and ${\widetilde{\eta }}^* = {\widetilde{\eta }} +\sum _{i:\mathbb {1}_{il}= 1,A_i=1}$$\exp (\widetilde{R}_{ih})\sum _{s:i\in {\mathcal {R}}(s)}\lambda _s\delta _s$; and

$$\begin{aligned}{} & {} {} p\left( \widetilde{\mathcal {M}}_h^{(M)}|\varvec{\widetilde{R}}_h^{(M)}, \widetilde{\mathcal {T}}_h, \sigma \right) \propto p\left( \varvec{\widetilde{R}}_h^{(M)}|\widetilde{\mathcal {M}}_h^{(M)},\widetilde{\mathcal {T}}_h, \sigma \right) p\left( \widetilde{\mathcal {M}}_h^{(M)}|\widetilde{\mathcal {T}}_h\right) \nonumber \\{}{} & {} {} \propto \prod _{l=1}^{b_h}\prod _{i:\mathbb {1}_{il}=1}\exp \left\{ -\frac{1}{2\sigma ^2}\left( \widetilde{R}_{ih}^{(M)} -\widetilde{\mu }_{hl}^{(M)}\right) ^2\right\} \exp \left\{ -\frac{1}{2\sigma _{\widetilde{\mu }}^2}\left( \widetilde{\mu }_{hl}^{(M)}\right) ^2\right\} \nonumber \\{}{} & {} {} \overset{d}{=}\ \prod _{l=1}^{b_h}N({\widetilde{a}}_l,{\widetilde{\sigma }}^2_l), \end{aligned}$$

(A.7)

with ${\widetilde{\sigma }}_l^2 = \frac{\sigma ^2\sigma _{\widetilde{\mu }}^2}{\sigma ^2 + \sigma _{\widetilde{\mu }}^2\sum _i\mathbb {1}_{il}}$, ${\widetilde{a}}_l = \frac{\sigma _{\widetilde{\mu }}^2\sum _{i:\mathbb {1}_{il} = 1}\widetilde{R}_{ih}^{(M)}}{\sigma ^2 + \sigma _{\widetilde{\mu }}^2\sum _i\mathbb {1}_{il}}$, and $\sigma ^2_{\widetilde{\mu }} = \frac{c^2_2}{4Hk_2^2}$.

Using the conjugate normal and log-gamma priors once again leads to closed-form marginal posteriors for $\widetilde{\mathcal {T}}_h$, which is integrated over $\{\widetilde{\mathcal {M}}_h,\widetilde{\mathcal {M}}_h^{(M)}\}$ and expressed as

$$\begin{aligned}{} & {} p\left( \widetilde{\mathcal {T}}_h|\varvec{\widetilde{R}}_h, \varvec{\widetilde{R}}_h^{(M)},\varvec{\lambda },\sigma \right) \propto \nonumber \\{} & {} \quad p(\widetilde{\mathcal {T}}_h)\int p\left( \varvec{\widetilde{R}}_h|\widetilde{\mathcal {M}}_h, \widetilde{\mathcal {T}}_h, \varvec{\lambda }\right) p\left( \widetilde{\mathcal {M}}_h|\widetilde{\mathcal {T}}_h\right) d\widetilde{\mathcal {M}}_h\nonumber \\{} & {} \quad \int p\left( \varvec{\widetilde{R}}_h^{(M)}|\widetilde{\mathcal {M}}_h^{(M)}, \widetilde{\mathcal {T}}_h, \sigma \right) p\left( \widetilde{\mathcal {M}}_h^{(M)}|\widetilde{\mathcal {T}}_h\right) d\widetilde{\mathcal {M}}_h^{(M)}\nonumber \\{} & {} \quad \propto p(\widetilde{\mathcal {T}}_h)\prod _{l=1}^{b_h} \int \exp \Bigg [\sum _{i:\mathbb {1}_{il}= 1}\delta _i\widetilde{R}_{ih} + \widetilde{\mu }_{hl}\sum _{i:\mathbb {1}_{il}= 1}\delta _iA_i\nonumber \\{} & {} \quad - \exp (\widetilde{\mu }_{hl})\sum _{i:\mathbb {1}_{il}=A_i= 1}\exp (\widetilde{R}_{ih})\sum _{s:i\in {\mathcal {R}}(s)}\lambda _s\delta _s\nonumber \\{} & {} \quad - \sum _{i:\mathbb {1}_{il}= 1,A_i=0}\exp (\widetilde{R}_{ih})\sum _{s:i\in {\mathcal {R}}(s)}\lambda _s\delta _s\Bigg ]\nonumber \\{} & {} \quad \times \frac{{\widetilde{\eta }}^{{\widetilde{\zeta }}}}{\Gamma ({\widetilde{\zeta }})}\exp \left( \widetilde{\mu }_{hl} {\widetilde{\zeta }} - \exp (\widetilde{\mu }_{hl}){\widetilde{\eta }}\right) d\widetilde{\mu }_{hl}\nonumber \\{} & {} \quad \times \int \prod _{i:\mathbb {1}_{il}=1}\frac{1}{\sqrt{2\pi \sigma ^2}}\exp \left\{ -\frac{1}{2\sigma ^2}\left( \widetilde{R}_{ih}^{(M)} - \widetilde{\mu }_{hl}^{(M)}\right) ^2\right\} \nonumber \\{} & {} \quad \frac{1}{\sqrt{2\pi \sigma ^2_{\widetilde{\mu }}}}\exp \left\{ -\frac{\left( \widetilde{\mu }_{hl}^{(M)}\right) ^2}{2\sigma ^2_{\widetilde{\mu }}}\right\} d\widetilde{\mu }_{hl}^{(M)}\nonumber \\{} & {} \quad \propto p(\widetilde{\mathcal {T}}_h)\prod _{l=1}^{b_h}\exp \left( \sum _{i:\mathbb {1}_{il}= 1}\delta _i\widetilde{R}_{ih}\right. \nonumber \\{} & {} \left. \quad - \sum _{i:\mathbb {1}_{il}= 1,A_i=0}\exp (\widetilde{R}_{ih})\sum _{s:i\in {\mathcal {R}}(s)}\lambda _s\delta _s\right) \frac{{\widetilde{\eta }}^{{\widetilde{\zeta }}}\Gamma ({\widetilde{\zeta }}^*)}{\Gamma ({\widetilde{\zeta }}) ({\widetilde{\eta }}^*)^{{\widetilde{\zeta }}^*}}\nonumber \\{} & {} \quad \times \left( \frac{1}{\sqrt{2\pi \sigma ^2}}\right) ^{n}\exp \left\{ -\frac{\sum _{i=1}^n(\widetilde{R}_{ih}^{(M)})^2}{2\sigma ^2}\right\} \nonumber \\{} & {} \quad \sqrt{\frac{{\widetilde{\sigma }}_l^2}{\sigma ^2_{\widetilde{\mu }}}}\exp \left\{ \frac{\left( \sum _{i:\mathbb {1}_{il} = 1}\widetilde{R}_{ih}^{(M)}\right) ^2{\widetilde{\sigma }}_l^2}{2\sigma ^4}\right\} . \end{aligned}$$

(A.8)

Posterior draws of $\widetilde{\mathcal {T}}_h$ can thus be obtained using MH algorithm with the same proposal distribution described above.

(III) Full conditional distribution of $\{\breve{\mathcal {T}}_k, \breve{\mathcal {M}}_k\}_{k=1}^K$

Let $\breve{R}_{ik} = \sum _{l \ne k}g({\textbf {x}}_i, M_i;\breve{\mathcal {T}}_l, \breve{\mathcal {M}}_l) + \tau ({\textbf {x}}_i)A_i$ and $\varvec{\breve{R}}_k = (\breve{R}_{1k}, \dots , \breve{R}_{nk})^T$ denote the partial residuals in the PH model w.r.t $\{\breve{\mathcal {T}}_k,\breve{\mathcal {M}}_k\}$. The full conditionals of $\{\breve{\mathcal {T}}_k, \breve{\mathcal {M}}_k\}_{k=1}^K$ are derived as

$$\begin{aligned}{} & {} {} p\left( \{\breve{\mathcal {T}}_k,\breve{\mathcal {M}}_k\}_{k = 1}^{K}|{\mathcal {D}}, \varvec{\lambda }, \{\widetilde{\mathcal {T}}_h,\widetilde{\mathcal {M}}_h\}_{h=1}^H \right) \\{}{} & {} {} = \prod _{k=1}^{K}p\left( \breve{\mathcal {M}}_k|\varvec{\breve{R}}_k, \breve{\mathcal {T}}_k, \varvec{\lambda }\right) p\left( \breve{\mathcal {T}}_k|\varvec{\breve{R}}_k,\varvec{\lambda }\right) , \end{aligned}$$

where

$$\begin{aligned}{} & {} {} p\left( \breve{\mathcal {M}}_k|\varvec{\breve{R}}_k, \breve{\mathcal {T}}_k, \varvec{\lambda }\right) \propto p\left( \varvec{\breve{R}}_k|\breve{\mathcal {M}}_k, \breve{\mathcal {T}}_k, \varvec{\lambda }\right) p\left( \breve{\mathcal {M}}_k|\breve{\mathcal {T}}_k\right) \nonumber \\{}{} & {} {} \propto \prod _{l=1}^{b_k}\prod _{i:\mathbb {1}_{il}=1}\exp \left[ \delta _i\left\{ \breve{R}_{ik} + \breve{\mu }_{kl}\right\} \right. \nonumber \\{}{} & {} {} \quad \left. - \exp \left\{ \breve{R}_{ik} + \breve{\mu }_{kl}\right\} \sum _{s:i\in {\mathcal {R}}(s)}\lambda _s\delta _s\right] \exp \left[ \breve{\mu }_{kl}\breve{\zeta } - \exp (\breve{\mu }_{kl})\breve{\eta }\right] \nonumber \\{}{} & {} {} \propto \prod _{l=1}^{b_k} \exp \left[ \breve{\mu }_{kl}\left( \sum _{i:\mathbb {1}_{il}= 1}\delta _i + \breve{\zeta }\right) \right. \nonumber \\{}{} & {} {} \quad \left. - \exp (\breve{\mu }_{kl})\left\{ \breve{\eta } +\sum _{i:\mathbb {1}_{il}= 1}\exp (\breve{R}_{ik})\sum _{s:i\in {\mathcal {R}}(s)}\lambda _s\delta _s\right\} \right] \nonumber \\{}{} & {} {} \overset{d}{=}\prod _{l=1}^{b_k}\text{ log-gamma }(\breve{\zeta }^*, \breve{\eta }^*), \end{aligned}$$

(A.9)

with $\breve{\zeta }^* = \sum _{i:\mathbb {1}_{il}= 1}\delta _i + \breve{\zeta }$ and $\breve{\eta }^* = \breve{\eta } +\sum _{i:\mathbb {1}_{il}= 1}\exp (\breve{R}_{ik}) \sum _{s:i\in {\mathcal {R}}(s)}\lambda _s\delta _s$. The marginal posterior distribution of $\breve{\mathcal {T}}_k$ is then derived as

$$\begin{aligned}{} & {} {} p(\breve{\mathcal {T}}_k|\varvec{\breve{R}}_k, \varvec{\lambda }) \propto p(\breve{\mathcal {T}}_k)\int p\left( \varvec{\breve{R}}_k|\breve{\mathcal {M}}_k, \breve{\mathcal {T}}_k, \varvec{\lambda }\right) p\left( \breve{\mathcal {M}}_k|\breve{\mathcal {T}}_k\right) d\breve{\mathcal {M}}_k,\nonumber \\{}{} & {} {} \propto p(\breve{\mathcal {T}}_k)\prod _{l=1}^{b_k}\int \exp \left[ \sum _{i:\mathbb {1}_{il}= 1}\delta _i\breve{R}_{ik} + \breve{\mu }_{kl}\sum _{i:\mathbb {1}_{il}= 1}\delta _i \right. \nonumber \\{}{} & {} {} \quad \left. - \exp (\breve{\mu }_{kl})\sum _{i:\mathbb {1}_{il}= 1}\exp (\breve{R}_{ik})\sum _{s:i\in {\mathcal {R}}(s)}\lambda _s\delta _s\right] \nonumber \\{}{} & {} {} \quad \times \frac{\breve{\eta }^{\breve{\zeta }}}{\Gamma (\breve{\zeta })}\exp \left\{ \breve{\mu }_{kl}\breve{\zeta } - \exp (\breve{\mu }_{kl})\breve{\eta }\right\} d\breve{\mu }_{kl}\nonumber \\{}{} & {} {} \propto p(\breve{\mathcal {T}}_k)\prod _{l=1}^{b_k}\exp \left( \sum _{i:\mathbb {1}_{il}= 1}\delta _i\breve{R}_{ik}\right) \nonumber \\{}{} & {} {} \quad \int \frac{\breve{\eta }^{\breve{\zeta }}}{\Gamma (\breve{\zeta })}\exp \Bigg [ \breve{\mu }_{kl}\bigg (\sum _{i:\mathbb {1}_{il}= 1}\delta _i + \breve{\zeta }\bigg ) \nonumber \\{}{} & {} {} \quad - \exp (\breve{\mu }_{kl})\bigg \{\breve{\eta }+\sum _{i:\mathbb {1}_{il}= 1}\exp (\breve{R}_{ik})\sum _{s: i\in {\mathcal {R}}(s)}\lambda _s\delta _s\bigg \}\Bigg ]d\breve{\mu }_{kl}\nonumber \\{}{} & {} {} \propto p(\breve{\mathcal {T}}_k)\prod _{l=1}^{b_k}\exp \left( \sum _{i:\mathbb {1}_{il}= 1}\delta _i\breve{R}_{ik}\right) \frac{\breve{\eta }^{\breve{\zeta }}\Gamma (\breve{\zeta }^*)}{\Gamma (\breve{\zeta })(\breve{\eta }^*)^{\breve{\zeta }^*}}. \end{aligned}$$

(A.10)