## 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)
 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447 \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}