A technique for approximating transition rates from published survival analyses

Background Quality-adjusted-life-years (QALYs) are used to concurrently quantify morbidity and mortality within a single parameter. For this reason, QALYs can facilitate the discussion of risks and benefits during patient counseling regarding treatment options. QALYs are often calculated using partitioned-survival modelling. Alternatively, QALYs can be calculated using more flexible and informative state-transition models populated with transition rates estimated using multistate modelling (MSM) techniques. Unfortunately the latter approach is considered not possible when only progression-free survival (PFS) and overall survival (OS) analyses are reported. Methods We have developed a method that can be used to estimate approximate transition rates from published PFS and OS analyses (we will refer to transition rates estimated using full multistate methods as true transition rates). Results The approximation method is more accurate for estimating the transition rates out of health than the transition rate out of illness. The method tends to under-estimate true transition rates as censoring increases. Conclusions In this article we present the basis for and use of the transition rate approximation method. We then apply the method to a case study and evaluate the method in a simulation study. Electronic supplementary material The online version of this article (10.1186/s12962-019-0182-7) contains supplementary material, which is available to authorized users.

The transition rates between states can be organized in a transition rate matrix R where rows correspond to the "from" state, and the columns to the "to" state. States in the illness-death model are healthy (h), ill (i), and dead (d). Disallowed transitions are given a value of 0, and all rows must sum to 0. We model the MESCC disease process using time-homogenous transition rates: that is transition rates do not change over time and thus R is constant over time. Disallowed transitions are i to h, and d to either h or i. The transition rate matrix is   − (λ hi + λ hd ) λ hi λ hd 0 −λ id λ id 0 0 0   . [3] The transition rate matrix can be converted to a transition probability matrix P. This matrix records the probability of making transitions in the time interval from t = 0 to the time t. Using the transition rate matrix R, the transition probability matrix is computed as The fraction of patients in each state at time t is recorded in the state membership vector m (t). The initial state membership fractions are m (0) = [θ h , θ i , θ d ]. m (t) is computed as the product of m (0) and the transition probability matrix P (t), APPENDIX A. TRANSITION PARAMETERS FOR AN ILLNESS-DEATH MODEL 2 The first element of m (t) is H (t), the fraction of healthy patients The second element of m (t) is I (t), the fraction of ill patients The third element of m (t) is D (t), the fraction of dead patients

Prior Distributions
Non-informative prior distributions were used for multistate model estimation log λ ∼ Normal (mean = 0, variance = 10) , β ∼ Normal (mean = 0, variance = 10) . [9] An upper truncation limit of 10 was applied to all prior distributions. [9] Implementation The model was estimated using the Bayesian modelling language Stan [12] run through the statistical programming language R. [10] Four Markov chains were implemented with different random initial values for each chain.
Initially all chains were run for 5000 iterations without thinning; 2500 burn-in iterations were discarded. Post-burn-in traceplots for each parameter were examined to ensure that all four chains reached a similar mean and sampled similar regions of the distribution. [6] Iterations were increased until (1) the potential scale reduction factor,R, was < 1.1; and (2) histograms of the posterior distribution were smooth and without gaps. [2,1,4] To mitigate autocorrelation, thinning was increased so that effective sample size, N ef f ,was similar for all parameters, and the effective sample size as ≥ 20% of the number of sampling iterations. [6] Post-burn in iterations were increased until Monte Carlo (MC) error for each parameter was < 5% of the sample standard deviation. [11] The final simulation conditions were termed the short conditions.
To ensure that global convergence was achieved (rather than local convergence), models was re-run with iterations doubled. [4] Simulation conditions were adjusted to ensure convergence and mitigate autocorrelation. These simulation conditions were termed the long conditions. Percent relative deviation for the estimates from the long conditions to the short conditions was computed. Global convergence was indicated if the percent relative deviation was less than ±2% for each parameter. [4]

Results
To meet pre-specified benchmarks for autocorrelation, thinning was increased to 10 to ensure that effective sample size as ≥ 20% of the number of sampling iterations (Table B.2).
Traceplots for λ, and β showed a stable mean and variance under both short (Figure B.1)and long (Figure B.2) conditions.R was < 1.1 for all parameters under all conditions; and percent relative deviation was < ±2% for all parameters (Table B.2). We conclude that global convergence was reached.
Histograms of the posterior distribution under both short (Figure B.3) and long (Figure B.4) simulation conditions show smooth histograms without empty bins. Under both short and long simulation conditions, Monte Carlo error was < 5% of the sample standard deviation for each parameter (Table B.2). We conclude that both the short and long simulation conditions were sufficient to characterize the shape of the posterior distribution.