## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  message =FALSE,
  warning =FALSE,
  fig.width = 7,
  comment = "#>"
)
if (capabilities(("cairo"))) {
  knitr::opts_chunk$set(dev.args = list(png = list(type = "cairo")))
}

library(coveffectsplot)
library(ggplot2)
library(dplyr)
library(tidyr)
library(mrgsolve)
library(ggridges)
library(ggstance)
library(patchwork)
library(ggdist)
library(ggrepel)
library(Rcpp)
library(egg)
theme_set(theme_bw())
nsim <- 100 # for vignette to make it run faster otherwise increase to 1000
#utility function to simulate varying one covariate at a time keeping the rest at the reference
expand.modelframe <- function(..., rv, covcol="covname") {
  args <- list(...)
  df <- lapply(args, function(x) x[[1]])
  df[names(rv)] <- rv
  res <- lapply(seq_along(rv), function(i) {
    df[[covcol]] <- names(rv)[i]
    df[[names(rv)[i]]] <- args[[names(rv)[i]]]
    as.data.frame(df)
  })
  do.call(rbind, res)
}

#useful for logit/inverselogit transforms
expit <- function(x) exp(x)/ (1 + exp(x) )
 

## ----exprespmodel, collapse=TRUE----------------------------------------------
# the typical probability from the model parameters will be :
TypicalProb<- 1/(1+exp(-(log(0.1/(1-0.1)) + (5*75/10/(7.5+75/10)))))
MaxProb<- 1/(1+exp(-(log(0.1/(1-0.1)) + (5*750/10/(7.5+750/10)))))
MinProb<- 1/(1+exp(-(log(0.1/(1-0.1)) + (5*0/10/(7.5+0/10)))))
#
# R::runif(0,1)< P1 ? 1 : 0
exprespmodel <- '
$PLUGIN Rcpp
$PARAM @annotated
TVCL    : 10  : Clearance CL (L/h)
WTCL    : 0.75: Weight on CL (ref. 70 kg)
TVEMAX  : 5   : Maximum Drug Effect
SEVEMAX : 2.5   : Severity Reduction of Drug Effect
AUC50   : 7.5 : Area Under the Curve providing half maximal response
TVBASEP : 0.1 : Baseline Probability of Response 

$PARAM @annotated // reference values for covariate
WT      : 70  : Weight (kg)
SEV     : 0   : Sex (0=Female, 1=Male)
DOSE    : 75  : Dose (mg)
$OMEGA @annotated @block
nCL     : 0.16 : ETA on CL

$OMEGA @annotated @block
nInt    : 0.05 : ETA on CL

$PRED
double CL = TVCL *
    pow((WT/70.0), WTCL)*exp(ETA(1)); 
double EMAX = TVEMAX - SEVEMAX*(SEV == 1) ; 
double BASEP = TVBASEP *exp(ETA(2)); 
double Intercept = log(BASEP/(1-BASEP));
capture CLi = CL;
capture AUC = DOSE/CL;
capture LGST = Intercept + (EMAX*AUC/(AUC50+AUC));
capture P1 = 1/(1+exp(-LGST));
capture DV = R::runif(0,1)< P1 ? 1 : 0;
'

modexprespsim <- mcode("exprespmodel", exprespmodel)
simdata <-  expand.idata(SEV=c(0),
               DOSE = c(0,75),
               ID = 1:1000) %>% 
  dplyr::mutate(WT = 70) #exp(rnorm(n(),log(70),0.3)
set.seed(466548)
simout <- modexprespsim %>%
  data_set(simdata) %>%
  carry.out(WT, DOSE, SEV) %>%
  mrgsim()%>%
  as.data.frame


## ----exprespmodeplotl, collapse=TRUE, fig.height= 7, fig.width= 9-------------
WT_names <- c(
  '70'="Weight: 70 kg"
)
SEV_names <- c(
  '0'="Severity: Not Severe"
)
simout$SEV_cat <- "Not Severe"

vlines <- quantile(simout$AUC[simout$AUC>0],probs = c(0,1/4,0.5,3/4,1),na.rm = TRUE)


