Single_arm_binomial_demo
2024-05-17
my_density_threshold = function(den,threshold,PP,
direction = "greater",
Criteria = "Go",title0,title1){
# Criteria = "Go"
# direction = "less"
# direction = "greater"
# title0 = "Go Rule:PP(Value>=0.2)>=0.8, n=25"
#
# lwd = 1.5
# n0 = 25
# x0 = seq(0,1,0.001)
# mu0 = 0.3
# sd0 = sqrt(mu0*(1-mu0)/n0)
# y0 = dnorm(x0,mu0,sd0)
# den = cbind(x0,y0)
# threshold = 0.2
# threshold = c(0.2,0.15)
# PP = 1-pnorm(threshold,mu0,sd0)
PP = round(PP*100,1)
text_height = den[ which.min(abs(x0 - threshold[1])) ,2]
df<-data.frame('x'=den[,1],'y'=den[,2])
tmp = subset(df,y>max(df$y)*0.1)
text_x = mean(tmp$x * tmp$y) / mean(tmp$y)
tmp = subset(df,y>max(df$y)*0.05)
color0 = ifelse(Criteria == "Go","green4","red3")
p = ggplot(df, aes(x=x, y=y)) + xlim(c(min(tmp$x),max(tmp$x)))+
ylab("Density") + xlab("Value")+geom_line()
if(direction == "greater" & length(threshold) == 1){
p = p + geom_ribbon(data=subset(df,x> threshold),aes(x=x,ymax=y),ymin=0,alpha=0.3,fill= color0)+
annotate("text", x = text_x , y = text_height/2,
label = paste("PP(Value >=", threshold, ") = ",PP,"%",sep=""),
size = 5, fontface = "bold.italic", colour = color0)
}else if(direction == "less" & length(threshold) == 1){
p = p + geom_ribbon(data=subset(df,x< threshold),aes(x=x,ymax=y),ymin=0,alpha=0.3,fill= color0)+
annotate("text", x = text_x , y = text_height/2,
label = paste("PP(Value <=", threshold, ") = ",PP,"%",sep=""),
size = 5, fontface = "bold.italic", colour = color0)
}else if(direction == "greater" & length(threshold) == 2){
p = p + geom_ribbon(data=subset(df,x> threshold[1]),aes(x=x,ymax=y),ymin=0,alpha=0.3,fill= color0)+
geom_ribbon(data=subset(df,x> threshold[2]),aes(x=x,ymax=y),ymin=0,alpha=0.3,fill= color0)+
annotate("text", x = text_x , y = text_height/2,
label = paste("PP(Value >=", threshold[1], ") = ",PP[1],"%,\n",
"PP(Value >=", threshold[2], ") = ",PP[2],"%.",sep=""),
size = 5, fontface = "bold.italic", colour = color0)
}else if(direction == "less" & length(threshold) == 2){
p = p + geom_ribbon(data=subset(df,x< threshold[1]),aes(x=x,ymax=y),ymin=0,alpha=0.3,fill= color0)+
geom_ribbon(data=subset(df,x< threshold[2]),aes(x=x,ymax=y),ymin=0,alpha=0.3,fill= color0)+
annotate("text", x = text_x , y = text_height/2,
label = paste("PP(Value <=", threshold[1], ") = ",PP[1],"%,\n",
"PP(Value <=", threshold[2], ") = ",PP[2],"%.",sep=""),
size = 5, fontface = "bold.italic", colour = color0)
}
p + theme_classic() + theme(legend.position = "bottom",
text = element_text(size = 15),
plot.title = element_text(hjust=0),
aspect.ratio=0.8) +
labs(title = title0, subtitle = title1) +
theme(
plot.title = element_text( face = "bold"),
plot.subtitle = element_text(size = 12)
)
}Goahead = function(Go_rule_1 = FALSE, Go_value_1 = NA,Go_PP_1=80,
Go_rule_2 = FALSE, Go_value_2 = NA,Go_PP_2=80,
NoGo_rule_1 = FALSE, NoGo_value_1 = NA,NoGo_PP_1=10,
NoGo_rule_2 = FALSE, NoGo_value_2 = NA,NoGo_PP_2=10,
n, accuracy = 0.001,Go_Combine_rule = "and",NoGo_Combine_rule = "and",
dominated_rule = "Go",direction_comparision = "greater",
size_true = 30, value_true = c(0.2,0.3,0.4,0.5),
method = c("Frequentist:exact","Bayesian"),
alpha = 1,beta = 1){
# method = "midp"
n = size_true
# For debug
# ---------------------------------------------------------------
# alpha = 1
# beta = 1
# n = size_true = 25
# accuracy = 0.0001
# dominated_rule = "Go"
# direction_comparision = "greater"
# # direction_comparision = "less"
#
# Go_rule_1 = TRUE
# Go_rule_2 = FALSE
# Go_Combine_rule = "and"
#
# Go_value_1 = 0.2
# Go_PP_1 = 80
# Go_value_2 = 0.25
# Go_PP_2 = 75
#
#
# NoGo_rule_1 = TRUE
# NoGo_rule_2 = FALSE
# NoGo_Combine_rule = "and"
#
# NoGo_value_1 = 0.3
# NoGo_PP_1 = 10
# NoGo_value_2 = 0.35
# NoGo_PP_2 = 5
#
# ---------------------------------------------------------
Go_PP_1 = Go_PP_1 / 100
Go_PP_2 = Go_PP_2 / 100
NoGo_PP_1 = NoGo_PP_1 / 100
NoGo_PP_2 = NoGo_PP_2 / 100
Go_criteria_1 = NA
Go_criteria_2 = NA
NoGo_criteria_1 = NA
NoGo_criteria_2 = NA
### Define functions
vec_theta = seq(0,n,1)
my_exact = function(x,n,alpha,sides = "left"){
n0 = length(x)
tmp = matrix(0,n0,2)
for(ii in 1:n0){
if(sides == "left"){
tmp[ii,1] = (qbeta(1-alpha,x[ii],n-x[ii]+1))
tmp[ii,2] = 1
}else{
tmp[ii,2] = (qbeta(alpha,x[ii]+1,n-x[ii]))
}
}
return( tmp)
}
my_beta = function(x,n,alpha,beta,CI,sides="left"){
n0 = length(x)
tmp = matrix(0,n0,2)
for(ii in 1:n0){
if(sides == "left"){
tmp[ii,1] = (qbeta(1-CI,x[ii]+alpha,n-x[ii]+beta))
tmp[ii,2] = 1
}else{
tmp[ii,2] = (qbeta(CI,x[ii]+alpha,n-x[ii]+beta))
}
}
return( tmp)
}
### Search for go criteria
if(Go_rule_1){
if(direction_comparision == "greater"){
# vec_Go_criteria_1 = BinomCI(vec_theta, n, conf.level = Go_PP_1, sides ="left",
# method =method,
# rand = 123, tol = 1e-05, std_est = TRUE)
# chosen = which(vec_Go_criteria_1[,2] >= Go_value_1)
if(method=="Frequentist:exact"){
vec_Go_criteria_1 = my_exact(vec_theta,n,Go_PP_1,sides = "left")
}else if(method=="Bayesian"){
vec_Go_criteria_1 = my_beta(vec_theta,n,alpha,beta,Go_PP_1,sides = "left")
}
chosen = which( !(vec_Go_criteria_1[,1] <= Go_value_1&
vec_Go_criteria_1[,2] >= Go_value_1))
Go_criteria_1 = min(vec_theta[chosen])
}else{
if(method=="Frequentist:exact"){
vec_Go_criteria_1 = my_exact(vec_theta,n,Go_PP_1,sides = "right")
}else if(method=="Bayesian"){
vec_Go_criteria_1 = my_beta(vec_theta,n,alpha,beta,Go_PP_1,sides = "right")
}
chosen = which( !(vec_Go_criteria_1[,1] <= Go_value_1&
vec_Go_criteria_1[,2] >= Go_value_1))
Go_criteria_1 = max(vec_theta[chosen])
}
}
if(Go_rule_2){
if(direction_comparision == "greater"){
if(method=="Frequentist:exact"){
vec_Go_criteria_2 = my_exact(vec_theta,n,Go_PP_2,sides = "left")
}else if(method=="Bayesian"){
vec_Go_criteria_2 = my_beta(vec_theta,n,alpha,beta,Go_PP_2,sides = "left")
}
chosen = which( !(vec_Go_criteria_2[,1] <= Go_value_2&
vec_Go_criteria_2[,2] >= Go_value_2))
Go_criteria_2 = min(vec_theta[chosen])
}else{
if(method=="Frequentist:exact"){
vec_Go_criteria_2 = my_exact(vec_theta,n,Go_PP_2,sides = "right")
}else if(method=="Bayesian"){
vec_Go_criteria_2 = my_beta(vec_theta,n,alpha,beta,Go_PP_2,sides = "right")
}
chosen = which( !(vec_Go_criteria_2[,1] <= Go_value_2&
vec_Go_criteria_2[,2] >= Go_value_2))
Go_criteria_2 = max(vec_theta[chosen])
}
}
### Search for stop criteria
# ------NoGo value 1------PP 20%----------------
if(NoGo_rule_1){
if(direction_comparision == "greater"){
if(method=="Frequentist:exact"){
vec_NoGo_criteria_1 = my_exact(vec_theta,n,NoGo_PP_1,sides = "left")
}else if(method=="Bayesian"){
vec_NoGo_criteria_1 = my_beta(vec_theta,n,alpha,beta,NoGo_PP_1,sides = "left")
}
chosen = which( (vec_NoGo_criteria_1[,1] <= NoGo_value_1&
vec_NoGo_criteria_1[,2] >= NoGo_value_1))
NoGo_criteria_1 = max(vec_theta[chosen])
}else{
if(method=="Frequentist:exact"){
vec_NoGo_criteria_1 = my_exact(vec_theta,n,NoGo_PP_1,sides = "right")
}else if(method=="Bayesian"){
vec_NoGo_criteria_1 = my_beta(vec_theta,n,alpha,beta,NoGo_PP_1,sides = "right")
}
chosen = which( (vec_NoGo_criteria_1[,1] <= NoGo_value_1&
vec_NoGo_criteria_1[,2] >= NoGo_value_1))
NoGo_criteria_1 = min(vec_theta[chosen])
}
}
if(NoGo_rule_2){
if(direction_comparision == "greater"){
if(method=="Frequentist:exact"){
vec_NoGo_criteria_2 = my_exact(vec_theta,n,NoGo_PP_2,sides = "left")
}else if(method=="Bayesian"){
vec_NoGo_criteria_2 = my_beta(vec_theta,n,alpha,beta,NoGo_PP_2,sides = "left")
}
chosen = which( (vec_NoGo_criteria_1[,1] <= NoGo_value_1&
vec_NoGo_criteria_1[,2] >= NoGo_value_1))
NoGo_criteria_2 = max(vec_theta[chosen])
}else{
if(method=="Frequentist:exact"){
vec_NoGo_criteria_2 = my_exact(vec_theta,n,NoGo_PP_2,sides = "right")
}else if(method=="Bayesian"){
vec_NoGo_criteria_2 = my_beta(vec_theta,n,alpha,beta,NoGo_PP_2,sides = "right")
}
chosen = which( (vec_NoGo_criteria_1[,1] <= NoGo_value_1&
vec_NoGo_criteria_1[,2] >= NoGo_value_1))
NoGo_criteria_2 = min(vec_theta[chosen])
}
}
## combine multiple criteria, case when direction is "greater"
if(direction_comparision == "greater"){
if(Go_rule_1 && Go_rule_2){
Value_in_Go = paste(Go_value_1,Go_value_2,sep = ",")
if(Go_Combine_rule =="and"){
Final_Go = max(Go_criteria_1,Go_criteria_2)
Decision_rule_for_Go = paste("PP(value >= ",Go_value_1,") >= ", Go_PP_1*100,
"% and\n PP(value >= ",Go_value_2,") >= ", Go_PP_2*100,"%",sep="")
}else{
Final_Go = min(Go_criteria_1,Go_criteria_2)
Decision_rule_for_Go = paste("PP(value >= ",Go_value_1,") >= ", Go_PP_1*100,
"% or\n PP(value >= ",Go_value_2,") >= ", Go_PP_2*100,"%",sep="")
}
}else if(Go_rule_1){
Value_in_Go = Go_value_1
Decision_rule_for_Go = paste("PP(value >= ",Go_value_1,") >= ", Go_PP_1*100,"%",sep="")
Final_Go = Go_criteria_1
}else if(Go_rule_2){
Value_in_Go = Go_value_2
Decision_rule_for_Go = paste("PP(value >= ",Go_value_2,") >= ", Go_PP_2*100,"%",sep="")
Final_Go = Go_criteria_2
}else{
Value_in_Go = NA
Decision_rule_for_Go = NA
cat("You must select a Go Rule!")
}
if(NoGo_rule_1 && NoGo_rule_2){
Value_in_NoGo = paste(NoGo_value_1,NoGo_value_2,sep = ",")
if(NoGo_Combine_rule =="and"){
Final_NoGo = max(NoGo_criteria_1,NoGo_criteria_2)
Decision_rule_for_NoGo = paste("PP(value >= ",NoGo_value_1,") < ", NoGo_PP_1*100,
"% and\n PP(value >= ",NoGo_value_2,") < ", NoGo_PP_2*100,"%",sep="")
}else{
Final_NoGo = min(NoGo_criteria_1,NoGo_criteria_2)
Decision_rule_for_NoGo = paste("PP(value >= ",NoGo_value_1,") < ", NoGo_PP_1*100,
"% or\n PP(value >= ",NoGo_value_2,") < ", NoGo_PP_2*100,"%",sep="")
}
}else if(NoGo_rule_1){
Value_in_NoGo = NoGo_value_1
Decision_rule_for_NoGo = paste("PP(value >= ",NoGo_value_1,") < ", NoGo_PP_1*100,"%",sep="")
Final_NoGo = NoGo_criteria_1
}else if(NoGo_rule_2){
Value_in_NoGo = NoGo_value_2
Decision_rule_for_NoGo = paste("PP(value >= ",NoGo_value_2,") < ", NoGo_PP_2*100,"%",sep="")
Final_NoGo = NoGo_criteria_2
}else{
Value_in_NoGo = NA
Decision_rule_for_NoGo = NA
cat("You must select a NoGo Rule!")
}
}
## combine multiple criteria, case when direction is "less"
if(direction_comparision == "less"){
if(Go_rule_1 && Go_rule_2){
Value_in_Go = paste(Go_value_1,Go_value_2,sep = ",")
if(Go_Combine_rule =="and"){
Final_Go = min(Go_criteria_1,Go_criteria_2)
Decision_rule_for_Go = paste("PP(value <= ",Go_value_1,") >= ", Go_PP_1*100,
"% and \n PP(value <= ",Go_value_2,") >= ", Go_PP_2*100,"%",sep="")
}else{
Final_Go = max(Go_criteria_1,Go_criteria_2)
Decision_rule_for_Go = paste("PP(value <= ",Go_value_1,") >= ", Go_PP_1*100,
"% or\n PP(value <= ",Go_value_2,") >= ", Go_PP_2*100,"%",sep="")
}
}else if(Go_rule_1){
Value_in_Go = Go_value_1
Decision_rule_for_Go = paste("PP(value <= ",Go_value_1,") >= ", Go_PP_1*100,"%",sep="")
Final_Go = Go_criteria_1
}else if(Go_rule_2){
Value_in_Go = Go_value_2
Decision_rule_for_Go = paste("PP(value <= ",Go_value_2,") >= ", Go_PP_2*100,"%",sep="")
Final_Go = Go_criteria_2
}else{
Value_in_Go = NA
Decision_rule_for_Go = NA
cat("You must select a Go Rule!")
}
if(NoGo_rule_1 && NoGo_rule_2){
Value_in_NoGo = paste(NoGo_value_1,NoGo_value_2,sep = ",")
if(NoGo_Combine_rule =="and"){
Final_NoGo = min(NoGo_criteria_1,NoGo_criteria_2)
Decision_rule_for_NoGo = paste("PP(value <= ",NoGo_value_1,") < ", NoGo_PP_1*100,
"% and\n PP(value <= ",NoGo_value_2,") < ", NoGo_PP_2*100,"%",sep="")
}else{
Final_NoGo = max(NoGo_criteria_1,NoGo_criteria_2)
Decision_rule_for_NoGo = paste("PP(value <= ",NoGo_value_1,") < ", NoGo_PP_1*100,
"% or\n PP(value <= ",NoGo_value_2,") < ", NoGo_PP_2*100,"%",sep="")
}
}else if(NoGo_rule_1){
Value_in_NoGo = NoGo_value_1
Decision_rule_for_NoGo = paste("PP(value <= ",NoGo_value_1,") < ", NoGo_PP_1*100,"%",sep="")
Final_NoGo = NoGo_criteria_1
}else if(NoGo_rule_2){
Value_in_NoGo = NoGo_value_2
Decision_rule_for_NoGo = paste("PP(value <= ",NoGo_value_2,") < ", NoGo_PP_2*100,"%",sep="")
Final_NoGo = NoGo_criteria_2
}else{
Value_in_NoGo = NA
Decision_rule_for_NoGo = NA
cat("You must select a NoGo Rule!")
}
}
## doninated rules
if(direction_comparision == "greater"){
if(dominated_rule == "Go"){
Final_NoGo = min(Final_NoGo,Final_Go)
}
if(dominated_rule == "NoGo"){
Final_Go = max(Final_NoGo,Final_Go)
}
}
if(direction_comparision == "less"){
if(dominated_rule == "Go"){
Final_NoGo = max(Final_NoGo,Final_Go)
}
if(dominated_rule == "NoGo"){
Final_Go = min(Final_NoGo,Final_Go)
}
}
### Specify sample size and true value
Go_zone = Final_Go
NoGo_zone = Final_NoGo
Cutoff_for_Go = paste(">=",Go_zone)
Cutoff_for_NoGo = paste("<=",NoGo_zone)
pbinom(NoGo_zone,size_true,value_true)
pbinom(Go_zone-1,size_true,value_true,lower.tail = FALSE)
1 - pbinom(NoGo_zone,size_true,value_true) - pbinom(Go_zone-1,size_true,value_true,lower.tail = FALSE)
n0 = length(value_true)
df.output = data.frame(x1 = rep(size_true,n0),
x2 = value_true,
x3 = pbinom(Go_zone-1,size_true,value_true,lower.tail = FALSE),
x4 = pbinom(NoGo_zone,size_true,value_true),
x5 = 1 - pbinom(NoGo_zone,size_true,value_true) -
pbinom(Go_zone-1,size_true,value_true,lower.tail = FALSE),
x6 = paste(Cutoff_for_Go,Cutoff_for_NoGo,sep = "/"),
x7 = Cutoff_for_Go,
x8 = Cutoff_for_NoGo)
df.output$x3 = round(df.output$x3,3)
df.output$x4 = round(df.output$x4,3)
df.output$x5 = round(df.output$x5,3)
colnames(df.output) = c("Sample Size","True Value","Pr(Go)","Pr(NoGo)",
"Pr(inconclusive)","Go/NoGo Zone","Cut off for Go","Cut off for NoGo")
Criteria = c(Go_criteria_1,Go_criteria_2,NoGo_criteria_1,NoGo_criteria_2)
names(Criteria) = c("Go1","Go2","NoGo1","NoGo2")
return(list(df.output = df.output,Criteria = Criteria,
Value_in_Go = Value_in_Go,
Value_in_NoGo = Value_in_NoGo,
Decision_rule_for_Go = Decision_rule_for_Go,
Decision_rule_for_NoGo = Decision_rule_for_NoGo,
Final_Go = Final_Go,
Final_NoGo = Final_NoGo))
}
Goahead( Go_rule_1 = TRUE, Go_value_1 = 0.3,Go_PP_1=80,
Go_rule_2 = FALSE, Go_value_2 = 0.25,Go_PP_2=80,
NoGo_rule_1 = TRUE, NoGo_value_1 = 0.4,NoGo_PP_1=10,
NoGo_rule_2 = FALSE, NoGo_value_2 = 0.35,NoGo_PP_2=5,
n, accuracy = 0.001,Go_Combine_rule = "and",NoGo_Combine_rule = "and",
dominated_rule = "Go", direction_comparision = "greater",
size_true = 41, value_true = c(0.2,0.3,0.4,0.5),method="Frequentist:exact")## $df.output
## Sample Size True Value Pr(Go) Pr(NoGo) Pr(inconclusive) Go/NoGo Zone
## 1 41 0.2 0.004 0.948 0.048 >= 16/<= 12
## 2 41 0.3 0.138 0.536 0.325 >= 16/<= 12
## 3 41 0.4 0.609 0.105 0.286 >= 16/<= 12
## 4 41 0.5 0.941 0.006 0.053 >= 16/<= 12
## Cut off for Go Cut off for NoGo
## 1 >= 16 <= 12
## 2 >= 16 <= 12
## 3 >= 16 <= 12
## 4 >= 16 <= 12
##
## $Criteria
## Go1 Go2 NoGo1 NoGo2
## 16 NA 12 NA
##
## $Value_in_Go
## [1] 0.3
##
## $Value_in_NoGo
## [1] 0.4
##
## $Decision_rule_for_Go
## [1] "PP(value >= 0.3) >= 80%"
##
## $Decision_rule_for_NoGo
## [1] "PP(value >= 0.4) < 10%"
##
## $Final_Go
## [1] 16
##
## $Final_NoGo
## [1] 12Goahead( Go_rule_1 = TRUE, Go_value_1 = 0.3,Go_PP_1=80,
Go_rule_2 = TRUE, Go_value_2 = 0.25,Go_PP_2=90,
NoGo_rule_1 = TRUE, NoGo_value_1 = 0.4,NoGo_PP_1=10,
NoGo_rule_2 = TRUE, NoGo_value_2 = 0.45,NoGo_PP_2=5,
n, accuracy = 0.001,Go_Combine_rule = "and",NoGo_Combine_rule = "and",
dominated_rule = "Go", direction_comparision = "greater",
size_true = 41, value_true = c(0.2,0.3,0.4,0.5),method="Bayesian")## $df.output
## Sample Size True Value Pr(Go) Pr(NoGo) Pr(inconclusive) Go/NoGo Zone
## 1 41 0.2 0.010 0.948 0.042 >= 15/<= 12
## 2 41 0.3 0.224 0.536 0.240 >= 15/<= 12
## 3 41 0.4 0.725 0.105 0.169 >= 15/<= 12
## 4 41 0.5 0.970 0.006 0.024 >= 15/<= 12
## Cut off for Go Cut off for NoGo
## 1 >= 15 <= 12
## 2 >= 15 <= 12
## 3 >= 15 <= 12
## 4 >= 15 <= 12
##
## $Criteria
## Go1 Go2 NoGo1 NoGo2
## 15 14 12 12
##
## $Value_in_Go
## [1] "0.3,0.25"
##
## $Value_in_NoGo
## [1] "0.4,0.45"
##
## $Decision_rule_for_Go
## [1] "PP(value >= 0.3) >= 80% and\n PP(value >= 0.25) >= 90%"
##
## $Decision_rule_for_NoGo
## [1] "PP(value >= 0.4) < 10% and\n PP(value >= 0.45) < 5%"
##
## $Final_Go
## [1] 15
##
## $Final_NoGo
## [1] 12method=c("Frequentist:exact","Bayesian")
Go_rule_1 = TRUE;Go_value_1 = 0.2;Go_PP_1=80
Go_rule_2 = FALSE;Go_value_2 = 0.25;Go_PP_2=80
NoGo_rule_1 = TRUE; NoGo_value_1 = 0.3;NoGo_PP_1=10;
NoGo_rule_2 = FALSE; NoGo_value_2 = 0.35;NoGo_PP_2=5;
n= 25;accuracy = 0.0001;Go_Combine_rule = "and";NoGo_Combine_rule = "and";
dominated_rule = "Go";direction_comparision = "greater";
size_true = 25;value_true = c(0.2,0.3,0.4,0.5)method = method[2]
res = Goahead( Go_rule_1 = TRUE, Go_value_1 = 0.2,Go_PP_1=80,
Go_rule_2 = FALSE, Go_value_2 = 0.25,Go_PP_2=80,
NoGo_rule_1 = TRUE, NoGo_value_1 = 0.3,NoGo_PP_1=10,
NoGo_rule_2 = FALSE, NoGo_value_2 = 0.35,NoGo_PP_2=5,
n, accuracy = 0.0001,Go_Combine_rule = "and",NoGo_Combine_rule = "and",
dominated_rule = "Go", direction_comparision = "greater",
size_true = 25, value_true = c(0.2,0.3,0.4,0.5),method=method)
knitr::kable(res$df.output)| Sample Size | True Value | Pr(Go) | Pr(NoGo) | Pr(inconclusive) | Go/NoGo Zone | Cut off for Go | Cut off for NoGo |
|---|---|---|---|---|---|---|---|
| 25 | 0.2 | 0.220 | 0.421 | 0.359 | >= 7/<= 4 | >= 7 | <= 4 |
| 25 | 0.3 | 0.659 | 0.090 | 0.250 | >= 7/<= 4 | >= 7 | <= 4 |
| 25 | 0.4 | 0.926 | 0.009 | 0.064 | >= 7/<= 4 | >= 7 | <= 4 |
| 25 | 0.5 | 0.993 | 0.000 | 0.007 | >= 7/<= 4 | >= 7 | <= 4 |
Fix True value, change sample size
value_true = 0.3
sample_size = seq(12,42,1)
record_criteria = matrix(0,length(sample_size),2)
method=c("Frequentist:exact","Bayesian")
method = method[2]
for(ii in 1:length(sample_size)){
res=Goahead( Go_rule_1 = TRUE, Go_value_1 = 0.2,Go_PP_1=80,
Go_rule_2 = FALSE, Go_value_2 = 0.25,Go_PP_2=80,
NoGo_rule_1 = TRUE, NoGo_value_1 = 0.3,NoGo_PP_1=10,
NoGo_rule_2 = FALSE, NoGo_value_2 = 0.35,NoGo_PP_2=5,
n, accuracy = 0.0001,Go_Combine_rule = "and",NoGo_Combine_rule = "and",
dominated_rule = "Go", direction_comparision = "greater",
size_true = sample_size[ii], value_true = value_true,
method=method)
record_criteria[ii,1] = res[[1]]$`Pr(NoGo)`
record_criteria[ii,2] = res[[1]]$`Pr(NoGo)` + res[[1]]$`Pr(inconclusive)`
}x0 = 1:10
df_vis = data.frame(x = sample_size,
y1 = record_criteria[,1]*100,
y2 = record_criteria[,2]*100)
breaks_x = c(12,15,20,25,30,36,42)
labels_x = c(12,15,20,25,30,36,42)
title0 = paste("True Value =", value_true)
ggplot(data = df_vis, mapping = aes(x = x)) +
xlab("Sample Size") +
ylab("Probability of Go/NoGo/Inconclusive (%)") +
theme(axis.title = element_text(size=20)) +
scale_x_continuous(limits = c(12,42), breaks=breaks_x,labels=labels_x) +
#scale_y_continuous(limits=c(0,1), sec.axis = sec_axis(~., name = "Absolute risk of LP without death (blue)")) +
scale_y_continuous(limits=c(0,100)) +
geom_ribbon(aes(ymin=0,ymax=y1), fill="red", alpha=1, show.legend = TRUE) +
geom_ribbon(aes(ymin=y1,ymax=y2), fill="#e4ce00", alpha=1, show.legend = TRUE) +
geom_ribbon(aes(ymin=y2,ymax=100), fill="#2E8B57", alpha=1, show.legend = TRUE) +
scale_fill_manual(labels=c("NoGo","go","Inconclusive"))+ theme_bw() +
theme(legend.position = "bottom",
text = element_text(size = 15),
plot.title = element_text(hjust=0.5),
aspect.ratio=0.8) + ggtitle(title0)
Fix sample size, change true value
size_true = 25
sample_value = seq(0.20,0.50,length.out=6)
record_prob = matrix(0,length(sample_value),3)
method=c("Frequentist:exact","Bayesian")
method = method[2]
for(ii in 1:length(sample_value)){
res=Goahead( Go_rule_1 = TRUE, Go_value_1 = 0.2,Go_PP_1=80,
Go_rule_2 = FALSE, Go_value_2 = 0.25,Go_PP_2=80,
NoGo_rule_1 = TRUE, NoGo_value_1 = 0.3,NoGo_PP_1=10,
NoGo_rule_2 = FALSE, NoGo_value_2 = 0.35,NoGo_PP_2=5,
n, accuracy = 0.0001,Go_Combine_rule = "and",NoGo_Combine_rule = "and",
dominated_rule = "Go", direction_comparision = "greater",
size_true = size_true, value_true = sample_value[ii],method=method)
record_prob[ii,1] = res[[1]]$`Pr(NoGo)`
record_prob[ii,2] = res[[1]]$`Pr(inconclusive)`
record_prob[ii,3] = res[[1]]$`Pr(Go)`
}
df_vis = data.frame(x = c(sample_value,sample_value,sample_value),
y = c(record_prob[,3],record_prob[,1],record_prob[,2]) * 100,
group = rep(c("Go","NoGo","Inconclusive"),each = length(sample_value)))
breaks_x = c(0.20,0.25,0.30,0.35,0.40,0.45,0.50)
labels_x = breaks_x
df_vis$group = as.factor(df_vis$group)
# word_go = paste(" <span style='color:#2E8B57;'>Go: Number of response>=</span>", 8)
# word_nogo = paste(" <span style='color:red;'>NoGo: Number of response<=</span>", 6)
# word_in = paste(" <span style='color:#e4ce00;'>Inconclusive: Anything else</span>")
# title0 = paste(word_go,word_nogo,word_in,sep="\n")
# title0 = paste("<span style='font-size:11pt'>",title0,"</span>")
title0 = paste("Sample Size =",size_true)https://stackoverflow.com/questions/39321483/multicolored-title-with-r
# library(ggtext)
p = ggplot(data = df_vis, mapping = aes(x = x,y= y,col = group,
linetype = group)) +
xlab("True Value") +
ylab("Probability of Go/NoGo/Inconclusive (%)") +
theme(axis.title = element_text(size=20), axis.text = element_text(size=18)) +
scale_x_continuous(limits = c(0.2,0.5), breaks=breaks_x,labels=labels_x) +
scale_y_continuous(limits=c(0,100)) +
geom_line(lwd=1.3)+geom_point(aes(shape = group),size = 2,stroke = 1.5)+
theme_classic()+theme(legend.position = "bottom")+
scale_shape_manual(values = c("Go"=16, "NoGo" = 4, "Inconclusive"=17))+
scale_linetype_manual(values=c("NoGo"="twodash",
"Go" = "solid", "Inconclusive" = "dashed"))+
scale_color_manual(values = c("NoGo"="red",
"Go" = "#2E8B57", "Inconclusive" = "#e4ce00"))+
theme_classic() +
scale_fill_discrete( labels = c("NoGo"="NoGo",
"Go" = "Go111", "Inconclusive" = "Inconclusive")) +
theme(legend.position = "bottom",
plot.title = element_text(hjust=0.5),
text = element_text(size = 14),aspect.ratio=0.8)+ggtitle(title0)
p$labels$colour =p$labels$shape =p$labels$linetype= "Decision"
p
# +
# ggtitle("**Fisher's *Iris* dataset**
# <span style='font-size:11pt'>Sepal width vs. sepal length for
# <span style='color:#0072B2;'>setosa</span>,
# <span style='color:#D55E00;'>versicolor</span>, and
# <span style='color:#009E73;'>virginica</span>
# </span>")
#
#
# p=p+tt("I spend more time with 'pickup' than\nwith 'other family members.'",
# "grey", 8.5, 0) +
# tt("I spend more time with 'pickup' than\nwith",
# "black", 8.5, 0)
#
# gt <- ggplot_gtable(ggplot_build(p))
# gt$layout$clip[gt$layout$name == "panel"] <- "off"
# grid.draw(gt)Cut Off Visualization
sample_value = seq(12,51,1)
record_rule = matrix(0,length(sample_value),2)
breaks_x = labels_x = seq(12,51,3)
for(ii in 1:length(sample_value)){
res=Goahead( Go_rule_1 = TRUE, Go_value_1 = 0.2,Go_PP_1=80,
Go_rule_2 = FALSE, Go_value_2 = 0.25,Go_PP_2=90,
NoGo_rule_1 = TRUE, NoGo_value_1 = 0.3,NoGo_PP_1=10,
NoGo_rule_2 = FALSE, NoGo_value_2 = 0.35,NoGo_PP_2=5,
n, accuracy = 0.0001,Go_Combine_rule = "and",NoGo_Combine_rule = "and",
dominated_rule = "Go", direction_comparision = "greater",
size_true = sample_value[ii], value_true = 0.5,method="Frequentist:exact")
record_rule[ii,1] = res$Final_Go
record_rule[ii,2] = (res$Final_NoGo)
}
df_vis = data.frame(value = as.numeric(record_rule)/sample_value,
x = rep(sample_value,2),
group = rep(c("Cut off of Go","Cut off of NoGo"),each = length(sample_value)))
p = ggplot(data = df_vis,aes(x = x,y = value,col = group)) + geom_line()+
xlab("Sample Size") +
ylab("Observed Ratio (%)") +
theme(axis.title = element_text(size=20), axis.text = element_text(size=18)) +
scale_x_continuous(limits = c(min(breaks_x),max(breaks_x)), breaks=breaks_x,labels=labels_x) +
scale_y_continuous() +
geom_line(lwd=1.3)+
theme_classic()+theme(legend.position = "bottom")+
scale_color_manual(values = c("Cut off of Go" = "#2E8B57",
"Cut off of NoGo"="red"))+
theme_classic() +
theme(legend.position = "bottom",
plot.title = element_text(hjust=0.5),
text = element_text(size = 14),aspect.ratio=0.8) + ggtitle("Cut Off / Sample Size")
p$labels$colour = "Rules"
p
### Cut off
breaks_y = labels_y = seq(min(record_rule),max(record_rule),2)
df_vis = data.frame(value = as.numeric(record_rule),
x = rep(sample_value,2),
group = rep(c("Cut off of Go","Cut off of NoGo"),each = length(sample_value)))
p2 = ggplot(data = df_vis,aes(x = x,y = value,col = group)) + geom_line()+
xlab("Sample Size") +
ylab("Observed Number") +
theme(axis.title = element_text(size=20), axis.text = element_text(size=18)) +
scale_x_continuous(limits = c(min(breaks_x),max(breaks_x)), breaks=breaks_x,labels=labels_x) +
scale_y_continuous(breaks=breaks_y,labels=labels_y) +
geom_line(lwd=1.3)+
theme_classic()+theme(legend.position = "bottom")+
scale_color_manual(values = c("Cut off of Go" = "#2E8B57",
"Cut off of NoGo"="red"))+
theme_classic() +
theme(legend.position = "bottom",
plot.title = element_text(hjust=0.5),
text = element_text(size = 14),aspect.ratio=0.8) + ggtitle("Cut Off")
p2$labels$colour = "Rules"
p2
Density Visualization
size_true = 25
method = c("Frequentist:exact","Bayesian")
method = method[2]
res=Goahead( Go_rule_1 = TRUE, Go_value_1 = 0.2,Go_PP_1=80,
Go_rule_2 = TRUE, Go_value_2 = 0.25,Go_PP_2=70,
NoGo_rule_1 = TRUE, NoGo_value_1 = 0.3,NoGo_PP_1=10,
NoGo_rule_2 = TRUE, NoGo_value_2 = 0.35,NoGo_PP_2=5,
n, accuracy = 0.0001,Go_Combine_rule = "or",NoGo_Combine_rule = "or",
dominated_rule = "Go", direction_comparision = "greater",
size_true = 25, value_true = 0.20,
method=method)
n0 = 25
x0 = seq(0,1,0.001)
alpha = 1
beta = 1
# For Go rule
cutoff = res$Final_Go
if(method=="Frequentist:exact"){
mu0 = cutoff / n0
sd0 = sqrt(mu0*(1-mu0)/n0)
y0 = dnorm(x0,mu0,sd0)
den = cbind(x0,y0)
threshold = c(0.2,0.25)
PP = 1-pnorm(threshold,mu0,sd0)
}else if(method == "Bayesian"){
y0 = dbeta(x0,cutoff + alpha,n0 - cutoff + beta)
den = cbind(x0,y0)
threshold = c(0.2,0.25)
PP = 1-pbeta(threshold,cutoff + alpha,n0 - cutoff + beta)
}
p1 = my_density_threshold(den,threshold,PP,direction = "greater",
Criteria = "Go",
title1=paste("When #Obs/n = ",res$Final_Go,"/",size_true,":",sep=""),
title0=paste("Go when #Obs", res$df.output$`Cut off for Go`,":\n ",res$Decision_rule_for_Go,sep="") )
# For NoGo rule
cutoff = res$Final_NoGo
if(method=="Frequentist:exact"){
mu0 = cutoff / n0
sd0 = sqrt(mu0*(1-mu0)/n0)
y0 = dnorm(x0,mu0,sd0)
den = cbind(x0,y0)
threshold = c(0.2,0.25)
PP = 1-pnorm(threshold,mu0,sd0)
}else if(method == "Bayesian"){
y0 = dbeta(x0,cutoff + alpha,n0 - cutoff + beta)
den = cbind(x0,y0)
threshold = c(0.3,0.35)
PP = 1-pbeta(threshold,cutoff + alpha,n0 - cutoff + beta)
}
p2 = my_density_threshold(den,threshold,PP,direction = "greater",
Criteria = "NoGo",
title1=paste("When #Obs/n = ",res$Final_NoGo,"/",size_true,":",sep=""),
title0=paste("NoGo when #Obs", res$df.output$`Cut off for NoGo`,":\n ",res$Decision_rule_for_NoGo,sep="") )
p1
p2