spower
Simulate Power of 2Sample Test for Survival under Complex Conditions
Given functions to generate random variables for survival times and
censoring times, spower
simulates the power of a usergiven
2sample test for censored data. By default, the logrank (Cox
2sample) test is used, and a logrank
function for comparing 2
groups is provided. Optionally a Cox model is fitted for each each
simulated dataset and the log hazard ratios are saved (this requires
the survival
package). A print
method prints various
measures from these. For composing R functions to generate random
survival times under complex conditions, the Quantile2
function
allows the user to specify the intervention:control hazard ratio as a
function of time, the probability of a control subject actually
receiving the intervention (dropin) as a function of time, and the
probability that an intervention subject receives only the control
agent as a function of time (noncompliance, dropout).
Quantile2
returns a function that generates either control or
intervention uncensored survival times subject to nonconstant
treatment effect, dropin, and dropout. There is a plot
method
for plotting the results of Quantile2
, which will aid in
understanding the effects of the two types of noncompliance and
nonconstant treatment effects. Quantile2
assumes that the
hazard function for either treatment group is a mixture of the control
and intervention hazard functions, with mixing proportions defined by
the dropin and dropout probabilities. It computes hazards and
survival distributions by numerical differentiation and integration
using a grid of (by default) 7500 equallyspaced time points.
The logrank
function is intended to be used with spower
but it can be used by itself. It returns the 1 degree of freedom
chisquare statistic, with the hazard ratio estimate as an attribute.
The Weibull2
function accepts as input two vectors, one
containing two times and one containing two survival probabilities, and
it solves for the scale and shape parameters of the Weibull distribution
($S(t) = exp(\alpha*t^\gamma)$)
which will yield
those estimates. It creates an R function to evaluate survival
probabilities from this Weibull distribution. Weibull2
is
useful in creating functions to pass as the first argument to
Quantile2
.
The Lognorm2
and Gompertz2
functions are similar to
Weibull2
except that they produce survival functions for the
lognormal and Gompertz distributions.
When cox=TRUE
is specified to spower
, the analyst may wish
to extract the two margins of error by using the print
method
for spower
objects (see example below) and take the maximum of
the two.
Usage
spower(rcontrol, rinterv, rcens, nc, ni, test=logrank, cox=FALSE, nsim=500, alpha=0.05, pr=TRUE)
"print"(x, conf.int=.95, ...)
Quantile2(scontrol, hratio, dropin=function(times)0, dropout=function(times)0, m=7500, tmax, qtmax=.001, mplot=200, pr=TRUE, ...)
"print"(x, ...)
"plot"(x, what=c("survival", "hazard", "both", "drop", "hratio", "all"), dropsep=FALSE, lty=1:4, col=1, xlim, ylim=NULL, label.curves=NULL, ...)
logrank(S, group)
Gompertz2(times, surv)
Lognorm2(times, surv)
Weibull2(times, surv)
Arguments
 rcontrol

a function of n which returns n random uncensored
failure times for the control group.
spower
assumes that noncompliance (dropin) has been taken into account by this function.  rinterv

similar to
rcontrol
but for the intervention group  rcens
 a function of n which returns n random censoring times. It is assumed that both treatment groups have the same censoring distribution.
 nc
 number of subjects in the control group
 ni
 number in the intervention group
 scontrol
 a function of a time vector which returns the survival probabilities for the control group at those times assuming that all patients are compliant.
 hratio
 a function of time which specifies the intervention:control hazard ratio (treatment effect)
 x

an object of class “Quantile2” created by
Quantile2
, or of class “spower” created byspower
 conf.int
 confidence level for determining foldchange margins of error in estimating the hazard ratio
 S

a
Surv
object or other twocolumn matrix for rightcensored survival times  group

group indicators have length equal to the number of rows in
S
argument.  times
 a vector of two times
 surv
 a vector of two survival probabilities
 test

any function of a
Surv
object and a grouping variable which computes a chisquare for a twosample censored data test. The default islogrank
.  cox

If true
TRUE
the two margins of error are available by using theprint
method forspower
objects (see example below) and taking the maximum of the two.  nsim
 number of simulations to perform (default=500)
 alpha
 type I error (default=.05)
 pr

If
FALSE
preventsspower
from printing progress notes for simulations. IfFALSE
preventsQuantile2
from printingtmax
when it calculatestmax
.  dropin
 a function of time specifying the probability that a control subject actually is treated with the new intervention at the corresponding time
 dropout

a function of time specifying the probability of an intervention
subject dropping out to control conditions. As a function of time,
dropout
specifies the probability that a patient is treated with the control therapy at time t.dropin
anddropout
form mixing proportions for control and intervention hazard functions.  m
 number of time points used for approximating functions (default is 7500)
 tmax

