Using joint models for longitudinal and time-to-event data to investigate the causal effect of salvage therapy after prostatectomy

Abstract

Prostate cancer patients who undergo prostatectomy are closely monitored for recurrence and metastasis using routine prostate-specific antigen measurements. When prostate-specific antigen levels rise, salvage therapies are recommended in order to decrease the risk of metastasis. However, due to the side effects of these therapies and to avoid over-treatment, it is important to understand which patients and when to initiate these salvage therapies. In this work, we use the University of Michigan Prostatectomy Registry Data to tackle this question. Due to the observational nature of this data, we face the challenge that prostate-specific antigen is simultaneously a time-varying confounder and an intermediate variable for salvage therapy. We define different causal salvage therapy effects defined conditionally on different specifications of the longitudinal prostate-specific antigen history. We then illustrate how these effects can be estimated using the framework of joint models for longitudinal and time-to-event data. All proposed methodology is implemented in the freely-available R package JMbayes2.

1. Introduction

Many prostate cancer (PCa) patients undergo surgical removal of the prostate gland (radical prostatectomy) as their initial treatment, after being diagnosed with prostate cancer. Even though this surgical procedure is generally successful, the risk of recurrence and metastasis remains. For this reason, urologists closely monitor the prostate-specific antigen (PSA) levels of these patients. PSA can be easily measured from a blood sample, and a persistent rise in the PSA values suggests the cancer may be regrowing, although it is generally not yet detectable on imaging. After the initial surgery, PSA levels drop to near zero; however, PSA may rise again for some patients, leading the treating physicians to recommend salvage therapy to reduce their risk of metastasis. Salvage therapy typically consists of radiotherapy with or without androgen deprivation therapy (ADT), although in some cases ADT alone may be utilized. After salvage therapy, PSA levels nearly always drop, sometimes substantially, but typically rise again if metastasis is going to occur. Because of the serious side effects of the aforementioned salvage therapies and to avoid over-treatment, it is critical to understand which patients are most likely to benefit from these treatments and when they should be initiated. To our knowledge, there are no randomized trials that are designed explicitly to address this question. However, there are clinical trials that compare different types of salvage therapy. ¹ In this article, we will not be considering that different salvage therapies have different effects.

We will use the University of Michigan Prostatectomy Registry (UMP) Data to tackle this question. This database includes 3634 PCa patients who underwent a radical prostatectomy during the period 1996–2013. Of those patients, 271 (15.6%) received salvage therapy, 102 (2.8%) developed metastasis, and 209 (5.8%) died without metastases. Of these 209 patients, 190 died before salvage therapy, and 19 died after salvage therapy without developing metastasis. For the patients who received salvage, the median PSA just prior to the initiation of the therapy was 0.7 ng/mL (min: 0.0009, Q1: 0.4, Q3: 2.06, max: 266 ng/mL), and 55 had a metastasis. For the patients who did not receive salvage therapy or experience metastases, the median PSA at their last measurement was 0.001 ng/mL (min: 0.0009, Q1: 0.001, Q3: 0.09, max: 11.7 ng/mL). Since death due to prostate cancer is extremely rare without prior metastases, these deaths are considered as due to other causes. Supplemental Figures 1 and 2 show the cumulative incidence functions for metastasis and death and initiation of salvage therapy and death. The median follow-up time was 6.4 years (min: 0.1, Q1: 3.3, Q3: 10.9, max: 17 years). This dataset has been previously described by Beesley et al. ² The detection of metastasis usually requires imaging. Although imaging is noninvasive, it is not routinely performed without an indication, and the decision to undergo imaging would usually be based on whether the patients exhibited symptoms or if the PSA levels were consistently rising and had attained higher values. There is variation in the schedule of when PSA measurements are taken after the initial surgery, but an average schedule might be every three months for the first year, every 6 months for the next 2 years, and annually thereafter. PSA measurements are also collected after salvage therapy is initiated. We aim to utilize the longitudinal PSA measurements and baseline information (Gleason score, T-stage of the tumors, age, race, and comorbidities) to quantify the effect of salvage therapy on reducing or delaying metastases and aid the decision-making process of urologists. Given the observational nature of the UMP data, we face several challenges in achieving our goal. First, the longitudinal PSA process affects both metastasis and the assignment of salvage therapy, and, therefore, it is a time-varying confounder. At the same time, the salvage therapy affects future PSA measurements, and hence the longitudinal PSA process will also be an intermediate variable. In addition, patients who received salvage therapy are more closely monitored, and they will tend to have more PSA measurements and more checks done to see if metastasis has occurred than patients who did not. Finally, death from other causes is a competing risk for metastasis and requires accounting for it appropriately.

Nowadays, there is a considerable body of literature on methods developed to estimate causal treatment effects in settings with time-varying treatments when there exists confounding by time-dependent covariates affected by earlier treatments. While matching methods^3,4 and joint modeling approaches^5,6 have been developed, most of this literature has focused on marginal structural models, dynamic treatment regimes, and related methods. Excellent overviews of these methods are given by Hernán and Robins ⁷ and Tsiatis et al. ⁸ The advantage of these methods is that they are primarily non-parametric, making very few assumptions for the time-varying confounders. However, most of them require specifying a single model for all patients that relates the decision to give treatment on past confounder values. This is especially challenging in our setting because doctors decide when to initiate therapy on different grounds. For example, one physician may use just the last observed value of PSA. In contrast, another may treat patients who showed a sudden increase/decrease in the biomarker in the previous two/three measurements. In situations like this, specifying a single model for the decision to start therapy conditional on past PSA values is challenging. An alternative approach to estimating a causal effect is through direct modeling of the outcome process.^5,9 We will adopt this approach and avoid having to model the treatment initiation process, and employ the framework of joint models for longitudinal and time-to-event data. ¹⁰ In particular, we postulate a linear mixed-effects model for the longitudinal PSA levels that explicitly accounts for the change in the subject-specific trajectories after initiating salvage therapy. For the instantaneous risk of metastasis, we specify a relative risk model that includes the time-varying salvage therapy and PSA effects, and for the hazard of death, we only include the time-varying salvage therapy effect. Based on the postulated joint model, and for each patient at risk at time $t$ , we predict the cumulative risk of metastasis under the two treatment scenarios, initiating salvage therapy at $t$ or continuing without it. Then, we obtain the causal salvage therapy effect by suitably averaging these predicted cumulative risks over appropriate groups of subjects. The benefit of our approach is that it does not require defining a model for initiating treatment. In this context, we have performed a simulation study to investigate the performance of joint models in settings with time-varying confounding.

The rest of the article is organized as follows. Section 2 presents the definition of causal salvage therapy effects we wish to estimate and the assumptions we make to identify these effects from observational data. Section 3 presents the modeling framework we use to analyze the data, and Section 4 shows the procedure to estimate the salvage therapy effects and derive their variance. In Section 5, we apply our proposed methods in the UMP dataset. Finally, Section 6 includes the results of a simulation study, and Section 7 concludes the article.

2. Salvage therapy effects

2.1. Definitions

We let $T^m$ denote the time to metastasis and $T^d$ the time to death. PSA measurements during follow-up are denoted as $Y(t_j)$ taken at times points $t_j$ ( $j=1,\dots ,J$ ), and $\mathcal{Y}(t)=\{Y(t_j)=y(t_j);0\le t_j\le t,j=1,\dots ,J\}$ denotes the available PSA measurements history up to time $t$ . We denote the time salvage therapy was initiated by $S$ , and $N(t)=I(t\ge S)$ , $\mathcal{N}(t)=\{N(t_j);0\le t_j<t,j=1,\dots ,J\}$ denotes the history of the salvage therapy process. Without the loss of generality, we presume that the urologists’ consideration to initiate salvage therapy took place at the same $J$ time points the patients provided PSA measurements. Baseline covariate information is denoted as $\mathcal{X}=\{X_q=x_q,q=1,\dots ,Q\}$ .

