DTR : An R Package for Estimation and Comparison of Survival Outcomes of Dynamic Treatment Regimes

Sequentially randomized designs, more recently known as sequential multiple assignment randomized trial (SMART) designs, are widely used in biomedical research, particularly in clinical trials, to assess and compare the eﬀects of various treatment sequences. In such designs, patients are initially randomized to one of the ﬁrst-stage therapies. Then patients meeting some criteria (e.g., no relapse of disease) participate in the second-stage randomization to one of the second-stage therapies. The advantage of such a design is that it allows the investigator to study various treatment sequences where the patients’ second-stage therapies can be adjusted based on their responses to the ﬁrst-stage therapies. In the past few years, substantial improvement has been made in the statistical methods for analyzing the data from SMARTs. Much of the proposed statistical approaches focus on estimating and comparing the survival outcomes of treatment sequences embedded in the SMART designs. In this article, we introduce the R package DTR , which provides a set of functions that can be used to estimate and compare the eﬀects of diﬀerent treatment sequences on survival outcomes using the newly proposed statistical approaches. The proposed package is also illustrated using simulated data from SMARTs. randomized ﬁrst-stage therapies. some no relapse of disease) participate in the second-stage randomization to one of the second-stage therapies. The second-stage therapy could be a rescue therapy if patients show primary resistance to the ﬁrst-stage therapy, or a maintenance therapy if favorable responses are observed. There are diﬀerent types of SMART designs with respect to the participants in the second-stage randomization: (i) SMART designs in which only the responders to one of the ﬁrst-stage therapies participate in the second-stage randomization to one of the maintenance therapies; (ii) SMART designs in which only the non-responders to one of the ﬁrst-stage therapies participate in the second-stage randomization to one of the rescue therapies; (iii) SMART designs in which the responders to all the ﬁrst-stage therapies participate in the second-stage randomization; (iv) SMART designs in which the non-responders to all the ﬁrst-stage therapies participate in the second-stage randomization; and (v) SMART designs in which all the responders and non-responders participate in the second-stage randomization (Lei, Nahum-Shani, Oslin, and Murphy 2012; Chakraborty and Moodie advantage of a design is that the investigator study various treatment sequences where the patients’ second-stage therapies can be adjusted based on their responses to the ﬁrst-stage therapies. Here is an combination myeloablative


Introduction
Sequentially randomized designs, more recently known as sequential multiple assignment randomized trial (SMART) designs, are widely used in biomedical research, particularly in clinical trials, to assess and compare the effects of various treatment sequences (Robins 1986(Robins , 1987Lavori and Dawson 2000;Murphy 2005;Bembom and van der Laan 2007;Chakraborty and Murphy 2013;Kidwell 2014). In such designs, patients are initially randomized to one of the first-stage therapies. Then patients meeting some criteria (e.g., no relapse of disease) participate in the second-stage randomization to one of the second-stage therapies. The secondstage therapy could be a rescue therapy if patients show primary resistance to the first-stage therapy, or a maintenance therapy if favorable responses are observed. There are different types of SMART designs with respect to the participants in the second-stage randomization: (i) SMART designs in which only the responders to one of the first-stage therapies participate in the second-stage randomization to one of the maintenance therapies; (ii) SMART designs in which only the non-responders to one of the first-stage therapies participate in the second-stage randomization to one of the rescue therapies; (iii) SMART designs in which the responders to all the first-stage therapies participate in the second-stage randomization; (iv) SMART designs in which the non-responders to all the first-stage therapies participate in the second-stage randomization; and (v) SMART designs in which all the responders and nonresponders participate in the second-stage randomization (Lei, Nahum-Shani, Lynch, Oslin, and Murphy 2012;Chakraborty and Moodie 2013). The advantage of such a design is that it allows the investigator to study various treatment sequences where the patients' secondstage therapies can be adjusted based on their responses to the first-stage therapies. Here is an example of a SMART design investigating the effect of a combination of myeloablative chemotherapy, total-body irradiation and transplantation of autologous bone marrow purged of cancer cells (ABMT) to a standard chemotherapy in treating children with high risk neuroblastoma (Matthay et al. 1999(Matthay et al. , 2009. All the children with high risk neuroblastoma were treated with the same initial regimen of chemotherapy, and those without disease progression participated in the first-stage randomization to either ABMT or three cycles of intensive chemotherapy. All the children who completed cytotoxic therapy without disease progression then participated in the second-stage randomization to either treatment with 13-cis-retinoic acid (cis-RA) for six months or no further therapy.
In the past few years, substantial improvement has been made regarding the statistical methods for analyzing data from SMART designs in terms of either continuous outcomes Lavori 2010, 2012), binary outcomes (Buyze and Goetghebeur 2013) or survival outcomes (Lunceford, Davidian, and Tsiatis 2002;Guo and Tsiatis 2005). This article focuses on the proposed statistical approaches for estimating and comparing the survival endpoints of treatment sequences from these trials. Those treatment sequences are usually referred to as dynamic treatment regimes (DTRs), adaptive treatment strategies, or treatment policies. A DTR is defined as a sequence of the first-stage therapy, an intermediate response to the first-stage therapy, and the second-stage therapy. The construction of DTRs is usually from SMART designs or longitudinal observational studies. However, this article only discusses the DTRs constructed from SMART designs. In the high risk neuroblastoma study described above, four DTRs can be constructed: (i) treat with ABMT followed by cis-RA if no disease progression; (ii) treat with ABMT followed by no further therapy if no disease progression; (iii) treat with chemotherapy followed by cis-RA if no disease progression; and (iv) treat with chemotherapy followed by no further therapy if no disease progression. Evaluating the effect of a sequence of therapies is more efficient and more clinically meaningful than looking at the first-and second-stage therapies separately (Chakraborty and Murphy 2013;Kidwell 2014). Lunceford et al. (2002) introduced the survival and mean restricted survival estimators for treatment regimes in a two-stage randomization design using two forms of inverse-probability weighting for second-stage randomization and censoring respectively. Guo and Tsiatis (2005) proposed a weighted risk set estimator (WRSE) for the cumulative hazard function. Tang and Wahed (2011) proposed a generalized Cox model for comparing any combination of treatment regimes after adjustment for auxiliary variables. Kidwell and Wahed (2013) proposed weighted log-rank test statistics to compare survival distributions of DTRs. Tang and Wahed (2015) introduced the cumulative hazard ratio (CHR) estimator between two DTRs, and a testing procedure to compare the effects of DTRs based on the CHR estimator.
Although we have seen considerable advancement in the development of statistical methods for estimating and comparing the effects of DTRs on survival outcomes, no statistical package has been developed to implement those newly-proposed statistical approaches. Investigators and statisticians hesitate to implement the newly-proposed statistical approaches, mainly because (i) the statistical formulae for estimating and comparing the survival distributions of DTRs are theoretically complex, and (ii) the newly-proposed statistical methodology has not been included in any widely-used statistical software. Therefore, we developed the R package DTR (Tang and Melguizo 2015) to implement the newly-proposed statistical approaches for estimating and comparing the survival outcomes of different DTRs. The current version of the package and documentation implemented the SMART designs in which the responders to all the first-stage therapies participate in the second-stage randomization. However, it can also be applied to SMART designs where the non-responders to all the first-stage therapies participate in the second-stage randomization by switching the responders with non-responders. In Section 2 we briefly introduce the notation and the formulae for calculating the survival quantities based on each statistical method included in this version of the package. In Section 3 the functionality of each function is described and illustrated using simulated data from SMARTs. A more complete overview of the functionality is given in the reference manual of the DTR package which is available from the Comprehensive R Archive Network (CRAN) at http://CRAN.R-project.org/package=DTR. Some concluding remarks are provided in Section 4.

Notation
In the SMART designs used for the DTR package, patients are initially randomized to J firststage therapies (A 1 , A 2 , . . . , A J ), then patients who responded to the first-stage therapies are then randomized to K second-stage therapies (B 1 , B 2 , . . . , B K ). The treatment regime A j B k , where j = 1, 2, . . . , J; k = 1, 2, . . . , K, is defined as "treat with A j , followed by B k if responds to A j ." The primary research interest is to estimate and compare the effects of different DTRs in terms of survival quantities (e.g., survival function, cumulative hazard function). The set of observed data from this design can be described as for j = 1, 2, . . . , J; k = 1, 2, . . . , K; and i = 1, 2, . . . , n; where X ji : the indicator for the jth first-stage therapy, X ji = 1 if the ith patient was assigned to first-stage therapy A j , and X ji = 0 if otherwise; R i : the indicator for response, R i = 1 if the ith patient responded, and R i = 0 if otherwise; T R i (optional): the time to response for the ith patient; Z ki : the indicator for kth second-stage therapy, Z ki = 1 if the ith patient was assigned to second-stage therapy B k , and R i = 1 and Z ki = 0 if otherwise; (note that R i = 0 and Z ki = 0 if the ith patient did not respond, and thus did not participate in the second-stage randomization); U i : the observed event/censoring time for the ith patient; ∆ i : the censoring indicator, ∆ i = 0 if the ith patient was censored, and ∆ i = 1 if the ith patient's event was observed; and V i (optional): the covariate vector for the ith patient.
When the response status is considered as a time-varying process, we define the time-varying response variable to be R i (t) = R i I(T R i < t). Based on the analytical framework of inverse-probability weighting (Lunceford et al. 2002;Guo and Tsiatis 2005;Tsiatis 2004, 2006;Tang and Wahed 2015;Kidwell and Wahed 2013), the weight functions used by various statistical approaches for treatment regime A j B k , j = 1, 2, . . . , J and k = 1, 2, . . . , K are defined in Table 1. Because patients randomized to A 1 , A 2 , . . ., and A J are independent samples, some statistical methods (Lunceford et al. 2002;Guo and Tsiatis 2005) estimate the survival quantities for treatment regime A j B k based on A j data only. Based on the counting process notation described in Fleming and Harrington (1991), the risk, event and censoring indicators used in this article are outlined in Table 2. The survival function, cumulative hazard function, and hazard function for treatment regime A j B k are denoted as S jk (t), Λ jk (t) and λ jk (t) respectively for j = 1, 2, . . . , J and k = 1, 2, . . . , K.

Inverse-probability weighting estimator
The survival estimator for treatment regime A j B k using the two forms of inverse-probability weighting (Lunceford et al. 2002) iŝ j = 1, 2, . . . , J and k = 1, 2, . . . , K whereK(t) is the Kaplan-Meier estimate of the censoring curve, and its variance estimator is written as where L is defined as the restricted lifetime, which is smaller than the maximum follow-up time of the SMART, and

Definition Estimation based Time-varying
Usage on A j data? response?

Yes
No Inverse-probability weighting estimator

Yes
Yes Weighted risk set estimator Table 2) Table 2) Table 1: Various weight functions for treatment regime A j B k for j = 1, 2, . . . , J and k = 1, 2, . . . , K, where π j = P(X ji = 1) and whereĜ The detailed steps for calculating the variance estimator are outlined in Appendix A. The inverse-probability weighting estimator is mostly used when the primary research question of a SMART is to estimate the survival probabilities of DTRs over time, and the time to response is not observed or it is known that the time to response does not affect the overall survival outcome.