simoutdataprobsdose <- simout %>% 
  group_by(DOSE) %>% 
  summarise(probs = mean(DV),
            n = n(),
            lower =  probs - qnorm(1-0.05/2)*sqrt((1/n)*probs*(1-probs)),
            upper =  probs + qnorm(1-0.05/2)*sqrt((1/n)*probs*(1-probs)),
            n0 = length(DV[DV==0]),
            n1= length(DV[DV==1]),
            medAUC= median(AUC))

simoutdataprobs <- simout %>% 
  mutate(AUC_bins=table1::eqcut(AUC,4,withhold = list(PLACEBO=AUC==0))) %>% 
  group_by(AUC_bins) %>% 
  summarise(probs = mean(DV),
            n = n(),
            lower =  probs - qnorm(1-0.05/2)*sqrt((1/n)*probs*(1-probs)),
            upper =  probs + qnorm(1-0.05/2)*sqrt((1/n)*probs*(1-probs)),
            n0 = length(DV[DV==0]),
            n1= length(DV[DV==1]),
            medAUC= median(AUC))

percentineachbreakcategory <- simout %>% 
  mutate(AUC_bins=table1::eqcut(AUC,4,withhold = list(PLACEBO=AUC==0))) %>% 
  group_by(DOSE) %>% 
  mutate(Ntot= n())%>% 
  group_by(DOSE,AUC_bins,WT,SEV) %>% 
  mutate(Ncat=n(),xmed=median(AUC))%>% 
  mutate(percentage=Ncat/Ntot)%>% 
  distinct(DOSE,xmed,AUC_bins,percentage,WT,SEV)

probplot <- ggplot(simout, aes(AUC,DV,linetype=factor(SEV_cat))) +
  facet_grid( WT~SEV,labeller=labeller(WT=WT_names,SEV=SEV_names))+
    geom_vline(xintercept = vlines)+
  geom_point(position=position_jitter(height=0.02,width=0.1),
             aes(color=factor(DOSE)),size=1,alpha=0.2)+
  geom_smooth(method = "glm", 
    method.args = list(family = "binomial"), 
    se = TRUE) +
  geom_pointrange(data=simoutdataprobs,aes(medAUC,probs,ymin=lower,ymax=upper),
                  inherit.aes = FALSE,alpha=0.2)+
  geom_text(data=simoutdataprobs, aes(x=medAUC,y=probs,
                label=paste(round(100*probs,1),"%")),
                  inherit.aes = FALSE,hjust=-0.15,vjust=-0.15,size=4)+
  geom_point(data=simoutdataprobsdose %>% filter(DOSE!=0),
             aes(x=medAUC,y=probs,color=factor(DOSE)),
             inherit.aes = FALSE,size=4, shape="diamond")+
  geom_text_repel(data=simoutdataprobs, aes(x=medAUC,y=Inf,
                label=paste(n1,n,sep = "/")),
                  inherit.aes = FALSE,direction="y")+
  labs(color="Dose (mg)",y="Probability of Response",
       linetype="Severity")+
  theme_bw(base_size = 16) + 
  theme(legend.position = "top")

exposureplot <- ggplot(simout,
         aes(y = as.factor(DOSE),x = AUC,
             fill = as.factor(DOSE))) +
    facet_grid( WT~SEV,labeller=labeller(WT=WT_names,SEV=SEV_names))+
    geom_vline(xintercept = vlines)+
   geom_text(data=data.frame(vlines), aes(x=vlines,y=-Inf,
                label=paste(round(vlines,1))),
                  inherit.aes = FALSE,vjust=-0.1,hjust=1)+
  stat_slab(scale = 1, alpha= 0.9,
            aes(fill_ramp =  after_stat(
              ifelse(x<= vlines[2],"",
                     ifelse(x>vlines[2] & x <=vlines[3],"50%",
                            ifelse(x>vlines[3] & x<=vlines[4],"50%","")))
              )
              )
  )+
  stat_pointinterval(.width = c(.5,1))+
    geom_label(data=percentineachbreakcategory,
                   aes(x=xmed, label=round(100*percentage,0) ),alpha=0.5)+
     theme_bw(base_size = 16)+
    theme(legend.position="none",
          strip.background.x = element_blank(),
          strip.text.x = element_blank(),
          plot.background = element_blank())+
   labs(x = "AUC", y = "DOSE (mg)")+
  scale_y_discrete(breaks= c(0,75),labels=c("PLB"," 75"))
