File: spower.Rd

package info (click to toggle)
hmisc 4.2-0-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 3,332 kB
  • sloc: asm: 27,116; fortran: 606; ansic: 411; xml: 160; makefile: 2
file content (447 lines) | stat: -rw-r--r-- 16,538 bytes parent folder | download | duplicates (4)
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
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
\name{spower}
\alias{spower}
\alias{print.spower}
\alias{Quantile2}
\alias{print.Quantile2}
\alias{plot.Quantile2}
\alias{logrank}
\alias{Gompertz2}
\alias{Lognorm2}
\alias{Weibull2}
\title{
  Simulate Power of 2-Sample Test for Survival under Complex Conditions
}
\description{
  Given functions to generate random variables for survival times and
  censoring times, \code{spower} simulates the power of a user-given
  2-sample test for censored data.  By default, the logrank (Cox
  2-sample) test is used, and a \code{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 \code{survival} package). A \code{print} method prints various
  measures from these.  For composing \R functions to generate random
  survival times under complex conditions, the \code{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 (non-compliance, dropout).
  \code{Quantile2} returns a function that generates either control or
  intervention uncensored survival times subject to non-constant
  treatment effect, dropin, and dropout.  There is a \code{plot} method
  for plotting the results of \code{Quantile2}, which will aid in
  understanding the effects of the two types of non-compliance and
  non-constant treatment effects.  \code{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 equally-spaced time points.
  
  The \code{logrank} function is intended to be used with \code{spower}
  but it can be used by itself.  It returns the 1 degree of freedom
  chi-square statistic, with the hazard ratio estimate as an attribute.

  The \code{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
  (\eqn{S(t) = e^{-\alpha {t}^{\gamma}}}{S(t) = exp(-\alpha*t^\gamma)})
  which will yield
  those estimates.  It creates an \R function to evaluate survival
  probabilities from this Weibull distribution.  \code{Weibull2} is
  useful in creating functions to pass as the first argument to
  \code{Quantile2}.

  The \code{Lognorm2} and \code{Gompertz2} functions are similar to
  \code{Weibull2} except that they produce survival functions for the
  log-normal and Gompertz distributions.

  When \code{cox=TRUE} is specified to \code{spower}, the analyst may wish
  to extract the two margins of error by using the \code{print} method
  for \code{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)

\method{print}{spower}(x, conf.int=.95, \dots)

Quantile2(scontrol, hratio, 
          dropin=function(times)0, dropout=function(times)0,
          m=7500, tmax, qtmax=.001, mplot=200, pr=TRUE, \dots)

\method{print}{Quantile2}(x, \dots)

\method{plot}{Quantile2}(x, 
     what=c("survival", "hazard", "both", "drop", "hratio", "all"),
     dropsep=FALSE, lty=1:4, col=1, xlim, ylim=NULL,
     label.curves=NULL, \dots)

logrank(S, group)

Gompertz2(times, surv)
Lognorm2(times, surv)
Weibull2(times, surv)
}
\arguments{
  \item{rcontrol}{
    a function of \var{n} which returns \var{n} random uncensored
    failure times for the control group.  \code{spower} assumes that
    non-compliance (dropin) has been taken into account by this
    function.
  }
  \item{rinterv}{
    similar to \code{rcontrol} but for the intervention group
  }
  \item{rcens}{
    a function of \var{n} which returns \var{n} random censoring times.
    It is assumed that both treatment groups have the same censoring
    distribution.
  }
  \item{nc}{
    number of subjects in the control group
  }
  \item{ni}{
    number in the intervention group
  }
  \item{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.
  }
  \item{hratio}{
    a function of time which specifies the intervention:control hazard
    ratio (treatment effect)
  }
  \item{x}{
    an object of class \dQuote{Quantile2} created by \code{Quantile2},
    or of class \dQuote{spower} created by \code{spower}
  }
  \item{conf.int}{
    confidence level for determining fold-change margins of error in
    estimating the hazard ratio
  }
  \item{S}{
    a \code{Surv} object or other two-column matrix for right-censored
    survival times 
  }
  \item{group}{
    group indicators have length equal to the number of rows in \code{S}
    argument.
  }
  \item{times}{
    a vector of two times
  }
  \item{surv}{
    a vector of two survival probabilities
  }
  \item{test}{
    any function of a \code{Surv} object and a grouping variable which
    computes a chi-square for a two-sample censored data test.  The
    default is \code{logrank}.
  }
  \item{cox}{
    If true \code{TRUE} the two margins of error are available by using
    the \code{print} method for \code{spower} objects (see example
    below) and taking the maximum of the two.
  }
  \item{nsim}{
    number of simulations to perform (default=500)
  }
  \item{alpha}{
    type I error (default=.05)
  }
  \item{pr}{
    If \code{FALSE} prevents \code{spower} from printing progress notes for
    simulations. 
    If \code{FALSE} prevents \code{Quantile2} from printing \code{tmax}
    when it calculates \code{tmax}.
  }
  \item{dropin}{
    a function of time specifying the probability that a control subject
    actually is treated with the new intervention at the corresponding
    time
  }
  \item{dropout}{
    a function of time specifying the probability of an intervention
    subject dropping out to control conditions.  As a function of time,
    \code{dropout} specifies the probability that a patient is treated
    with the control therapy at time \var{t}.  \code{dropin} and
    \code{dropout} form mixing proportions for control and intervention
    hazard functions.
  }
  \item{m}{
    number of time points used for approximating functions (default is
    7500)
  }
  \item{tmax}{
    maximum time point to use in the grid of \code{m} times.  Default is
    the time such that \code{scontrol(time)} is \code{qtmax}.
  }
  \item{qtmax}{
    survival probability corresponding to the last time point used for
    approximating survival and hazard functions.  Default is 0.001.  For
    \code{qtmax} of the time for which a simulated time is needed which
    corresponds to a survival probability of less than \code{qtmax}, the
    simulated value will be \code{tmax}.
  }
  \item{mplot}{
    number of points used for approximating functions for use in
    plotting (default is 200 equally spaced points)
  }
  \item{\dots}{
    optional arguments passed to the \code{scontrol} function when it's
    evaluated by \code{Quantile2}.  Unused for \code{print.spower}.
  }
  \item{what}{
    a single character constant (may be abbreviated) specifying which
    functions to plot.  The default is \samp{"both"} meaning both
    survival and hazard functions.  Specify \code{what="drop"} to just
    plot the dropin and dropout functions, \code{what="hratio"} to plot
    the hazard ratio functions, or \samp{"all"} to make 4 separate plots
    showing all functions (6 plots if \code{dropsep=TRUE}).
  }
  \item{dropsep}{
    If \code{TRUE} makes \code{plot.Quantile2} separate pure and
    contaminated functions onto separate plots
  }
  \item{lty}{
    vector of line types
  }
  \item{col}{
    vector of colors
  }
  \item{xlim}{
    optional x-axis limits
  }
  \item{ylim}{
    optional y-axis limits
  }
  \item{label.curves}{
    optional list which is passed as the \code{opts} argument to
    \code{\link{labcurve}}.
  }
}
\value{
  \code{spower} returns the power estimate (fraction of simulated
  chi-squares greater than the alpha-critical value).  If
  \code{cox=TRUE}, \code{spower} returns an object of class
  \dQuote{spower} containing the power and various other quantities.

  \code{Quantile2} returns an \R function of class \dQuote{Quantile2}
  with attributes that drive the \code{plot} method.  The major
  attribute is a list containing several lists.  Each of these sub-lists
  contains a \code{Time} vector along with one of the following:
  survival probabilities for either treatment group and with or without
  contamination caused by non-compliance, hazard rates in a similar way,
  intervention:control hazard ratio function with and without
  contamination, and dropin and dropout functions.

  \code{logrank} returns a single chi-square statistic.

  \code{Weibull2}, \code{Lognorm2} and \code{Gompertz2} return an \R
  function with three arguments, only the first of which (the vector of
  \code{times}) is intended to be specified by the user.
}
\section{Side Effects}{
  \code{spower} prints the interation number every 10 iterations if
  \code{pr=TRUE}.
}
\author{
  Frank Harrell
  \cr
  Department of Biostatistics
  \cr
  Vanderbilt University School of Medicine
  \cr
  \email{f.harrell@vanderbilt.edu}
}
\references{
  Lakatos E (1988): Sample sizes based on the log-rank statistic in complex
  clinical trials.  Biometrics 44:229--241 (Correction 44:923).

  Cuzick J, Edwards R, Segnan N (1997): Adjusting for non-compliance and 
  contamination in randomized clinical trials. Stat in Med 16:1017--1029.

  Cook, T (2003): Methods for mid-course 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 multi-arm survival trials allowing for non-uniform accrual,
  non-proportional hazards, loss to follow-up and cross-over.  Stat in Med
  25:2521--2542.
}
\seealso{
  \code{\link{cpower}}, \code{\link{ciapower}}, \code{\link{bpower}},
  \code{\link[rms]{cph}}, \code{\link[survival]{coxph}},
  \code{\link{labcurve}}
}
\examples{
# Simulate a simple 2-arm clinical trial with exponential survival so
# we can compare power simulations of logrank-Cox test with cpower()
# Hazard ratio is constant and patients enter the study uniformly
# with follow-up ranging from 1 to 3 years
# Drop-in probability is constant at .1 and drop-out probability is
# constant at .175.  Two-year survival of control patients in absence
# of drop-in 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

\dontrun{
# 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')]
}

# Simulate a 2-arm 5-year follow-up study for which the control group's
# survival distribution is Weibull with 1-year survival of .95 and
# 3-year 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 + (b-a) * u ^ (1/(d+1))


rcens <- function(n) 1 + (5-1) * (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 follow-up 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, ...)

\dontrun{
# 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)
}
}
\keyword{htest}
\keyword{survival}
\concept{power}
\concept{study design}