GoNoGo
#read historical placebo data
# installed.packages("metafor",repos="http://R-Forge.R-project.org")
# installed.packages("ggplot2",repos="http://R-Forge.R-project.org")
# install.packages("mvtnorm", repos="http://R-Forge.R-project.org")
library(metafor)## Loading required package: Matrix## Loading 'metafor' package (version 2.4-0). For an overview
## and introduction to the package please type: help(metafor).library(ggplot2)
library(mvtnorm)R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
Old Code
Function and Go Plot
#######################Bayesian decision rule, codes written by Wei Zhong
#######################GO: Prob(TRUE diff > B | OBSERVED diff) >= pupper
#######################NO-GO: Prob(TRUE diff > B | OBSERVED diff) < plower
#######################Gray Zone: plower < Prob(TRUE diff > B | OBSERVED diff) < puppe
#######################This function serves to obtain operating characteristics of a 2-arm trial
#######################Completely noninformatvie prior
#######################Assume equal known variance for both arms -- sigma^2
#######################B: delta0 cut off
require(tidyverse)## Loading required package: tidyverse## ── Attaching packages ───────────────────────────────────────────────────────────── tidyverse 1.3.0 ──## ✓ tibble 2.1.3 ✓ dplyr 0.8.5
## ✓ tidyr 1.0.2 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0
## ✓ purrr 0.3.3## ── Conflicts ──────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x tidyr::expand() masks Matrix::expand()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x tidyr::pack() masks Matrix::pack()
## x tidyr::unpack() masks Matrix::unpack()#require(scales)
Bayes.normal<-function(delta0=c(8, -4),sigma=c(10,6),pupper=c(0.7, 0.7),plower=c(0.2,0.2),
Npergroup=c(12,12), nsim=1000, delta_true_1=c(6, 10, 14, 18, 22, 26),
delta_true_2=c(-12.8, -10.8, -8.8, -6.8, -4.8)){
decision<-function(x) { ###x for true delta
#obs_diffvec<- rnorm(nsim, mean=x, sd=sqrt(2/Npergroup)*sigma)
obs_diffvec<- matrix(rnorm(nsim*2, mean=x, sd=sqrt(1/Npergroup)*sigma), byrow = T, ncol = 2)
dec_govec<-dec_nogovec<-dec_greyvec<-rep(NA,nsim)
for (i in 1:nsim) {
#postprob<-1-pnorm(delta0, mean=obs_diffvec[i], sd=sqrt(2/Npergroup)*sigma)
#for crossover:
postprob<-1-pnorm(delta0, mean=obs_diffvec[i, ], sd=sqrt(1/Npergroup)*sigma)
#futility
#postprob<-pnorm(delta0, mean=obs_diffvec[i], sd=sqrt(1/Npergroup)*sigma)
dec_govec[i]<- prod(ifelse( postprob >= pupper, 1, 0))
dec_nogovec[i]<- prod(ifelse( postprob < plower, 1,0))
dec_greyvec[i]<- 1- max(dec_govec[i], dec_nogovec[i])
}
prob_dec_go<-sum(dec_govec)/nsim
prob_dec_nogo<-sum(dec_nogovec)/nsim
prob_dec_grey<-sum(dec_greyvec)/nsim
return(c(prob_dec_go,prob_dec_nogo, prob_dec_grey))
}
TrueEff_mat = expand.grid(endpoint1 = delta_true_1, endpoint2 = delta_true_2)
result.mat<-t(apply(TrueEff_mat, FUN = decision, MARGIN = 1))
colnames(result.mat) = c("Go_prob", "Nogo_prob", "Grey_prob")
result<-cbind(TrueEff_mat, result.mat)
result = as_tibble(result)
result = result %>%
mutate(Go_prob_text = sprintf("%.1f %%", 100*Go_prob),
Nogo_prob_text = sprintf("%.1f %%", 100*Nogo_prob),
Grey_prob_text = sprintf("%.1f %%", 100*Grey_prob))
return(result)
} ##end of the function Bayes.normal
result = Bayes.normal(delta0=c(20, 4),sigma=c(10,8),pupper=c(0.7, 0.7),plower=c(0.1,0.1),
Npergroup=c(18,18), nsim=1000, delta_true_1=c(0,4,8,12,16,20,24,30),
delta_true_2=c(0,2,3,5,7,8,10,12))
Go_plot = ggplot(result, aes(endpoint1, endpoint2)) +
geom_tile(aes(fill = Go_prob)) +
scale_fill_gradient(low = "white", high = "green", guide = guide_colorbar(title = "Probility of Go"), limits = c(0, 1), breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) +
geom_text(aes(label = Go_prob_text)) +
theme_classic()+
scale_x_continuous(breaks=c(0,4,8,12,16,20,24,30))+
scale_y_continuous(breaks=c(0,2,3,5,7,8,10,12))+
ggtitle("Probability of Go")+
labs(x="MWT", y="ESS")
Go_plot
Nogo Plot
Nogo_plot = ggplot(result, aes(endpoint1, endpoint2)) +
geom_tile(aes(fill = Nogo_prob)) +
scale_fill_gradient(low = "white", high = "red", guide = guide_colorbar(title = "Probility of No-Go"), limits = c(0, 1), breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) +
geom_text(aes(label = Nogo_prob_text)) +
theme_classic()+
scale_x_continuous(breaks=c(0,4,8,12,16,20,24,30))+
scale_y_continuous(breaks=c(0,2,3,5,7,8,10,12))+
ggtitle("Probability of No-Go")+
labs(x="MWT", y="ESS")
Nogo_plot
Grey Plot
Grey_plot = ggplot(result, aes(endpoint1, endpoint2)) +
geom_tile(aes(fill = Grey_prob)) +
scale_fill_gradient(low = "white", high = "grey", guide = guide_colorbar(title = "Probility of Grey Zone"), limits = c(0, 1), breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) +
geom_text(aes(label = Grey_prob_text)) +
theme_classic()+
scale_x_continuous(breaks=c(0,4,8,12,16,20,24,30))+
scale_y_continuous(breaks=c(0,2,3,5,7,8,10,12))+
ggtitle("Probability of Grey Zone")+
labs(x="MWT", y="ESS")
Grey_plot
New Function
# install.packages("metafor")
# installed.packages("ggplot2",repos="http://R-Forge.R-project.org")
# install.packages("mvtnorm", repos="http://R-Forge.R-project.org")
# require(tidyverse)
# require(metafor)
# require(ggplot2)
# require(mvtnorm)
Bayes.normal.combo <- function(delta0MWT=13, sigmaMWT=11, MUpMWT=3,pupperMWT=0.7, plowerMWT=0.2,
delta_true_rangeMWT=c(3,33), nsegMWT=20,
muPRt_MWT=0, sigmaPRt_MWT=Inf, muPRp_MWT=0, sigmaPRp_MWT=Inf,
delta0ESS=-6, sigmaESS=7, MUpESS=0, pupperESS=0.7, plowerESS=0.1,
muPRt_ESS=0, sigmaPRt_ESS=Inf, muPRp_ESS=0, sigmaPRp_ESS=Inf,
TrueESS=-2,
r.seed=11235,nsim=1000, Nt=12, Np=12, rho,
myylab="Probability Go/No-Go/Further Assessment (%)",
myxlab="True Mean Difference (min)"){
set.seed(r.seed)
mycor <- matrix(rho, nrow=2, ncol=2)
diag(mycor) <- c(sigmaMWT, sigmaESS)
decision <- function(xMWT, xESS) { ###x for true delta for MWT
obs_plb <- rmvnorm(nsim, mean = c(MUpMWT, MUpESS), sigma = mycor*2/Np)
obs_trt <- rmvnorm(nsim, mean = c(MUpMWT+xMWT, MUpESS+xESS), sigma = mycor*2/Nt)
#obs_plb_MWT <- rnorm(nsim, mean=MUpMWT, sd=sqrt(1/Np)*sigmaMWT)
#obs_trt_MWT <- rnorm(nsim, mean=MUpMWT+xMWT, sd=sqrt(1/Nt)*sigmaMWT)
for(i in 1:length(obs_trt[,1])){if(obs_trt[i,1] > 40){obs_trt[i,1]=40}}
#obs_plb_ESS <- rnorm(nsim, mean=MUpESS, sd=sqrt(1/Np)*sigmaESS)
#obs_trt_ESS <- rnorm(nsim, mean=MUpESS+xESS, sd=sqrt(1/Nt)*sigmaESS)
dec_govec <- dec_nogovec <- dec_greyvec<-rep(NA,nsim)
post_meanMWT <- (obs_trt[,1]/((1/Nt)*sigmaMWT^2)+muPRt_MWT/sigmaPRt_MWT^2)/(1/sigmaPRt_MWT^2+1/((1/Nt)*sigmaMWT^2))-
(obs_plb[,1]/((1/Np)*sigmaMWT^2)+muPRp_MWT/sigmaPRp_MWT^2)/(1/sigmaPRp_MWT^2+1/((1/Np)*sigmaMWT^2))
post_sigma_MWT <- sqrt(1/(1/sigmaPRt_MWT^2+1/((1/Nt)*sigmaMWT^2))+1/(1/sigmaPRp_MWT^2+1/((1/Np)*sigmaMWT^2)))
post_meanESS <- (obs_trt[,2]/((1/Nt)*sigmaESS^2)+muPRt_ESS/sigmaPRt_ESS^2)/(1/sigmaPRt_ESS^2+1/((1/Nt)*sigmaESS^2))-
(obs_plb[,2]/((1/Np)*sigmaESS^2)+muPRp_ESS/sigmaPRp_ESS^2)/(1/sigmaPRp_ESS^2+1/((1/Np)*sigmaESS^2))
post_sigma_ESS <- sqrt(1/(1/sigmaPRt_ESS^2+1/((1/Nt)*sigmaESS^2))+1/(1/sigmaPRp_ESS^2+1/((1/Np)*sigmaESS^2)))
postprobMWT <- 1-pnorm(delta0MWT, mean=post_meanMWT, sd=post_sigma_MWT) #for MWT, superiority
postprobESS <- 1-pnorm(delta0ESS, mean=post_meanESS, sd=post_sigma_ESS) #for ESS, noninferiority
dec_govec <- ifelse( postprobMWT >= pupperMWT & postprobESS >=pupperESS, 1, 0) #strict go: Go/Go
dec_nogovec <- ifelse( postprobMWT < plowerMWT & postprobESS < plowerESS, 1,0) #strict no-go: No-Go/No-Go
dec_greyvec <- 1- pmax(dec_govec, dec_nogovec)
prob_dec_go <- sum(dec_govec)/nsim
prob_dec_nogo <- sum(dec_nogovec)/nsim
prob_dec_grey <- sum(dec_greyvec)/nsim
return(c(prob_dec_go,prob_dec_nogo, prob_dec_grey))
}
#TrueEffMWT <- seq(delta_true_rangeMWT[1],delta_true_rangeMWT[2], (delta_true_rangeMWT[2]-delta_true_rangeMWT[1])/nsegMWT)
#TrueEffESS <- rep(TrueESS, length(TrueEffMWT))
TrueEff_mat = expand.grid(endpoint1 = delta_true_rangeMWT, endpoint2 = TrueESS)
# result.mat <- mapply(FUN=decision, TrueEffMWT, TrueEffESS)
result.mat <- mapply(FUN=decision, TrueEff_mat$endpoint1, TrueEff_mat$endpoint2)
return(
list(result=data.frame(TrueEffMWT=TrueEff_mat$endpoint1, TrueEffESS=TrueEff_mat$endpoint2, ProbNoGo = result.mat[2,], ProbGray=result.mat[3,],ProbGo = result.mat[1,] ),
Nt=Nt, Np=Np)
)
}
delta_true_rangeMWT = c(0,4,8,12,16,20,24,30)
TrueESS = c(0,2,3,5,7,8,10,12)
output <- Bayes.normal.combo(delta0MWT=13, sigmaMWT=11, MUpMWT=3,pupperMWT=0.7, plowerMWT=0.1,
delta_true_rangeMWT=delta_true_rangeMWT, nsegMWT=20,
muPRt_MWT=0, sigmaPRt_MWT=Inf,
muPRp_MWT=2, sigmaPRp_MWT=1, ###using informative prior on placebo
TrueESS=TrueESS,
delta0ESS=-6, sigmaESS=7, pupperESS=0.7, plowerESS=0.1,
muPRt_ESS=0, sigmaPRt_ESS=Inf,
muPRp_ESS=1, sigmaPRp_ESS=1, #using informative prior on placebo
Nt=12, Np=12, rho=0.9)
result1 = output$result
TrueEffMWT = output$TrueEffMWT
TrueEffESS = output$TrueEffESS
colnames(result1) = c("endpoint1", "endpoint2", "Nogo_prob", "Grey_prob","Go_prob")
result1 = result1 %>%
mutate(Go_prob_text = sprintf("%.1f %%", 100*Go_prob),
Nogo_prob_text = sprintf("%.1f %%", 100*Nogo_prob),
Grey_prob_text = sprintf("%.1f %%", 100*Grey_prob))
Go_plot = ggplot(result1, aes(endpoint1, endpoint2)) +
geom_tile(aes(fill = Go_prob)) +
scale_fill_gradient(low = "white", high = "green", guide = guide_colorbar(title = "Probility of Go"), limits = c(0, 1), breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) +
geom_text(aes(label = Go_prob_text)) +
theme_classic()+
scale_x_continuous(breaks=delta_true_rangeMWT)+
scale_y_continuous(breaks=TrueESS)+
ggtitle("Probability of Go")+
labs(x="MWT", y="ESS")
Go_plot
## Compare of old and new Dataframe
result## # A tibble: 64 x 8
## endpoint1 endpoint2 Go_prob Nogo_prob Grey_prob Go_prob_text Nogo_prob_text
## <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 0 0 0 0.807 0.193 0.0 % 80.7 %
## 2 4 0 0 0.802 0.198 0.0 % 80.2 %
## 3 8 0 0 0.8 0.2 0.0 % 80.0 %
## 4 12 0 0 0.792 0.208 0.0 % 79.2 %
## 5 16 0 0 0.535 0.465 0.0 % 53.5 %
## 6 20 0 0 0.073 0.927 0.0 % 7.3 %
## 7 24 0 0 0.001 0.999 0.0 % 0.1 %
## 8 30 0 0.002 0 0.998 0.2 % 0.0 %
## 9 0 2 0 0.403 0.597 0.0 % 40.3 %
## 10 4 2 0 0.418 0.582 0.0 % 41.8 %
## # … with 54 more rows, and 1 more variable: Grey_prob_text <chr>result1## endpoint1 endpoint2 Nogo_prob Grey_prob Go_prob Go_prob_text Nogo_prob_text
## 1 0 0 0 1.000 0.000 0.0 % 0.0 %
## 2 4 0 0 1.000 0.000 0.0 % 0.0 %
## 3 8 0 0 1.000 0.000 0.0 % 0.0 %
## 4 12 0 0 0.909 0.091 9.1 % 0.0 %
## 5 16 0 0 0.056 0.944 94.4 % 0.0 %
## 6 20 0 0 0.000 1.000 100.0 % 0.0 %
## 7 24 0 0 0.000 1.000 100.0 % 0.0 %
## 8 30 0 0 0.000 1.000 100.0 % 0.0 %
## 9 0 2 0 1.000 0.000 0.0 % 0.0 %
## 10 4 2 0 1.000 0.000 0.0 % 0.0 %
## 11 8 2 0 1.000 0.000 0.0 % 0.0 %
## 12 12 2 0 0.917 0.083 8.3 % 0.0 %
## 13 16 2 0 0.054 0.946 94.6 % 0.0 %
## 14 20 2 0 0.000 1.000 100.0 % 0.0 %
## 15 24 2 0 0.000 1.000 100.0 % 0.0 %
## 16 30 2 0 0.000 1.000 100.0 % 0.0 %
## 17 0 3 0 1.000 0.000 0.0 % 0.0 %
## 18 4 3 0 1.000 0.000 0.0 % 0.0 %
## 19 8 3 0 1.000 0.000 0.0 % 0.0 %
## 20 12 3 0 0.900 0.100 10.0 % 0.0 %
## 21 16 3 0 0.056 0.944 94.4 % 0.0 %
## 22 20 3 0 0.000 1.000 100.0 % 0.0 %
## 23 24 3 0 0.000 1.000 100.0 % 0.0 %
## 24 30 3 0 0.000 1.000 100.0 % 0.0 %
## 25 0 5 0 1.000 0.000 0.0 % 0.0 %
## 26 4 5 0 1.000 0.000 0.0 % 0.0 %
## 27 8 5 0 1.000 0.000 0.0 % 0.0 %
## 28 12 5 0 0.921 0.079 7.9 % 0.0 %
## 29 16 5 0 0.060 0.940 94.0 % 0.0 %
## 30 20 5 0 0.000 1.000 100.0 % 0.0 %
## 31 24 5 0 0.000 1.000 100.0 % 0.0 %
## 32 30 5 0 0.000 1.000 100.0 % 0.0 %
## 33 0 7 0 1.000 0.000 0.0 % 0.0 %
## 34 4 7 0 1.000 0.000 0.0 % 0.0 %
## 35 8 7 0 1.000 0.000 0.0 % 0.0 %
## 36 12 7 0 0.899 0.101 10.1 % 0.0 %
## 37 16 7 0 0.038 0.962 96.2 % 0.0 %
## 38 20 7 0 0.000 1.000 100.0 % 0.0 %
## 39 24 7 0 0.000 1.000 100.0 % 0.0 %
## 40 30 7 0 0.000 1.000 100.0 % 0.0 %
## 41 0 8 0 1.000 0.000 0.0 % 0.0 %
## 42 4 8 0 1.000 0.000 0.0 % 0.0 %
## 43 8 8 0 1.000 0.000 0.0 % 0.0 %
## 44 12 8 0 0.923 0.077 7.7 % 0.0 %
## 45 16 8 0 0.050 0.950 95.0 % 0.0 %
## 46 20 8 0 0.000 1.000 100.0 % 0.0 %
## 47 24 8 0 0.000 1.000 100.0 % 0.0 %
## 48 30 8 0 0.000 1.000 100.0 % 0.0 %
## 49 0 10 0 1.000 0.000 0.0 % 0.0 %
## 50 4 10 0 1.000 0.000 0.0 % 0.0 %
## 51 8 10 0 1.000 0.000 0.0 % 0.0 %
## 52 12 10 0 0.915 0.085 8.5 % 0.0 %
## 53 16 10 0 0.064 0.936 93.6 % 0.0 %
## 54 20 10 0 0.000 1.000 100.0 % 0.0 %
## 55 24 10 0 0.000 1.000 100.0 % 0.0 %
## 56 30 10 0 0.000 1.000 100.0 % 0.0 %
## 57 0 12 0 1.000 0.000 0.0 % 0.0 %
## 58 4 12 0 1.000 0.000 0.0 % 0.0 %
## 59 8 12 0 1.000 0.000 0.0 % 0.0 %
## 60 12 12 0 0.910 0.090 9.0 % 0.0 %
## 61 16 12 0 0.051 0.949 94.9 % 0.0 %
## 62 20 12 0 0.000 1.000 100.0 % 0.0 %
## 63 24 12 0 0.000 1.000 100.0 % 0.0 %
## 64 30 12 0 0.000 1.000 100.0 % 0.0 %
## Grey_prob_text
## 1 100.0 %
## 2 100.0 %
## 3 100.0 %
## 4 90.9 %
## 5 5.6 %
## 6 0.0 %
## 7 0.0 %
## 8 0.0 %
## 9 100.0 %
## 10 100.0 %
## 11 100.0 %
## 12 91.7 %
## 13 5.4 %
## 14 0.0 %
## 15 0.0 %
## 16 0.0 %
## 17 100.0 %
## 18 100.0 %
## 19 100.0 %
## 20 90.0 %
## 21 5.6 %
## 22 0.0 %
## 23 0.0 %
## 24 0.0 %
## 25 100.0 %
## 26 100.0 %
## 27 100.0 %
## 28 92.1 %
## 29 6.0 %
## 30 0.0 %
## 31 0.0 %
## 32 0.0 %
## 33 100.0 %
## 34 100.0 %
## 35 100.0 %
## 36 89.9 %
## 37 3.8 %
## 38 0.0 %
## 39 0.0 %
## 40 0.0 %
## 41 100.0 %
## 42 100.0 %
## 43 100.0 %
## 44 92.3 %
## 45 5.0 %
## 46 0.0 %
## 47 0.0 %
## 48 0.0 %
## 49 100.0 %
## 50 100.0 %
## 51 100.0 %
## 52 91.5 %
## 53 6.4 %
## 54 0.0 %
## 55 0.0 %
## 56 0.0 %
## 57 100.0 %
## 58 100.0 %
## 59 100.0 %
## 60 91.0 %
## 61 5.1 %
## 62 0.0 %
## 63 0.0 %
## 64 0.0 %Nogo Plot
Nogo_plot = ggplot(result1, aes(endpoint1, endpoint2)) +
geom_tile(aes(fill = Nogo_prob)) +
scale_fill_gradient(low = "white", high = "red", guide = guide_colorbar(title = "Probility of No-Go"), limits = c(0, 1), breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) +
geom_text(aes(label = Nogo_prob_text)) +
theme_classic()+
scale_x_continuous(breaks=c(0,4,8,12,16,20,24,30))+
scale_y_continuous(breaks=c(0,2,3,5,7,8,10,12))+
ggtitle("Probability of No-Go")+
labs(x="MWT", y="ESS")
Nogo_plot
Grey Plot
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.