egg::ggarrange(
  (probplot+
  theme(axis.title.x.bottom = element_blank(),
        axis.text.x.bottom = element_blank(),
        axis.ticks.x =  element_blank(),
        plot.margin = unit(c(0,0,0,0), "cm"))) ,
  (exposureplot +
  theme(plot.margin = unit(c(0,0,0,0), "cm") )) ,heights = c(1,0.5))
   

## ----bsvrangeplot, collapse=TRUE----------------------------------------------
set.seed(466548)
simdata <-  expand.idata(SEV=c(0),
               DOSE = c(0,75),
               ID = 1:10000) %>% 
  dplyr::mutate(WT = 70) #exp(rnorm(n(),log(70),0.3)
simout <- modexprespsim %>%
  data_set(simdata) %>%
  carry.out(WT, DOSE, SEV) %>%
  mrgsim()%>%
  as.data.frame

simoutbsvref <- simout %>% 
  gather(paramname, paramvalue,P1)%>%
  group_by(paramname,DOSE)%>%
  dplyr::summarize(Pstdref = quantile(paramvalue, 0.5,na.rm =TRUE)) %>%
  ungroup() %>% 
  select(DOSE, Pstdref)
simout <- left_join(simout,simoutbsvref)
simoutbsvlong  <- simout  %>%
  group_by(DOSE)%>%
  mutate(P1std=P1/Pstdref) %>% 
  gather(paramname, paramvalue,P1std,P1)


yvar_names <- c(
  'P1std'="Standardized Probability",
  'P1'="Probability"
)

pbsvranges<-  ggplot(simoutbsvlong %>% 
                       filter(DOSE!=0), aes(
        x      = paramvalue,
        y      = paramname,
        fill   = factor(..quantile..),
        height = ..ndensity..)) +
  facet_wrap(paramname~DOSE , scales="free_y", ncol=1,
             labeller=labeller(paramname=yvar_names) ) +
  stat_density_ridges(
    geom="density_ridges_gradient", calc_ecdf=TRUE,
    quantile_lines=TRUE, rel_min_height=0.001, scale=0.9,
    quantiles=c(0.05, 0.25, 0.5, 0.75, 0.95)) +
  scale_fill_manual(
    name="Probability",
    values=c("white", "#FF000050", "#FF0000A0", "#FF0000A0", "#FF000050", "white"),
    labels = c("(0, 0.05]", "(0.05, 0.25]",
             "(0.25, 0.5]", "(0.5, 0.75]",
             "(0.75, 0.95]", "(0.95, 1]")) +
  theme_bw() +
  theme(
    legend.position = "none",
    axis.text.y     = element_blank(),
    axis.ticks.y    = element_blank(),
    axis.title.y    = element_blank()) +
  labs(x="Parameters", y="") +
  scale_x_log10(breaks=c(0.5,0.8,1,1.25,2)) +
  coord_cartesian(expand=FALSE,xlim=c(0.2,2))

simoutbsvranges <- simoutbsvlong %>%
  group_by(paramname)%>%
  dplyr::summarize(
   P05 = quantile(paramvalue, 0.05),
    P25 = quantile(paramvalue, 0.25),
    P50 = quantile(paramvalue, 0.5),
    P75 = quantile(paramvalue, 0.75),
    P95 = quantile(paramvalue, 0.95))
simoutbsvranges
pbsvranges

## ----collapse=TRUE------------------------------------------------------------
set.seed(678549)
thmeans <- c(10,0.75, #TVCL WTCL
             5,3, # TVEMAX  SEVEMAX
             7.5, # AUC50
              0.1) #BASEP
thvariances<- (thmeans*0.25)^2
thecorrelations <- matrix(ncol=length(thmeans),nrow=length(thmeans))
diag(thecorrelations)<- 1
thecorrelations[lower.tri(thecorrelations, diag = FALSE)]<- 0.2
thecorrelations[upper.tri(thecorrelations, diag = FALSE)]<- 0.2
thevarcovmatrix<- diag(sqrt(thvariances))%*%thecorrelations%*%diag(sqrt(thvariances))
sim_parameters <- MASS::mvrnorm(n = nsim, mu=as.numeric(thmeans),
                                Sigma=thevarcovmatrix, empirical = TRUE)
