demo/Simulation2026_discrete_summary.R

In bpcp: Beta Product Confidence Procedure for Right Censored Data

#' ---
#' title: "Discrete Simulation"
#' author: "Allyson Mateja and Megan Grieco and Mike Fay"
#' date: "`r Sys.Date()`"
#' output: word_document
#' ---
#' 
## ----setup, include=FALSE------------------------------------------------------------------------------
#knitr::opts_chunk$set(echo = FALSE, message=F)
library(dplyr)
library(tidyr)
library(data.table)
library(abind)
library(ggplot2)

#' 
#' 
#' ## Simple Mixture of Discrete  Distribution
#' 
#' There are three types of individuals:
#' 
#' * always fail: have event at $t=1$
#' * effected: if treated have event at $t=4$, otherwise have event at $t=2$.
#' * never fail: have event  after $t=6$ (study end)
#' 
#' Let the proportions in the three populations be $p_A$, $p_E$, and $p_N$, respectively.
#' We have $p_A+p_E+p_N=1$, so we only need to define $p_A$ and $p_E$. 
#' 
#' There is independent censoring: where $\pi_C$ proportion of the study population
#' are censored at $t=3$ and the rest are censored at the study end, $t=6$. 
#' 
#' We test for a difference between survival distributions at $t=5$. This is a 
#' test of the type I error rate. 
#' 
#' The true survival curves are equal after $t=4$. Simulated data will have n=300 in each arm, $p_A=.5$ and $\pi_c=0.9$. We will vary $p_E$ from 0 to 0.4 for a total of 41 scenarios, each with 10,000 reps.
#' 
#' ## Simulation Results
#' 
#' 
## ------------------------------------------------------------------------------------------------------
#sim_files <- list.files("C:/Users/matejaam/OneDrive - National Institutes of Health/_HDrive/Projects/bpcp/simulations/discrete/sim_discrete", pattern=".rds", full.names = T) %>% purrr::map(readRDS)
#sim_files_diff <- list.files("C:/Users/matejaam/OneDrive - National Institutes of Health/_HDrive/Projects/bpcp/simulations/discrete/sim_discrete_diff", pattern=".rds", full.names = T) %>% purrr::map(readRDS)

#' 
## ------------------------------------------------------------------------------------------------------
summarizeCI<-function(ci,nullBeta=0){
  N<- nrow(ci)
  N.lo<- sum(!is.na(ci[,1]))
  N.hi<- sum(!is.na(ci[,2]))
  prop.Res.lo<- N.lo/N
  prop.Res.hi<- N.hi/N
  reject.lo<- sum(!is.na(ci[,1]) & ci[,1]>nullBeta)/N.lo
  reject.hi<- sum(!is.na(ci[,1]) & ci[,2]<nullBeta)/N.hi
  coverage.lo <- sum(!is.na(ci[,1]) & ci[,1]<nullBeta)/N.lo
  coverage.hi <- sum(!is.na(ci[,2]) & ci[,2]>nullBeta)/N.lo
  rejectTotal.lo<- sum(!is.na(ci[,1]) & ci[,1]>nullBeta)/N
  rejectTotal.hi<- sum(!is.na(ci[,1]) & ci[,2]<nullBeta)/N
  out<-c(N=N,N.lo=N.lo,N.hi=N.hi,prop.Res.lo=prop.Res.lo,
         prop.Res.hi=prop.Res.hi,reject.lo=reject.lo,
         reject.hi=reject.hi,
         coverage.lo, coverage.hi,
         rejectTotal.lo=rejectTotal.lo,
         rejectTotal.hi=rejectTotal.hi)
  out
}


#summarizeCI(out[,5:6,1])
#' 
#' 
## ------------------------------------------------------------------------------------------------------
pE <- seq(0, 0.4, 0.01)

#' 
#' 
## ------------------------------------------------------------------------------------------------------
#combined <- abind(sim_files, along=1)
#combined_diff <- abind(sim_files_diff, along=1)
combined<- readRDS("./simBPCP/simResults/efflogs2026_discrete_nsim1000.rds")
#' 
#' 
## ------------------------------------------------------------------------------------------------------
out<- array(NA,c(4,11,length(pE)),
            dimnames=list(c("meld","meld.midp","delta.zo","delta.adjhy"),
                          c("N","N.lo","N.hi","prop.Res.lo","prop.Res.hi","reject.lo","reject.hi", "coverage.lo", "coverage.hi","rejectTotal.lo","rejectTotal.hi"),
                          c(paste0("Scenario ",1:length(pE))))
)

