Plot data with data simulated and sensitivity of Biomarker varying

Hi everyone,

I am simulating a clinical trial to study the impact of a imperfect biomarker on the sample size and on the power. So I would like to have a graph with the sample size ( determined with calculation) and the power depending on the imperfection of the biomarker. I would like to know first how can I extract the power. Here is my code


# p_pos_bm is proportion of the population that will test positive for the biomarker
# NPV_bm is the negative predictive value of the biomarker
# PPV_bm is the positive predictive value of the biomarker
# hr is the hazard ratio of the treatment on the G+ (we suppose 1 for the G-)
# hr0 is the HR of the BK+ in the population (with or without tmt); 1 by default
# lambda is the constant risk of survival 
# npop is the size of the population 

#Input variables
p_pos_bm <- 0.4 
NPV_bm <-  seq (0.1, 1, by=0.1) 
PPV_bm <- seq(0.3, 1.0, by=0.1) 
npop <- 1e+06
lambda <- 0.067
n<- 200

##Construction of the Bk stratefy design 
simdat <- function(npop=1e+06, p_pos_bm=0.4, NPV_bm=0.9, PPV_bm=0.8, hr=0.7, hr0=1, lambda=0.067)
  dat <- data.frame(id=1:npop, bk=rbinom(npop, 1, p_pos_bm))
  dat$gen[dat$bk==1] <- rbinom(sum(dat$bk==1), 1, PPV_bm)
  dat$gen[dat$bk==0] <- rbinom(sum(dat$bk==0), 1, 1-NPV_bm)
  dat$T0 <- rexp(nrow(dat),0.067*hr0^dat$gen)
  dat$T1 <- rexp(nrow(dat),0.067*(hr0*hr)^dat$gen)

dat <- simdat()

Ddesign <- function(dat, x=1:n, a=12, b=36)
  temp <- dat[x,]
  n <- length(x)
  temp$arm <- rbinom(n,1,0.5)
  temp$treat <- ifelse(temp$arm==0,0,temp$bk)
  temp$temps <- temp$T1*temp$treat + temp$T0*(1-temp$treat) 
  # censure
  Cens <- runif(n, a, b)
  temp$survie <- pmin(Cens, temp$temps)
  temp$dc <- as.numeric(temp$temps <= Cens)

test <- Ddesign(dat=dat)
- coxph(Surv(survie,dc)~arm, data=dat)
  zfit <- cox.zph(fit)

##Simulation of essays 

n <- 500
res2 <- matrix(nrow=nsim, ncol=4)
for(i in 1:nsim)
  test <- Ddesign(dat=dat, x=seq(n*(i-1)+1,n*(i-1)+n,by=1))
  res2[i,] <- analysetrial(test)

# Extract Power ?? 

## Need of the sample size ( ici mon HR est constant, psi aussi)

samplesize <- function(HR=0.7, p=0.5, psi=0.4, alpha=0.05, beta=0.1)
  n <- (qnorm(1-alpha/2)+qnorm(1-beta))^2/(log(HR)^2*p*(1-p)*psi)

cox.zph does not return a value with a power statistic. See

  • The CR package proposes power calculation for weighted Log-Rank tests in cure rate models.
  • The NPHMC permits to calculate sample size based on proportional hazards mixture cure models.
  • The powerSurvEpi package provides power and sample size calculation for survival analysis (with a focus towards epidemiological studies).
  • Power analysis and sample size calculation for SNP association studies with time-to-event outcomes can be done using the survSNP package.

from the task view

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.