rpsftm : An R Package for Rank Preserving Structural Failure Time Models

Treatment switching in a randomised controlled trial occurs when participants change from their randomised treatment to the other trial treatment during the study. Failure to account for treatment switching in the analysis (i.e. by performing a standard intention-to-treat analysis) can lead to biased estimates of treatment efficacy. The rank preserving structural failure time model (RPSFTM) is a method used to adjust for treatment switching in trials with survival outcomes. The RPSFTM is due to Robins and Tsiatis (1991) and has been developed by White et al. (1997, 1999). The method is randomisation based and uses only the randomised treatment group, observed event times, and treatment history in order to estimate a causal treatment effect. The treatment effect, ψ, is estimated by balancing counter-factual event times (that would be observed if no treatment were received) between treatment groups. G-estimation is used to find the value of ψ such that a test statistic Z (ψ) = 0. This is usually the test statistic used in the intention-to-treat analysis, for example, the log rank test statistic. We present an R package, rpsftm, that implements the method.


Introduction
In a two-arm randomised controlled trial, participants are randomly allocated to receive one of two treatments or interventions.Ideally, all participants would fully receive their allocated treatment and no other treatment.However, in a recent review, 98 of 100 trials published in four high quality general medical journals reported some form of departure from randomised treatment (Dodd et al., 2012).When treatment is 'all-or-nothing' (for example, a surgical procedure), possible departures include not receiving the allocated treatment, receiving the other trial treatment, or receiving a non-trial treatment.When treatment is given over time (for example, a drug for HIV treatment), departures also occur over time, and often include starting a new treatment in response to a disease-related event.
This paper specifically focuses on the case of treatment switching over time, where participants may switch to receive the other trial treatment during the trial (Latimer et al., 2014).A common example is in a trial of active treatment versus placebo where placebo arm participants may start the active treatment in response to disease progression.
A randomised controlled trial with departures from randomised treatment is commonly analysed by intention-to-treat, in which departures from randomised treatment are ignored and all randomised participants are compared in the groups to which they were randomised (Higgins et al., 2011).This provides an unbiased comparison of two treatment policies, which accept that treatment may be changed, but does not compare the efficacies of the treatments themselves (White et al., 1999).
In order to compare the efficacies of the treatments, analysts often use per-protocol analysis, which censors participants when they depart from randomised treatment.However, this loses the advantage that randomisation produces comparable groups and instead introduces selection bias.Another possibility is to include treatment as a time-dependent variable in a Cox regression model.The issue with this method is that the estimate of treatment effect may not have a causal interpretation when switchers are prognostically different to non-switchers (Robins and Tsiatis, 1991).
The IPCW method is an extension of the per-protocol censoring approach, whereby the bias associated with censoring participants that depart from randomised treatment is removed by weighting the remaining non-switchers (Latimer et al., 2016).A model is formed for the probability of switching using baseline and time-dependent covariates.Time-varying weights are then obtained for each participant based on the inverse probability of not switching until a given time (Watkins et al., 2013).The method relies on data on prognostic factors for mortality that also affect the probability of switching being collected at both baseline and over time.This is called the 'no unmeasured confounders' assumption.Since this assumption is untestable, the method relies on good planning at the study design stage in order to collect information on all possible confounders.The method also has two potential issues: a) if many prognostic factors are included in the model to calculate weights the model may fail to converge, and b) if only a few participants switch then large weights are assigned to the remaining participants, which may result in biased analyses.The method has already been implemented in R in the package ipw (van der Wal and Geskus, 2011) and so is not considered further here.
Instead, this paper focuses on a popular causal method -the rank preserving structural failure time model (RPSFTM) (Robins and Tsiatis, 1991).In contrast to the IPCW method which requires potential confounders to be collected over time, the RPSFTM is randomisation based and only requires information on the randomised treatment group, observed event times, and treatment history in order to estimate a causal treatment effect.The method has been used, for example, in submissions to the UK's National Institute for Health and Clinical Excellence whose Decision Support Unit commissioned a guidance document (Latimer and Abrams, 2014).The RPSFTM has been implemented in Stata (White et al., 2002) but not in R.
This paper describes the theory of the RPSFTM and presents a new implementation of the method with the R package rpsftm (Bond and Allison, 2017), using data based on a trial in HIV infection.

