



ORIGINAL ARTICLE 

Year : 2018  Volume
: 1
 Issue : 1  Page : 5157 

A Bayesian approach for dynamic treatment regimes in the presence of competing risk analysis
Atanu Bhattacharjee^{1}, Gajendra K Vishwakarma^{2}, Souvik Banerjee^{2}
^{1} Centre for Cancer Epidemiology, Tata Memorial Centre, Navi Mumbai, Maharashtra, India ^{2} Department of Applied Mathematics, Indian Institute of Technology (ISM), Dhanbad, Jharkhand, India
Date of Web Publication  12Dec2018 
Correspondence Address: Dr. Atanu Bhattacharjee Centre for Cancer Epidemiology, Tata Memorial Centre, Navi Mumbai, Maharashtra India
Source of Support: None, Conflict of Interest: None
DOI: 10.4103/CRST.CRST_6_18
Background: A sequencing rule is considered to formulate the dynamic treatment regime (DTR). This sequence rule is based on the clinical relevance, prior evidence about the best performing therapy, and the requirement to treat a patient in a specific scenario. The challenge occurs when the study offers a concluding remark about best effective therapy among all possible combinations of treatment management schedules treated with a sequence rule. The timetoevent data analysis is the only available method to figure out the best effective treatment in the context of oncology research. However, the presence of the competing risk event of death in the timetoevent analysis is unavoidable, and it becomes challenging to decide regarding the best effective treatment strategy. Methodology: In this article, we describe the statistical methodology to handle the competing risk timetoevent data analysis in DTR. The analysis is performed with the Bayesian approach to help determine the best effective treatment strategy. Results: We introduce the OpenBUGS function, which provides the comparison and estimation of different treatment sequences in timetoevent competing risk data analysis adopting the newly proposed statistical approach. Conclusion: This method is efficient to guide the personalized medicine in oncology setup through the supportive decision rule.
Keywords: Bayesian, competing risk, dynamic treatment regime, personalized medicine, sequential rule
How to cite this article: Bhattacharjee A, Vishwakarma GK, Banerjee S. A Bayesian approach for dynamic treatment regimes in the presence of competing risk analysis. Cancer Res Stat Treat 2018;1:517 
How to cite this URL: Bhattacharjee A, Vishwakarma GK, Banerjee S. A Bayesian approach for dynamic treatment regimes in the presence of competing risk analysis. Cancer Res Stat Treat [serial online] 2018 [cited 2019 Jun 19];1:517. Available from: http://www.crstonline.com/text.asp?2018/1/1/51/247331 
Introduction   
The incidence rates of various cancers are increasing,^{[1]} and therapeutic progress is still limited, with survival improvement ranging from only a few weeks to few months.^{[2]} Occasionally, the priority shifts toward cost and quality of life^{[3],[4]} and away from primary survival benefit, when there is no possible way to prolong the duration of survival. Even after several decades of cancer research, the treatment strategy is still decided based on the medical history, clinical status, disease symptoms, pathology, and prognostic status as determined by the pathologist and decided by the clinician. At present, the diagnostic procedures are well developed. The epigenetic, biochemical, genetic, imaging, metabolomic, and proteomic data are combined to help decide the best possible treatment.^{[5]} This is called personalized medicine.^{[6]} Personalized medicine provides us an opportunity to select the best effective treatment for a particular patient by considering his/her different features that may impact the treatment outcome. The most commonly considered contributing factors are the proven genetic biomarkers. The expression of these genetic biomarkers is often dynamic; they sometimes lead to enhanced responsiveness to chemotherapy and sometimes to chemoresistance. The ideal way to manage this in clinical practice is to continue the same regimen or to change to an alternative regimen, based on the changes in the sensitive genetic biomarker values.
This frequent change of decision about the treatment regimen based on the values of drugsensitive biomarkers generates a new treatment management protocol known as dynamic treatment regimes (DTRs) module.^{[7],[8],[9]} The closely observed timevarying covariates and intermediate response are a set of rules to help guide the treatment strategy. The DTR is beneficial in terms of costeffectiveness and patient compliance by avoiding both undertreatment and overtreatment. The suitable design to adopt the DTR approach is known as a sequential multipleassignment randomized trial (SMART).^{[10],[11],[12]} The literature about SMART and connecting approaches with DTR is well described.^{[12],[13]} In this manuscript, we are not going to further describe this SMART design. In our opinion, DTR approaches compatible only with SMART design are not openended enough to prescribe the best possible treatment option for a patient since the SMART design restricts a clinician only the study drug. The aim of this current manuscript is to overlook the SMART design for acceptability of DTR in general clinical practice. The schematic diagram of a different method of allocation of treatment is provided in [Figure 1]. The detection of the mean response of a DTR is developed by marginal modeling and optimal response modeling.^{[14]} The patient selection criteria for DTR are explored.^{[15]} The Qlearning is another development to perform DTR.^{[16],[17],[18]} The required optimum time for shifting of treatment in secondline therapy is also presented.^{[19]} The DTR estimation is determined by penalized least squares method.^{[13]} The DTR in the context of survival analysis particularly in censored data worked by the adjusted Qlearning algorithm.^{[20]} The development of DTR in survival context is very limited. However, inverse weighting method and the marginal mean model for survival analysis are documented.^{[21],[22],[23]} In nonparametric setup, the extension over Nelson–Aalen estimator is proposed by weighted risk set estimator.^{[24]}  Figure 1: A schematic diagram of different way allocation of A_{1}, B_{1}, A_{2}, and B_{2} treatments, and thereafter, patients final status defined as censored, died, or competing risk
Click here to view 
In DTR, an individual is expected to get differing combinations of treatment as per requirement during the course of therapy. The decision about the best effective treatment in DTR scheme is required to update the clinical practice evidence for the future patient care. In oncology, it becomes more interesting where the final or permanent outcome is death. This irreversible outcome is contributed to by differently scheduled combined therapy. However, it becomes challenging to take a decision about the best effective treatment in DTR setup. The concept of DTR is relatively new, and there is very limited literature with DTR.^{[2],[3],[9],[25],[26]} There has been a recent attempt toward the development of DTR in the context of survival analysis.^{[27],[28]}
The cause of death in a cancer patient is not necessarily always disease related. Death could be due to other reasons; these are defined as competing risk events.^{[29],[30],[31],[32]} Evaluation of the competing risk becomes, especially important when the survival benefit observed in the treatment group is limited. Death comparison through competing risk approach enhances the statistical inference and helps to decide the actual treatment effect. In the context of timetoevent analysis, apart from death, the relapsefree survival can also serve as the event of interest. The relapsefree survival is useful to capture the intermediate effect in DTR. The basic documentation on competing risk analysis is well described.^{[33],[34]} The data monitoring procedures for competing risk^{[35]} and competing risk with system failures have been recently attempted.^{[36]} Competing risk with missing data technique has also been explored.^{[32]} There are details about prior selection in competing risk analysis.^{[37]} The competing risk with interval censored through the Weibull parametric distribution has been recently described.^{[38]} There are several procedures available to handle competing risk procedure.^{[39],[40],[41]} These methods are compatible with opensource software, like R packages. However, there is no any available method for the development of competing risk in DTR. In this article, we have developed and described the statistical methodology to handle the competing risk timetoevent data analysis in DTR.
Methodology   
Kaplan–Meier estimator in dynamic treatment regime
The sample size of a study is denoted as n. The survival function is defined as S(T_{ijk}). It is defined as the duration of survival of the i^{th} individual at the j^{th} treatment effect in order k. It is formulated as:
Any patient may opt for a maximum of six lines of therapy (j) in his/her course of treatment, and corresponding order (k) of treatment is maximized with 6.
The estimates of S (T_{i}) are generated as follows:
Further, the variance is estimated as:
Competing risk in dynamic treatment regime
The covariate is defined as x. The cause of death is defined as m. If m = 1, then we consider that the death is due to cancer, if m = 2, then it is due to competing risk. If patient is alive, then it is defined as m = 0.
The hazard function is defined as:
The conditional probability that a subject with covariates x, dies in the interval (t_{jk}, t_{jk} + dt_{jk}), is calculated with m^{th} cause of death, given that the subject is alive before the time t_{jk}. The conversion of probability into a rate is done by dividing by dt and then take the limit as dt → 0.
Competing risk density in dynamic treatment regime
The density function of death or due to competing risk at time t is defined as:
It represents the unconditional risk of death of a patient at time t of cause m. The overall density of death can be represented with:
Competing risk in dynamic treatment regime with regression
The likelihood is defined as:
The probability of dead or alive will be taken care by S(t_{i}, x_{i}).
Now,
The survival function is presented as:
The time of observations are ordered as 0 <t_{1} <t_{2} <....< t_{N}. This order represents the distinct point of time T. Now, the risk set is defined as R (t_{p}) for the time point t_{p}. The number of individuals at risk at the time point t_{p} is defined as n_{p}. The number of occurrence of event is d_{p}. The time spent by the individual is assumed to be discrete. The probability of occurrence of event at t_{p} is:
Now, the term λ(t_{p}) can be represented as
The probability of surviving at time t is:
It shows that As n_{ip} increases, the number of event time increases as well and Kaplan–Meier starts to become continuous. The covariates of interest (Z) are specified as:
The β stands for the regression coefficient for the covariates of interest and λ_{0} is the baseline hazard. The β is estimated as:
Now, the cause of failure D is computed by:
Suppose the time of failure due to cause m is T_{m}. The cause specific function can be defined as:
The estimated value is defined as:
The proportion of subjects at risk that fail from cause m is:
The baseline hazard function is defined as λ_{m,0}(t) for cause of death m. The covariate of interest is Z and regression coefficient for m cause is β_{m}. Now, we consider m = 2 and this is defined as:
Results   
There is no consensus on the most appropriate procedure to perform DTR in a clinical research setting. We are approaching a future when all oncology patients will be offered some form of personalized medicine. However, due to nonavailability of the appropriate statistical tool, the suitable data are not available for reanalysis. We explored different possible sources but failed to get DTR data for reanalysis. The recently available R package dedicated to DTR also generated data (Tang and Melguizo). To illustrate the process, we perform a simulation study to assess the DTR in the oncology setup. Datasets were generated to resemble the skewed distributions seen in a headandneck cancer trial example. A total of six lines of therapy are generated from binomial distribution with a sample size of 300, and these are named as A_{1}, A_{2}, A_{3}, A_{4}, A_{5}, and A_{6} treatment arms, respectively. The random variables A_{1} to A_{6} are regenerated separately. The intention to generate treatment arms (A_{1} to A_{6}) is to make use of this methodology in general clinical practice rather than restricting it to a randomized controlled trial setting alone. The data are generated considering that any of the treatment regimens from A_{1} to A_{6} could be preferred by the clinician at any point of time for an individual patient. For example, there are 66 patients in whom is considered as first choice over. If any of these patients received this treatment, then it is coded as 1; if not, then it is coded as 0. However, the conditions are maintained as . A total of 160 patients are assigned to A_{1} as the first choice of therapy. The number of patients assigned A_{1}, A_{2}, A_{3}, A_{4}, A_{5}, and A_{6} is 160, 158, 159, 162, 155, and 162, respectively. The duration of survival is generated from normal distribution with a median duration of 24 months. The corresponding death status is obtained from rbinom (300, 2, 0.3). As a result, 161 individuals are observed with censored observations, 117 individuals who have died, and 22 with competing risk. [Figure 2]a and [Figure 2]b shows the Kaplan–Meier estimates on deaths and competing risk separately for different treatment groups. The preliminary competing risk analysis is performed with “CumIncidence” function.^{[40]} The median followup duration for different treatment groups with censored, death, and competing status is detailed in [Table 1]. The corresponding estimates obtained are detailed in [Table 2]. The treatmentwise competing risk for death and death due cancer is coded as 2 and 1, respectively. The survival function defined in equations (1) to (4) is used to handle the DTR problem. In this survival function, the survival duration is multiplied with a corresponding weight of each treatment (A_{i}, i = 1 to 6). It is assumed that each patient spent each unit of his/her unit of survival/followup time within a particular arm, i.e., A_{1} to A6. Now, the treatment that the patients spent the maximum time receiving is the maximum weight, followed by other treatments according to their time spend. If any of the treatments, A_{1} to A6 is not assigned to patients, then the corresponding weight would be zero for survival function. The unique feature of DTR is that one patient is expected to have been exposed to different treatments sequentially, and all the treatments can be allocated randomly in any sequence. However, the treatment that is assigned last and is followed by death should be penalized. There is a limitation of the existing statistical methodology. However, this proposed survival function can remedy that. This random sequence allocation is considered in the context of random variable data generation purpose only. In a reallife setting, it is expected that treatment allocation in any sequence will be based on the clinical evidence and the experience of a clinician, and it will not be based on blind randomization sequentially. This proposed method is generalized in this context and not bound by a blind randomized sequence. Once, we consider the modified survival duration as proposed and plots with [Figure 2]a and [Figure 2]b separately, and then, we consider this distribution for competing risk. The competing risk analyses are performed at two levels [Figure 3]. Initially, the classical competing risk is attempted and results presented in [Table 1] and [Table 2]. In the next step, the Bayesian counterpart is also extended to provide scope data exploration with covariates in competing risk. The likelihood defined in equation 10 is considered to perform competing risk with covariates. In this analysis, the covariate x_{i} is defined as treatment. This analysis is performed through opensource software OpenBUGS (www.openbugs.net). A total of 20,000 iterations with Markov chain Monte Carlo (MCMC) procedure are performed with 100 refresh in time point. This inference is based on the past 5000 measured iterations. The generated posterior mean and standard deviation with 95% credible interval is presented in [Table 3]. There are six arms and their regression coefficients are defined as β_{i1}, i = 1−6 toward magnitude of death due to cancer. The similar regression coefficients for competing risk are presented with β_{i1}, i = 1to6. The MCMC output provides the posterior mean estimate for each treatment with β_{i2}, i = 1to6 he significant level of contribution toward hazard for any treatment can be detected through 95% credible interval for each category.  Figure 2: (a) Treatment effect comparison with the duration of survival between different way allocation of A_{1}, A_{2}, A_{3}, A_{4}, A_{5}, and A_{6} treatments, and thereafter, patients' final status defined as censored or died. (b) Treatment effect comparison with the duration of survival between different way allocation of A_{1}, A_{2}, A_{3}, A_{4}, A_{5}, and A_{6} treatments, and thereafter, patients' final status defined as censored or competing risk
Click here to view 
 Table 1: Descriptive representation of median duration of followup among different treatment groups
