#
# This is the user-interface definition of a Shiny web application. You can
# run the application by clicking 'Run App' above.
#
# Find out more about building applications with Shiny here:
#
# http://shiny.rstudio.com/
#
library(shiny)
#library(clinfun)
#library(Biobase)
library(knitr)
library(rmarkdown)
library(DT)
library(lubridate)
library(janitor)
library(shinyjs)
library(ggplot2)
library(shinythemes)
library(survival)
library(tidyverse)
#---one-sample survival analysis
power_one_sample_logrank <- function(medH0=30, medH1=42, n=120, time.accrual=24, time.followup=24 ,n.sim=10000, con.CI=.9)
{
# require(survminer)
hazard.rate0<-log(2)/medH0
med0<-median(rexp(10000,rate=hazard.rate0))
hazard.rate1<-log(2)/medH1
med1<-median(rexp(10000,rate=hazard.rate1))
time.total<-time.accrual+time.followup
event.rate.H0=1-1/((time.followup)*hazard.rate0)*(exp(-time.accrual*hazard.rate0)-exp(-time.total*hazard.rate0))
event.rate.H1=1-1/((time.followup)*hazard.rate1)*(exp(-time.accrual*hazard.rate1)-exp(-time.total*hazard.rate1))
p.tmp<-numeric(0)
for(i in 1:n.sim)
{
x<-rexp(n,rate=hazard.rate1)
c1=runif(n,time.accrual,time.total)
data1<-data.frame(S_time=x,C_time=c1,Obs_time=ifelse(x>c1,c1,x),censor=ifelse(x>c1,0,1))
#L1=(sum(data1$censor)-event.rate.H0*n)/sqrt(event.rate.H0*n)
s1=Surv(data1$Obs_time,data1$censor)
fit1<-survfit(s1~1,conf.int=con.CI)
#surv.med<-summary(fit1)$table[c("median" , "0.9LCL" , "0.9UCL")]
surv.med<-summary(fit1)$table[-(1:6)]
s0=1-pexp(data1$Obs_time,rate=hazard.rate0)
d11=survdiff(s1~offset(s0))$chisq
Event.E0<-sum(-log(s0))
Event.O0<-sum(data1$censor)
d1<-(Event.O0-Event.E0)/sqrt(Event.E0)
p.tmp<- rbind(p.tmp,c(Event.O0,d1,1-pchisq(d1^2,df=1),Event.E0,d11,1-pchisq(d11,df=1),surv.med))
}
p.tmp
}
#---one-sample survival analysis with two-stage design for futility
power_one_sample_logrank_interim <- function(medH0=30, medH1=42, n1=60, n2=60, time.accrual=24, time.followup=24, n.sim=10000)
{
tt.all<-numeric(0)
p.list<-list()
med.list<-c(medH1,medH0)
for(j in 1:length(med.list))
{
hazard.rate0<-log(2)/medH0
med0<-median(rexp(10000,rate=hazard.rate0))
hazard.rate1<-log(2)/med.list[j]
med1<-median(rexp(10000,rate=hazard.rate1))
n=n1+n2
time.total<-time.accrual+time.followup
event.rate.H0=1-1/((time.total-time.accrual)*hazard.rate0)*(exp(-time.accrual*hazard.rate0)-exp(-time.total*hazard.rate0))
event.rate.H1=1-1/((time.total-time.accrual)*hazard.rate1)*(exp(-time.accrual*hazard.rate1)-exp(-time.total*hazard.rate1))
p0.tmp<-p.tmp<-matrix(numeric(),n.sim,4)
for(i in 1:n.sim)
{
k1<-ceiling(n1/(time.accrual/2))
#--1st stage: accrual k1 (e.g., 3) subjects per time unit (per month)--
x1<-rexp(n1,rate=hazard.rate1)
c10<-time.accrual/2-rep((0:(ceiling(n1/k1)-1)),each=k1)
c10<-c10[1:n1]
data10<-data.frame(S_time=x1,C_time=c10,Obs_time=ifelse(x1>c10,c10,x1),censor=ifelse(x1>c10,0,1))
# s10=Surv(data10$Obs_time,data10$censor)
s00=1-pexp(data10$Obs_time,rate=hazard.rate0)
Event.E0<-sum(-log(s00))
Event.O0<-sum(data10$censor)
d1<-(Event.O0-Event.E0)/sqrt(Event.E0)
p0.tmp[i,]<- c(Event.O0,d1,1-pchisq(d1^2,df=1),Event.E0)
#p0.tmp<- rbind(p0.tmp,c(Event.O0,d1,1-pchisq(d1^2,df=1),Event.E0))
#--2nd stage at end of study: full analysis--
c1_all=runif(n,time.accrual,time.total)
x2<-rexp(n-n1,rate=hazard.rate1)
x_all<-c(x1,x2)
data_all<-data.frame(S_time=x_all,C_time=c1_all,Obs_time=ifelse(x_all>c1_all,c1_all,x_all),censor=ifelse(x_all>c1_all,0,1))
# s1_all=Surv(data_all$Obs_time,data_all$censor)
s00_all=1-pexp(data_all$Obs_time,rate=hazard.rate0)
Event.E<-sum(-log(s00_all))
Event.O<-sum(data_all$censor)
d1_all<-(Event.O-Event.E)/sqrt(Event.E)
p.tmp[i,]<-c(Event.O,d1_all,1-pchisq(d1_all^2,df=1),Event.E)
#p.tmp<- rbind(p.tmp,c(Event.O,d1_all,1-pchisq(d1_all^2,df=1),Event.E))
}
p.list[[j]]<-cbind(p.tmp,p0.tmp)
tt1<-numeric(0)
for(p1 in seq(.05,.15,by=0.01)) # p1 as p value to deterimne significance
for(p2 in seq(-0.1,.5,by=0.01)) # p2 as a threshold to stop at end of 1st stage (stop if > threshold)
{
tt0<-as.vector(table((p.tmp[,3]<p1)&(p.tmp[,2]<0),p0.tmp[,2]>p2))/n.sim
tt1<-rbind(tt1,c(p1,p2,tt0))
}
tt.all<-cbind(tt.all,tt1)
}
names(p.list)<-med.list
dimnames(tt.all)[[2]]<-c('p.2ndStage','boundary.1st','H1_NonSig_Pass','H1_Sign_Pass','H1_NonSig_Stop','H1_Sig_Stop','p.2ndStage','boundary.1st','H0_NonSig_Pass','H0_Sign_Pass','H0_NonSig_Stop','H0_Sig_Stop')
return(list(p.list = p.list, summary = tt.all, n.accrual.per.time.unit = k1))
}
# Define UI for application that draws a histogram ####
ui <- shinyUI(fluidPage(
# Application title
titlePanel("One Sample Interim Analysis"),
# Sidebar with a slider input for number of bins
sidebarLayout(
sidebarPanel(
numericInput("medH0",
"Median Survival time (months) in Control group:",
min = 1,
max = 120,
step = 1,
value = 30),
numericInput("medH1",
"Median Survival time (months) in Treatment group:",
min = 1,
max = 120,
step = 1,
value = 42),
numericInput("n.sim",
"number of simulations:",
min = 1,
max = 10000,
step = 100,
value = 2000),
numericInput("n1",
"number of patients enrolled initially:",
min = 1,
max = 1000,
step = 1,
value = 120),
numericInput("n2",
"number of additional patients enrolled:",
min = 1,
max = 1000,
step = 1,
value = 120),
numericInput("time.accrual",
"Accrual time (months):",
min = 1,
max = 120,
step = 1,
value = 24),
numericInput("time.followup",
"Follow up time (months):",
min = 1,
max = 120,
step = 1,
value = 24),
# radioButtons("con.CI", "Confidence interval:",
# choiceNames = list("95%" , "90%" ),
# choiceValues = list( 0.95 , 0.90 ) ),
hr(),
textInput("endpointShort", "Short Endpoint (eg. OS, DFS, PFS)", "DFS"),
textInput("endpointLong", "Long Endpoint (eg. Overall Survival,
disease-free survival)", "disease-free survival"),
textInput("endpointDescription", "Describe the endpoint", "the time from randomization to the earliest event which includes local, regional, or distant recurrence, new lung cancer, or death"),
hr(),
actionButton("go", "Submit"),
hr(),
radioButtons('format', 'Document format', c('Word', 'PDF'),
inline = TRUE),
downloadButton('downloadReport', "Download Report"),
),
# Show a plot of the generated distribution
mainPanel(
textOutput("writing"),
br(),
textOutput("writing2"),
br(),
plotOutput("scatterplot1"),
verbatimTextOutput("text"),
dataTableOutput(outputId = "table2") ,
verbatimTextOutput("text2"),
dataTableOutput(outputId = "table1")
)
)
))
# Define server logic ####
server <- shinyServer(function(input, output) {
# output$text <- renderText( {"Best case simulation results" })
# output$text2 <- renderText( {"All simulation results" })
#
text <- eventReactive(input$go, {"Best case simulation results" })
text2 <- eventReactive(input$go, {"All simulation results" })
output$text <- renderText({ text() })
output$text2 <- renderText({ text2() })
writing <- eventReactive(input$go, {paste("Statistical Plan:
Sample size justification: A total of ", input$n1+input$n2," patients will be enrolled in the first ", input$time.accrual," months with
an estimated accrual of ", ceiling((input$n1+input$n2)/input$time.accrual)," subjects per month. All the patients will be followed up
for at least ", input$time.followup," months. We hypothesize a median ", input$endpointLong," (", input$endpointShort,") of ", input$medH1,"
months (H1: alternative hypothesis) in the treatment group compared to ", input$medH0," months (H0: null hypothesis) in the historic control
group (Hazard ratio (HR)=", input$medH0,"/", input$medH1,"=", round(input$medH0/input$medH1,2),"). ", input$endpointShort," is defined as ", input$endpointDescription,".
An interim futility analysis will be performed after the first ", round(input$time.accrual/2)," months (", input$n1," patients enrolled).
A statistic, W, utilizes the expected events (Event.E) under the null hypothesis to compare to the observed events (Event.O) and defines
the stopping boundary where W = (Event.O-Event.E)/sqrt(Event.E). If W is higher than ", somevalues()$cut1," (i.e., the number of observed events
is higher than expected, meaning a shorter ", input$endpointShort,"), the treatment will be considered ineffective and the trial will be stopped.
Otherwise, additional ", input$n2," patients will be enrolled in the second ", round(input$time.accrual/2)," months. At end of the study, if W is
less than ", somevalues()$cut2," (i.e., the number of events is lower than expected, meaning a longer ", input$endpointShort,"), the treatment will be
considered promising. Simulation analysis show that the design has ", somevalues()$pow*100,"% power to detect the effect size of HR=", round(input$medH0/input$medH1,2),"
controlled at one-sided ", somevalues()$typI*100,"% type I error. Probability of early termination (PET) is ", somevalues()$PETH0*100,"% if the true median ", input$endpointShort,"
is ", input$medH0," months under the null hypothesis with an expected sample size of ", round( somevalues()$PETH0*input$n1+(input$n1+input$n2)*(1 - somevalues()$PETH0)),"
(", somevalues()$PETH0," * ", input$n1," + ", input$n1+input$n2," * ", 1- somevalues()$PETH0,"). When the true median ", input$endpointShort," is ", input$medH1," months (H1),
PET is ", somevalues()$PETH1*100," % and the estimated sample size is
", round( somevalues()$PETH1*input$n1+(input$n1+input$n2)*(1- somevalues()$PETH1))," (", somevalues()$PETH1," * ", input$n1," + ", input$n1+input$n2," * "
, 1- somevalues()$PETH1,"). \n")})
writing2 <- eventReactive(input$go, {paste("Data Analysis:
Kaplan-Meier curves of estimated ", input$endpointShort," will be generated. One-sided log-rank test will be used to test if the treatment
group yields a longer ", input$endpointShort," compared to the historical control (median ", input$endpointShort," of ", input$medH0," months).
Median ", input$endpointShort," times with 95% confidence intervals will also be determined. As defined in the protocol, the survival analysis
will be based on the intent-to-treat population, which will include all eligible patients enrolled in this study. The intent-to-treat survival
analyses using multivariable Cox regression model by incorporating appropriate covariates in the model. Data will be summarized overall using
descriptive statistics. For example, continuous data will be summarized with number of patients (n), mean, median, minimum, maximum, standard
deviation, coefficient of variation, and geometric mean (where applicable). Categorical data will be summarized using frequency counts and
percentages.") })
output$writing <- renderText({ writing() })
output$writing2 <- renderText({ writing2() })
get.result <- eventReactive(input$go,{
power_one_sample_logrank_interim(medH0 = input$medH0,
medH1 = input$medH1,
n1 = input$n1,
n2 = input$n2,
time.accrual = input$time.accrual,
time.followup = input$time.followup,
n.sim = input$n.sim)
})
output$table1 <- renderDataTable({
x <- as.data.frame(get.result()$summary);
return(x)} , rownames = TRUE)
somevalues <- eventReactive(input$go,{
test <- as.data.frame(get.result()$summary[,c(-1,-2)] ) %>% filter(H0_Sign_Pass < 0.05)
type1lt05 <- as.data.frame(test ) %>% filter(H0_Sign_Pass < 0.05) %>% arrange(desc(H0_Sign_Pass), desc(boundary.1st), desc(H1_Sign_Pass))
type1lt05 <- type1lt05[1, ]
# list( pow = round(type1lt05[ ,2], 2),
# typI = round(type1lt05[ ,8], 2),
# cut1 = round(type1lt05[ ,6], 2),
# cut2 = -round(sqrt(qchisq(1-type1lt05[ ,6],df=1)),2))
return(list( pow = round(type1lt05$H1_Sign_Pass, 3),
typI = round(type1lt05$H0_Sign_Pass, 3),
cut1 = round(type1lt05$boundary.1st, 2),
cut2 = -round(sqrt(qchisq(1 - type1lt05$p.2ndStage, df = 1)), 2),
PETH0 = round(type1lt05$H0_NonSig_Stop + type1lt05$H0_Sig_Stop, 2),
PETH1 = round(type1lt05$H1_NonSig_Stop + type1lt05$H1_Sig_Stop, 2))
)
})
observe({
print("DOES THIS WORK:")
print(somevalues())
print(somevalues()$pow)
print(somevalues()$typI)
print(somevalues()$cut1)
print(somevalues()$cut2)
plotdata <- get.result()$summary[,c(-1,-2)]
print(str(plotdata))
# test<-as.data.frame(plotdata ) %>% filter(H0_Sign_Pass < 0.05)
# print(length(test$H0_Sign_Pass[test$H0_Sign_Pass == max(test$H0_Sign_Pass)]))
#
# print( NROW(type1lt05 <- as.data.frame(plotdata) %>% filter(H0_Sign_Pass < 0.05) %>% arrange(H0_Sign_Pass, boundary.1st, H1_Sign_Pass) ))
# print(summary(type1lt05$H1_Sign_Pass))
# print( type1lt05)
})
output$table2 <- renderDataTable({
plotdata <- get.result()$summary[,c(-1,-2)]
test <- as.data.frame(plotdata ) %>% filter(H0_Sign_Pass < 0.05)
type1lt05 <- as.data.frame(plotdata) %>% filter(H0_Sign_Pass < 0.05) %>% arrange(desc(H0_Sign_Pass), desc(boundary.1st), desc(H1_Sign_Pass))
type1lt05 <- head(type1lt05, length(test$H0_Sign_Pass[test$H0_Sign_Pass == max(test$H0_Sign_Pass)]))
return(type1lt05)} , rownames = TRUE)
output$scatterplot1 <- renderPlot({
x <- as.data.frame(get.result()$summary[, c(-1,-2)]);
print(str(x))
ggplot(x, aes(x = H1_Sign_Pass, y = H0_Sign_Pass)) + geom_point() +
geom_hline(yintercept = 0.05, linetype = "dashed", color = "red")
# Change the point size, and shape
#geom_point(size=2, shape=23)
})
#---output report----
output$downloadReport <- downloadHandler(
filename = function() {
paste('InterimReport', sep = '.',
switch(input$format, PDF = 'pdf', HTML = 'html', Word = 'docx'))
},
content = function(file) {
out <- rmarkdown::render(input = 'F:\\myGitRepo\\clinical_trial_design\\clinical_trial_design\\OneSampleInterimAnalysis\\InterimReport.Rmd',
output_format =
switch(input$format,
PDF = rmarkdown::pdf_document(),
HTML = rmarkdown::html_document(),
Word = rmarkdown::word_document()
),
params = list(set_title = input$project_title, set_author = input$author_input)
)
#file.rename(out, file)
file.copy(out,file)
}
)
})
# END O'script
shinyApp(ui=ui,server=server)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.