Weighted risk set estimator
The weighted risk set estimator (WRSE) of the survival function for treatment regime A j B k (Guo and Tsiatis 2005) iŝ and its variance estimator is written as The detailed steps for calculatingσ 2 are outlined in Appendix B. The weighted risk set estimator is mostly used when the primary research question of a SMART is to estimate the survival probabilities of DTRs over time, and it is suspected that the time to response would impact the overall survival outcome.

Cumulative hazard ratio estimator
Based on stratified proportional hazards with treatment regimes as strata (Tang and Wahed 2015), the hazard function for treatment regime A j B k , j = 1, . . . , J and k = 1, . . . , K can be written as where λ jk0 (t) is the baseline hazard function for treatment regime A j B k , and β is a vector of coefficients corresponding to baseline covariates V . The coefficient estimateβ can be obtained by solving the pseudo-score equation The cumulative baseline hazard for treatment regime A j B k can be obtained aŝ Based on the estimated cumulative baseline hazards, the estimated cumulative hazard ratio (CHR) for comparing treatment regimes A j B k to A j B k can be obtained aŝ The cumulative hazard ratio estimatorθ jkj k (t) follows a Gaussian process with mean θ jkj k (t) and variance function σ 2 jkj k (t). The detailed steps for computing the variance estimator σ 2 jkj k (t) forθ jkj k (t) are described in Section 3 and Appendix B of Tang and Wahed (2015). The cumulative hazard ratio estimator is mostly used when the primary research interest is to compare the survival curves of different DTRs after adjusting for other potential risk factors, and it is suspected that the proportional hazards assumption is violated.