colnames(sim_parameters) <- colnames(thevarcovmatrix) <- c("TVCL","WTCL",
                                                           "TVEMAX","SEVEMAX","AUC50",
                                                           "BASEP")
sim_parameters<- as.data.frame(sim_parameters)

reference.values <- data.frame(WT = 70, DOSE = 75, SEV = 0 )   
covcomb <- expand.modelframe(
  WT  = c(50,60,70,80,90),
  DOSE = c(0,25,50,75,100,125,150),
  SEV = c(0,1),
  rv = reference.values)
covcomb <- covcomb[!duplicated(
  paste(covcomb$WT,covcomb$WT,covcomb$DOSE,covcomb$SEV)),]
covcomb$ID <- 1:nrow(covcomb)

iter_sims <- NULL
for(i in 1:nsim) {
  idata <- as.data.frame(covcomb)
  idata$covname<- NULL
  data.all <- idata
  data.all$TVCL <- as.numeric(sim_parameters[i,1])
  data.all$WTCL <- as.numeric(sim_parameters[i,2])
  data.all$TVEMAX   <- as.numeric(sim_parameters[i,3])
  data.all$SEVEMAX      <- as.numeric(sim_parameters[i,4])
  data.all$AUC50 <- as.numeric(sim_parameters[i,5])
  data.all$BASEP <- as.numeric(sim_parameters[i,6])
  out <- modexprespsim %>%
    data_set(data.all) %>%
    carry.out(CL,WT, DOSE, SEV, AUC) %>%
    zero_re() %>% 
    mrgsim()
  dfsimunc <- as.data.frame(out%>% mutate(rep = i) )
  iter_sims <- rbind(iter_sims,dfsimunc)
}

## ----collapse=TRUE------------------------------------------------------------

wt.labs <- c("weight: 50 kg","weight: 60 kg","weight: 70 kg","weight: 80 kg","weight: 90 kg","(all)")
names(wt.labs) <- c("50","60","70","80","90","(all)")

SEV_names <- c(
  '0'="Not Severe",
  '1'="Severe"

)

stdprobplot<- ggplot(iter_sims, aes(DOSE,P1,col=factor(SEV) ) )+
  geom_point(aes(group=interaction(ID,rep)),alpha=0.5,size=3)+
  geom_hline(yintercept=TypicalProb)+
  facet_grid(SEV~ WT,labeller = labeller(WT=wt.labs, SEV = SEV_names))+
  labs(y="Probability of Response", colour="Severity")
stdprobplot

## ----fig.height= 7, collapse=TRUE---------------------------------------------
TypicalProb <- iter_sims %>%
  filter(ID==3) %>% 
  summarize(medP1=median(P1)) %>% 
  pull(medP1)

iter_sims <- iter_sims %>%
  mutate(P1std=P1/TypicalProb)%>%
  gather(paramname,paramvalue,P1std)%>% 
  ungroup() %>% 
  dplyr::mutate( covname = case_when(
    ID== 1 ~ "Weight",
    ID== 2 ~ "Weight",
    ID== 3 ~ "REF",
    ID== 4 ~ "Weight",
    ID== 5 ~ "Weight",
    ID== 6 ~ "DOSE",
    ID== 7 ~ "DOSE",
    ID== 8 ~ "DOSE",
    ID== 9 ~ "DOSE",
    ID== 10 ~ "DOSE",
    ID== 11 ~ "DOSE",
    ID== 12 ~ "SEV"
  ),
  covvalue =case_when(
    ID== 1 ~ paste(WT,"kg"), 
    ID== 2 ~ paste(WT,"kg"),
    ID== 3 ~ "70 kg\nNot Severe\n75 mg",
    ID== 4 ~ paste(WT,"kg"),
    ID== 5 ~ paste(WT,"kg"),
    ID== 6 ~ paste(DOSE,"mg"),
    ID== 7 ~ paste(DOSE,"mg"),
    ID== 8 ~ paste(DOSE,"mg"),
    ID== 9 ~ paste(DOSE,"mg"),
    ID== 10 ~ paste(DOSE,"mg"),
    ID== 11 ~ paste(DOSE,"mg"),
    ID== 12 ~ "Severe"
  ) )
