## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(survobj)
library(survival)

## ----simulation1, fig.align='center', fig.width= 7, fig.height=5--------------

# Number of simulations
nsim = 1000

# Participants in each group
nsubjects = 250

# Vaccine efficacy
ve = 40

# Hazard ratio
hr = 1-ve/100

# Follow-up time
ftime <- 12

# Fail events in controls 
fail_control = 0.4

# Define Object with exponential distribution for events in controls
s_events <- s_exponential(fail = fail_control, t = ftime)


## ----simulation, eval= FALSE--------------------------------------------------
# set.seed(12345)
# 
# # Define the group for the subjects
# group = c(rep(0, nsubjects), rep(1, nsubjects))
# 
# # Define the hazard ratio according to the group
# hr_vector <- ifelse(group ==0,1,hr)
# 
# # Loop
# sim <- lapply(
#   1:nsim,
#   function(x){
#     # Simulate survival times for event
#     sim_time_event <- s_events$rsurvhr(hr_vector)
# 
#     # Censor events at end of follow-up.
#     cevent <- censor_event(censor_time = ftime, time = sim_time_event, event = 1)
#     ctime <- censor_time(censor_time = ftime, time = sim_time_event)
# 
#     # Analyze the data using cox regression
#     reg <- summary(coxph(Surv(ctime, cevent)~ group))
# 
#     # Collect the information
#     pval = reg$coefficients["group","Pr(>|z|)"]
#     ve = (1- exp(reg$coefficients["group","coef"]))*100
#     nevents = reg$nevent
# 
#     # return values
#     return(data.frame(simid = x, pval,ve, nevents))
#   }
# )
# 
# # Join all the simulations in a single data frame
# sim_df <- do.call(rbind, sim)

## ----loadsimul, include=FALSE-------------------------------------------------
 # The simulation takes too much time to be included in CRAN
 # Load a previous simulation
 load("sim_df.rda")

## ----analyze------------------------------------------------------------------
empirical_power = binom.test(sum(sim_df$pval <= 0.05), length(sim_df$pval))
empirical_power$estimate
empirical_power$conf.int

# Distribution of the simulated VEs
summary(sim_df$ve)

# Distribution of the simulated number of events
summary(sim_df$nevents)