Wald-type test
Wald-type tests can be used for comparing the survival quantities of different treatment regimes based on either the inverse-probability weighting estimator (Lunceford et al. 2002), WRSE (Guo and Tsiatis 2005) or the logarithm of the CHR estimator (Tang and Wahed 2015). For example, let us denote the estimates for the survival functions to bê and their estimated variance covariance matrix to beΣ. For testing the null hypothesis , and it follows a chi-square distribution with d degrees of freedom. The Wald-type test is mostly used when the primary research interest is to compare the survival rates of different DTRs at certain times (e.g., 3-year survival, and 5-year survival).
2.6. Weighted log-rank test Feng and Wahed (2008) introduced a supremum weighted log-rank test statistic for testing for the survival differences between two DTRs. However, the supremum weighted log-rank tests can only be applied for comparing two separate-path DTRs. Kidwell and Wahed (2013) extended Feng and Wahed (2008)'s work and proposed the weighted log-rank tests for comparing the survival distributions of shared-path DTRs. The current version of the package incorporated the weighted log-rank tests proposed by Kidwell and Wahed (2013). The standardized weighted log-rank test statistic for comparing treatment regime A j B k and A j B k can be written as where and The detailed steps for calculating Z W jkjk (t) andσ 2 (t) are outlined in Appendix C. The weighted log-rank test is mostly used when the primary research interest is to compare the survival curves of different DTRs, and there is no other risk factors to be considered.