iter_sims$covname <-factor(as.factor(iter_sims$covname ),
                          levels =  c("Weight","DOSE","SEV","REF"))
iter_sims$covvalue <- factor(as.factor(iter_sims$covvalue),
                          levels =  c("0 mg","25 mg","50 mg",
                                      "100 mg","125 mg","150 mg",
                          "50 kg","60 kg","80 kg", "90 kg",
                          "70 kg\nNot Severe\n75 mg",  "Severe"))


ggplot(iter_sims,aes(x=paramvalue,y=covvalue))+
  stat_density_ridges(aes(fill=factor(..quantile..),height=..ndensity..),
    geom = "density_ridges_gradient", calc_ecdf = TRUE,
    quantile_lines = TRUE, rel_min_height = 0.001,scale=0.9,
    quantiles = c(0.025,0.5, 0.975))+
    facet_grid(covname~paramname,scales="free",switch="both",
             labeller = labeller(paramname=yvar_names))+ 
  scale_fill_manual(
    name = "Probability", values = c("white","#0000FFA0", "#0000FFA0","white"),
    labels = c("(0, 0.025]","(0.025, 0.5]","(0.5, 0.975]","(0.975, 1]")
  )+
  theme_bw()+
  theme(axis.title = element_blank(),strip.placement = "outside")

## ----plot3, collapse=TRUE, fig.height=7, fig.width=9--------------------------
coveffectsdatacovrep <- iter_sims %>%
  dplyr::group_by(paramname,ID,WT,DOSE,SEV,covname,covvalue) %>% 
  dplyr::summarize(
    mid= median(paramvalue),
    lower= quantile(paramvalue,0.025),
    upper = quantile(paramvalue,0.975))%>% 
  dplyr::filter(!is.na(mid))
simoutbsvranges<-simoutbsvranges[simoutbsvranges$paramname=="P1std",]

coveffectsdatacovrepbsv <- coveffectsdatacovrep[coveffectsdatacovrep$covname=="REF",]
coveffectsdatacovrepbsv$covname <- "BSV"
coveffectsdatacovrepbsv$covvalue <- "90% of patients"
coveffectsdatacovrepbsv$label <-    "90% of patients"
coveffectsdatacovrepbsv$lower <- simoutbsvranges$P05
coveffectsdatacovrepbsv$upper <- simoutbsvranges$P95
coveffectsdatacovrepbsv2 <- coveffectsdatacovrep[coveffectsdatacovrep$covname=="REF",]
coveffectsdatacovrepbsv2$covname <- "BSV"
coveffectsdatacovrepbsv2$covvalue <- "50% of patients"
coveffectsdatacovrepbsv2$label <-    "50% of patients"
coveffectsdatacovrepbsv2$lower <- simoutbsvranges$P25
coveffectsdatacovrepbsv2$upper <- simoutbsvranges$P75

coveffectsdatacovrepbsv<- rbind(coveffectsdatacovrep,coveffectsdatacovrepbsv2,
                                coveffectsdatacovrepbsv)

coveffectsdatacovrepbsv <- coveffectsdatacovrepbsv %>% 
  mutate(
    label= covvalue,
    LABEL = paste0(format(round(mid,2), nsmall = 2),
                   " [", format(round(lower,2), nsmall = 2), "-",
                   format(round(upper,2), nsmall = 2), "]"))
coveffectsdatacovrepbsv<- as.data.frame(coveffectsdatacovrepbsv)

coveffectsdatacovrepbsv$label <-factor(as.factor(coveffectsdatacovrepbsv$label ),
            levels = c("All Subjects","90% of patients","50% of patients",
              "50 kg","60 kg","80 kg","90 kg",
              "0 mg","25 mg","50 mg","100 mg","125 mg","150 mg",
              "Severe","70 kg\nNot Severe\n75 mg"
              ))
coveffectsdatacovrepbsv$covname <-factor(as.factor(coveffectsdatacovrepbsv$covname ),
                    levels = c("Weight","DOSE","SEV","REF","BSV"))