for (i in 1:length(pE)) {
  
  CImat <- combined[,,i]
  sout<-matrix(NA,ncol(combined)/2,11)

  s1<-summarizeCI(CImat[,1:2])
  sout[1,]<- s1 
  sout[2,]<- summarizeCI(CImat[,3:4])
  sout[3,]<- summarizeCI(CImat[,5:6])
  sout[4,]<- summarizeCI(CImat[,7:8])
    
  dimnames(sout) <- list(c("meld","meld.midp","delta.zo","delta.adjhy"),
      names(s1))
  
  out[,,i]<- sout
  
}

#out[,c(1,6,7),]

#' 
#' 
#' 
## ------------------------------------------------------------------------------------------------------
# make 3D output into data frame
out_df <- as.data.frame.table(out[,c(1,6:9),])
names(out_df) <- c("method","calc","Scenario","res")

# calculated expected # at risk - n*(1-pA-pE)*(1-pic)
scenario_calcs <- data.frame(
  pE=pE,
  pA=0.5,
  n=300,
  pi.c=0.9,
  Scenario=paste0("Scenario ",seq(1:length(pE)))
) %>% mutate(n.risk=n*(1-pA-pE)*(1-pi.c))

# combine with full results data
out_df_2 <- out_df %>% left_join(scenario_calcs)

#out_df_2$method <- factor(out_df_2$method, levels = c("meld", "meld.midp", "delta.zo", "delta.adjhy"), 
#                          labels = c("Melding", "Melding (mid-p)", "Delta (Greenwood)", "Delta (Borkowf)"))

out_df_2$method <- factor(out_df_2$method, levels = c("meld", "meld.midp", "delta.zo", "delta.adjhy"), 
                          labels = c("meld BPCP", "meld mid-p BPCP", "Delta (standard, Z-O)", "Delta (adj hybrid, Z-O)"))

#' 
#' The following plots show the simulation results for each method and by lower or upper-tailed tests. A red dashed line is drawn at 0.025 to correspond to a two-sided 95% confidence interval. All melded Type I error rates remained below 2.5%.
#' 
#' The first plot shows data by pE. The second plot shows data by the expected number at risk in each arm just before the test time, $n*(1-p_A-p_E)*(1-\pi_c)$. As to be expected, $p_E$ and expected number of risk have an inverse relationship in terms of Type I error rate. 
#' 
## ----fig.height=6--------------------------------------------------------------------------------------
# Option 1: plot by pE (facet by method and low/high)
#ggplot(out_df_2 %>% filter(calc!="N" & method !="delta.none" & !(grepl("coverage", calc))) %>% mutate(calc=ifelse(calc=="reject.lo","Lower","Upper")), aes(x=pE, y=res)) + geom_hline(yintercept=0.025, linetype="dashed",color="red") + geom_point() + facet_grid(method~calc, ) + theme_bw() + theme(strip.background = element_blank()) + labs(y="Type I Error Rate", x="Proportion of Effected Individuals")
#ggsave("C:/R/work/simBPCP/simResults/pE_plot.pdf", height=6, width=5)

#' 
## ----fig.height=6--------------------------------------------------------------------------------------
# Option 2: plot by number at risk (facet by method and low/high)
ggplot(out_df_2 %>% filter(calc!="N" & method !="delta.none" & !(grepl("coverage", calc))) %>% mutate(calc=ifelse(calc=="reject.lo","Lower","Upper")), aes(x=n.risk, y=res)) + geom_hline(yintercept=0.025, linetype="dashed",color="red") + geom_point() + facet_grid(method~calc) + theme_bw() + theme(strip.background = element_blank()) + labs(y="Type I Error Rate", x="Expected Number at Risk")
ggsave("C:/R/work/simBPCP/simResults/num_risk_plot.pdf", height=6, width=5)