Theory RPSFTM method and assumptions
The RPSFTM uses a causal model to produce counter-factual event times (the time that would be observed if no treatment were received) in order to estimate a causal treatment effect.Let T i = T off i + T on i be the observed event time for participant i, where T off i and T on i represent the time spent off and on treatment, respectively.The T i are related to the counter-factual event times, U i , via the following causal model: where exp (−ψ 0 ) is the acceleration factor associated with treatment and ψ 0 is the true causal parame- ter.
A grid search (g-estimation) procedure is used to estimate the treatment effect that balances the counter-factual event times across randomised treatment groups.To estimate the causal treatment effect, ψ, we assume that the U i are independent of randomised treatment group R i , i.e. if the groups are similar with respect to all other characteristics except treatment, the average event times should be the same in each group if no participant were treated.In general, this assumption is plausible in a randomised controlled trial.
A g-estimation procedure is used to find the value of ψ such that U i is independent of R i .For a range of ψs, the hypothesis ψ 0 = ψ is tested by computing U i (ψ) , subject to censoring, and calculating the test statistic Z (ψ) .This is usually the same test statistic as for the intention-to-treat analysis, for example, the log rank test statistic to compare survival distributions between groups.In the rpsftm() function, the possible test options are the log rank, and the Wald test from a Cox or Weibull regression model.For the parametric Weibull test, the point estimate ( ψ) is the value of ψ for which Z (ψ) = 0.For the non-parametric tests (log rank, Cox), ψ is the value of ψ for which Z (ψ) crosses 0, since Z (ψ) is a step function.Confidence intervals are similarly found with the 100 (1 After finding ψ, an adjusted hazard ratio can be calculated.For example, by comparing counterfactual event times at ψ in the control group to the observed event times in the experimental group we can estimate the treatment effect that would have been observed in the absence of switching (in the case where only the control group switch).
As well as assuming that the only difference between randomised groups is the treatment received, the RPSFTM also assumes a 'common treatment effect'.The common treatment effect assumption states that the treatment effect is the same for all participants (with respect to time spent on treatment) regardless of when treatment is received, which is implicit in equation (1).

Re-censoring
We assume that censoring occurs if participants survive to a specific calendar date (Robins and Tsiatis, 1991).Thus the potential censoring time C i is known for all participants.The censoring indicators of the observed event times are initially carried over to the counter-factual event times.However, the uninformative censoring on the T i scale may be informative on the U i scale.Suppose we have two participants with the same U i , one of whom receives the superior treatment.The participant receiving The R Journal Vol.9/2, December 2017 ISSN 2073-4859 the superior treatment has their U i extended so that they are censored whilst the other participant may observe the event.
In detail, C i , the censoring times for the T i , are transformed to C i (1 − P i + P i exp (ψ)) when we work on the U i scale, where P i is a random variable, the proportion of time on treatment, . This transformation of censoring times occurs by representing the censoring status in two equivalent ways: {C i < T i }, and re-scaling both sides {C i (1 We cannot assume the variable P i is independent of U i .For example, adherence to a protocolised treatment may depend on the underlying prognosis.
We can overcome this induced dependence by considering the sample space for P i marginally over U i and, optionally, conditioning on R i .If we replace P i with a function of its sample space then any dependency on U i is thus removed.Any alternative transformed censoring times must be smaller than the original censoring times, else we may have to impute uncensored U i values for censored T i observations, which would be impossible.Hence taking the minimum value from the sample space for P i meets both desiderata.If the sample space differs between arms, say switching is impossible in one arm, then this can be utilised to potentially observe more events with gains in efficiency.
Operationally, let C i be the potential censoring time for participant i.A participant is then recensored at the minimum possible censoring time: , then U i is replaced by D * i and the censoring indicator is replaced by 0. For treatment arms where switching does not occur, there can be no informative censoring and so re-censoring is not applied.

Sensitivity analysis
As previously mentioned, the RPSFTM has two assumptions: 1.The only difference between randomised groups is the treatment received.
2. The treatment effect is the same for all participants regardless of when treatment is received.
Whilst the first assumption is plausible in a randomised controlled trial, the latter may be unlikely to hold.For example, if control group participants can only switch at disease progression then the treatment benefit may be different in these participants compared to those randomised to the experimental treatment.The rpsftm() function allows for investigation of deviations from the common treatment effect assumption by featuring a treatment-effect modifier variable which means the treatment effect can be varied across participants.The assumption that defines the counter-factual treatment-free event times U i is thus modified by multiplying ψ by some factor k i > 0: The value taken by k i is derived from observed data and is left to the user to assign.The package assumes that k i is determined at baseline, and re-censoring is undertaken in a similar way by re-censoring at the minimum possible censoring time: , then U i is replaced by D * i and the censoring indicator is replaced by 0. The case where k i is observed post randomisation is equally valid in terms of defining U i , but more complicated re-censoring may be required, taking the minimum value of D * i over the sample space of k i conditional on arm.This is not implemented in the package.

Using package rpsftm
The main function in the package rpsftm is of the same name, rpsftm(), which returns an object that has print, summary, and plot methods, and that inherits from the class used to define the test statistic (survdiff, coxph or survreg) .The arguments to rpsftm() are given in table 1 below.

Example
The rpsftm function will be illustrated using a simulated dataset immdef taken from (White et al., 2002), based on a randomized controlled trial (Concorde Coordinating Committee, 1994) The number of points between hi_psi and low_psi at which to evaluate the Z-statistics in the estimating equation.Default is 100.

Intention-to-treat analysis
First, we estimate the effect of giving zidovudine ignoring any treatment changes during the trial.This is found by fitting an accelerated failure time model using the eha package (Broström, 2017): > library(eha) > itt_fit <-aftreg (Surv(progyrs, prog)  The intention-to-treat estimate is −0.147 which means that lifetime is used up exp (−0.147) = 0.863 times as fast when on zidovudine as when off zidovudine (or zidovudine extends lifetime by a factor of exp (0.147) = 1.16).

Fitting the RPSFTM
We now show how to use rpsftm with the immdef data.First, a variable rx for the proportion of time spent on treatment must be created: rx <-with(immdef, 1 -xoyrs / progyrs) This sets rx to 1 in the immediate treatment arm (since no participants could switch to the deferred arm, and xoyrs = 0 in such participants), 0 in the deferred arm participants that did not receive treatment (as xoyrs = progyrs in such participants) and 1 -xoyrs / progyrs in the deferred arm participants that did receive treatment.Using the default options, the fitted model is rpsftm_fit_lr <-rpsftm(formula = Surv(progyrs, prog) ~rand(imm, rx), data = immdef, censor_time = censyrs) The above formula fits an RPSFTM where progyrs is the observed event time, prog is the indicator of disease progression, imm is the randomised treatment group indicator, rx is the proportion of time spent on treatment and censyrs is the censoring time.The log rank test is used in finding the point estimate of ψ.Re-censoring is performed since the censor_time parameter is specified; if not specified then re-censoring would not be performed.After finding ψ, rpsftm refits the model at ψ and produces a survdiff object of the counter-factual event times to be used in plotting Kaplan-Meier curves.The list of objects that the function returns is given in table 3.
The point estimate and 95% confidence interval (CI) can be returned using rpsftm_fit_lr$psi and rpsftm_fit_lr$CI which gives ψ = −0.181(−0.350, 0.002) , slightly larger than in the intention- to-treat analysis.This means that lifetime is used up exp (−0.181) = 0.834 times as fast when on zidovudine as when off zidovudine, i.e. the time to progression to AIDS, or CDC group IV disease, or death is longer when on zidovudine (treatment is beneficial).However, the confidence interval contains zero, which suggests the treatment effect is non-significant.The function plot() produces Kaplan-Meier curves of the counter-factual event times in each group and can be used to check that the distributions are indeed the same at ψ as shown in figure 1 We now provide examples of using the Cox regression model and the Weibull model in place of the log rank test.To use the Wald test from a Cox regression model, we specify test = coxph in the function parameters.Covariates can also be included in the estimation procedure by adding them to the right hand side of the formula.For example, baseline covariates that are included in the intention-to-treat analysis may also be incorporated into the estimation procedure of the RPSFTM.In the following example we add entry time as a covariate and use summary() to find the value of ψ and its 95% CI.
> rpsftm_fit_cph <-rpsftm(formula = Surv(progyrs, prog)   As we can see from the output for the Cox and Weibull models, the estimates of ψ are similar to the estimate obtained from using the log rank test.Both fitted models could be used as inputs to the plot() function, which produce figures very similar to figure 1.
As a sensitivity analysis we can investigate what would happen to the estimate of ψ if the treatment effect in the deferred treatment group was half of that in the immediate treatment group by setting k i = 1 for participants in the latter group and k i = 0.5 for participants in the former group.
The limits of the confidence intervals may not exist at all if Z (•) reaches an asymptote before it hits the required quantile of the standard normal distribution, in which case the limits should be reported as ±∞.This scenario would be illustrated by a similar plot to figure 2.
Another possibility is for the coxph function to fail to converge.This occurs when the maximum likelihood estimate of a coefficient is infinity, e.g. if one of the treatment groups has no events.The coxph documentation states that the Wald statistic should be ignored in this case and therefore the rpsftm output should be taken with caution.

Limitations
As well as the potential computational issues highlighted in the previous section, the RPSFTM itself has some limitations.The method relies on the assumption that the treatment effect is the same for all participants regardless of when treatment is received.We have allowed for investigation of deviations from this assumption within the rpsftm() function by adding a treatment-effect modifier variable.Another possible limitation of the model is its requirement for only the total amount of time spent on/off treatment.In instances where participants can switch back and forth between treatments the model may be inefficient.

Conclusion
The intention is to fill a void in the methodology available to R users by providing this package and thus facilitate the adoption of newer methods in application to real data in future.The sensitivity analyses, which generalise the definition of the counter-factual treatment-free event times, are an original development.
Further developments to the package may include expanding the functionality to allow multiarmed studies with three or more arms.Such an extension would have at its core a set of g-estimating equations to solve, one for each element of the vector of ψ parameters.Defining the syntax to capture the multi-dimensional metric of time on treatments is challenging, as will be the numerical computational details when the g-estimating equations are step functions.
The R Journal Vol.9/2, December 2017 ISSN 2073-4859 Further terms can be added to the right hand side to adjust for covariates.data an optional data frame containing the variables censor_time variable or constant giving the time at which censoring would, or has occurred.This should be provided for all observations unlike standard Kaplan-Meier or Cox regression where it is only given for censored observations.If no value is given then re-censoring is not applied.subset an expression indicating which subset of the rows of data should be used in the fit.This can be a logical vector, a numeric vector indicating which observation numbers are to be included, or a character vector of row names to be included.All observations are included by default.na.action a missing-data filter function.This is applied to the model.frameafter any subset argument has been used.Default is options()$na.action.test one of survdiff, coxph or survreg.Describes the test to be used in the estimating equation.Default is survdiff.low_psi the lower limit of the range to search for the causal parameter.

Table 1 :
Arguments for the rpsftm() two policies (immediate or deferred treatment) of zidovudine treatment in symptom free participants infected with HIV.The immediate treatment arm received treatment at randomisation whilst the deferred arm received treatment either at onset of AIDS related complex or AIDS (CDC group IV disease) or development of persistently low CD4 count.The endpoint considered here was time from study entry to progression to AIDS, or CDC group IV disease, or death.

Table 2 :
Description of simulated data set Surv() data using the estimate value of psi to give counterfactual untreated failure times.rand the rand() object used to specify the allocated and observed amount of treatment.ansthevalues from uniroot.all used to solve the estimating equation, but embedded within a list as per uniroot, with an extra element root_all, a vector of all roots found in the case of multiple solutions.The first element of root_all is subsequently used.eval_z a data frame with the Z-statistics from the estimating equation evaluated at a sequence of values of psi.Used to plot and check if the range of values to search for solution and limits of confidence intervals need to be modified.Further elements corresponding to either a survdiff, coxph, or survreg object.This will always include:call the R call object formula a formula representing any adjustments, strata or clusters -used for the update() function terms a more detailed representation of the model formula

Table 3 :
Outputs from the rpsftm()For the Weibull model we specify test = survreg in the function parameters.: function Figure 1: Output from plot(rpsftm_fit_lr)