Title

strbee -- Estimate a treatment effect adjusting for treatment switches

Introduction

strbee analyzes a two-arm clinical trial with a survival outcome, in which
some subjects may "crossover" or "switch" to receive the treatment of the other arm.

The model is the Rank-Preserving Structural Nested Failure Time Model
(RPSFTM) of Robins and Tsiatis (1991).
This model uses an accelerated life model to relate the observed event time T to
an underlying event time U that would have been observed in the absence of
treatment.
The effect of treatment is expressed through a parameter psi
which represents the log of the "acceleration factor",
the factor by which life is accelerated by treatment.
Negative values of psi mean that treatment improves survival.

Psi is estimated as the value at which U is balanced between the treatment arms
(on a user-specified test applied to all patients whether or not they switch).
Point estimation ("G-estimation") is performed by searching over a range of values of psi for a value where the test statistic Z(psi) is close to 0.
Confidence limits are similarly found where the test statistic Z(psi) is close to its critical values.

strbee was originally described in White et al (2002). This version has new features described below.

Syntax

Estimation syntax:

strbee treatvar [if exp] [in range]
[,
xo0(timevar eventvar)
xo1(timevar eventvar)
endstudy(varname)
psimult(exp)
test(logrank|wilcoxon|cox|weibull|exponential)
adjvars(varlist)
strata(varlist)
ipe
ipecens
psistep(#)
tol(#)
noci
trace
psistart(#)
maxiter(#)
savedta(filename[,append|replace])
ipestore(name)
common_options
]

where treatvar is the treatment arm, which must have values 0 and 1;
timevar and eventvar are time to switch and an indicator of switching.

This is the standard syntax. Results for psi and Z(psi) may be
stored or appended to an existing results file.

Replay syntax:

strbee [using filename]
[,
common_options
]

This uses previously stored results for psi and Z(psi).
If using filename is omitted and strbee was the last r-class command run, then the results from the last run of strbee are used.

common_options are

list
psimin(#)
psimax(#)
zgraph[(graph_options)]
level(#)
hr
kmgraph[(graph_options)]
gen(newvarname)

Description

strbee is an st command: data must be stset before running strbee.
Each subject must have a single record starting at time 0.

Censored data present an extra problem for strbee, since noninformative
censoring of T implies informative censoring of U. To avoid bias, users must
specify the potential censoring time for all subjects: strbee then computes
a recensoring time, which may be earlier than the actual censoring time, and
recensors the data.

Options to specify the model (estimation syntax)

Simple switching patterns may be specified using xo0() and/or xo1().
More complex switching patterns may be specified using ton() and/or toff().
These options cannot be used together.

xo0(timevar eventvar) specifies time to switch and a switch
indicator in arm 0. If not specified, it is assumed that there is no
switch from arm 0 to arm 1.

xo1(timevar eventvar) does the same for arm 1. If not specified, it is
assumed that there is no switch from arm 1 to arm 0.

ton(varname[,suboptions]) is an alternative way to specify times to switch.
It is suitable for more complex switching patterns.
varname must be a variable containing the total time on treatment for each individual.

The suboptions are min0(exp), min1(exp), max0(exp)
and max1(exp),
which allow the user to make the recensoring procedure more efficient by specifying logical limits on
the lengths of time spent on treatment. For example, max0(1) max1(1) would be appropriate if it was not possible to receive treatment, in either arm, for more than 1 year (and assuming year is the unit of time).

toff(varname[,suboptions]) specifies the total time off treatment for each individual.
The suboptions are the same as for ton().
toff() may be used as an alternative to ton() or as a complement.
Continuing the example given under ton(), if it were also known that no control arm individual could switch in the first 6 months, then one might code ton(tonvar, max0(1) max1(1)) toff(, min(0.5)).

endstudy(varname) specifies the time of the end of study (the potential
censoring time). U values are then recensored at the minimum possible
potential censoring time on the U scale, where the minimization is carried
out over all treatment profiles possible for the subject's randomized arm.
The potential censoring time must be specified both for censored
and uncensored subjects.

Censoring due to random events (e.g. competing risk or loss to follow-up)
should ideally be treated differently from censoring at end of study.
However, if small amounts of random censoring are present, a reasonable
approximation is to set the potential censoring time for subjects who are
randomly censored equal to the actual time of random censoring.

psimult(exp) allows the acceleration factor to vary between individuals.
Specifying psimult(k), for example, specifies that the log acceleration factor
is psi*k where psi is to be estimated.
If k is a variable then this allows the acceleration factor to vary between individuals.
This option is typically used to assess sensitivity to departures from the common treatment effect assumption:
for example, to fit the model with the effect of treatment 30% smaller in the control arm,
you would specify psimult(k) where k is 1 in the experimental arm and 0.7 in the control arm.

Options to specify the estimator (estimation syntax)

test(sts_test_option|stcox|cox|streg_dist)
specifies the test used
to compare the underlying event time U between the two randomized arms.
test(sts_test_option) uses sts test with the given option: for example, the commonest choice test(logrank) uses sts test treatvar, logrank.
test(streg_dist) uses the Wald test from streg with the specified distribution: for example, test(weibull) uses streg treatvar, dist(weibull).
test(stcox|cox) uses the Wald test from stcox.
The default is the logrank test (or the weibull test if ipe is specified).

adjvars(varlist) modifies
regression-based tests (that is, using test(streg_dist) or
test(stcox|cox) by
adjusting for the specified variables.

strata(varlist) requests the tests to be stratified by the specified variables.

ipe invokes Branson and Whitehead's Iterative Parameter Estimation (IPE).
This uses a parametric accelerated life model and estimates the parameter psi by iteratively comparing the counterfactual survival times of arm 1 if completely treated with those of arm 0 if completely untreated.
test(weibull|exponential) is required.
Note that the standard error and confidence interval from this method do not allow for all the uncertainty, and bootstrap methods are preferable.

ipecens modifies IPE to use Branson and Whitehead's procedure for recensoring. The standard Robins-Tsiatis recensoring procedure is preferable (White, 2006) and is used by default.

Options to specify the search procedure (estimation syntax)

psimin(#) and psimax(#): with the estimation syntax, these specify the extreme permitted values of the parameter psi.
Defaults are -10 and 10.
With the replay syntax, they instead affect the zgraph - see options controlling output.

psistep(#) with #>0 specifies the step size between the extreme values for
a grid search, while
psistep(0) (the default) specifies an interval-bisection approach for psi.

tol(#) specifies the convergence criterion for interval-bisection
estimation. strbee searches until bounds for each solution differ by
less than 10(-tol). tol() also defines the number of decimal places
reported. The default is tol(3).

noci suppresses searching for the confidence limits. strbee still reports
what limited information it has about them.

level(#) specifies searching the confidence interval, in percent, for
confidence intervals. The default is level(95) or as set by set
level; see level.

trace gives details on recensoring and prints psi and the test statistic at
each step. It also saves a file rbeetrace.dta containing variables u,
du, z0, dz0, z1, dz1, and recens evaluated at the last value
of psi used. Note that with grid search this is the value of psimax(),
and with interval-bisection search it is the upper confidence limit
without noci or the point estimate with noci.

psistart(#) specifies a starting value for the IPE method. Default is 0.

maxiter(#) limits the number of iterations for the IPE method. Default is 100.

Options controlling storage of results (estimation syntax)

savedta(filename[,append|replace]) directs the values of psi and the test statistic to the specified file which can later be used with the replay syntax.
If append is specified, strbee checks that the current test statistics
and the stored test statistics were computed using the same model. In this
case, the results reported only use test statistics computed in this run of
strbee: to get full results, run strbee using filename.
If savedta() is omitted then results are stored in _strbee_results.dta.

ipestore(name) stores the last model fit in the IPE algorithm. It can later be retrieved using estimates store name.

Options controlling output (both syntaxes)

list lists the values of psi and the test statistic.

zgraph[(graph_options)]
graphs the test statistic Z against psi.
graph_options are most of the options allowed with graph, twoway; see graph.

psimin(#) and psimax(#): with the replay syntax only, these control the values plotted on the graph, if zgraph is specified.

hr outputs the estimated hazard ratio comparing arm 1 if always
treated with arm 0 if never treated.
A test-based confidence interval is given.
For details, see White et al (1999).
If adjvars(varlist) is also specified
then the hazard ratio is adjusted for varlist.

kmgraph[(suboptions)]
draws, for each arm, the Kaplan-Meier graph for the observed event times (labelled "as observed")
and also the counterfactual event times
if that arm had never received treatment (labelled "if untreated").
The counterfactual untreated graph is not drawn for arm 0 if there are no switches in that arm.

Possible suboptions are:

most of the options allowed with sts graph; see sts graph.

showall causes the Kaplan-Meier graph for the counterfactual event times
if that arm had always received treatment (labelled "if fully treated").
However the fully-treated graph for arm 1 is not drawn if there are no switches in that arm.

untreated causes only the two graphs for the counterfactual untreated event times to be drawn.

lpatterns(pattern1 pattern2) lists the line patterns
for the control arm (default: dash) and then the treatment arm (default: solid).

lcolors(color1 color2 color3) lists the line colours
for the observed data (default: black),
then the counterfactual untreated data (default: orange),
and then the counterfactual fully treated data (default: blue).

show(#....) specifies exactly which graphs to show:
0, 10, 20 for arm 0 as observed, if untreated and if fully-treated;
1, 11, 21 for arm 1 as observed, if untreated and if fully-treated.

gen(newvarname) generates
a new variable newvarname containing the values of the counterfactual untreated outcome;
a new variable dnewvarname containing the event indicator for newvarname;
and if endstudy() is specified, a new variable cnewvarname containing the censoring time of newvarname.

With the IPE method, newvarname is instead the ideal outcome, defined as the counterfactual untreated outcome in arm 0 and the counterfactual fully-treated outcome in arm 1.

Changes from version 1.2 to version 1.5

Main new options:
hr
adjvars()
kmgraph()
gen().

ipe, formerly an undocumented option, is now fully documented.

psimin() and psimax() are now unlikely to be needed,
since by default the search for suitable parameter values starts with the range -1 to 1,
but is extended as far as necessary towards -10 and 10.

The replay syntax is new.

Some bugs have been fixed, including:
switch times greater than event times are now handled correctly;
test(wilcoxon) now works correctly.

graph has been renamed zgraph, but the former syntax is still allowed.

Graphs use the improved graphics introduced in Stata 8.

Changes from version 1.5 to version 1.8

Main new options:
strata()
ton()
toff().

More tests available: almost everything available with sts test or streg.

New options for kmgraph.

Estimation by interval bisection

1. An initial interval is found by searching for two values of psi with Z(psi) above the upper critical value (1.96 with 5% significance level) and below the lower critical value.

2. A value of psi solving Z(psi)=0 is found by evaluating Z(psi) at the midpoint of the interval, narrowing psi down to the appropriate half of the interval, and repeating until the desired accuracy is achieved.

3. Value of psi solving Z(psi)= the lower and upper critical values are found similarly.

Interval bisection may give wrong answers if the test statistic is
nondecreasing in psi, and should always be checked using the zgraph
option.

Examples

Example using simulated data of immediate vs. deferred treatment:

. use http://www.homepages.ucl.ac.uk/~rmjwiww/stata/noncomp/immdef.dta, clear
. stset progyrs prog

Intention-to-treat analysis

. strbee imm

RBEE analysis without recensoring

. strbee imm, xo0(xoyrs xo)

RBEE analysis with recensoring

. strbee imm, xo0(xoyrs xo) endstudy(censyrs) savedta(recens)

Grid search to check for multiple axis-crossings

. strbee imm, xo0(xoyrs xo) endstudy(censyrs) savedta(recens,append) psimin(-0.5) psimax(0.1) psistep(0.02) zgraph(title(RBEE with recensoring))

Hazard-based interpretation of results

. strbee using recens, hr kmgraph

The same using the ton() syntax

. gen ton = progyrs if imm
. replace ton = progyrs - xoyrs if !imm & xo
. replace ton = 0 if !imm & !xo
. strbee imm, ton(ton) endstudy(censyrs)

Troubleshooting

If you run into trouble, try the following:

1. Use the trace option to see what values of psi are being considered and how much recensoring is occurring.

2. Use the zgraph option to inspect the graph of Z(psi) against psi. Ideally it should start above 1.96 (or other critical value if used) and descend to below -1.96 as psi increases.

3. Use grid search instead of interval bisection.

4. Try updating the software (see below).

If all else fails, feel free to email me, ian.white@ucl.ac.uk, sending me a log file showing your problem.
The log file should include at least which strbee, a summary of your data and an strbee call with the trace option.
It's very helpful if you can also send me some data to illustrate your problem.

Technical notes

strbee is an r-class command.
Typed without options (or with only output options), it should redisplay the latest results provided no later r-class command has been run.

Results for psi and Z(psi) are saved in the file specified by savedta() or else in _strbee_results.dta.
This file contains details of the original data set as characteristics.
If results are appended to this file then strbee checks that the characteristics match.

Limitations

. Censored switches occurring before the event are assumed to represent no switch.

. Switches occurring after the event are censored at the event time.

. Only two-arm trials are allowed.

. Treatment at any time must be yes/no, so that it can be summarised by the total time on or off treatment. Thus varying doses can not be allowed for in strbee, though the RPSFTM does allow for them.

References

Branson M, Whitehead J (2002). Estimating a treatment effect in survival studies in which patients switch treatment. Statistics in Medicine 21: 2449-2463.

Robins JM, Tsiatis AA (1991). Correcting for non-compliance in randomized trials using rank preserving structural failure time models. Communications in Statistics - Theory and Methods 20: 2609-2631.

*White IR, Walker S, Babiker A (2002). strbee: Randomisation-based efficacy estimator. Stata Journal 2: 140-150.

White IR, Babiker AG, Walker S, Darbyshire JH (1999). Randomisation-based methods for correcting for treatment changes: examples from the Concorde trial. Statistics in Medicine 18: 2617-2634.

White IR (2006).
Estimating treatment effects in randomised trials with treatment switching.
Statistics in Medicine 25: 1619-1622.

* Please use this reference to cite this program.

Author and updates

Ian White, MRC Clinical Trials Unit at UCL, London, UK.
Email ian.white@ucl.ac.uk.

I put updates from time to time on my website.
You can install them using
net install http://www.homepages.ucl.ac.uk/~rmjwiww/stata/noncomp/strbee, replace.

You can get the latest version of all my Stata software using
net from http://www.homepages.ucl.ac.uk/~rmjwiww/stata/.

The original code was co-authored with
Sarah Walker and Abdel Babiker,
MRC Clinical Trials Unit at UCL, London, UK.