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
#Definition
# 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)
return(dat)
}
set.seed(1234)
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)
return(temp[,c("bk","arm","treat","survie","dc")])
}
test <- Ddesign(dat=dat)
- coxph(Surv(survie,dc)~arm, data=dat)
zfit <- cox.zph(fit)
return(c(summary(fit)$coefficients[c(1,3,5)],zfit$table[3]))
}
##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)
}
summary(res2)
# 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)
return(n)
}