maximum time point to use in the grid of
m
times. Default is the time such thatscontrol(time)
isqtmax
.  qtmax

survival probability corresponding to the last time point used for
approximating survival and hazard functions. Default is 0.001. For
qtmax
of the time for which a simulated time is needed which corresponds to a survival probability of less thanqtmax
, the simulated value will betmax
.  mplot
 number of points used for approximating functions for use in plotting (default is 200 equally spaced points)
 ...

optional arguments passed to the
scontrol
function when it's evaluated byQuantile2
. Unused forprint.spower
.  what

a single character constant (may be abbreviated) specifying which
functions to plot. The default is "both" meaning both
survival and hazard functions. Specify
what="drop"
to just plot the dropin and dropout functions,what="hratio"
to plot the hazard ratio functions, or "all" to make 4 separate plots showing all functions (6 plots ifdropsep=TRUE
).  dropsep

If
TRUE
makesplot.Quantile2
separate pure and contaminated functions onto separate plots  lty
 vector of line types
 col
 vector of colors
 xlim
 optional xaxis limits
 ylim
 optional yaxis limits
 label.curves

optional list which is passed as the
opts
argument tolabcurve
.
Value
spower
returns the power estimate (fraction of simulated
chisquares greater than the alphacritical value). If
cox=TRUE
, spower
returns an object of class
“spower” containing the power and various other quantities.Quantile2
returns an R function of class “Quantile2”
with attributes that drive the plot
method. The major
attribute is a list containing several lists. Each of these sublists
contains a Time
vector along with one of the following:
survival probabilities for either treatment group and with or without
contamination caused by noncompliance, hazard rates in a similar way,
intervention:control hazard ratio function with and without
contamination, and dropin and dropout functions.logrank
returns a single chisquare statistic.Weibull2
, Lognorm2
and Gompertz2
return an R
function with three arguments, only the first of which (the vector of
times
) is intended to be specified by the user.
Side Effects
spower
prints the interation number every 10 iterations if
pr=TRUE
.
References
Lakatos E (1988): Sample sizes based on the logrank statistic in complex clinical trials. Biometrics 44:229241 (Correction 44:923).
Cuzick J, Edwards R, Segnan N (1997): Adjusting for noncompliance and contamination in randomized clinical trials. Stat in Med 16:10171029.
Cook, T (2003): Methods for midcourse corrections in clinical trials with survival outcomes. Stat in Med 22:34313447.
Barthel FMS, Babiker A et al (2006): Evaluation of sample size and power for multiarm survival trials allowing for nonuniform accrual, nonproportional hazards, loss to followup and crossover. Stat in Med 25:25212542.
See Also
Examples
# Simulate a simple 2arm clinical trial with exponential survival so
# we can compare power simulations of logrankCox test with cpower()
# Hazard ratio is constant and patients enter the study uniformly
# with followup ranging from 1 to 3 years
# Dropin probability is constant at .1 and dropout probability is
# constant at .175. Twoyear survival of control patients in absence
# of dropin is .8 (mortality=.2). Note that hazard rate is log(.8)/2
# Total sample size (both groups combined) is 1000
# % mortality reduction by intervention (if no dropin or dropout) is 25
# This corresponds to a hazard ratio of 0.7283 (computed by cpower)
cpower(2, 1000, .2, 25, accrual=2, tmin=1,
noncomp.c=10, noncomp.i=17.5)
ranfun < Quantile2(function(x)exp(log(.8)/2*x),
hratio=function(x)0.7283156,
dropin=function(x).1,
dropout=function(x).175)
rcontrol < function(n) ranfun(n, what='control')
rinterv < function(n) ranfun(n, what='int')
rcens < function(n) runif(n, 1, 3)
set.seed(11) # So can reproduce results
spower(rcontrol, rinterv, rcens, nc=500, ni=500,
test=logrank, nsim=50) # normally use nsim=500 or 1000
## Not run:
# # Run the same simulation but fit the Cox model for each one to
# # get log hazard ratios for the purpose of assessing the tightness
# # confidence intervals that are likely to result
#
# set.seed(11)
# u < spower(rcontrol, rinterv, rcens, nc=500, ni=500,
# test=logrank, nsim=50, cox=TRUE)
# u
# v < print(u)
# v[c('MOElower','MOEupper','SE')]
# ## End(Not run)
# Simulate a 2arm 5year followup study for which the control group's
# survival distribution is Weibull with 1year survival of .95 and
# 3year survival of .7. All subjects are followed at least one year,
# and patients enter the study with linearly increasing probability after that
# Assume there is no chance of dropin for the first 6 months, then the
# probability increases linearly up to .15 at 5 years
# Assume there is a linearly increasing chance of dropout up to .3 at 5 years
# Assume that the treatment has no effect for the first 9 months, then
# it has a constant effect (hazard ratio of .75)
# First find the right Weibull distribution for compliant control patients
sc < Weibull2(c(1,3), c(.95,.7))
sc
# Inverse cumulative distribution for case where all subjects are followed
# at least a years and then between a and b years the density rises
# as (time  a) ^ d is a + (ba) * u ^ (1/(d+1))
rcens < function(n) 1 + (51) * (runif(n) ^ .5)
# To check this, type hist(rcens(10000), nclass=50)
# Put it all together
f < Quantile2(sc,
hratio=function(x)ifelse(x<=.75, 1, .75),
dropin=function(x)ifelse(x<=.5, 0, .15*(x.5)/(5.5)),
dropout=function(x).3*x/5)
par(mfrow=c(2,2))
# par(mfrow=c(1,1)) to make legends fit
plot(f, 'all', label.curves=list(keys='lines'))
rcontrol < function(n) f(n, 'control')
rinterv < function(n) f(n, 'intervention')
set.seed(211)
spower(rcontrol, rinterv, rcens, nc=350, ni=350,
test=logrank, nsim=50) # normally nsim=500 or more
par(mfrow=c(1,1))
# Compose a censoring time generator function such that at 1 year
# 5% of subjects are accrued, at 3 years 70% are accured, and at 10
# years 100% are accrued. The trial proceeds two years past the last
# accrual for a total of 12 years of followup for the first subject.
# Use linear interporation between these 3 points
rcens < function(n)
{
times < c(0,1,3,10)
accrued < c(0,.05,.7,1)
# Compute inverse of accrued function at U(0,1) random variables
accrual.times < approx(accrued, times, xout=runif(n))$y
censor.times < 12  accrual.times
censor.times
}
censor.times < rcens(500)
# hist(censor.times, nclass=20)
accrual.times < 12  censor.times
# Ecdf(accrual.times)
# lines(c(0,1,3,10), c(0,.05,.7,1), col='red')
# spower(..., rcens=rcens, ...)
## Not run:
# # To define a control survival curve from a fitted survival curve
# # with coordinates (tt, surv) with tt[1]=0, surv[1]=1:
#
# Scontrol < function(times, tt, surv) approx(tt, surv, xout=times)$y
# tt < 0:6
# surv < c(1, .9, .8, .75, .7, .65, .64)
# formals(Scontrol) < list(times=NULL, tt=tt, surv=surv)
#
# # To use a mixture of two survival curves, with e.g. mixing proportions
# # of .2 and .8, use the following as a guide:
# #
# # Scontrol < function(times, t1, s1, t2, s2)
# # .2*approx(t1, s1, xout=times)$y + .8*approx(t2, s2, xout=times)$y
# # t1 < ...; s1 < ...; t2 < ...; s2 < ...;
# # formals(Scontrol) < list(times=NULL, t1=t1, s1=s1, t2=t2, s2=s2)
#
# # Check that spower can detect a situation where generated censoring times
# # are later than all failure times
#
# rcens < function(n) runif(n, 0, 7)
# f < Quantile2(scontrol=Scontrol, hratio=function(x).8, tmax=6)
# cont < function(n) f(n, what='control')
# int < function(n) f(n, what='intervention')
# spower(rcontrol=cont, rinterv=int, rcens=rcens, nc=300, ni=300, nsim=20)
#
# # Do an unstratified logrank test
# library(survival)
# # From SAS/STAT PROC LIFETEST manual, p. 1801
# days < c(179,256,262,256,255,224,225,287,319,264,237,156,270,257,242,
# 157,249,180,226,268,378,355,319,256,171,325,325,217,255,256,
# 291,323,253,206,206,237,211,229,234,209)
# status < c(1,1,1,1,1,0,1,1,1,1,0,1,1,1,1,1,1,1,1,0,
# 0,rep(1,19))
# treatment < c(rep(1,10), rep(2,10), rep(1,10), rep(2,10))
# sex < Cs(F,F,M,F,M,F,F,M,M,M,F,F,M,M,M,F,M,F,F,M,
# M,M,M,M,F,M,M,F,F,F,M,M,M,F,F,M,F,F,F,F)
# data.frame(days, status, treatment, sex)
# table(treatment, status)
# logrank(Surv(days, status), treatment) # agrees with p. 1807
# # For stratified tests the picture is puzzling.
# # survdiff(Surv(days,status) ~ treatment + strata(sex))$chisq
# # is 7.246562, which does not agree with SAS (7.1609)
# # But summary(coxph(Surv(days,status) ~ treatment + strata(sex)))
# # yields 7.16 whereas summary(coxph(Surv(days,status) ~ treatment))
# # yields 5.21 as the score test, not agreeing with SAS or logrank() (5.6485)
# ## End(Not run)