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(α*t^γ))
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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24  spower(rcontrol, rinterv, rcens, nc, ni,
test=logrank, cox=FALSE, nsim=500, alpha=0.05, pr=TRUE)
## S3 method for class 'spower'
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, ...)
## S3 method for class 'Quantile2'
print(x, ...)
## S3 method for class 'Quantile2'
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)

rcontrol 
a function of n which returns n random uncensored
failure times for the control group. 
rinterv 
similar to 
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 
conf.int 
confidence level for determining foldchange margins of error in estimating the hazard ratio 
S 
a 
group 
group indicators have length equal to the number of rows in 
times 
a vector of two times 
surv 
a vector of two survival probabilities 
test 
any function of a 
cox 
If true 
nsim 
number of simulations to perform (default=500) 
alpha 
type I error (default=.05) 
pr 
If 
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,

m 
number of time points used for approximating functions (default is 7500) 
tmax 
maximum time point to use in the grid of 
qtmax 
survival probability corresponding to the last time point used for
approximating survival and hazard functions. Default is 0.001. For

mplot 
number of points used for approximating functions for use in plotting (default is 200 equally spaced points) 
... 
optional arguments passed to the 
what 
a single character constant (may be abbreviated) specifying which
functions to plot. The default is "both" meaning both
survival and hazard functions. Specify 
dropsep 
If 
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 
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.
spower
prints the interation number every 10 iterations if
pr=TRUE
.
Frank Harrell
Department of Biostatistics
Vanderbilt University School of Medicine
f.harrell@vanderbilt.edu
Lakatos E (1988): Sample sizes based on the logrank statistic in complex clinical trials. Biometrics 44:229–241 (Correction 44:923).
Cuzick J, Edwards R, Segnan N (1997): Adjusting for noncompliance and contamination in randomized clinical trials. Stat in Med 16:1017–1029.
Cook, T (2003): Methods for midcourse corrections in clinical trials with survival outcomes. Stat in Med 22:3431–3447.
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:2521–2542.
cpower
, ciapower
, bpower
,
cph
, coxph
,
labcurve
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164  # 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)

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
All documentation is copyright its authors; we didn't write any of that.