ref_legend_text      <- "Reference (vertical line)"
#png("./Figure_8_4.png",width =9 ,height = 7,units = "in",res=72)
fpp <- forest_plot(coveffectsdatacovrepbsv,
                            strip_placement = "outside",
                            show_ref_area = FALSE,
                            show_ref_value=TRUE,
                            ref_legend_text = ref_legend_text,
                            plot_table_ratio = 2,
                            interval_shape = "circle small",
                            bsv_shape = "triangle",
                            legend_order = c("pointinterval","shape","area","ref"),
                            combine_interval_shape_legend = TRUE,
                            base_size = 12,
                            table_text_size = 4,
                            y_label_text_size      = 12,
                            xlabel= " ",
                            facet_formula = "covname~paramname",
                            facet_labeller      = labeller(paramname=yvar_names),
                            facet_scales           = "free",
                            facet_space ="free",
                            logxscale              = TRUE,
                            major_x_ticks          = c(0.1,0.25, 0.5,1,1.5),
                            x_range                = c(0.1, 1.5))

#dev.off()
#![Covariate Effects Plot.](./Figure_8_4.png)

## ----plot4, collapse=TRUE, fig.height= 6, fig.width= 12-----------------------
covcombdr <- expand.grid(
  WT  = c(50,70,90),
  DOSE = c(0,25,50,75,100,150),
  SEV = c(0,1)
)

covcombdr <- covcombdr[!duplicated(
  paste(covcombdr$WT,covcombdr$WT,covcombdr$DOSE,covcombdr$SEV)),]
covcombdr$ID <- 1:nrow(covcombdr)

iter_sims <- NULL
for(i in 1:nsim) {
  idata <- as.data.frame(covcombdr)
  idata$covname<- NULL
  data.all <- idata
  data.all$TVCL     <- as.numeric(sim_parameters[i,1])
  data.all$WTCL     <- as.numeric(sim_parameters[i,2])
  data.all$TVEMAX   <- as.numeric(sim_parameters[i,3])
  data.all$SEVEMAX  <- as.numeric(sim_parameters[i,4])
  data.all$AUC50    <- as.numeric(sim_parameters[i,5])
  data.all$BASEP    <- as.numeric(sim_parameters[i,6])
  out <- modexprespsim %>%
    data_set(data.all) %>%
    carry.out(CL,WT, DOSE, SEV, AUC) %>%
    zero_re() %>% 
    mrgsim()
  dfsimunc <- as.data.frame(out%>% mutate(rep = i) )
  iter_sims <- rbind(iter_sims,dfsimunc)
}
iter_sims <- data.table(iter_sims)

## ----plot4b, collapse=TRUE, fig.height= 6, fig.width= 12----------------------
summary_P1 <- function(P1) {
  x <- c(
    P1med = median(P1),
    P1low = quantile(P1, probs = 0.05),
    P1up =  quantile(P1, probs = 0.95)
  )
  data.table(paramname=names(x), paramvalue=x)
}

iter_sims_sum <- iter_sims[, summary_P1(P1), by=.( DOSE, WT, SEV)]
iter_sims_sum <- spread(iter_sims_sum,paramname,paramvalue)
iter_sims_sum <- as.data.frame(iter_sims_sum)
wt.labs <- c("weight: 50 kg","weight: 60 kg","weight: 70 kg","weight: 80 kg","weight: 90 kg","(all)")
names(wt.labs) <- c("50","60","70","80","90","(all)")

SEV_names <- c(
  '0'="Not Severe",
  '1'="Severe"

)

iter_sims_sum_ref <- iter_sims_sum %>% 
  filter(WT==70,SEV==0)
iter_sims_sum_ref$WT <- NULL
iter_sims_sum_ref$SEV <- NULL

