A multistate transition model for survival estimation in randomized trials with treatment switching and a cured subgroup | BMC Medical Research Methodology

Multistate transition model

We considered four states transitions, including cured, uncured, PD and death, in cancer clinical trials as shown in Fig. 1. Dotted lines represent the patient transitioning to a specific state with a certain probability (indicated by the characters on the dotted line). Initially, the study population consists of two heterogeneous subgroups as uncured and latent cured states with the cure probability of π. For a specific patient, his cure status can only be one of the two types, which will not change thereafter nor exist in an intermediate state. Patients who are cured will be no longer at risk of PD and subsequent treatment switching. They face a death hazard (shown as ({h}_{14})) similar to the general healthy population. For the uncured subgroup, they are exposed to a mortality risk as ({h}_{24}), or progress to PD (({h}_{23})) and then to death (({h}_{34})). The patients randomized to the control group can switch to the experimental group with a probability of p, and those who experience later PD are assumed with a lower probability of switching [23]. When estimating the long-term treatment effect (i.e., OS), we adopted hypothetical strategy within estimand framework and proposed a novel multistate transition model which is a continuous-time discrete-state Markov process in essence [24].

Specifically, let (c=1) or 0 representing the latent cured patient or uncured patient, respectively. The probability of being cured (denoted by π) is expressed as

where (g=1) or 0 denotes the experimental treatment or the control treatment.

For the cured subgroup, their baseline hazard aligns with the general healthy population, while situation of the uncured subgroup is more complicated. Suppose ({M}_{PD},{M}_{S},{M}_{D}) to be the state of PD, treatment switching and death when these indicators = 1, and ({t}_{PD},{t}_{D},{t}_{PPD}) be the time from randomization to PD, to death, and the time from PD to death, respectively. For patients who have experienced PD, there is ({t}_{D}={t}_{PD}+{t}_{PPD}). We convert ({t}_{PD}) to an ordinal categorical variable ({Q}_{{t}_{PD}}) according to the four quantiles of ({t}_{PD}) for convenience. A logistic model can be employed to establish the probability of treatment switching (denoted by p) as

Let (gamma) denote the patient-specific shared frail and follow a gamma distribution with mean one and variance θ, i.e., (Galeft({theta }^{-1},{theta }^{-1}right)) [25], which reflects the correlation between PD and the death hazards for each patient. That is, patients with higher PD risk tend to face greater death hazard, and those who experience later PD would have longer survival time. Based on similar studies [14, 21], we elucidate the transition hazards between states 1–4 as follows.

Therefore, the first equation ({h}_{14}) expresses just the baseline hazard similar to the general healthy population. The second and third equations describe the risk of PD (({h}_{23})) and death for the uncured subgroup (({h}_{24})), respectively. While the last equation ({h}_{34}) indicates the hazard of death after PD for the uncured subgroup with the coefficient β_34,1 and β_34,2 for the treatment effects of non-switchers and switchers, respectively. The main parameter estimates and their interpretation are listed as Table 1.

Maximum likelihood estimation (MLE) of the treatment effect

Suppose a trial with N patients, the possible observed status scenarios for each patient can be listed as Table 2. Patients with observed PD are undeniably uncured as in statuses A-D ((c=0)), while patients without observed PD could belong to either the latent cured or the uncured subgroup. Since patients in the control group are allowed to switch to the experimental group only after the occurrence of PD, the switching indicator is not applicable in statuses E and F.

Based on the observed data, the MLE method can be employed to estimate the parameters. Firstly, the likelihood functions ({f}_{1}sim {f}_{6}) are the marginal distributions with respect to the distribution of shared frailty ((gamma)), i.e., ({f}_{k}left(tright)=int {f}_{k}left(t|gamma right)dgamma,k=1,…,6). Moreover, the likelihood functions ({f}_{5}sim {f}_{6}) should be marginalized with respect to the distribution of c. Let ({P}_{(c)}) and ({P}_{({M}_{S})}) denote the certain state probability of cured and treatment switching, and (Ileft(gright)) be the indicator function for the treatment group or control group, where the values and meanings are consistent with the definitions above. For example,({P}_{left(c=0right)})​ represents the uncured probability, ({{P}_{({M}_{S}=1)}}^{I(g=0)}) represents the probability of the uncured control group patients switching to the treatment group. Consequently, the corresponding conditional distributions can be calculated as follows,

where ({S}_{24}left(tright)=mathit{exp}left[-{int }_{0}^{t}{h}_{24}left(u|gamma right)duright]), ({S}_{23}left(tright)=mathit{exp}left[-{int }_{0}^{t}{h}_{23}left(u|gamma right)duright]),({S}_{34}left(tright)=mathit{exp}left[-{int }_{0}^{t}{h}_{34}left(u|gamma right)duright]), and ({S}_{14}left(tright)=mathit{exp}left[-{int }_{0}^{t}{h}_{14}left(u|gamma right)duright]).

Assume the transitional baseline hazards (lambda) are constant, with marginal distributions of the conditional distributions above, the observed data likelihood for (Phi =left({a}_{0},{a}_{1},{b}_{0},{b}_{1},{lambda }_{14},{lambda }_{23},{beta }_{23},{lambda }_{24},{beta }_{24},{lambda }_{34},{beta }_{text{34,1}},{beta }_{text{34,2}},theta right)) is

To maximize (Lleft(Phi right)) in Eq. (4), let

Then, the goal is to find (Phi) that maximize (Z(Phi )) in Eq. (5).

Particle swarm optimization (PSO) algorithm

Although Laplace transformations-based approaches has be utilized by some researchers [26], we employed a simple and easily implemented PSO method in this study. The PSO algorithm is a population-based stochastic optimization technique developed by Kennedy and Eberhart [27] with the operational principle shown as Fig. 2. The globally optimal likelihood value is obtained via stochastic searching and optimization of candidate parameter estimate (denoted by particles). Once maximum number of iterations is reached or the change of the estimated likelihood falls below a certain threshold, the algorithm converges. Therefore, it is easily applied without complex formula derivations and widely applicable for optimization problems without the requirement of the concavity property of objective functions [28]. Given the simple iterations with low computational complexity, the PSO algorithm is theoretically fast -converging to meet the practical requirement. The step-by-step descriptions of PSO and more details of the initial particle setup are provided in Appendix 1.

Flow chart of the PSO algorithm [29]