Generalized Cox model
For comparing the effects of different DTRs after adjustment for auxiliary variables, Tang and Wahed (2011) proposed to use the following version of the Cox model based on two-stage randomization designs where J = K = 2: where λ(t) is the general hazard function at time t; λ 0 (t) denotes the baseline hazard function (when all the covariates are equal to 0); R(t) denotes the time-varying measurement of response and consent as defined before; and β is the vector of coefficients denoted as 1 , β 1 , β 11 , γ . Under model (9), comparisons of treatment regimes in terms of their hazards can be interpreted based on the coefficient vector β. For example, for comparing all four treatment regimes in a simple two-stage randomization design, the null hypothesis H 0 : The generalized Cox model is mostly used when the primary research interest is to compare the survival curves of different DTRs after adjusting for other potential risk factors under the assumption of proportional hazards across regimes.

The DTR package
The DTR package for the R environment for statistical computing and graphics (R Core Team 2015) can be downloaded from CRAN at http://CRAN.R-project.org/package=DTR. This package is designed for the estimation and comparison of survival distributions of DTRs from SMARTs in which the responders to all the first-stage therapies participate in the secondstage randomization. In SMART designs, there could be more than two therapies available at each stage. For simplicity, and to maintain the similarity to the most common SMARTs, a two-stage randomization design allowing for two treatment options at each stage is used in the current version of the package. The basic structure of a SMART design used for this package is depicted in Figure 1. In detail, patients are initially randomized to either A 1 or A 2 as first-stage therapy, then patients who responded to the first-stage therapies are randomized to either B 1 or B 2 as second-stage therapy. Four DTRs are embedded in this design: A 2 B 1 : treat with A 2 , followed by B 1 if responds to A 2 ; and A 2 B 2 : treat with A 2 , followed by B 2 if responds to A 2 .
The current version of the DTR package includes the functions for simulating data from various two-stage randomization designs, estimating the survival quantities (e.g., survival function, restricted mean survival, cumulative hazard ratio) of treatment regimes, plotting the survival estimates over time, and comparing the survival distributions of different regimes. The statistical approaches proposed in Lunceford et al. (2002); Guo and Tsiatis (2005); Tang and Wahed (2011Wahed ( , 2015; and Kidwell and Wahed (2013) are implemented in the current version of the DTR package. A brief introduction of each statistical method is given in Sections 2.2-2.7. For the convenience of the reader, we also summarize the functionality of the DTR package in Table 3.

Output/input data
The data from a two-stage randomization design shown in Figure 1 is stored in an R data frame consisting of the following arguments: delta: Censoring indicator, delta = 1 for observed events, and delta = 0 for censored.
V: Covariate(s) to be adjusted. V may include one column for one covariate, or more than one column for multiple covariates.
Because (i) only the data from the A j arm was described in the simulation studies of Lunceford et al. (2002) and Guo and Tsiatis (2005), (ii) the approaches described in Guo and Tsiatis (2005) and Kidwell and Wahed (2013) treated response status as a time-varying process, (iii) Tang and Wahed (2011) and Tang and Wahed (2015) took into account the information from auxiliary covariates, different sets of the above arguments are generated by different data simulation functions, and required as the input data for different estimation and comparison functions. For clarification, we summarize the output/input data for each function in Table 4.

Inverse-probability weighting estimator functions
The function LDTestimate() computes the inverse-probability weighting survival estimates and their estimated standard errors for DTRs at observed event times based on Equations 1 and 2. The technical details are briefly described in Section 2.2 (see more details in Lunceford et al. 2002). The examples in this section will be illustrated using a simulated dataset called LDTdata which is included in the DTR package. The LDTdata dataset was simulated as  described in the simulation study of Lunceford et al. (2002) using the simLDTdata() function as follows.
After installing package DTR, load the package:
R> dim(LDTdata)     There was a significant difference in the 1-year survival among four treatment regimes (overall p = 0.048), and patients following regime A 2 B 2 appeared to have a better 1-year survival compared to patients following regime A 2 B 1 (A 2 B 1 vs. A 2 B 2 : p = 0.04; Figure 2).

Weighted risk set estimator functions
The function WRSEestimate() calculates the weighted risk set estimates of the survival functions and their estimated standard errors based on Equations 3 and 4. The technical details are briefly described in Section 2.3 (see more details in Guo and Tsiatis 2005). The examples in this section will be illustrated using a simulated dataset called WRSEdata which is included in the DTR package. The WRSEdata dataset was simulated as described in the simulation study of Guo and Tsiatis (2005) using the simWRSEdata() function.
Load the WRSEdata dataset into the workspace: R> data("WRSEdata", package = "DTR") The WRSEdata dataset is a data frame with 200 rows corresponding to patients and 6 columns corresponding to variables. In this dataset, 200 eligible patients were equally assigned to first-stage therapy A 1 or A 2 . All the patients were followed from the time of randomization until they responded/did not respond to the first-stage therapy, and the time to response varied across patients. 98 (49%) of the patients responded to the first-stage therapy, and were further assigned to one of the second-stage therapies B 1 or B 2 at a 1 : 1 ratio. The maximum follow-up time was 3.5 years.
R> dim(WRSEdata) The function WRSEestimate() creates an object of the same class 'DTR' as described above.  A 1 B 1 , A 1 B 2 , A 2 B 1 , and A 2 B 2 respectively. The functions summary() and plot() can be used to display and plot the survival estimates at observed event times respectively. The plot of the survival estimates and their 95% confidence bands across study period is shown in Figure 3. Comparisons of survival estimates can be carried out between different DTRs by testing the equality of the weighted risk set estimates at a specific time point by calling the function contrast_wald().

Cumulative hazard ratio estimator functions
Based on the statistical methods proposed in Tang and Wahed (2015), we developed the CHRestimate() function for calculating the CHR estimates between pairwise DTRs across times. A brief introduction of the statistical method can be found in Section 2.4. The examples in this section will be illustrated using a simulated dataset called CHRdata which is included in the DTR package. The CHRdata dataset was simulated as described in the simulation study of Tang and Wahed (2015) using the simCHRdata() function.
Load the CHRdata dataset into the workspace: R> data("CHRdata", package = "DTR") The CHRdata dataset is a data frame with 200 rows corresponding to patients and 7 columns corresponding to variables. In this dataset, a total of 200 eligible patients participated in the first-stage randomization, and were equally assigned to either A 1 or A 2 as first-stage therapy. Among the 200 patients, 121 (60.5%) responded to the first-stage therapies they were assigned to, and then participated in the second-stage randomization. At the secondstage randomization 121 patients were equally assigned to either B 1 or B 2 as second-stage therapy. Approximately 25% of the patients were censored during the 5-year study period. The investigator suspected that gender (column V1, 1 for males and 0 for females) may have an impact on the survival outcomes, and thus would like to adjust for the gender effect when assessing the effects of the four DTRs in terms of survival.
R> est <-CHRestimate(data = CHRdata, covar = "V1") R> est$coefficients The gender coefficient was estimated to be 0.37. The hazards of death among males were estimated to be 1.45 times the hazards among females.

R> est
Call: CHRestimate(data = CHRdata, covar = "V1") Only the first few rows of the summary() output are shown above. The pairwise log CHR estimates between any two different DTRs and their 95% confidence bands can be plotted by calling the plot() function with log.CHR = TRUE and confidence.interval = TRUE (Figure 4). The comparisons of different treatment regimes are carried out by performing the Wald-type tests based on the log CHR estimates by calling the contrast_chr() function.  As shown in Figure 4, there was no significant difference in the 3-year survival between treatment regimes A 1 B 2 and A 1 B 1 (p = 0.75). Other than that, all other treatment regimes were significantly different from each other at 3 years (overall p < 0.001).

The weighted log-rank test function
The contrast_logrank() function was developed with regard to the statistical methods proposed in Kidwell and Wahed (2013). We briefly introduced the methods in Section 2.6. The examples in this section will be illustrated using a simulated dataset called LRdata which is included in the DTR package. The LRdata dataset was simulated as described in the simulation study of Kidwell and Wahed (2013) using the simLRdata() function.
Load the LRdata dataset into the workspace: R> data("LRdata", package = "DTR") The LRdata dataset is a data frame with 100 rows corresponding to patients and 6 columns corresponding to variables. In this dataset, 100 eligible patients were equally assigned to one of the first-stage therapies A 1 or A 2 . They were followed until they responded to or did not respond to the first-stage therapies. 61 (61%) responders were then equally assigned to one of the second-stage therapies B 1 or B 2 . The maximum follow-up time was 12 years, and 24% of the patients were lost to follow-up during the study.
R> dim(LRdata) The contrast_logrank() function compares the survival distributions of different DTRs using the weighted log-rank tests (Kidwell and Wahed 2013 We did not find any significant difference in the survival distributions of four DTRs (overall p = 0.52).

Generalized Cox model functions
The PHfit() and contrast_ph() functions were developed with regard to the statistical methods proposed in Tang and Wahed (2011). We briefly introduced the methods in Section 2.7. The examples in this section will be illustrated using a simulated dataset called PHdata which is included in the DTR package. The PHdata dataset was simulated as described in the simulation study of Tang and Wahed (2011) using the simPHdata() function.
Load the PHdata dataset into the workspace: R> data("PHdata", package = "DTR") The PHdata dataset is a data frame with 400 rows corresponding to patients and 7 columns corresponding to variables. In this dataset, 400 patients participated in the first-stage randomization to either A 1 or A 2 , and then were followed until their response status was observed at different times. Among 400 patients, 227 (57%) patients responded, and then participated in the second-stage randomization to either B 1 or B 2 . The censoring rate was 53% throughout the study with a maximum follow-up time of 14 years.
R> dim(PHdata) The generalized Cox model can be fitted by calling the PHfit() function. The function returns an object of class 'coxph' that is introduced in the survival package (Therneau and Grambsch 2000;Therneau 2015). The Cox proportional hazard model output can be printed by calling the summary() function of the 'coxph' object.
R> fit <-PHfit(data = PHdata, covar = "V") R> summary(fit) There was a significant difference in the survival distribution among the four DTRs (overall p < 0.001). Patients assigned to A 1 at the first stage had significantly different survival distribution compared to those assigned to A 2 . Among patients who were assigned to A 1 at the first stage, whether they were assigned to B 1 or B 2 significantly affected their survival outcomes. However, if patients were assigned to A 2 at the first stage, there was no significant difference in the survival irrespective of their assignment at the second stage.

Discussion
In this article, we illustrated the DTR package designed for estimating and comparing the effects of treatment regimes from SMARTs in the context of survival endpoints. The package implements the statistical estimating and testing procedures proposed in Lunceford et al. (2002); Guo and Tsiatis (2005); Tang and Wahed (2011); Tang and Wahed (2015); and Kidwell and Wahed (2013). Each statistical method was briefly introduced, and the functionality of the DTR package was demonstrated using the simulated data from various SMARTs as described in aforementioned papers. The simulated datasets were generated for illustration purpose, and may not reflect the dataset in a real-life situation. For example, in the WRSEdata and LRdata datasets, some response times were observed after the censoring times, which may not be realistic in SMARTs. The package allows users to perform the generation of artificial data from SMARTs, point estimation of the survival quantities, estimates of their standard errors, overall and pairwise comparisons of the effects of different DTRs, and a visualization of the survival curves for DTRs over the study period. We hope that this package can be useful for researchers in several areas where SMART designs with survival outcomes are applicable. Although it is still common to separately evaluate the first-and second-stage therapies when analyzing the survival data from SMARTs, we hope that the analyses could be substantially improved by using the functions provided by this package, as there has been a growing interest in studying the effects of different treatment sequences in biomedical research.

A. Computation of Equation 2
The detailed steps for computing VAR{Ŝ jk (t)} using Equation 2 are as follows: Step 1: Calculate the survival estimate regardless of treatment regimes: Step 2: Calculate: Step 3: Calculate: Step 4: Calculate:

B. Computation of Equation 4
The detailed steps for computingσ 2 using Equation 4 are as follows:

C. Computation of Equations 7 and 8
The detailed steps for computing Z W jkjk (t) using Equation 7 are as follows: where Y * jk (U i ) = n p=1 W jkp (t)Y p (U i ) and Y * jk (U i ) = n p=1 W jk p (t)Y p (U i ).