To quantify the causal effect of salvage therapy, we will use the framework of counterfactual outcomes. Namely, we would like to decide if salvage therapy should be initiated at the follow-up time $t$ by comparing the counterfactual cumulative risks of metastasis under the two regimes in the medically relevant time interval $(t,t+\mathrm{\Delta }t]$ conditional on survival, no prior salvage therapy and no metastasis up to $t$ . In particular, we let $[T^{m(1)},{\mathcal{F}}^{(1)}(v,t)\mid T^m>t,T^d>t,\mathcal{Y}(t),\mathcal{N}(t)=0,\mathcal{X}]$ denote the joint conditional distribution of $T^{m(1)}$ and ${\mathcal{F}}^{(1)}(v,t)$ for all $v>t$ , where ${\mathcal{F}}^{(1)}(v,t)=\{Y(t_j);t<t_j\le v,j=1,\dots ,J^{″}\}$ denotes the $J^{″}$ future PSA measurements after $t$ and up to the horizon time $v$ , if salvage therapy was initiated at $t$ given the available PSA information up to $t$ and the fact that neither metastasis nor death has occurred up to this time. Analogously, we let $[T^{m(0)},{\mathcal{F}}^{(0)}(v,t)\mid T^m>t,T^d>t,\mathcal{Y}(t),\mathcal{N}(t)=0,\mathcal{X}]$ denote the joint conditional distribution of $T^{m(0)}$ and ${\mathcal{F}}^{(0)}(v,t)$ for all $v>t$ , if salvage therapy will not be initiated in the interval $(t,v]$ .

The marginal effect over all possible longitudinal histories up to time $t$ can be defined as

$$\begin{array}{rl}& {\text{ST}}^M(t+\mathrm{\Delta }t,t)\\ & =\int [Pr\{T^{m(1)}\le t+\mathrm{\Delta }t\mid T^m>t,T^d>t,\mathcal{Y}(t),\mathcal{X}\}\\ & -Pr\{T^{m(0)}\le t+\mathrm{\Delta }t\mid T^m>t,T^d>t,\mathcal{Y}(t),\mathcal{X}\}]\\ & \prod _{j=1}^{J^{\prime }}\mathrm{d}G\{Y(t_j),\mathcal{X}\mid \mathcal{Y}(t_{j-1})\}\\ & =Pr\{T^{m(1)}\le t+\mathrm{\Delta }t\mid T^m>t,T^d>t\}\\ & -Pr\{T^{m(0)}\le t+\mathrm{\Delta }t\mid T^m>t,T^d>t\}\end{array}$$

(1)

where $G(\cdot )$ denotes the cumulative joint distribution function for $Y$ and $\mathcal{X}$ , and $J^{\prime }$ denotes the number of longitudinal measurements up to $t$ . In this and the following expressions, and for simplicity of exposition, we have omitted that we condition on salvage not being initiated up to $t$ , that is, $\mathcal{N}(t)=0$ . Also, we have not included the term $\prod _{j=1}^{J^{\prime }}p\{N(t_j)\mid {\mathcal{Y}}_i(t_j),\mathcal{N}(t_{j-1}),\mathcal{X}\}$ because, as we will explain in Section 3.2, this can be ignored under our modeling approach. Even though this type of marginal effect is often used in the context of causal inference, it may be of less interest to urologists because they would typically decide to initiate salvage therapy for patients with elevated PSA. A more clinically interesting setting is to restrict to specific PSA values, that is, to condition on the longitudinal history ${\mathcal{Y}}_i(t)$ of a specific subject $i$ :

$$\begin{array}{rl}& {\text{ST}}^C(t+\mathrm{\Delta }t,t)\\ & =Pr\{T_i^{m(1)}\le t+\mathrm{\Delta }t\mid T_i^m>t,T_i^d>t,{\mathcal{Y}}_i(t),{\mathcal{X}}_i\}\\ & -Pr\{T_i^{m(0)}\le t+\mathrm{\Delta }t\mid T_i^m>t,T_i^d>t,{\mathcal{Y}}_i(t),{\mathcal{X}}_i\}\end{array}$$

(2)

This effect is conditional/individualized in the sense that it is for the group of patients with baseline variables ${\mathcal{X}}_i$ and longitudinal history ${\mathcal{Y}}_i(t)$ , that is, the same PSA values as for the patient the doctor wants to decide on initiating salvage therapy. Following the discussion by Taylor et al., ⁶ ${\text{ST}}^C(t+\mathrm{\Delta }t,t)$ is a conditional causal effect relevant for subject-specific treatment decisions. The difference between the two salvage therapy effects presented above is a bias versus variance tradeoff. In particular, the marginal effect (1) marginalizes over a bigger group of patients and will have a smaller variance than the other effects; however, as explained above, it will also be less relevant for the practicing urologist. The conditional effect (2) is the one most relevant to the doctor, but it will have a larger variance because it is based on very little data. A compromise between (1) and (2) is to quantify the salvage therapy effect for patients who had PSA levels above a threshold value $c$ at their last visit, that is, ${\mathcal{Y}}^{*}(t)=\{Y(t):Y(t)>c\}$ :

$$\begin{array}{rl}& {\text{ST}}^{MC}(t+\mathrm{\Delta }t,t)\\ & =\int [Pr\{T^{m(1)}\le t+\mathrm{\Delta }t\mid T^m>t,T^d>t,{\mathcal{Y}}^{*}(t)\}\\ & -Pr\{T^{m(0)}\le t+\mathrm{\Delta }t\mid T^m>t,T^d>t,{\mathcal{Y}}^{*}(t)\}]\\ & \prod _{j=1}^{J^{\prime }}\mathrm{d}\tilde{G}\{Y(t_j),\mathcal{X}\mid \mathcal{Y}(t_{j-1})\}\end{array}$$

(3)

where $\tilde{G}(\cdot )$ denotes the cumulative distribution function for the truncated distribution of $Y$ and $\mathcal{X}$ satisfying the condition ${\mathcal{Y}}^{*}(t)$ . This effect marginalizes over all possible longitudinal histories of PSA that end up having PSA greater than $c$ at time $t$ . The effect (3) offers a compromise by considering a more relevant group of patients, but smaller than the group considered in (1) (i.e. it will have a bigger variance than (1)). We should note that by changing the specification of ${\mathcal{Y}}^{*}(t)$ , alternative ST effects may be defined. For example, we could define the ST effect for patients with elevated PSA levels in the last $K$ months (e.g. $K=6$ ).

2.2. Assumptions

Because of the observational nature of the University of Michigan Prostatectomy Data, we will need a set of assumptions to identify and unbiasedly estimate the salvage therapy effects defined in the previous section. In particular, we will make the standard assumptions for causal inference with observational data, time-varying confounding, and intermediate variables.

-

Consistency: The observed outcomes equal the counterfactual outcomes for the actually assigned treatment.

-

Sequential exchangeability: The counterfactual outcomes are independent of the assigned treatment conditionally on the history of PSA measurements, the history of salvage therapy treatments, and baseline covariates, that is,

$$T^{m(a)}⫫N(t)\mid \mathcal{Y}(t),\mathcal{N}(t),\mathcal{X}$$

and

$${\mathcal{F}}^{(a)}(v)⫫N(t)\mid \mathcal{Y}(t),\mathcal{N}(t),\mathcal{X},v>t$$

where $a=\{0,1\}$ . Note that we assume here that the decision to start salvage therapy at time $t$ can only causally affect future PSA measurements $Y(v)$ , with $v>t$ but not $Y(t)$ .

The different salvage therapies are sufficiently well-defined in our motivating dataset not to cause any issues with the consistency assumption. The sequential exchangeability assumption also seems to be satisfied in our setting because, as mentioned earlier, urologists typically decide to initiate salvage therapy based on the history of PSA values and possibly other factors such as age and comorbidities (captured in $\mathcal{X}$ ). Because of the parametric nature of our modeling approach (presented in the following section), the positivity assumption is not required to identify and estimate the three causal effects we introduced above. However, for the marginal effect (1), we extrapolate beyond the support of the data because the urologists will not prescribe salvage medication for patients with a zero PSA value. Nonetheless, because this effect is of little practical interest, we are not concerned about this extrapolation.

3. Modeling

3.1. Sub-models specification

We will use the framework of joint models to associate the risk of metastasis with the longitudinal PSA measurements and account for the risk of death. Joint models will account for the endogenous nature of the PSA process, and under the set of assumptions presented in Section 2.2 will provide valid predictions of the risk of metastasis. We will use the subscript $i$ ( $i=1,\dots ,n$ ) to denote the subject in all the random variables introduced above. Also, we will use $T_i=min(T_i^m,T_i^d,C_i)$ to denote the observed event times with $C_i$ denoting the censoring time, and ${\delta }_i\in \{0,1,2\}$ , where “0” is for censoring, “1” for metastasis, and “2” for death. The deaths we consider as competing events are the ones that occurred before metastasis; deaths and PSA measurements after metastasis are ignored. As explained in the introduction, we expect a drop in PSA levels after the initiation of salvage therapy. Following previous literature in modeling PSA profiles for prostate cancer patients, we specify a linear mixed-effects model for the logarithm of PSA with a change point in the subject-specific trajectories after the initiation of salvage therapy:

$$\begin{array}{r}{y}_i(t)=\{\begin{array}{ll}{\eta }_i(t)+{\epsilon }_i(t)={\mathit{x}}_i{(t)}^{\mathrm{\top }}\mathit{\beta }+{\mathit{z}}_i{(t)}^{\mathrm{\top }}{\mathit{b}}_i+{\epsilon }_i(t),& t<S_i\\ {\tilde{\eta }}_i(t)+{\epsilon }_i(t)={\eta }_i(t)+\{{\tilde{\mathit{x}}}_i(\tilde{t}{)}^{\mathrm{\top }}\tilde{\mathit{\beta }}+{\tilde{\mathit{z}}}_i(\tilde{t}{)}^{\mathrm{\top }}{\tilde{\mathit{b}}}_i\}+{\epsilon }_i(t),& t\ge S_i\end{array}\end{array}$$

(4)

where $y_i(t)=\log \{{\text{PSA}}_i(t)+1\}$ the design vectors ${\mathit{x}}_i$ and ${\mathit{z}}_i$ for the fixed $\mathit{\beta }$ and random effects ${\mathit{b}}_i$ , respectively, describe the subject-specific PSA evolutions before salvage therapy. Analogously, the design vectors ${\tilde{\mathit{x}}}_i$ and ${\tilde{\mathit{z}}}_i$ for the fixed $\tilde{\mathit{\beta }}$ and random effects ${\tilde{\mathit{b}}}_i$ , respectively, describe the change in the subject-specific PSA evolutions after salvage therapy. The latter design vectors use the relative time variable $\tilde{t}=t-S_i$ . The covariates present in the four design vectors defined above are assumed to be a subset of ${\mathcal{X}}_i$ . The random effects is the model component that allows for different PSA profiles before and after salvage per patient. The distributional assumptions for the random effects and the error terms are ${\mathit{u}}_i=({\mathit{b}}_i^{\mathrm{\top }},{\tilde{\mathit{b}}}_i^{\mathrm{\top }}{)}^{\mathrm{\top }}\sim \mathcal{N}(\mathbf{0},\mathbf{\Omega })$ , ${\epsilon }_i(t)\sim \mathcal{N}(0,{\sigma }^2)$ , and $\text{cov}\{{\mathit{u}}_i,{\epsilon }_i(t)\}=0$ . The variance–covariance matrix $\mathbf{\Omega }$ is assumed completely unstructured.

Metastasis and death are considered competing risks, and we specify a different hazard model for each one. For metastasis, we postulate the model:

$$\begin{array}{rl}{h}_i^m(t)& =\underset{ϵ\to 0}{lim}{ϵ}^{-1}Pr\{t\le T_i^m<t+ϵ\mid T_i^m>t,T_i^d>t,{\mathcal{H}}_i(t),{\mathcal{N}}_i(t),{\mathcal{X}}_i\}\\ & =h_0^m(t)\exp ({\mathit{\psi }}_m^{\mathrm{\top }}{\mathit{w}}_i+{\gamma }_m{N}_i(t)+{\mathit{\alpha }}_m^{\mathrm{\top }}[\{1-N_i(t)\}\times f\{{\mathcal{H}}_i(t)\}]\\ & +{\mathit{\xi }}_m^{\mathrm{\top }}[N_i(t)\times g\{{\mathcal{H}}_i(t)\}])\end{array}$$

(5)

where $h_i^m(t)$ is the metastasis-specific hazard function for patient $i$ , and $h_0^m(t)$ is the baseline hazard. In our model, the logarithm of $h_0^m(t)$ is modeled using penalized B-splines, that is,

$$\log h_0^m(t)={\psi }_{h_m,0}+\sum _{q=1}^Q{\psi }_{h_m,q}{B}_q(t)$$

where $B_q(t)$ denotes the $q$ -th basis function of a B-spline with knots $v_1,\dots ,v_{Q+1}$ and ${\mathit{\psi }}_{h_m}$ the vector of spline coefficients. The design vector ${\mathit{w}}_i$ with the corresponding coefficients vector $\mathit{\psi }$ is for the baseline covariates (subset of ${\mathcal{X}}_i$ ) relevant to the risk of metastasis. The term ${\mathcal{H}}_i(t)=\{{\eta }_i(s);0\le s<min(S_i,t)\}\bigcup \{{\tilde{\eta }}_i(s);S_i\le s<t\}$ denotes the history of the subject-specific linear predictor and ${\mathcal{N}}_i(t)$ is defined as in Section 2.1. The vector functions $f(\cdot )$ and $g(\cdot )$ determine the functional form for the dependence of the hazard on the PSA evolutions before and after ST, for example, $f\{{\eta }_i(t)\}=[{\eta }_i(t),d{\eta }_i(t)/dt{]}^{\mathrm{\top }}$ or $f\{{\eta }_i(t)\}=[{\eta }_i(t),{\eta }_i(t)-{\eta }_i(t-1){]}^{\mathrm{\top }}$ . From this model, we can compute the hazard ratio for salvage therapy for subject $i$ , that is,

$$\text{HR}(t)=\exp ({\gamma }_m+{\mathit{\xi }}_m^{\mathrm{\top }}g\{{\tilde{\eta }}_i(t)\}-{\mathit{\alpha }}_m^{\mathrm{\top }}f\{{\eta }_i(t)\}),t>S_i$$

(6)

where in this expression ${\eta }_i(t)$ represents the expected value of $y_i(t)$ as if subject $i$ had not received salvage therapy at time $S_i$ . We note that the metastasis model (5) could be extended in several ways, including to include interaction terms between ${\mathit{w}}_i$ and $N_i(t)$ or by allowing ${\mathit{\psi }}_m$ to depend on $t$ .

We presume that the PSA is not associated with the hazard of death:

$$\begin{array}{rl}{h}_i^d(t)& =\underset{ϵ\to 0}{lim}{ϵ}^{-1}Pr\{t\le T_i^d<t+ϵ\mid T_i^m>t,T_i^d>t,{\mathcal{N}}_i(t),{\mathcal{X}}_i\}\\ & =h_0^d(t)\exp \{{\mathit{\psi }}_d^{\mathrm{\top }}{\mathit{w}}_i+{\gamma }_d{N}_i(t)\}\end{array}$$

where $h_i^d(t)$ is the death-specific hazard function for patient $i$ , and $h_0^d(t)$ is the baseline hazard that is again modeled using penalized B-splines with an associated spline coefficients vector ${\mathit{\psi }}_{h_d}$ . Likewise, interaction terms between ${\mathit{w}}_i$ and $N_i(t)$ could also be included.

3.2. Estimation

The longitudinal and event time processes are linked via the random effects to define their joint distribution. We will use a Bayesian approach to fit the postulated joint model. Inference proceeds via the posterior distribution of the parameters $\{{\mathit{u}}_i,\mathit{\theta };i=1,\dots ,n\}$ given the data $\mathcal{D}=\{T_i,{\delta }_i,{\mathit{Y}}_i;i=1,\dots ,n\}$ , where $\mathit{\theta }$ denotes the vector of model parameters. We start by formulating the joint posterior of $\{{\mathit{u}}_i,\mathit{\theta };i=1,\dots ,n\}$ and ${\mathit{\theta }}_N$ , with ${\mathit{\theta }}_N$ denoting the parameter vector for distribution of the salvage therapy assignment process. Using telescoping, we get

$$\begin{array}{rl}& p(\mathit{\theta },\mathit{u},{\mathit{\theta }}_N\mid \mathit{T},\mathit{\delta },\mathit{Y},\mathit{N})\\ & \propto \prod _{i=1}^n\prod _{j=1}^{n_i}p\{Y_i(t_{ij}),T_i,{\delta }_i\mid {\mathcal{Y}}_i(t_{i,j-1}),{\mathcal{N}}_i(t_{i,j-1}),{\mathcal{X}}_i,{\mathit{u}}_i,\mathit{\theta }\}\\ & \times \prod _{j=1}^{n_i}p\{N_i(t_{ij})\mid {\mathcal{Y}}_i(t_{i,j-1}),{\mathcal{N}}_i(t_{i,j-1}),Y_i(t_{ij}),T_i,{\delta }_i,{\mathcal{X}}_i,{\mathit{u}}_i,{\mathit{\theta }}_N\}\\ & \times p({\mathit{u}}_i\mid \mathit{\theta })\times p(\mathit{\theta })\times p({\mathit{\theta }}_N)\end{array}$$

where $\mathit{T}$ , $\mathit{\delta }$ , $\mathit{Y}$ , and $\mathit{N}$ denote the vectors of the event times, event indicators, longitudinal measurements, and treatments decisions, respectively, ${\mathcal{Y}}_i(t_{i,0})={\mathcal{N}}_i(t_{i,0})=\mathrm{\varnothing }$ , and $p(\cdot )$ denotes an appropriate probability density or probability mass function. Under sequential exchangeability, we have that

$$\begin{array}{rl}& p\{N_i(t_{ij})\mid {\mathcal{Y}}_i(t_{i,j-1}),{\mathcal{N}}_i(t_{i,j-1}),{\mathcal{F}}_i^{(a)}(v_{ij},t_{ij}),T_i^{(a)},{\delta }_i^{(a)},{\mathcal{X}}_i,{\mathit{u}}_i,{\mathit{\theta }}_N\}\\ & =p\{N_i(t_{ij})\mid {\mathcal{Y}}_i(t_{i,j-1}),{\mathcal{N}}_i(t_{i,j-1}),{\mathcal{X}}_i,{\mathit{\theta }}_N\}\end{array}$$

where ${\mathcal{F}}_i^{(a)}(v_{ij},t_{ij})$ denotes the future counterfactual PSA measurements for $v_{ij}>t_{ij}$ and $\{T_i^{(a)},{\delta }_i^{(a)}\}$ is the counterfactual event times. Hence, assuming that the parameters $\{{\mathit{u}}_i,\mathit{\theta };i=1,\dots ,n\}$ and ${\mathit{\theta }}_N$ are functionally independent, inference for the posterior distribution of the counterfactual outcomes $\{\mathit{\theta },\mathit{u}\mid {\mathit{T}}^{(a)},{\mathit{\delta }}^{(a)},{\mathit{Y}}^{(a)},{\mathit{N}}^{(a)}\}$ can be obtained using the first term (i.e. the observed data model) and ignore the second term. We should note that this decomposition allows the salvage therapy distribution to depend on a set of instrumental variables ${\mathcal{I}\mathcal{V}}_i$ that determine who gets therapy but are unrelated to the outcomes and the random effects (i.e. the $\mathcal{I}\mathcal{V}$ s are not a part of $\mathcal{X}$ ). An example could be the urologist’s preferences on when to start treatment. As noted in Section 1, different physicians may use in other manners the PSA history when deciding to whom they will give the therapy. Under this setting, we obtain:

$$\begin{array}{rl}& p(\mathit{\theta },\mathit{u}\mid \mathit{T},\mathit{\delta },\mathit{Y},\mathit{N})\\ & \propto \prod _{i=1}^n\{\prod _j\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp (-\frac{1}{2{\sigma }^2}\{y_{ij}-{\mu }_{ij}(\mathit{\theta },{\mathit{u}}_i,{\mathcal{N}}_i){\}}^2)\}\\ & \times \text{det}(2\pi \mathbf{\Omega }{)}^{-1/2}\exp (-\frac{1}{2}{\mathit{u}}_i^{\mathrm{\top }}{\mathbf{\Omega }}^{-1}{\mathit{u}}_i)\\ & \times \{h_i^m(T_i;\mathit{\theta },{\mathit{u}}_i,{\mathcal{N}}_i){\}}^{I({\delta }_i=1)}\{h_i^d(T_i;\mathit{\theta },{\mathcal{N}}_i){\}}^{I({\delta }_i=2)}\\ & \times \exp (-{\displaystyle {\int }_0^{T_i}\{h_i^m(s;\mathit{\theta },{\mathit{u}}_i,{\mathcal{N}}_i)+h_i^d(s;\mathit{\theta },{\mathcal{N}}_i)\}\mathrm{d}s)}\\ & \times p(\mathit{\theta })\end{array}$$

where ${\mu }_{ij}(\mathit{\theta },{\mathit{u}}_i,{\mathcal{N}}_i)$ denotes the mean of the linear mixed model (4), and $\text{det}(A)$ is the determinant of matrix $A$ . We use standard priors for $\mathit{\theta }$ , that is, normal priors for all regression coefficients $(\mathit{\beta },\tilde{\mathit{\beta }},{\mathit{\psi }}_{h_m},{\mathit{\psi }}_m,{\gamma }_m,{\mathit{\alpha }}_m,{\mathit{\xi }}_m,{\mathit{\psi }}_{h_d},{\mathit{\psi }}_d,{\gamma }_d)$ , inverse-Gamma priors for ${\sigma }^2$ and the diagonal elements of $\mathbf{\Omega }$ , and the LKJ prior for the correlation matrix of the random effects. ¹¹ To ensure smoothness of the baseline hazard functions $h_0^m(t)$ and $h_0^d(t)$ , we postulate a “penalized” prior distribution for the regression coefficients ${\mathit{\psi }}_{h_m}$ and ${\mathit{\psi }}_{h_d}$ (we only show the formulation for the former):

$$p({\mathit{\psi }}_{h_m}\mid {\tau }_m)\propto {\tau }_m^{\rho (K)/2}\exp (-\frac{{\tau }_m}{2}{\mathit{\psi }}_{h_m}^{\mathrm{\top }}\mathit{K}{\mathit{\psi }}_{h_m})$$

where ${\tau }_m$ is the smoothing parameter that takes a $\text{Gamma}(5,0.05)$ hyper-prior in order to ensure a proper posterior for ${\mathit{\psi }}_{h_m}$ , $\mathit{K}={\mathrm{\Delta }}_r^{\mathrm{\top }}{\mathrm{\Delta }}_r$ , where ${\mathrm{\Delta }}_r$ denotes $r$ -th difference penalty matrix, and $\rho (\mathit{K})$ denotes the rank of $\mathit{K}$ .

We use a Markov chain Monte Carlo (MCMC) approach to obtain samples from the posterior distribution for all model parameters and the random effects. This algorithm is implemented in the R package JMbayes2 ¹² that we used to fit the model and calculate the salvage therapy effects.

4. Salvage therapy effects estimation

4.1. Estimates

We start with the estimation of the conditional effect (2) that we calculate for a specific patient with longitudinal history ${\mathcal{Y}}_i(t)$ and covariates ${\mathcal{X}}_i$ . Both terms in the definition of this effect are posterior predictive distributions, which under the postulated joint model are written as

$$\begin{array}{rl}& Pr\{T_i^{(a)}\le t+\mathrm{\Delta }t,{\delta }_i^{(a)}=1\mid T_i>t,{\mathcal{Y}}_i(t),{\mathcal{X}}_i\}\\ & ={\displaystyle \int {\displaystyle \int Pr\{T_i^{(a)}\le t+\mathrm{\Delta }t,{\delta }_i^{(a)}=1\mid T_i>t,{\mathit{u}}_i,{\mathcal{X}}_i,\mathit{\theta }\}}}\\ & \times p\{{\mathit{u}}_i\mid T_i>t,{\mathcal{Y}}_i(t),{\mathcal{X}}_i,\mathit{\theta }\}p(\mathit{\theta }\mid \mathcal{D})\mathrm{d}{\mathit{u}}_i\mathrm{d}\mathit{\theta }\end{array}$$

(7)

where $a=\{0,1\}$ , and the term $p(\mathit{\theta }\mid \mathcal{D})$ denotes the posterior distribution for the model parameters $\mathit{\theta }$ . Again for simplicity of exposition, we have omitted the conditioning on $\mathcal{N}(t)=0$ . We should note that these risk predictions are marginalized over both the parameters and the random effects. The first term in the integrand is written as:

$$Pr\{T_i^{(a)}\le t+\mathrm{\Delta }t,{\delta }_i^{(a)}=1\mid T_i>t,{\mathit{u}}_i,{\mathcal{X}}_i,\mathit{\theta }\}=\mathcal{A}/\mathcal{B}$$

where

$$\begin{array}{rl}& \mathcal{A}={\displaystyle {\int }_t^{t+\mathrm{\Delta }t}{h}_i^{m(a)}(v)\exp (-{\displaystyle {\int }_t^v\{h_i^{m(a)}(s)+h_i^{d(a)}(s)\}\mathrm{d}s}}\\ & -{\displaystyle {\int }_0^t\{h_i^{m(0)}(s)+h_i^{d(0)}(s)\}\mathrm{d}s)\mathrm{d}v}\end{array}$$

and

$$\mathcal{B}={\displaystyle \exp (-{\displaystyle {\int }_0^t\{h_i^{m(0)}(s)+h_i^{d(0)}(s)\}\mathrm{d}s)}}$$

with $h_i^{m(1)}(t)=h_0^m(t)\exp ({\mathit{\psi }}_m^{\mathrm{\top }}{\mathit{w}}_i+{\gamma }_m+{\mathit{\xi }}_m^{\mathrm{\top }}g\{{\tilde{\eta }}_i(t)\})$ and $h_i^{m(0)}(t)=h_0^m(t)\exp ({\mathit{\psi }}_m^{\mathrm{\top }}{\mathit{w}}_i+{\mathit{\alpha }}_m^{\mathrm{\top }}f\{{\eta }_i(t)\})$ . The counterfactual hazard functions for death $h_i^{d(1)}(t)$ and $h_i^{d(0)}(t)$ are defined analogously. In particular, $h_i^{d(1)}(t)=h_0^d(t)\exp ({\mathit{\psi }}_d^{\mathrm{\top }}{\mathit{w}}_i+{\gamma }_d)$ and $h_i^{d(0)}(t)=h_0^d(t)\exp ({\mathit{\psi }}_d^{\mathrm{\top }}{\mathit{w}}_i)$ . In the specification of the conditional distribution of the random effects given the observed information $\{T_i>t,{\mathcal{Y}}_i(t),{\mathcal{X}}_i\}$ , it is assumed that salvage therapy has not been initiated by time $t$ . Combining the integral equations presented above, we can estimate (7) using the following Monte-Carlo scheme:

• S1

  We sample ${\breve{\mathit{\theta }}}^{(l)}$ from the MCMC sample of the posterior distribution $[\mathit{\theta }\mid \mathcal{D}]$ .

• S2

  We sample ${\breve{\mathit{u}}}_i^{(l)}$ from the posterior distribution $[{\mathit{u}}_i\mid T_i>t,{\mathcal{Y}}_i(t),{\mathcal{X}}_i,{\breve{\mathit{\theta }}}^{(l)}]$ . Because this distribution cannot be written in closed-form, we sample from it using the Metropolis-Hastings algorithm.

• S3

  We calculate the term ${\pi }_i^{(l)}(t+\mathrm{\Delta }t\mid t,a)=Pr\{T_i^{(a)}\le t+\mathrm{\Delta }t,{\delta }_i^{(a)}=1\mid T_i>t,{\breve{\mathit{u}}}_i^{(l)},{\mathcal{X}}_i,{\breve{\mathit{\theta }}}^{(l)}\}$ . The integrals in the definition of the overall survival functions are approximated numerically using the 15-point Gauss-Kronrod quadrature rule.

We repeat Steps S1–S3, $L$ times, and we take as an estimate of ${\text{ST}}_i^C(t+\mathrm{\Delta }t,t)$ , the mean over the Monte-Carlo samples, that is,

$${\widehat{\text{ST}}}_i^C(t+\mathrm{\Delta }t,t)=\frac{1}{L}\sum _{l=1}^L{\pi }_i^{(l)}(t+\mathrm{\Delta }t\mid t,a=1)-{\pi }_i^{(l)}(t+\mathrm{\Delta }t\mid t,a=0)$$

(8)

The estimation of the marginal effects (1) and (3) proceeds by averaging the conditional effects over the respective groups of patients in the sample. In particular, for (1), we define $\mathcal{R}(t)$ to denote the subset of patients at risk at time $t$ . For each patient in $\mathcal{R}(t)$ , we calculate ${\widehat{\text{ST}}}_i^C(t+\mathrm{\Delta }t,t)$ . Then, we obtain the estimator

$${\widehat{\text{ST}}}^M(t+\mathrm{\Delta }t,t)=n_r^{-1}\sum _{i:i\in R(t)}{\widehat{\text{ST}}}_i^C(t+\mathrm{\Delta }t,t)$$

where $n_r$ denotes the number of subjects in $\mathcal{R}(t)$ . For (3), the summation is over ${\mathcal{R}}^{*}(t)$ that denotes the subset of patients at risk at time $t$ and had longitudinal histories that satisfy the definition of ${\mathcal{Y}}^{*}(t)$ ( $n_{r^{*}}$ is defined analogously).

4.2. Variance

To derive the variance of ${\widehat{\text{ST}}}^M(t+\mathrm{\Delta }t,t)$ and ${\widehat{\text{ST}}}^{MC}(t+\mathrm{\Delta }t,t)$ , we need to take into account that they are a function of both the parameters $\mathit{\theta }$ and the data $\mathcal{D}$ (in the definition of $\mathcal{D}$ here we also include $\{{\mathcal{X}}_i,i=1,\dots ,n_r\}$ ). The variance of both estimators can be derived using similar arguments, and here we will only show how to calculate the variance of the former. To make the dependence on the data $\mathcal{D}$ more explicit, we write the estimator of the marginal salvage therapy effect as

$${\widehat{\text{ST}}}^M(t+\mathrm{\Delta }t,t;\mathcal{D})=E_{\mathit{\theta }\mid \mathcal{D}}\{{\text{ST}}^M(t+\mathrm{\Delta }t,t;\mathit{\theta },\mathcal{D})\}$$

To account for both the variability in the parameters and the sampling variability, our target variance is

$$\begin{array}{rl}& {\text{var}}_{\mathcal{D}}\{{\widehat{\text{ST}}}^M(t+\mathrm{\Delta }t,t;\mathit{\theta },\mathcal{D})\}\\ & ={\text{var}}_{\mathcal{D}}[E_{\mathit{\theta }\mid \mathcal{D}}\{{\text{ST}}^M(t+\mathrm{\Delta }t,t;\mathit{\theta },\mathcal{D})\}]\end{array}$$

(9)

We will adapt the approach of Antonelli et al. ¹³ to obtain an estimate of (9). More specifically, we let $\{{\mathcal{D}}^{(1)},\dots ,{\mathcal{D}}^{(M)}\}$ to denote $M$ datasets sampled with replacement from $\mathcal{D}$ . For each of these datasets, we calculate

$$\begin{array}{rl}& {\widehat{\text{ST}}}^M(t+\mathrm{\Delta }t,t;\mathit{\theta },{\mathcal{D}}^{(m)})\\ & =E_{\mathit{\theta }\mid \mathcal{D}}\{{\text{ST}}^M(t+\mathrm{\Delta }t,t;\mathit{\theta },{\mathcal{D}}^{(m)})\},m=1,\dots ,M\end{array}$$

Note that the expectation is taken with respect to the posterior distribution of the parameters using the original data $\mathcal{D}$ (i.e. we do not refit the model for each dataset ${\mathcal{D}}^{(m)}$ ). Calculating the sample variance of these $M$ estimates we obtain

$${\text{var}}_{{\mathcal{D}}^{(m)}}[E_{\mathit{\theta }\mid \mathcal{D}}\{{\text{ST}}^M(t+\mathrm{\Delta }t,t;\mathit{\theta },{\mathcal{D}}^{(m)})\}]$$

(10)

Even though (10) resembles our target variance, it ignores the variability in the posterior due to the different samples $\{{\mathcal{D}}^{(1)},\dots ,{\mathcal{D}}^{(M)}\}$ . Hence, to get correct inferences, we use the correction term:

$$\begin{array}{rl}& {\widehat{\text{var}}}_{\mathcal{D}}[E_{\mathit{\theta }\mid \mathcal{D}}\{{\text{ST}}^M(t+\mathrm{\Delta }t,t;\mathit{\theta },\mathcal{D})\}]\\ & ={\text{var}}_{{\mathcal{D}}^{(m)}}[E_{\mathit{\theta }\mid \mathcal{D}}\{{\text{ST}}^M(t+\mathrm{\Delta }t,t;\mathit{\theta },{\mathcal{D}}^{(m)})\}]\\ & +{\text{var}}_{\mathit{\theta }\mid \mathcal{D}}\{{\text{ST}}^M(t+\mathrm{\Delta }t,t;\mathit{\theta },\mathcal{D})\}\end{array}$$

(11)

To derive the variance of ${\widehat{\text{ST}}}_i^C(t+\mathrm{\Delta }t,t)$ we need to account for the sampling variability in ${\mathcal{Y}}_i(t)$ , that is, the PSA values patient $i$ would have shown if we “cloned” him. To account for this variability, we cannot use the same idea as in ${\widehat{\text{ST}}}^M(t+\mathrm{\Delta }t,t)$ above because we cannot obtain samples with replacement from ${\mathcal{Y}}_i(t)$ . As an alternative, we employ a parametric Bootstrap approach,that is, we create different versions of the history $\{{\mathcal{Y}}_i^{(1)}(t),\dots ,{\mathcal{Y}}_i^{(M)}(t)\}$ using the Monte-Carlo scheme

• Q1

  We sample ${\ddot{\mathit{\theta }}}^{(l)}$ from the MCMC sample of the posterior distribution $[\mathit{\theta }\mid \mathcal{D}]$ .

• Q2

  We sample ${\ddot{\mathit{u}}}_i^{(l)}$ from the posterior distribution $[{\mathit{u}}_i\mid T_i>t,{\mathcal{Y}}_i(t),{\mathcal{X}}_i,{\ddot{\mathit{\theta }}}^{(l)}]$ .

• Q3

  We simulate ${\mathcal{Y}}_i^{(m)}(t)$ using independent draws from $[y_i\mid {\ddot{\mathit{u}}}_i,{\mathcal{X}}_i,{\ddot{\mathit{\theta }}}^{(l)}]$ , that is from the linear mixed model (4) with $t<S_i$ .

Subsequently, we use the same procedure as for the variance of ${\widehat{\text{ST}}}^M(t+\mathrm{\Delta }t,t;\mathcal{D})$ . Namely, along the lines of (11) we use

$$\begin{array}{rl}& {\widehat{\text{var}}}_{{\mathcal{Y}}_i}[E_{\mathit{\theta }\mid {\mathcal{Y}}_i}\{{\text{ST}}^C(t+\mathrm{\Delta }t,t;\mathit{\theta },{\mathcal{Y}}_i(t))\}]\\ & ={\text{var}}_{{\mathcal{Y}}_i^{(m)}}[E_{\mathit{\theta }\mid {\mathcal{Y}}_i}\{{\text{ST}}^C(t+\mathrm{\Delta }t,t;\mathit{\theta },{\mathcal{Y}}_i^{(m)}(t))\}]\\ & +{\text{var}}_{\mathit{\theta }\mid {\mathcal{Y}}_i}\{{\text{ST}}^C(t+\mathrm{\Delta }t,t;\mathit{\theta },{\mathcal{Y}}_i(t))\}\end{array}$$

Again, the first term accounts for the variability in ${\mathcal{Y}}_i(t)$ and the second for the parameters’ variability. The calculation of the first term entails obtaining the Monte Carlo estimate (8) for each realization ${\mathcal{Y}}_i^{(m)}(t)$ (i.e. we have nested Monte-Carlo schemes). Because of its parametric nature, this estimator relies more heavily on the model’s formulation.

An example on how package JMbayes2 can be used to estimate these causal effects and their variances is given in the following URL: https://drizopoulos.github.io/JMbayes2/articles/Causal_Effects.html. The running time for this example in a laptop with an Intel(R) Core(TM) i9-10885H CPU @ 2.40 GHz with eight physical cores and 32.0 GB RAM with a Microsoft Windows 11 Pro operating system is 40 s.

5. University of Michigan Prostatectomy Data – Analysis

We return to the UMP dataset, which we will use to estimate the salvage therapy effects introduced in Section 2. More information on the dataset and some pre-processing we applied is given in Supplemental Section 1.1. In the original version of the dataset, some patients received salvage therapy multiple times. For those patients, we have only considered the first time they received salvage therapy and used in the analysis only the PSA measurements taken up to 1.5 years after salvage therapy.

The joint model we fitted to the final dataset had the following specification. For the PSA trajectories, we use a modified version of the generic model (4), that is,

$$\log \{{\text{PSA}}_i(t_{ij})+1\}=\{\begin{array}{l}{\eta }_i(t_{ij})+{\epsilon }_i(t_{ij})={\beta }_{0i}+\sum _{k=1}^8{\beta }_{ki}{\mathcal{B}}_k(t_{ij},v)+{\mathit{x}}_{base,i}^{\mathrm{\top }}\mathit{\lambda }+{\epsilon }_i(t_{ij}),t_{ij}\le S_i\\ \widetilde{{\eta }_i(t_{ij})}+{\epsilon }_i(t_{ij})={\eta }_i(t_{ij})+{\tilde{\beta }}_{0i}+{\tilde{\beta }}_{1i}(t_{ij}-S_i)+{\mathit{x}}_{base,i}^{\mathrm{\top }}\tilde{\mathit{\lambda }}+{\epsilon }_i(t_{ij}),t_{ij}>S_i\end{array}$$

The subject-specific coefficients, ${\beta }_{0i},\dots ,{\beta }_{8i}$ , ${\tilde{\beta }}_{0i}$ , and ${\tilde{\beta }}_{1i}$ are decomposed into a fixed and random part. The combined random-effects vector ${\mathit{b}}_i=(b_{0i},\dots ,b_{8i},{\tilde{b}}_{0i},{\tilde{b}}_{1i}{)}^{\mathrm{\top }}$ is assumed to follow a multivariate normal distribution with mean zero and covariance matrix $\mathbf{\Omega }$ , and is independent of the error terms ${\epsilon }_{ij}$ that follow a normal distribution with mean zero and variance ${\sigma }^2$ . In both branches of the model, and via the terms ${\mathit{x}}_{base,i}^{\mathrm{\top }}\mathit{\lambda }$ and ${\mathit{x}}_{base,i}^{\mathrm{\top }}\tilde{\mathit{\lambda }}$ we control for the baseline covariates age, baseline PSA, Gleason score, and the Charlson comorbidity index. After prostatectomy, most patients will exhibit close to zero PSA levels; however, for some patients, at some point, PSA will start to rapidly increase, triggering the urologists to initiate salvage therapy. The time PSA will start increasing differs per patient, requiring a flexible model to capture these evolutions. Hence, we postulate a nonlinear model with a cubic B-spline for the time variable, with internal knots placed at $v=\{0.5,1,3,5,7,9\}$ years after prostatectomy and boundary knots placed at 0 and 17.1 years after prostatectomy (i.e. the terms ${\mathcal{B}}_k(t_{ij},v)$ ). After salvage, we assume that the PSA will drop by ${\tilde{\beta }}_{0i}$ and then have a modified linear increase before metastasis. We allow for a nonlinear PSA profile during follow-up with a change in the linear slope after salvage because of the limited available information in the data for $t>S_i$ . Before salvage, patients had a median of six PSA measurements (IQR = 6), whereas, after salvage, the median of available PSA measurements was three (IQR = 0). We should note that our model does not account for undetectable PSA levels; in the dataset these were set to zero. A better approach would be to account for the limit of detection in the formulation of the likelihood of the longitudinal submodel by using the cumulative distribution function instead of the probability density function for these values.

For the hazard of metastasis, we postulate the relative risk models

$$h_i^m(t)=h_0^m(t)\exp \{{\mathit{\psi }}_m^{\mathrm{\top }}{\mathit{w}}_{mi}+\mathcal{A}(t)\}$$

where for the time-varying component $\mathcal{A}(t)$ we consider the following versions

$$\begin{array}{rl}{M}_1:\mathcal{A}(t)& =\{\begin{array}{ll}{\alpha }_{m1}{\eta }_i(t),& t\le S_i\\ {\gamma }_{m1}(t-S_i)+{\gamma }_{m2}\{(t-S_i)\times {\mathtt{basePSA}}_i\}+{\xi }_{m1}{\tilde{\eta }}_i(t),& t>S_i\end{array}\\ M_2:\mathcal{A}(t)& =\{\begin{array}{ll}{\alpha }_{m1}{\eta }_i(t)+{\alpha }_{m2}{\displaystyle \frac{\mathrm{d}{\eta }_i(t)}{\mathrm{d}t},}& t\le S_i\\ {\gamma }_{m1}(t-S_i)+{\gamma }_{m2}\{(t-S_i)\times {\mathtt{basePSA}}_i\}+{\xi }_{m1}{\tilde{\eta }}_i(t),& t>S_i\end{array}\\ M_3:\mathcal{A}(t)& =\{\begin{array}{ll}{\alpha }_{m1}{\eta }_i(t)+{\alpha }_{m2}{\displaystyle \frac{\mathrm{d}{\eta }_i(t)}{\mathrm{d}t},}& t\le S_i\\ {\gamma }_{m1}(t-S_i)+{\gamma }_{m2}\{(t-S_i)\times {\mathtt{basePSA}}_i\}+{\xi }_{m2}{\tilde{\beta }}_{0i}+{\xi }_{m3}{\tilde{\beta }}_{1i},& t>S_i\end{array}\\ M_4:\mathcal{A}(t)& =\{\begin{array}{ll}{\alpha }_{m1}{\eta }_i(t)+{\alpha }_{m3}\{{\int }_0^t{\eta }_i(v)dv/t\},& t\le S_i\\ {\gamma }_{m1}(t-S_i)+{\gamma }_{m2}\{(t-S_i)\times {\mathtt{basePSA}}_i\}+{\xi }_{m2}{\tilde{\beta }}_{0i}+{\xi }_{m3}{\tilde{\beta }}_{1i},& t>S_i\end{array}\end{array}$$

Coefficient ${\alpha }_{m1}$ , ${\alpha }_{m2}$ , and ${\alpha }_{m3}$ quantify the association between the current value, current slope, and average $\log (\mathtt{PSA}+1)$ in the period before salvage and the hazard of metastasis, respectively. Coefficient ${\xi }_{m1}$ quantifies the association of the current value of $\log (\mathtt{PSA}+1)$ after salvage and the instantaneous risk of metastasis; coefficients ${\xi }_{m2}$ and ${\xi }_{m3}$ quantify the association between the drop of PSA just after salvage, and the change in the linear slope of $\log (\mathtt{PSA}+1)$ after salvage and the instantaneous risk of metastasis, respectively. Also, the coefficient ${\gamma }_{m1}$ quantifies the association between the length of the period the patient was on salvage therapy and the hazard of metastasis, and coefficient ${\gamma }_{m2}$ the interaction between the time a patient has been on salvage and baseline PSA. The baseline covariates part ${\mathit{\psi }}_m^{\mathrm{\top }}{\mathit{w}}_{mi}$ includes the effects of age, baseline PSA, the Gleason score, and the Charlson comorbidity index. The baseline hazard is approximated with B-splines, as explained earlier. For the hazard of death, we assume the model

$$h_i^d(t)=h_0^d(t)\exp \{{\psi }_{d1}(\mathtt{Age}\mathtt{-}\mathtt{50})+{\gamma }_{d1}{N}_i(t)\}$$

The prior distributions we assumed can be found in the Supplemental Material. Samples from the posterior distribution of the model parameters from these four joint models have been obtained from the MCMC algorithm provided in the R package JMbayes2 using three parallel chains started from different initial values. Each chain run for 55,000 iterations, discarding 5000 iterations as burn-in, and applying a thinning of 10 iterations. Hence, each chain contained 5000 realization of all model parameters. The priors used for the various parameters were the same as in Section 3.2. Trace-plots and the potential scale reduction factor ( $\hat{R}$ ) showed satisfactory convergence of the MCMC algorithm. In particular, for all model parameters the $\hat{R}$ -values were smaller than 1.09. We have also performed a sensitivity analysis for the choice of the prior distributions. We found that for almost all parameters the choice of the prior had little influence in the results. For the only parameters for which the impact was greater was for ${\tau }_m$ and ${\tau }_d$ that control the shape of the baseline hazard functions for metastasis and death, respectively.

We first evaluate the model fit before focusing on the results. All fitted joint models showed a similar fit to the longitudinal data; therefore, we only show the results from model $M_1$ . Supplemental Figures 3 to 5 show the observed data and fitted longitudinal trajectories for 48 selected patients from the dataset. These patients have been selected to showcase different profiles seen in the data, including patients who did and did not initiate salvage therapy. We observe from these figures that the postulated longitudinal submodel provides a good fit to the observed data. As previously observed, ¹⁴ this is an important aspect in joint models. In addition, using residuals plots we assessed the normality and homoscedastiscity assumptions for the error terms; these plots did not indicate any violations of these assumptions. Table 1 shows the Deviance information criterion (DIC) and the Watanabe-Akaike information criterion (WAIC) for the four fitted joint models.

Table 1.

Deviance information criterion (DIC), the Watanabe-Akaike information criterion (WAIC), and the log pseudo marginal likelihood (LPML) for the fitted joint models $M_1-M_4$ .

	DIC	WAIC	LPML
$M_4$	23505.49	12902244.27	−130702.53
$M_3$	25920.73	17954385.81	−177836.90
$M_1$	43554.78	33007334.39	−263741.71
$M_2$	91596.92	44668135.08	−270099.58

According to both WAIC and DIC, Model $M_4$ provides better predictions compared to the other three models.

For the event process, we show in Supplemental Table 2 the posterior means and the corresponding 95% credible intervals for the coefficients of the relative risk models for metastasis and death. All models suggest that before initiating salvage, high PSA levels at time $t$ , high PSA velocity at $t$ , and high average PSA from baseline to $t$ translate to a greater risk of metastasis at $t$ . Models $M_1$ and $M_2$ indicate that after salvage has been initiated, high PSA levels and high PSA velocity at time $t$ translate to a higher risk of metastasis. Models $M_3$ and $M_4$ suggest that the greater the drop in PSA levels after salvage initiation, the lower the risk of metastasis afterward, and the greater the PSA velocity after salvage, the higher the risk of metastasis. For the longitudinal process, the corresponding coefficients are shown in Supplemental Table 3. We observe that the baseline PSA, the Gleason score, and the Charlson comorbidity index seem to be associated with the level of post-surgery PSA. There is no indication that the effect of these confounders change after the initiation of salvage therapy.

We now turn our focus on estimating causal salvage therapy effects from the fitted joint models. To contrast the different effect types, we show the conditional effect (2) for Patients 490 and 327, the marginal-conditional effect (3) that averages over the group of patients who at their last visit has a PSA value > 0.5 ng/mL, and the marginal effect (1) that averages over all patients irrespective of their PSA values. These effects have been calculated at $t=5,7,9$ , and $13$ years after the initial surgery and for a $\mathrm{\Delta }t=2$ -year medically relevant time window. The 95% confidence intervals are calculated using the variance estimate presented in Section 4.2. The number of patients based on which the marginal and marginal-conditional effects have been calculated is shown in Supplemental Table 1. Figures 1 to 4 present these results.

Figure 1.

Salvage therapy effects for follow-up times $t=5,7,9$ , and $13$ years and $\mathrm{\Delta }t=2$ under model $M_1$ . For Patients 490 and 327, conditional causal effects are shown. The marginal-conditional causal effect is for patients who had a prostate-specific antigen (PSA) value greater or equal to 0.5 ng/mL at their last visit. The marginal effect is for all patients at risk at the corresponding $t$ .

Figure 2.

Salvage therapy effects for follow-up times $t=5,7,9$ , and $13$ years and $\mathrm{\Delta }t=2$ under model $M_2$ . For Patients 490 and 327, conditional causal effects are shown. The marginal-conditional causal effect is for patients who, at their last visit, had a prostate-specific antigen (PSA) value greater or equal to 0.5 ng/mL. The marginal effect is for all patients at risk at the corresponding $t$ .

Figure 3.

Salvage therapy effects for follow-up times $t=5,7,9$ , and $13$ years and $\mathrm{\Delta }t=2$ under model $M_3$ . For Patients 490 and 327, conditional causal effects are shown. The marginal-conditional causal effect is for patients who had a prostate-specific antigen (PSA) value greater or equal to 0.5 ng/mL at their last visit. The marginal effect is for all patients at risk at the corresponding $t$ .

Figure 4.

Salvage therapy effects for follow-up times $t=5,7,9$ , and $13$ years and $\mathrm{\Delta }t=2$ under model $M_4$ . For Patients 490 and 327, conditional causal effects are shown. The marginal-conditional causal effect is for patients who had a prostate-specific antigen (PSA) value greater or equal to 0.5 ng/mL at their last visit. The marginal effect is for all patients at risk at the corresponding $t$ .

The estimated values of the conditional and marginal-conditional salvage therapy effects tend to be negative, as expected, since salvage therapy is considered to be beneficial but are also small in magnitude. They are small because the probability of metastasis in the next two years is itself small, even in the absence of salvage therapy, so the room for improvement is also small. The results for the salvage therapy effects and across models align with the points raised in Section 2.1. In particular, we observe that the marginal effects are the smallest in magnitude, followed by the marginal-conditional and conditional effects. The reverse is observed regarding the variance of these effects, with the marginal effects having the smallest variance and the conditional effects the largest. Also, we observe that the conditional effects are more adaptable to the shape of the PSA profile. For example, while the marginal and marginal-conditional effects become smaller in magnitude at later follow-up times, the conditional effect for Patient 490 increases in size because this subject shows a steeply increasing PSA trajectory. This also happens for Patient 327 but to a lesser degree because his PSA profile is less steep. The variance of the conditional effects also increases with increasing PSA values. This reflects the fact that because in the sample, the majority of patients showed stable PSA profiles close to zero, the model is less “certain” for the shape of increasing PSA trajectories.

In this section, we estimated the marginal-conditional effect at time $t$ , by averaging the causal effects of individuals who had their last PSA > 0.5. A different marginal-conditional effect would have been obtained if we had used different criteria for who to average over. For example, we could have used a range of PSA, say 0.5–4.0. The criteria could have also included requirements such as Gleason score was at least 8 and age was < 75. We might also exclude from the set any patient whose last PSA was not current, for example, was more than 2 years ago. Restrictions such as these could make the causal estimates more targeted to an individual patient, but at the expense of greater variance.

6. Simulation

We performed a simulation study to assess the performance of the approaches presented in the previous sections. Note that because the causal effects (1)–(3) entail most of the model parameters in their definition, it is very challenging to simulate data with specific values for the causal effects. This prevents us from setting up a simulation to assess whether specific values for the causal effects can be unbiasedly estimated. However, the key assumption behind our approach is that joint models can be unbiasedly estimated under time-varying confounding and do not require to include a model for the treatment initiation process. This property of joint models has not been established before to our knowledge. Therefore, we focus here on validating the finite sample performance of joint models when the decision to initiate salvage therapy heavily depends on past PSA values. By setting the decision to perform salvage therapy on previous PSA values, we thus, mimic a causal relationship between the decision to perform salvage therapy and the history of PSA. Furthermore, to reduce the computational burden, we opted for a simplified setting with no additional covariates and linear time evolutions. Also, when we simulate PSA values we have not imposed the constraint that they should be positive.

More specifically, we assumed 1000 subjects and then randomly selected follow-up visits, $t_{ij}$ from a uniform distribution between 0 and 20. To mimic the causal relationship between salvage therapy and PSA history, the timing of salvage therapy was assumed to depend on the value of PSA. Specifically, if PSA was smaller than 2 ng/mL, the probability of receiving salvage at that visit time is 0.01; when PSA was in the interval (2 and 4 ng/mL), the probability of receiving treatment was 0.5, and for PSA > 4 ng/mL, the probability was set to 0.9. If salvage therapy was not given according to the model described previously, the subject was assumed not to undergo salvage therapy but still be at risk for metastasis or death. More details for the settings of the simulation study are given in Supplemental Section 2. We then simulated 300 datasets under three different scenarios for the association structure between features of the longitudinal PSA values and the hazard of metastasis. In Scenario 1, we considered an association with the current value of PSA; in Scenario 2, the current value and the current slope before salvage therapy and only the current value after salvage therapy; and in Scenario 3, the value and cumulative effect before salvage therapy and the current value after salvage therapy. We assumed no association between the PSA history and the hazard of death. The detailed model specification and parameter values are provided in the Supplemental Material, along with visualizations of the simulated data in Supplemental Figures 6 to 8.

The simulation study results are presented in Table 2.

Table 2.

Mean, bias, and root mean square error (RMSE) for the regression coefficients from the three different simulation scenarios.

	Scenario 1			Scenario 2			Scenario 3
Parameter	Mean	Bias	RMSE	Mean	Bias	RMSE	Mean	Bias	RMSE
${\beta }_0$	−1.745	−0.005	−0.036	−1.745	−0.005	−0.036	−1.745	−0.005	−0.036
${\beta }_1$	−0.038	−0.004	−0.012	−0.039	−0.005	−0.012	−0.038	−0.005	−0.012
${\tilde{\beta }}_0$	−5.840	−0.000	−0.065	−5.838	−0.002	−0.066	−5.838	−0.002	−0.065
${\tilde{\beta }}_1$	−0.194	−0.012	−0.016	−0.194	−0.013	−0.017	−0.194	−0.013	−0.017
${\gamma }_{d1}$	−0.250	−0.000	−0.223	−0.246	−0.004	−0.233	−0.238	−0.012	−0.265
${\gamma }_{m1}$	−0.266	−0.016	−0.150	−0.264	−0.014	−0.150	−0.262	−0.012	−0.150
${\alpha }_{m1}$	−0.810	−0.035	−0.096	−0.808	−0.037	−0.127	−0.780	−0.065	−0.165
${\alpha }_{m2}$	–	–	–	−1.048	−0.152	−1.816	–	–	–
${\alpha }_{m3}$	–	–	–	–	–	–	−0.558	−0.058	−0.267
${\xi }_{m1}$	−0.603	−0.008	−0.027	−0.605	−0.010	−0.027	−0.604	−0.009	−0.026

A detailed discussion of these results is given in Supplemental Section 2.2. In general, the results suggest that joint models can unbiasedly estimate the parameters of the joint distribution of the longitudinal and event time outcomes in the presence of time-varying confounding. Hence, we expect that well-specified joint models can be used to estimate causal effects for time-varying treatments.

7. Discussion

In this article, we have showcased how causal effects for time-varying treatments (or exposures) can be estimated using the framework of joint models for longitudinal and time-to-event data. These models will account for time-varying confounding without requiring an explicit specification of a model for the probability of receiving treatment conditional on the history of longitudinal confounders and past treatments. The causal effects (1)–(3) are in the flavor of the parametric G-formula and, by conditioning on different specifications of the longitudinal PSA history correspond to different targets of inference. In our model, we specified a linear mixed-effects model to link the longitudinal PSA measurements and the risk of metastasis. An alternative approach is to postulate a mixed model with nonlinear ordinary differential equations and link the PSA kinetics to the risk of metastasis. ¹⁵

Our approach relies on parametric assumptions and is expected to use the available data efficiently. However, the estimated causal effects will be biased if these assumptions are seriously violated. To minimize the chance of biased estimated effects, performing a thorough check of the model’s assumption using residual plots and evaluating the model’s fit with figures such as Supplemental Figures 3 to 5 is advisable. As we have seen in Section 3.2, an advantage of our full likelihood approach is that we do not need to model the salvage therapy initiation process. This property of our approach also holds for all other processes that may depend on the observed PSA history. For example, if the treating urologists decide to change the visiting process (i.e. when patients should come back for their next PSA test), and if this decision is solely based on the past observed PSA history of the patient (and possibly covariates), then our modeling approach will still provide valid results without requiring modeling this process. The same consideration also holds for the censoring process, that is, our approach allows censoring to depend in a complex manner on past longitudinal measurements and does not require a model to derive censoring weights. In particular, we feel that there are two general routes to follow for deriving causal effects from observational data; either to make no assumptions for the measurements process but require to derive weighting models for the other competing processes or to make stronger assumptions for the measurement model, and no assumptions for the other competing mechanisms. We have selected here the latter approach.

In this article, we have focused on the effect of salvage therapy at fixed time $t$ for the individual patient or for a group of patients with similar PSA histories. In the former case, we envisage using our approach in the context of shared decision-making for an individual patient. Namely, the doctor and the patient can assess the current risk of metastasis and how much this can be lowered by starting salvage therapy. This should weigh against the potential side effects of salvage therapy. In the context of a group of patients with similar PSA histories, our approach can be used to estimate the causal effect of different policies. For example, one policy might start salvage therapy when PSA first goes above 1.0 ng/mL. In contrast, another policy could delay the start of salvage therapy until PSA first goes above 4.0 ng/mL. A micro-simulation approach could then be used in which the parameter estimates from the joint model allow us to simulate realistic data under both these scenarios. Then, the causal effect would be defined as the difference in the incidence rate of metastasis.

As mentioned in Section 5, we have elected only to consider the first time patients received salvage therapy. However, we should stress that this is not a limitation of our proposed modeling framework. Namely, the model could be adjusted to specify the longitudinal subject-specific PSA profiles under multiple salvage therapy interventions. However, we have not done this in our analysis because few patients had received salvage therapy more than once. Hence, there was insufficient information in the data to estimate the changes in the PSA profiles after these interventions. In addition, our model could also be extended to account for deaths after metastasis by specifying a multi-state process model with a third hazard equation for the transition from metastasis to death.

Supplemental Material

sj-pdf-1-smm-10.1177_09622802241239003 - Supplemental material for Using joint models for longitudinal and time-to-event data to investigate the causal effect of salvage therapy after prostatectomy

Supplemental material, sj-pdf-1-smm-10.1177_09622802241239003 for Using joint models for longitudinal and time-to-event data to investigate the causal effect of salvage therapy after prostatectomy by Dimitris Rizopoulos, Jeremy MG Taylor, Grigorios Papageorgiou and Todd M Morgan in Statistical Methods in Medical Research