plot1 <- ggplot(iter_sims_sum, aes(DOSE,P1med,col=factor(SEV) ) )+
  geom_ribbon(data=iter_sims_sum_ref,aes(x=DOSE,
                         ymin =`P1low.5%`,
                         ymax =`P1up.95%`),inherit.aes = FALSE,
              fill = "gray", alpha=0.2, color="transparent")+
    geom_line(data=iter_sims_sum_ref,aes(x=DOSE,
                         y =P1med),inherit.aes = FALSE,
              color = "black", alpha=1)+
  geom_pointrange(aes(group=interaction(WT,SEV,DOSE),
                         ymin =`P1low.5%`,
                         ymax =`P1up.95%`,
                      shape = as.factor(WT)),alpha=0.5,
                  position = position_dodge(width=10))+
  geom_hline(yintercept=TypicalProb)+
  facet_grid(SEV~ WT,labeller = labeller(WT=wt.labs, SEV = SEV_names))+
  theme_bw(base_size = 14)+
  labs(color="Severity",y = "Probability of Response", x = "Dose (mg)",
       shape = "Weight (kg)")+
  scale_x_continuous(breaks=unique(iter_sims_sum$DOSE),guide = guide_axis(n.dodge = 2))+
  theme(legend.position = "none")




plot2 <- plot1 +
  facet_grid(SEV~ .,labeller = labeller(WT=wt.labs, SEV = SEV_names))+
  theme(axis.title.y.left = element_blank(),
        axis.ticks.y.left = element_blank(),
        axis.text.y.left = element_blank(),legend.position = "right")

egg::ggarrange( plot1 ,plot2 ,widths = c(1,0.3),
                bottom = "The black line and 95%CI gray area reference curve for Weight = 70 kg and Not Severe\nis kept in the background of all panels" )


## ----plot5, collapse=TRUE, fig.height= 6, fig.width= 12-----------------------

vlines <- quantile(iter_sims$AUC[iter_sims$AUC>0],probs = c(0,1/3,2/3,1),na.rm = TRUE)
iter_sims$PLACEBO <- ifelse(iter_sims$DOSE>0,"Active","PLACEBO")
simoutdataprobs <- iter_sims %>% 
  mutate(AUC_bins=table1::eqcut(AUC,3,withhold = list(PLACEBO=AUC==0))) %>% 
  group_by(PLACEBO, WT,SEV,AUC_bins) %>% 
  summarise(probs = mean(DV),
            n = n(),
            lower =  probs - qnorm(1-0.05/2)*sqrt((1/n)*probs*(1-probs)),
            upper =  probs + qnorm(1-0.05/2)*sqrt((1/n)*probs*(1-probs)),
            n0 = length(DV[DV==0]),
            n1= length(DV[DV==1]),
            medAUC= median(AUC))

simoutdataprobsdose <- iter_sims %>% 
  group_by(PLACEBO, WT,SEV,DOSE) %>% 
  summarise(probs = mean(DV),
            n = n(),
            lower =  probs - qnorm(1-0.05/2)*sqrt((1/n)*probs*(1-probs)),
            upper =  probs + qnorm(1-0.05/2)*sqrt((1/n)*probs*(1-probs)),
            n0 = length(DV[DV==0]),
            n1= length(DV[DV==1]),
            medAUC= median(AUC))


SEV_names <- c(
  '0'="Not Severe",
  '1'="Severe"

)
wt.labs <- c("weight: 50 kg","weight: 70 kg","weight: 90 kg")
names(wt.labs) <- c("50","70","90")



a <- ggplot(data=data.frame(iter_sims),aes(AUC, P1,
                                           linetype = as.factor(SEV),
                                           color = as.factor(SEV),
                                           fill = as.factor(SEV),
                                           group=interaction(WT,SEV)))+
    facet_grid(.~ WT,labeller = labeller(WT=wt.labs, SEV = SEV_names))+
  geom_vline(xintercept = vlines)+
  geom_smooth(method = "glm", method.args = list(family = "binomial")
              ) +
   theme_bw(base_size = 14)+
  theme(axis.title.x.bottom = element_blank(),
        axis.text.x.bottom = element_blank(),
         axis.ticks.x =  element_blank(),
        legend.position = "none")+
  labs(y="Probability of response",linetype="Severity",
       color = "Weight (kg)", fill = "Weight (kg)")
 b <-  ggplot(iter_sims,
         aes(y = as.factor(DOSE),x = AUC,fill = as.factor(SEV))) +
       facet_grid(.~ WT,labeller = labeller(WT=wt.labs, SEV = SEV_names))+
    geom_vline(xintercept = vlines)+
   stat_slab(scale = 4, alpha=0.5,position = "dodgejust")+
   stat_pointinterval(.width = c(.5,1),position = "dodgejust")+
     theme_bw(base_size = 14)+
    theme(legend.position="none",
          strip.background.x = element_blank(),
          strip.text.x = element_blank())+
   labs(x = "AUC", y = "DOSE (mg)")
