1spower <- function(rcontrol, rinterv, rcens, nc, ni, 2 test=logrank, cox=FALSE, nsim=500, alpha=.05, pr=TRUE) 3{ 4 crit <- qchisq(1-alpha, 1) 5 group <- c(rep(1,nc), rep(2,ni)) 6 nexceed <- 0 7 if(cox) 8 { 9 require(survival) 10 beta <- numeric(nsim) 11 } 12 13 maxfail <- 0; maxcens <- 0 14 for(i in 1:nsim) { 15 if(pr && i %% 10 == 0) cat(i,'\r') 16 17 yc <- rcontrol(nc) 18 yi <- rinterv(ni) 19 cens <- rcens(nc+ni) 20 y <- c(yc, yi) 21 maxfail <- max(maxfail, max(y)) 22 maxcens <- max(maxcens, max(cens)) 23 S <- cbind(pmin(y,cens), 1*(y <= cens)) 24 nexceed <- nexceed + (test(S, group) > crit) 25 if(cox) 26 { 27 fit <- coxph.fit(as.matrix(group), S, strata=NULL, 28 offset=NULL, init=NULL, 29 control=coxph.control(iter.max=10, eps=.0001), 30 method="efron", rownames=NULL) 31 beta[i] <- fit$coefficients 32 } 33 } 34 cat('\n') 35 if(maxfail < 0.99*maxcens) 36 stop(paste('Censoring time distribution defined at later times than\nsurvival time distribution. There will likely be uncensored failure times\nstacked at the maximum allowed survival time.\nMaximum simulated failure time:', max(y),'\nMaximum simulated censoring time:', max(cens))) 37 38 power <- nexceed/nsim 39 if(cox) structure(list(power=power, betas=beta, nc=nc, ni=ni, 40 alpha=alpha, nsim=nsim), class='spower') else power 41} 42 43print.spower <- function(x, conf.int=.95, ...) 44 { 45 b <- x$betas 46 hr <- exp(b) 47 pp <- (1+conf.int)/2 48 cl <- quantile(hr, c((1-conf.int)/2, pp)) 49 meanbeta <- mean(b) 50 medbeta <- median(b) 51 hrmean <- exp(meanbeta) 52 hrmed <- exp(medbeta) 53 moehi <- cl[2]/hrmed 54 moelo <- hrmed/cl[1] 55 g <- function(w) round(w, 4) 56 mmoe <- max(moehi, moelo) 57 cat('\nTwo-Group Event Time Comparison Simulation\n\n', 58 x$nsim,' simulations\talpha: ', x$alpha, '\tpower: ', x$power, 59 '\t', conf.int, ' confidence interval\n', 60 '\nHazard ratio from mean beta : ', g(hrmean), 61 '\nHazard ratio from median beta : ', g(hrmed), 62 '\nStandard error of log hazard ratio: ', g(sd(b)), 63 '\nConfidence limits for hazard ratio: ', g(cl[1]), ', ', g(cl[2]), 64 '\nFold-change margin of error high : ', g(moehi), 65 '\t(upper CL/median HR)', 66 '\nFold-change margin of error low : ', g(moelo), 67 '\t(median HR/lower CL)', 68 '\nMax fold-change margin of error : ', g(mmoe),'\n\n') 69 70 cat('The fold change margin of error of', g(mmoe), 71 'represents the margin of error\n', 72 'the study is likely to achieve in estimating the intervention:control\n', 73 'hazard ratio. It is the ratio of a', conf.int, 'confidence limit on the\n', 74 'hazard ratio to the median hazard ratio obtained over the', x$nsim, 'simulations.\n', 75 'The confidence limit was obtained by computing the', pp, 'quantile of the\n', 76 x$nsim, 'observed hazard ratios. The standard error is the standard deviation\n', 77 'of the', x$nsim, 'simulated log hazard ratios.\n\n') 78 79 res <- c(cl, hrmean, hrmed, sd(b), moelo, moehi, x$power) 80 names(res) <- c('CLlower','CLupper','HRmean','HRmedian','SE', 81 'MOElower','MOEupper','Power') 82 invisible(res) 83 } 84 85Quantile2 <- function(scontrol, hratio, 86 dropin=function(times)0, 87 dropout=function(times)0, 88 m=7500, tmax, qtmax=.001, mplot=200, pr=TRUE, 89 ...) 90{ 91 ## Solve for tmax such that scontrol(t)=qtmax 92 dlist <- list(...) 93 k <- length(dlist) && !is.null(dlist) 94 f <- if(k) function(x, scontrol, qt, ...) scontrol(x, ...) - qt 95 else function(x, scontrol, qt) scontrol(x) - qt 96 97 if(missing(tmax)) { 98 if(k) tmax <- uniroot(f, c(0,1e9), scontrol=scontrol, qt=qtmax, ...)$root 99 else tmax <- uniroot(f, c(0,1e9), scontrol=scontrol, qt=qtmax)$root 100 } 101 102 if(pr) 103 cat('\nInterval of time for evaluating functions:[0,', 104 format(tmax),']\n\n') 105 106 ## Generate sequence of times to use in all approximations and sequence 107 ## to use for plot method 108 109 times <- seq(0, tmax, length=m) 110 tim <- seq(0, tmax, length=mplot) 111 tinc <- times[2] 112 113 ## Approximate hazard function for control group 114 sc <- scontrol(times, ...) 115 hc <- diff(-logb(sc)) 116 hc <- c(hc, hc[m-1])/tinc ## to make length=m 117 118 ## hazard function for intervention group 119 hr <- rep(hratio(times), length=m) 120 hi <- hc*hr 121 122 ## hazard for control group with dropin 123 di <- rep(dropin(times),length=m) 124 hc2 <- (1-di)*hc + di*hi 125 126 ## hazard for intervention group with dropout 127 do <- rep(dropout(times),length=m) 128 hi2 <- (1-do)*hi + do*hc 129 130 ## survival for intervention group 131 si <- exp(-tinc*cumsum(hi)) 132 133 ## Compute contaminated survival function for control and intervention 134 sc2 <- if(any(di>0))exp(-tinc*cumsum(hc2)) 135 else sc 136 137 si2 <- exp(-tinc*cumsum(hi2)) 138 139 140 ## Store all functions evaluated at shorter times vector (tim), for 141 ## plotting 142 asing <- if(.R.)function(x)x 143 else as.single 144 145 sc.p <- asing(approx(times, sc, xout=tim)$y) 146 hc.p <- asing(approx(times, hc, xout=tim)$y) 147 sc2.p <- asing(approx(times, sc2, xout=tim)$y) 148 hc2.p <- asing(approx(times, hc2, xout=tim)$y) 149 150 si.p <- asing(approx(times, si, xout=tim)$y) 151 hi.p <- asing(approx(times, hi, xout=tim)$y) 152 si2.p <- asing(approx(times, si2, xout=tim)$y) 153 hi2.p <- asing(approx(times, hi2, xout=tim)$y) 154 155 dropin.p <- asing(approx(times, di, xout=tim)$y) 156 dropout.p <- asing(approx(times, do, xout=tim)$y) 157 hratio.p <- asing(approx(times, hr, xout=tim)$y) 158 hratio2.p <- hi2.p/hc2.p 159 160 tim <- asing(tim) 161 162 plot.info <- list("C Survival" =list(Time=tim,Survival=sc.p), 163 "I Survival" =list(Time=tim,Survival=si.p), 164 "C Survival w/Dropin" =list(Time=tim,Survival=sc2.p), 165 "I Survival w/Dropout" =list(Time=tim,Survival=si2.p), 166 "C Hazard" =list(Time=tim,Hazard=hc.p), 167 "I Hazard" =list(Time=tim,Hazard=hi.p), 168 "C Hazard w/Dropin" =list(Time=tim,Hazard=hc2.p), 169 "I Hazard w/Dropout" =list(Time=tim,Hazard=hi2.p), 170 "Dropin" =list(Time=tim,Probability=dropin.p), 171 "Dropout" =list(Time=tim,Probability=dropout.p), 172 "Hazard Ratio" =list(Time=tim,Ratio=hratio.p), 173 "Hazard Ratio w/Dropin+Dropout"=list(Time=tim,Ratio=hratio2.p)) 174 175 ## Create S-Plus functions for computing random failure times for 176 ## control and intervention subject to dropin, dropout, and hratio 177 178 r <- function(n, what=c('control','intervention'), 179 times, csurvival, isurvival) 180 { 181 what <- match.arg(what) 182 approx(if(what=='control')csurvival 183 else isurvival, 184 times, xout=runif(n), rule=2)$y 185 } 186 187 asing <- if(.R.)function(x)x 188 else as.single 189 190 formals(r) <- list(n=integer(0), 191 what=c('control','intervention'), 192 times=asing(times), csurvival=asing(sc2), 193 isurvival=asing(si2)) 194 195 structure(r, plot.info=plot.info, 196 dropin=any(di>0), dropout=any(do>0), 197 class='Quantile2') 198} 199 200 201print.Quantile2 <- function(x, ...) 202{ 203 attributes(x) <- NULL 204 print(x) 205 invisible() 206} 207 208plot.Quantile2 <- function(x, 209 what=c('survival','hazard','both','drop','hratio', 210 'all'), dropsep=FALSE, 211 lty=1:4, col=1, xlim, ylim=NULL, 212 label.curves=NULL, ...) 213{ 214 what <- match.arg(what) 215 pi <- attr(x, 'plot.info') 216 if(missing(xlim)) 217 xlim <- c(0,max(pi[[1]][[1]])) 218 219 dropin <- attr(x, 'dropin') 220 dropout <- attr(x, 'dropout') 221 i <- c(1,2, 222 if(dropin)3, 223 if(dropout)4) 224 225 if(what %in% c('survival','both','all')) { 226 if(dropsep && (dropin|dropout)) { 227 labcurve(pi[1:2], pl=TRUE, lty=lty, col.=col, xlim=xlim, ylim=ylim, 228 opts=label.curves) 229 labcurve(pi[i[-(1:2)]], pl=TRUE, lty=lty, col.=col, xlim=xlim, ylim=ylim, 230 opts=label.curves) 231 } else 232 labcurve(pi[i], pl=TRUE, lty=lty, col.=col, xlim=xlim, ylim=ylim, 233 opts=label.curves) 234 } 235 236 if(what %in% c('hazard','both','all')) { 237 if(dropsep && (dropin|dropout)) { 238 labcurve(pi[5:6], pl=TRUE, lty=lty, col.=col, xlim=xlim, ylim=ylim, 239 opts=label.curves) 240 labcurve(pi[4+i[-(1:2)]], pl=TRUE, lty=lty, col.=col, xlim=xlim, 241 ylim=ylim, opts=label.curves) 242 } else 243 labcurve(pi[4+i], pl=TRUE, lty=lty, col.=col, xlim=xlim, ylim=ylim, 244 opts=label.curves) 245 } 246 247 if(what=='drop' || (what=='all' && (dropin | dropout))) { 248 i <- c(if(dropin)9, 249 if(dropout)10) 250 251 if(length(i)==0) 252 i <- 10 253 254 labcurve(pi[i], pl=TRUE, lty=lty, col.=col, xlim=xlim, ylim=ylim, 255 opts=label.curves) 256 } 257 258 if(what %in% c('hratio','all')) { 259 i <- c(11, 260 if(dropin|dropout) 12) 261 262 labcurve(pi[i], pl=TRUE, lty=lty, col.=col, xlim=xlim, ylim=ylim, 263 opts=label.curves) 264 } 265 266 invisible() 267} 268 269logrank <- function(S, group) 270{ 271 group <- as.factor(group) 272 i <- is.na(S) | is.na(group) 273 if(any(i)) 274 { 275 i <- !i 276 S <- S[i,,drop=FALSE] 277 group <- group[i] 278 } 279 group <- as.integer(group) 280 y <- S[,1] 281 event <- S[,2] 282 i <- order(-y) 283 y <- y[i] 284 event <- event[i] 285 group <- group[i] 286 x <- cbind(group==1, group==2, (group==1)*event, (group==2)*event) 287 if(TRUE) 288 { 289 s <- rowsum(x, y, FALSE) 290 nr1 <- cumsum(s[,1]) 291 nr2 <- cumsum(s[,2]) 292 d1 <- s[,3] 293 d2 <- s[,4] 294 rd <- d1+d2 295 rs <- nr1+nr2-rd 296 n <- nr1+nr2 297 oecum <- d1 - rd*nr1/n 298 vcum <- rd * rs * nr1 * nr2 / n / n / (n-1) 299 chisq <- sum(oecum)^2 / sum(vcum,na.rm=TRUE) 300 hr <- sum(d1*(nr1-d1)/n)/sum(d2*(nr2-d2)/n) 301 } 302 else 303 { # non-working code; trying to get stratification to work 304 OE <- v <- hrn <- hrd <- 0 305 for(strat in unique(strata)) 306 { 307 j <- strata==strat 308 s <- rowsum(x[j,], y[j], FALSE) 309 nr1 <- cumsum(s[,1]) 310 nr2 <- cumsum(s[,2]) 311 d1 <- s[,3] 312 d2 <- s[,4] 313 rd <- d1+d2 314 rs <- nr1+nr2-rd 315 n <- nr1+nr2 316 oecum <- d1 - rd*nr1/n 317 vcum <- rd * rs * nr1 * nr2 / n / n / (n-1) 318 OE <- OE + sum(oecum) 319 v <- v + sum(vcum, na.rm=TRUE) 320 hrn <- hrn + sum(d1*(nr1-d1)/n) 321 hrd <- hrd + sum(d2*(nr2-d2)/n) 322 } 323 chisq <- OE^2 / v 324 hr <- hrn/hrd 325 } 326 structure(chisq, hr=hr) 327} 328 329 330Weibull2 <- function(times, surv) 331{ 332 z1 <- -logb(surv[1]) 333 z2 <- -logb(surv[2]) 334 t1 <- times[1] 335 t2 <- times[2] 336 gamma <- logb(z2/z1)/logb(t2/t1) 337 alpha <- z1/(t1^gamma) 338 339 g <- function(times, alpha, gamma) 340 { 341 exp(-alpha*(times^gamma)) 342 } 343 344 formals(g) <- list(times=NULL, alpha=alpha, gamma=gamma) 345 g 346} 347 348 349## Function to fit a Gompertz survival distribution to two points 350## The function is S(t) = exp[-(1/b)exp(a+bt)] 351## Returns a list with components a and b, and a function for 352## generating S(t) for a vector of times 353Gompertz2 <- function(times, surv) 354{ 355 z1 <- logb(-logb(surv[1])) 356 z2 <- logb(-logb(surv[2])) 357 t1 <- times[1] 358 t2 <- times[2] 359 b <- (z2-z1)/(t2-t1) 360 a <- z1 + logb(b)-b*t1 361 362 g <- function(times, a, b) { 363 exp(-exp(a+b*times)/b) 364 } 365 366 formals(g) <- list(times=NULL, a=a, b=b) 367 g 368} 369 370 371Lognorm2 <- function(times, surv) 372{ 373 z1 <- qnorm(1-surv[1]) 374 z2 <- qnorm(1-surv[2]) 375 sigma <- logb(times[2]/times[1])/(z2-z1) 376 mu <- logb(times[1]) - sigma*z1 377 378 g <- function(times, mu, sigma) { 379 1 - pnorm((logb(times)-mu)/sigma) 380 } 381 382 formals(g) <- list(times=NULL, mu=mu, sigma=sigma) 383 g 384} 385