In Knusprikus/BSSRed: Blinded Sample Size Reestimation For Event Driven Trials

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%"
)

BSSRed (Blinded Sample Size Reestimation with Evnt-Driven Data.)

The goal of BSSRed is to provide the methods published by T. Friede et al. [@friede2010blinded], which main aspect lies on reestimation procedures.

Installation

You can install the released version of BSSRed from GITHUB with:

library("devtools")
install_git(url="https://github.com/ChristophAnten/BSSRed.git",
            upgrade=FALSE)
library("BSSRed")

Study Exapmle

This is a basic example to evaluate and varify a study plan with event driven data.

For this example it is assumed, that we have some prior information from other studies on witch we base our assumptions on the hazard and drop out rates.

We assume, that 30% of the patients in the control group have an event after 24 month and 20% of the patients across both group drop out of the study within 24 month. Assuming a 30% reduction of hazard in the treatment group.

library(BSSRed)

# assumed reduction of hazard by 30% 
# defining the hazard ratio 'theta'
theta = 1-.3

# hazard rate of control group is 30% after 24 month
lambdaC <- get_hazardRate(.3,24)
lambdaC

# the treatment group has a reduction of 30% after 24 month
lambdaT <- lambdaC*theta
lambdaT

# the overall dropout proportion is 20% after 24 months
gamma <- get_hazardRate(.2,24)
gamma

With a nominal alpha niveau of 5% (two-sided), the study aims for a power of 90%. The treatment allocation is planned to be 2:1 (treatment:control). With these additional assumptions we can calculate the required number of events to achive the targetet power as follows.

# define the alpha and beta error
alpha = .05
beta = .1

# set the alternative 
alternative = "two-sided"

# define the allocation parameter
k=2

# calculating the required number of events
nEvents <- nschoenfeld(theta = theta, k=k, alpha=alpha, beta=beta, alternative = "two-sided")
nEvents

nEvents <- ceiling(nEvents)
nEvents

The recruitment is assumed to be linear with with a total of 1071 patients and an treatment allocation of 2:1 (tratment:control). The recruitment is assumend to be available at the start of each month.

# define the recruitment table
N_recruit <- set_recruitment(targetN=1071, k=2, timePoints=0:23)
head(N_recruit,5)

colSums(N_recruit[,c("T","C")])

# short overview of recruitments
barplot(names.arg =N_recruit$time,height=(t(cumsum(N_recruit[,c("T","C")]))),
        main = "Cummulative number of recruitments",col=c("grey33","grey88"))
legend("topleft",legend = c("T","C"),fill = c("grey33","grey88"),bty="n")

The study is planned to run for 48 months. With these information we can check if our recruitments will suffice to achieve the calculated number of events \code{nEvent} before 48 months.

# define the administrative censoring at the end of the study
L = 48

# calculate expected number of events
calc_expEvents(N_recruit[,1:2],theta,L,N_recruit[,3],lambdaC,gamma = gamma)

Thes are not enough events to hold the assumptions. With \code{rEstimate} we can calculate the additional recruitements we need.

# that is not enough we need
rEstimate(N=N_recruit[,1:2], tn=N_recruit[,3], theta=theta, nEvents=nEvents, L=L,
                  lambda=lambdaC, gamma = gamma, distC="exponential",distS="exponential")

So we need to recruit r sum(rEstimate(N=N_recruit[,1:2], tn=N_recruit[,3], theta=theta, nEvents=nEvents, L=L,lambda=lambdaC, gamma = gamma, distC="exponential",distS="exponential")$rec_total_number) patients or we need r sum(rEstimate(N=N_recruit[,1:2], tn=N_recruit[,3], theta=theta, nEvents=nEvents, L=L,lambda=lambdaC, gamma = gamma, distC="exponential",distS="exponential")$rec_batch) batches of recruits. The last recruitment arrives at timepoint r sum(rEstimate(N=N_recruit[,1:2], tn=N_recruit[,3], theta=theta, nEvents=nEvents, L=L,lambda=lambdaC, gamma = gamma, distC="exponential",distS="exponential")$rec_time).

AS the final step we can simulate the study and observe the impact of misspecification of the assumed hazard ratio theta on the study duration.

load("./vignette/res.RData")

# simulating the study
res <- sim_BSSRed(N=N_recruit[,1:2],tn=N_recruit[,3],M=as.matrix(N_recruit[24,1:2,drop=F]),
                  tm=25,addRecT=c(2,4,6),
                  sim.lambdaP = get_hazardRate(seq(.2,.4,length.out = 9),24), sigma=1,
                  theta=theta,gamma=gamma,kappa=1,distS="exponential",distC="exponential",nEvents=nEvents,
                  L=L,nSim=1000, par.est=FALSE, fixed.min.studytime=TRUE, seed=235711)

Show the results ...

library(plyr)
library(dplyr)
library(ggplot2)
res$results %>%
  group_by(lambda,type) %>%
  summarise_all(mean) %>%
  ungroup() %>%
  mutate(lambda = 1-exp(-lambda*24)) %>%
  ggplot(aes(x=lambda,y=study.end,col=factor(type))) +
  geom_point() +
  geom_hline(yintercept = L,lty=2,col="darkgrey") +
  geom_line() +
  ggtitle(sprintf("Mean Time to End of Study (nEvent=%i)",ceiling(nEvents))) +
  theme_bw() +
  ylab("Mean study duration") +
  xlab("24 month conifrmed progression probability")# +
# scale_y_continuous(breaks = seq(30,70,by=10),limits = c(30,70))
res$results %>%
  group_by(lambda,type) %>%
  summarise_all(mean) %>%
  ungroup() %>%
  mutate(lambda = 1-exp(-lambda*24)) %>%
  ggplot(aes(x=lambda,y=nPop,col=factor(type))) +
  geom_point() +
  geom_line() +
  ggtitle(sprintf("Mean number of Additional Patients (nEvent=%i)",ceiling(nEvents))) +
  ylab("Mean number of additional patients") +
  xlab("24 month conifrmed progression probability") +
  theme_bw()# +
# scale_y_continuous(breaks = seq(0,700,by=100),limits = c(0,700))
res$results %>%
  group_by(lambda,type) %>%
  mutate(power = BSSRed::pschoenfeld(theta = theta,
                                     nEvent = events,
                                     k = 2,
                                     alpha = .05,
                                     alternative = "two-sided")) %>%
  summarise_all(mean) %>%
  ungroup() %>%
  mutate(lambda = 1-exp(-lambda*24)) %>%
  ggplot(aes(x=lambda,y=power,col=factor(type))) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 1-.1,lty=2,col="darkgrey") +
  ggtitle(sprintf("Power at End of Study (nEvent=%i)",ceiling(nEvents))) +
  ylab("Mean power") +
  xlab("24 month conifrmed progression probability") +
  theme_bw()