Click here to view 
 Table 2: Statistical inference on cumulative incidence function estimated from competing risks data
Click here to view 
 Figure 3: Cumulative incidence plots on competing risks and death for each treatment groups
Click here to view 
 Table 3: Posterior estimates with treatment arms through competing risk analysis
Click here to view 
Discussion   
In this paper, we have developed a Bayesian approach to the analysis of competing risk data. Our Bayesian justification is the first such development in the DTR context with competing risk data. The outcome of the method presented combines the classical available method with competing risk. The occurrence of death and the cause of death are both equally important in the setup of competing risk. The causes of death may be independent or dependent in nature in competing risk assumption. The influence of covariates cannot be ignored in competing risk analysis. The assumption test criteria about failure time distribution are well documented.^{[33],[42]} There is also existing methodology toward shape and scale testing properties of survival duration.^{[43]} However, those approaches are limited with competing risk analysis and not explored in the presence of DTR. However, this work is also not attempted toward assumption testing in competing risk. It is applicable only in a nonparametric setup. This is the limitation of our study, but the assumption testing work can be extended in the future. This work is illustrated with basic application of competing risk modeling. If parametric approach is needed to apply in this setup, then it requires different sensitive assumption testing. We would prefer Fine and Gray method for causespecific hazard function modeling.^{[44]} As a counterpart, the Bayesian method is found to be equally effective.^{[32]} The SAS macro is available for Fine and Gray testing.^{[45]} It is also feasible with R through mstate and timereg packages. There is no available Bayesian code to perform competing risk. The OpenBUGS code is available on request by email to the corresponding author. Moreover, this proposed method can also be used for an exploratory study with sequential treatment allocation. However, more care is necessary in providing the next stage treatment arm allocation probabilities.
Conclusion   
There is a need for competing risks approaches in oncology research. Simultaneously, there is progress in oncology therapy toward the establishment of personalized medicine. The DTR is the way to conduct inference about personalized medicine in the future. However, there is a gap between the availability and requirement of statistical methods to cover DTR approaches. The challenge with competing risk analysis is one of them. In this study, Bayesian analysis for competingrisk models is derived for DTR data analysis. This method provides flexibility in modeling, without creating the stratum in the dataset for competing risk and death separately. There is also the provision to consider covariates in this modeling. Furthermore, it gives an assessment of the influencing covariates with hazard rate in DTR data analysis. Our method is feasible based on the outcome of this work.
Financial support and sponsorship
Nil.
Conflicts of interest
There are no conflicts of interest.
References   
1.  Thakkar JP, McCarthy BJ, Villano JL. Agespecific cancer incidence rates increase through the oldest age groups. Am J Med Sci 2014;348:6570. 
2.  Thippeswamy R, Noronha V, Krishna V, Joshi A, Bal MM, Purandare N, et al. Stage IV lung cancer: Is cure possible? Indian J Med Paediatr Oncol 2013;34:1215. [ PUBMED] [Full text] 
3.  Joshi A, Patil VM, Noronha V, Ghosh J, Bhattacharjee A, Prabhash K, et al. A crosssectional observation study regarding patients and their physician willingness to wait for driver mutation report in nonsmallcell lung cancer. Indian J Med Paediatr Oncol 2016;37:748. [ PUBMED] [Full text] 
4.  Noronha V, Joshi A, Patil VM, Agarwal J, GhoshLaskar S, Budrukkar A, et al. Onceaweek versus onceevery3weeks cisplatin chemoradiation for locally advanced head and neck cancer: A phase III randomized noninferiority trial. J Clin Oncol 2018;36:106472. 
5.  Verma M. Personalized medicine and cancer. J Pers Med 2012;2:14. 
6.  van de Kerkhof PC. Stratified or personalised medicine in the treatment of psoriasis? J Dermatol Treat 2017;28:683. 
7.  Shortreed SM, Moodie EE. Estimating the optimal dynamic antipsychotic treatment regime: Evidence from the sequential multiple assignment randomized CATIE schizophrenia study. J R Stat Soc Ser C Appl Stat 2012;61:57799. 
8.  Chakraborty B, Moodie EE. Statistical Methods for Dynamic Treatment Regimes. New York: Springer; 2013. p. 3152. 
9.  Lizotte DJ, Tahmasebi A. Prediction and tolerance intervals for dynamic treatment regimes. Stat Methods Med Res 2017;26:161129. 
10.  Lavori PW, Dawson R. Dynamic treatment regimes: Practical design considerations. Clin Trials 2004;1:920. 
11.  Murphy SA. A generalization error for Qlearning. J Mach Learn Res 2005;6:107397. 
12.  Murphy SA, Lynch KG, Oslin D, McKay JR, TenHave T. Developing adaptive treatment strategies in substance abuse research. Drug Alcohol Depend 2007;88:S2430. 
13.  Qian M, Murphy SA. Performance guarantees for individualized treatment rules. Ann Stat 2011;39:1180210. 
14.  Murphy SA. Optimal dynamic treatment regimes. J R Stat Soc 2003;65:33155. 
15.  Murphy SA, Bingham D. Screening experiments for developing dynamic treatment regimes. J Am Stat Assoc 2009;104:391408. 
16.  Watkins CJ. Learning from Delayed Rewards. Doctoral Dissertation. Cambridge: King's College; 1989. 
17.  Murphy SA, van der Laan MJ, Robins JM; Conduct Problems Prevention Research Group. Marginal mean models for dynamic regimes. J Am Stat Assoc 2001;96:141023. 
18.  Zhao YQ, Laber EB. Estimation of optimal dynamic treatment regimes. Clin Trials 2014;11:4007. 
19.  Schulte PJ, Tsiatis AA, Laber EB, Davidian M. Q and Alearning methods for estimating optimal dynamic treatment regimes. Stat Sci 2014;29:64061. 
20.  Goldberg Y, Kosorok MR. Qlearning with censored data. Ann Stat 2012;40:52960. 
21.  Robins JM, Rotnitzky A, Zhao LP. Estimation of regression coefficients when some regressors are not always observed. J Am Stat Assoc 1994;89:84666. 
22.  Murphy SA. An experimental design for the development of adaptive treatment strategies. Stat Med 2005;24:145581. 
23.  Lunceford JK, Davidian M, Tsiatis AA. Estimation of survival distributions of treatment policies in twostage randomization designs in clinical trials. Biometrics 2002;58:4857. 
24.  Guo X, Tsiatis A. A weighted risk set estimator for survival distributions in twostage randomization designs with censored survival data. Int J Biostat 2005;1:115. 
25.  Moodie EE, Richardson TS, Stephens DA. Demystifying optimal dynamic treatment regimes. Biometrics 2007;63:44755. 
26.  Zhao YQ, Zeng D, Laber EB, Kosorok MR. New statistical learning methods for estimating optimal dynamic treatment regimes. J Am Stat Assoc 2015;110:58398. 
27.  Tang X, Melguizo M. DTR: An R package for estimation and comparison of survival outcomes of dynamic treatment regimes. J Stat Softw 2015;65:128. 
28.  Wallace MP, Moodie EE, Stephens DA. Dynamic treatment regimen estimation via regressionbased techniques: Introducing package DTRreg. J Stat Softw 2017;80:120. 
29.  Satagopan JM, BenPorat L, Berwick M, Robson M, Kutler D, Auerbach AD, et al. A note on competing risks in survival data analysis. Br J Cancer 2004;91:122935. 
30.  Pintilie M. Competing Risks: A Practical Perspective. Vol. 58. The Atrium, Southern Gate, England: John Wiley and Sons; 2006. 
31.  Bhattacharyya T, Bhattacharjee A. Competing risk: An illustration with aspiration pneumonia in head and neck cancer patients undergoing radical radiotherapy: A biostatistician's perspective. Indian J Cancer 2014;51:4069. [ PUBMED] [Full text] 
32.  Bhattacharjee A. Bayesian competing risks model: An application to breast cancer clinical trial with incomplete observations. J Stat Manag Syst 2015;18:381404. 
33.  Crowder MJ. Classical Competing Risks. Boca Raton: Chapman and Hall/CRC; 2001. 
34.  BandeenRoche K, Liang KY. Modelling multivariate failure time associations in the presence of a competing risk. Biometrika 2002;89:299314. 
35.  Steiner SH, MacKay RJ. Monitoring processes with data censored owing to competing risks by using exponentially weighted moving average control charts. J R Stat Soc 2001;50:293302. 
36.  Sen A, Basu S, Banerjee M. Analysis of masked failure data under competing risks. Handbook of Statistics. Vol. 20. Ch. 19. Elsevier North Holland; 2001. p. 52340. 
37.  Lin DK, Usher JS, Guess FM. Bayes estimation of componentreliability from masked systemlife data. IEEE Trans Reliab 1996;45:2337. 
38.  Kundu D, Pradhan B. Bayesian analysis of progressively censored competing risks data. Sankhya B 2011;73:27696. 
39.  Putter H, Fiocco M, Geskus RB. Tutorial in biostatistics: Competing risks and multistate models. Stat Med 2007;26:2389430. 
40.  Scrucca L, Santucci A, Aversa F. Competing risk analysis using R: An easy guide for clinicians. Bone Marrow Transplant 2007;40:3817. 
41.  de Wreede LC, Fiocco M, Putter H. Mstate: An R package for the analysis of competing risks and multistate models. J Stat Softw 2011;38:130. 
42.  Kalbfleisch JD, Prentice RL. The Statistical Analysis of Failure Time Data. Vol. 360. John Wiley & Sons; 2011. 
43.  Rao BR, Talwalker S, Kundu D. Confidence intervals for the relative risk ratio parameter from survival data under a random censorship model in biomedical and epidemiologic studies (by simulation). Biom J 1991;33:95984. 
44.  Fine JP, Gray RJ. A proportional hazards model for the subdistribution of a competing risk. J Am Stat Assoc 1999;94:496509. 
45.  Rosthøj S, Andersen PK, Abildstrom SZ. SAS macros for estimation of the cumulative incidence functions based on a cox regression model for competing risks survival data. Comput Methods Programs Biomed 2004;74:6975. 
[Figure 1], [Figure 2], [Figure 3]
[Table 1], [Table 2], [Table 3]