egg::ggarrange(a,b,heights=c(1,0.6))


## ----plot667, collapse=TRUE, include=TRUE, fig.height= 16, fig.width= 12------

percentineachbreakcategory <- iter_sims %>% 
  mutate(AUC_bins=table1::eqcut(AUC,3,withhold = list(PLACEBO=AUC==0))) %>% 
  group_by(DOSE, WT,SEV) %>% 
  mutate(Ntot= n())%>% 
  group_by(DOSE,WT,SEV,AUC_bins) %>% 
  mutate(Ncat=n(),xmed=median(AUC))%>% 
  mutate(percentage=Ncat/Ntot)%>% 
  distinct(DOSE,xmed,AUC_bins,percentage,WT,SEV)


probplot<- ggplot(iter_sims , aes(AUC,DV)) +
  facet_grid(SEV ~WT,labeller=labeller(WT=wt.labs,SEV=SEV_names))+
  geom_smooth(data=iter_sims %>%
                filter(WT==70,SEV ==0) %>% 
                mutate(SEV=NULL,WT=NULL),method = "glm", 
              method.args = list(family = "binomial"), 
              se = TRUE,color="white",fill="lightgray",alpha=0.3) +
  geom_vline(xintercept = vlines)+
  geom_point(position=position_jitter(height=0.02,width=0.1),size=1,alpha=0.2,
             aes(colour=factor(DOSE)))+
  geom_smooth(method = "glm", 
              method.args = list(family = "binomial"), 
              se = TRUE,color="black",fill="blue",alpha=0.5,
              aes(linetype=as.factor(SEV))) +
  geom_pointrange(data=simoutdataprobs,aes(medAUC,probs,ymin=lower,ymax=upper),
                  inherit.aes = TRUE,alpha=0.5)+
  geom_text(data=simoutdataprobs, aes(x=medAUC,y=probs,
                                      label=paste(paste(round(100*probs,1),"%")
                           #, paste(n1,n,sep = "/"),sep = "\n"
                                            )),
            inherit.aes = TRUE,hjust=-0.15,vjust=-0.15,size=4)+
  geom_point(data=simoutdataprobsdose, aes(x=medAUC,y=probs,colour=factor(DOSE)),
            inherit.aes = FALSE,shape="diamond",size=4,alpha=0.3)+
  theme_bw(base_size = 16)+ 
  theme(legend.position = "top")+
  labs(color="Dose (mg)",y="Probability of Response",
       linetype="Severity")


explot<- ggplot(iter_sims,aes(y = as.factor(DOSE),x = AUC,
                           fill = as.factor(DOSE))) +
  facet_grid(SEV ~WT,labeller=labeller(WT=wt.labs,SEV=SEV_names))+
  geom_vline(xintercept = vlines)+
  geom_text(data=data.frame(vlines), aes(x=vlines,y=-Inf,
                                         label=paste(round(vlines,1))),
            inherit.aes = FALSE,vjust=-0.1,hjust=1)+
  stat_slab( aes(fill_ramp =  after_stat(
              ifelse(x<= vlines[2],"lower",
                     ifelse(x>vlines[2] & x <=vlines[3],"middle",
                            ifelse(x>vlines[3],"upper","")))
            )
            ),scale=1,alpha=0.8,normalize ="groups"
  )+
  stat_pointinterval(.width = c(.5,1))+
  geom_label_repel(data=percentineachbreakcategory,
             aes(x=xmed, label=round(100*percentage,0) ),alpha=0.5,direction = "y")+
  theme_bw(base_size = 16)+
  theme(legend.position="none",
        strip.background.x = element_blank(),
        strip.text.x = element_blank(),
        plot.background = element_blank())+
  labs(x = "AUC", y = "DOSE (mg)")

egg::ggarrange(probplot +
    theme(axis.title.x.bottom = element_blank(),
          axis.text.x.bottom = element_blank(),
          axis.ticks.x =  element_blank() ) , explot,
  heights=c(1,2))