spower <- function(rcontrol, rinterv, rcens, nc, ni,
test=logrank, nsim=500, alpha=.05, pr=TRUE)
{
crit <- qchisq(1-alpha, 1)
group <- c(rep(1,nc), rep(2,ni))
nexceed <- 0
for(i in 1:nsim) {
if(pr && i %% 10 == 0)
cat(i,'')
yc <- rcontrol(nc)
yi <- rinterv(ni)
cens <- rcens(nc+ni)
y <- c(yc, yi)
S <- cbind(pmin(y,cens), 1*(y <= cens))
nexceed <- nexceed + (test(S, group) > crit)
}
nexceed/nsim
}
Quantile2 <- function(scontrol, hratio,
dropin=function(times)0,
dropout=function(times)0,
m=7500, tmax, qtmax=.001, mplot=200, pr=TRUE,
...)
{
## Solve for tmax such that scontrol(t)=qtmax
dlist <- list(...)
k <- length(dlist) && !is.null(dlist)
f <- if(k) function(x, scontrol, qt, ...) scontrol(x, ...) - qt
else function(x, scontrol, qt) scontrol(x) - qt
if(missing(tmax)) {
if(k) tmax <- uniroot(f, c(0,1e9), scontrol=scontrol, qt=qtmax, ...)$root
else tmax <- uniroot(f, c(0,1e9), scontrol=scontrol, qt=qtmax)$root
}
if(pr)
cat('\nInterval of time for evaluating functions:[0,',
format(tmax),']\n\n')
## Generate sequence of times to use in all approximations and sequence
## to use for plot method
times <- seq(0, tmax, length=m)
tim <- seq(0, tmax, length=mplot)
tinc <- times[2]
## Approximate hazard function for control group
sc <- scontrol(times, ...)
hc <- diff(-logb(sc))
hc <- c(hc, hc[m-1])/tinc ## to make length=m
## hazard function for intervention group
hr <- rep(hratio(times), length=m)
hi <- hc*hr
## hazard for control group with dropin
di <- rep(dropin(times),length=m)
hc2 <- (1-di)*hc + di*hi
## hazard for intervention group with dropout
do <- rep(dropout(times),length=m)
hi2 <- (1-do)*hi + do*hc
## survival for intervention group
si <- exp(-tinc*cumsum(hi))
## Compute contaminated survival function for control and intervention
sc2 <- if(any(di>0))exp(-tinc*cumsum(hc2))
else sc
si2 <- exp(-tinc*cumsum(hi2))
## Store all functions evaluated at shorter times vector (tim), for
## plotting
asing <- if(.R.)function(x)x
else as.single
sc.p <- asing(approx(times, sc, xout=tim)$y)
hc.p <- asing(approx(times, hc, xout=tim)$y)
sc2.p <- asing(approx(times, sc2, xout=tim)$y)
hc2.p <- asing(approx(times, hc2, xout=tim)$y)
si.p <- asing(approx(times, si, xout=tim)$y)
hi.p <- asing(approx(times, hi, xout=tim)$y)
si2.p <- asing(approx(times, si2, xout=tim)$y)
hi2.p <- asing(approx(times, hi2, xout=tim)$y)
dropin.p <- asing(approx(times, di, xout=tim)$y)
dropout.p <- asing(approx(times, do, xout=tim)$y)
hratio.p <- asing(approx(times, hr, xout=tim)$y)
hratio2.p <- hi2.p/hc2.p
tim <- asing(tim)
plot.info <- list("C Survival" =list(Time=tim,Survival=sc.p),
"I Survival" =list(Time=tim,Survival=si.p),
"C Survival w/Dropin" =list(Time=tim,Survival=sc2.p),
"I Survival w/Dropout" =list(Time=tim,Survival=si2.p),
"C Hazard" =list(Time=tim,Hazard=hc.p),
"I Hazard" =list(Time=tim,Hazard=hi.p),
"C Hazard w/Dropin" =list(Time=tim,Hazard=hc2.p),
"I Hazard w/Dropout" =list(Time=tim,Hazard=hi2.p),
"Dropin" =list(Time=tim,Probability=dropin.p),
"Dropout" =list(Time=tim,Probability=dropout.p),
"Hazard Ratio" =list(Time=tim,Ratio=hratio.p),
"Hazard Ratio w/Dropin+Dropout"=list(Time=tim,Ratio=hratio2.p))
## Create S-Plus functions for computing random failure times for
## control and intervention subject to dropin, dropout, and hratio
r <- function(n, what=c('control','intervention'),
times, csurvival, isurvival)
{
what <- match.arg(what)
approx(if(what=='control')csurvival
else isurvival,
times, xout=runif(n), rule=2)$y
}
asing <- if(.R.)function(x)x
else as.single
formals(r) <- list(n=integer(0),
what=c('control','intervention'),
times=asing(times), csurvival=asing(sc2),
isurvival=asing(si2))
structure(r, plot.info=plot.info,
dropin=any(di>0), dropout=any(do>0),
class='Quantile2')
}
print.Quantile2 <- function(x, ...)
{
attributes(x) <- NULL
print(x)
invisible()
}
plot.Quantile2 <- function(x,
what=c('survival','hazard','both','drop','hratio',
'all'), dropsep=FALSE,
lty=1:4, col=1, xlim, ylim=NULL,
label.curves=NULL, ...)
{
what <- match.arg(what)
pi <- attr(x, 'plot.info')
if(missing(xlim))
xlim <- c(0,max(pi[[1]][[1]]))
dropin <- attr(x, 'dropin')
dropout <- attr(x, 'dropout')
i <- c(1,2,
if(dropin)3,
if(dropout)4)
if(what %in% c('survival','both','all')) {
if(dropsep && (dropin|dropout)) {
labcurve(pi[1:2], pl=TRUE, lty=lty, col=col, xlim=xlim, ylim=ylim,
opts=label.curves)
labcurve(pi[i[-(1:2)]], pl=TRUE, lty=lty, col=col, xlim=xlim, ylim=ylim,
opts=label.curves)
} else
labcurve(pi[i], pl=TRUE, lty=lty, col=col, xlim=xlim, ylim=ylim,
opts=label.curves)
}
if(what %in% c('hazard','both','all')) {
if(dropsep && (dropin|dropout)) {
labcurve(pi[5:6], pl=TRUE, lty=lty, col=col, xlim=xlim, ylim=ylim,
opts=label.curves)
labcurve(pi[4+i[-(1:2)]], pl=TRUE, lty=lty, col=col, xlim=xlim, ylim=ylim,
opts=label.curves)
} else
labcurve(pi[4+i], pl=TRUE, lty=lty, col=col, xlim=xlim, ylim=ylim,
opts=label.curves)
}
if(what=='drop' || (what=='all' && (dropin | dropout))) {
i <- c(if(dropin)9,
if(dropout)10)
if(length(i)==0)
i <- 10
labcurve(pi[i], pl=TRUE, lty=lty, col=col, xlim=xlim, ylim=ylim,
opts=label.curves)
}
if(what %in% c('hratio','all')) {
i <- c(11,
if(dropin|dropout) 12)
labcurve(pi[i], pl=TRUE, lty=lty, col=col, xlim=xlim, ylim=ylim,
opts=label.curves)
}
invisible()
}
logrank <- function(S, group)
{
y <- S[,1]
event <- S[,2]
i <- order(-y)
y <- y[i]
event <- event[i]
group <- group[i]
x <- cbind(group==1, group==2, (group==1)*event, (group==2)*event)
s <- rowsumFast(x, y, FALSE)
nr1 <- cumsum(s[,1])
nr2 <- cumsum(s[,2])
d1 <- s[,3]
d2 <- s[,4]
rd <- d1+d2
rs <- nr1+nr2-rd
n <- nr1+nr2
oecum <- d1 - rd*nr1/n
vcum <- rd * rs * nr1 * nr2 / n / n / (n-1)
sum(oecum)^2 / sum(vcum,na.rm=TRUE)
}
Weibull2 <- function(times, surv)
{
z1 <- -logb(surv[1])
z2 <- -logb(surv[2])
t1 <- times[1]
t2 <- times[2]
gamma <- logb(z2/z1)/logb(t2/t1)
alpha <- z1/(t1^gamma)
g <- function(times, alpha, gamma)
{
exp(-alpha*(times^gamma))
}
formals(g) <- list(times=NULL, alpha=alpha, gamma=gamma)
g
}
## Function to fit a Gompertz survival distribution to two points
## The function is S(t) = exp[-(1/b)exp(a+bt)]
## Returns a list with components a and b, and a function for
## generating S(t) for a vector of times
Gompertz2 <- function(times, surv)
{
z1 <- logb(-logb(surv[1]))
z2 <- logb(-logb(surv[2]))
t1 <- times[1]
t2 <- times[2]
b <- (z2-z1)/(t2-t1)
a <- z1 + logb(b)-b*t1
g <- function(times, a, b) {
exp(-exp(a+b*times)/b)
}
formals(g) <- list(times=NULL, a=a, b=b)
g
}
Lognorm2 <- function(times, surv)
{
z1 <- qnorm(1-surv[1])
z2 <- qnorm(1-surv[2])
sigma <- logb(times[2]/times[1])/(z2-z1)
mu <- logb(times[1]) - sigma*z1
g <- function(times, mu, sigma) {
1 - pnorm((logb(times)-mu)/sigma)
}
formals(g) <- list(times=NULL, mu=mu, sigma=sigma)
g
}