library(survival) # censored exponential sampling: samp.expo <- function(n=100,rate.surv = 0.2, rate.cens = 0.1, end.of.study=Inf) { t.vec <- rexp(n,rate.surv) c.vec <- rexp(n,rate.cens) c.vec[c.vec> end.of.study] <- end.of.study delta.vec <- t.vec <= c.vec y.vec <- ifelse(delta.vec,t.vec,c.vec) return(y.vec,delta.vec) } # 1. Effect of effect-size on power. p.mat <- array(dim=c(3,100)) haz.fact <- c(1,1.3,1.5) for(j in 1:3) for(i in 1:100) { a1 <- samp.expo(n=100,rate.surv = 0.2, rate.cens = 0.1, end.of.study=20) a2 <- samp.expo(n=100,rate.surv = 0.2*haz.fact[j], rate.cens = 0.15, end.of.study=20) surv.obj <- Surv(c(a1$y.vec,a2$y.vec),c(a1$delta.vec,a2$delta.vec)) b1 <- survreg(surv.obj ~ c(rep(1,100),rep(2,100)),dist = "expo") p.mat[j,i] <- summary(b1)$table[2,4] } plot(sort(p.mat[1,])) points(sort(p.mat[2,]),pch=2) points(sort(p.mat[3,]),pch=3) apply(p.mat<0.05,1,mean) # 2. Sampling distribution of exponential log(theta.hat) lambda.hat <- rep(NA,10000) d.vec <- rep(NA,10000) for(i in 1:10000) { a <- samp.expo(n=200,rate.surv = 0.2, rate.cens = 0.1, end.of.study=20) d.vec[i] <- sum(a$delta.vec) lambda.hat[i] <- d.vec[i]/sum(a$y.vec) } hist(log(lambda.hat)) mean(log(lambda.hat)) var(log(lambda.hat)) mean(d.vec) log(0.2) 1/mean(d.vec) # 3. Analytical power computations p.1 <- 0.2 /(0.1+0.2) *(1-exp(-(0.1 +0.2) *6)) ; p.1 p.2 <- 0.25/(0.15+0.25)*(1-exp(-(0.15+0.25)*6)) ; p.2 (qnorm(.975) + qnorm(0.90))^2 / (log(.2/.25))^2 * ( 1/p.1 + 1/p.2) p.vec <- rep(NA,1000) d1.vec <- rep(NA,1000) d2.vec <- rep(NA,1000) for(i in 1:1000) { a1 <- samp.expo(n=751,rate.surv = 0.2, rate.cens = 0.1, end.of.study=6) a2 <- samp.expo(n=751,rate.surv = 0.25,rate.cens = 0.15,end.of.study=6) surv.obj <- Surv(c(a1$y.vec,a2$y.vec),c(a1$delta.vec,a2$delta.vec)) b1 <- survreg(surv.obj ~ c(rep(1,751),rep(2,751)),dist = "expo") p.vec[i] <- summary(b1)$table[2,4] d1.vec[i] <- sum(a1$delta.vec) d2.vec[i] <- sum(a2$delta.vec) } mean(p.vec < 0.05) mean(d1.vec) / 751 ; p.1 mean(d2.vec) / 751 ; p.2 a1 <- samp.expo(n=100,rate.surv = 0.2, rate.cens = 0.1, end.of.study=20) a2 <- samp.expo(n=100,rate.surv = 0.2*haz.fact[j], rate.cens = 0.15, end.of.study=20) # 4. Relative efficiency of non-parametric tests p.expo <- rep(NA,1000) p.logr <- rep(NA,1000) p.peto <- rep(NA,1000) p.cox <- rep(NA,1000) for(i in 1:1000) { n1 <- 150 n2 <- 200 a1 <- samp.expo(n=n1,rate.surv = 0.2, rate.cens = 0.1, end.of.study=20) a2 <- samp.expo(n=n2,rate.surv = 0.3, rate.cens = 0.15, end.of.study=20) surv.obj <- Surv(c(a1$y.vec,a2$y.vec),c(a1$delta.vec,a2$delta.vec)) b1 <- survreg(surv.obj ~ c(rep(1,n1),rep(2,n2)),dist = "expo") b2 <- survdiff(surv.obj ~ c(rep(1,n1),rep(2,n2)),rho=0) b3 <- survdiff(surv.obj ~ c(rep(1,n1),rep(2,n2)),rho=1) b4 <- coxph(surv.obj ~ c(rep(1,n1),rep(2,n2))) p.expo[i] <- summary(b1)$table[2,4] p.logr[i] <- 1-pchisq(b2$chisq,1) p.peto[i] <- 1-pchisq(b3$chisq,1) p.cox[i] <- summary(b4)$logtest[3] } summary(b4) ; b2 mean(p.expo < 0.05) mean(p.logr < 0.05) mean(p.peto < 0.05) mean(p.cox < 0.05) par(mfcol=c(1,3)) plot(sort(p.expo),type="l") lines(sort(p.logr),lty=2) lines(sort(p.peto),lty=3) lines(sort(p.cox), lty=4) plot(sort(p.expo),log="y",type="l") lines(sort(p.logr),lty=2) lines(sort(p.peto),lty=3) lines(sort(p.cox), lty=4) plot(200:400,sort(p.expo) [200:400],log="y",type="l",xlab="Index") lines(200:400,sort(p.logr)[200:400],lty=2) lines(200:400,sort(p.peto)[200:400],lty=3) lines(200:400,sort(p.cox) [200:400], lty=4) # 5. Power profile simulation pow.mat <- array(0,dim=c(8,11)) haz.fact <- c(1,1.05,1.1,1.15,1.2,1.25,1.3,1.35,1.4,1.45,1.5) sam.size <- c(50,100,200,350,550,800,1100,1400) for(j in 1:8) for(k in 1:11) for(i in 1:1000) { a1 <- samp.expo(n=sam.size[j],rate.surv = 0.2, rate.cens = 0.1, end.of.study=20) a2 <- samp.expo(n=sam.size[j],rate.surv = 0.2*haz.fact[k], rate.cens = 0.15, end.of.study=20) surv.obj <- Surv(c(a1$y.vec,a2$y.vec),c(a1$delta.vec,a2$delta.vec)) b2 <- survdiff(surv.obj ~ c(rep(1,sam.size[j]),rep(2,sam.size[j])),rho=0) p.value <- 1-pchisq(b2$chisq,1) if(p.value < 0.05) pow.mat[j,k] <- pow.mat[j,k]+1/1000 } contour(sam.size,haz.fact,pow.mat,levels=c(.1,.5,.8,.95,.99), xlab="Sample size",ylab="Hazard ratio")