## A good introduction to sample size /power calculation @ https://bvag.com.vn/wp-content/uploads/2013/01/k2_attachments_CT-Lecture-23.-Sample-size---part-2.pdf 
   # slide 44 onwards: Estimation of sample size using R
   # slide 58 onwards: Sample size for time to an event
## R Package gsDesign::nSurvival() is used to calculate the sample size for a clinical trial with a time-to-event endpoint 
## nSurvival(): Time-to-event sample size calculation @ https://rdrr.io/cran/gsDesign/man/nSurvival.html (including examples)

if (!require("gsDesign")) install.packages("gsDesign")   # will install "gsDesign" if it is not already installed, and then load it.

# Generic function to calculate the sample size: 
# nSurvival(lambda1 = 1/12, lambda2 = 1/24,	# lambda1 & 2: event hazard rate for control/placebo and intervention/treatment group, respectively
#           Ts = 24, Tr = 12,               	# Ts: maximum study duration; Tr: accrual (recruitment) duration 
#           eta = 0, 				# eta: (equal) dropout hazard rate for the two groups (0 to 1)
#           ratio = 1,				# ratio: randomization ratio between placebo and treatment group. Default is balanced design, i.e., randomization ratio is 1
#           alpha = 0.025, 			# alpha: type I error rate. Default is 0.025 since 1-sided testing is default (0.05 if two-sided test is used: sided = 2)
#           beta = 0.10, 			# beta: type II error rate. (1 - beta) = statistical power (0.90 in this case)
#           sided = 1, 				# one or two-sided test? Default is one-sided test (sided = 1). For two-sided test: sided = 2
#           approx = FALSE,			# logical (TRUE or FALSE). If TRUE, the approximation sample size formula for risk difference is used.
#           type = c("rr", "rd"), 		# type of sample size calculation: risk ratio ("rr") or risk difference ("rd"). Choose one as: type = "rr" or type = "rd"
#           entry = c("unif", "expo"), 		# patient entry type: uniform entry ("unif") or exponential entry ("expo"). Choose one as: entry = "unif" or entry = "expo" 
#           gamma = NA 				# rate parameter for exponential entry. NA if entry type is "unif" (uniform). A non-zero value is supplied if entry type is "expo" (exponential)
#	    )                			# close parantheses

# An example:
library(gsDesign)
nSurvival(lambda1 = 1/24, lambda2 = 1/48, Ts = 12, Tr = 6, eta = 0, ratio = 1,
      alpha = 0.05, beta = 0.20, sided = 2, approx = FALSE,
      type = "rr", entry = "unif", gamma = NA)
# See also the examples @ https://rdrr.io/cran/gsDesign/man/nSurvival.html

# Result and interpretation: 
# Fixed design, two-arm trial with time-to-event outcome (Lachin and Foulkes, 1986)
Study duration (fixed):          Ts=12
Accrual duration (fixed):        Tr=6
Uniform accrual:              entry="unif", therefore, gamma = NA
Control median:      log(2)/lambda1=16.6
Experimental median: log(2)/lambda2=33.3
Censoring only at study end (eta=0)
Control failure rate:       lambda1=0.042
Experimental failure rate:  lambda2=0.021
Censoring rate:                 eta=0
Power:                 100*(1-beta)=80%
Type I error (2-sided):   100*alpha=5%
Equal randomization:          ratio=1
Sample size based on hazard ratio=0.5 (type="rr")
Sample size (computed):           n=278             # sample size required to have the desired power (0.90 in this example)
Events required (computed): nEvents=67              # Number of events required to achieve the power 

## To double-check the result, use the R package "epiR":
install.packages("epiR")
library(epiR)
# Function: epi.sscomps() replacing epi.studysize() in the most recent (2021) version of epiR. It is used to detect the required number of events for a given power level. 
  # See: https://cran.r-project.org/web/packages/epiR/epiR.pdf  &   https://rdrr.io/cran/epiR/man/epi.sscomps.html 
  # Note that: The argument "treat" is the proportion of treated subjects that will have not experienced the event of interest 
    # at the end of the study period and "control" is the proportion of control subjects that
    # will have not experienced the event of interest at the end of the study period. 
# Generic function: 
  # epi.sscomps(treat, 			# expected survival probability of patients receiving the intervention (treatment group)
		control, 		# expected survival probability of patients receiving a standard treatment (control group)
		n = NA, 		# for sample size determination, leave it as n = NA
		power = 0.80, 		# desired level of statistical power (usually 0.80)
		r = 1, 			# ratio of the sample size 
		design = 1, 
		sided.test = 2, 
		nfractional = FALSE, 
		conf.level = 0.95
		)
# Value (output)
  # n.crude:  the crude estimated total number of events required for the specified level of confidence and power
  # n.total:  the total estimated number of events required for the specified level of confidence and power, respecting the requirement for r times as many events in the treatment group compared with the control group
  # hazard:  the minimum detectable hazard ratio >1 and the maximum detectable hazard ratio <1 
  # power:  the power of the study given the number of events, the expected hazard ratio and level of confidence

epi.sscomps(0.30, 0.20, n = NA, 0.80, r = 1, design = 1, sided.test = 2, nfractional = FALSE, conf.level = 0.95)
## Result: 
   # n.crude
   # [1] 373
   # n.total
   # [1] 374
   # n.treat
   # [1] 187	# meaning a minimum of 187 events are required in the treatment group
   # n.control
   # [1] 187	# meaning a minimum of 187 events are required in the control group

## Examples from epiR reference manual @ https://rdrr.io/cran/epiR/man/epi.sscomps.html
 ## EXAMPLE 1 (from Therneau and Grambsch 2000 p. 63):
   ## The 5-year survival probability of patients receiving a standard treatment is 0.30 and we anticipate that a new treatment will increase it to 0.45.
   ## Assume that a study will use a two-sided test at the 0.05 level with 0.90 power to detect this difference. How many events are required?
epi.sscomps(treat = 0.45, control = 0.30, n = NA, power = 0.90, r = 1, design = 1, sided.test = 2, nfractional = FALSE, conf.level = 0.95)
   ## Result: A total of 250 events are required. Assuming one event per individual assign 125 individuals to the treatment group and 125 to the control group.

 ## EXAMPLE 2 (from Therneau and Grambsch 2000 p. 63):
   ## What is the minimum detectable hazard in a study involving 500 subjects where the treatment to control ratio is 1:1, assuming a power of 0.90 and a 2-sided test at the 0.05 level?
epi.sscomps(treat = NA, control = NA, n = 500, power = 0.90, r = 1, design = 1, sided.test = 2, nfractional = FALSE, conf.level = 0.95)
   ## Result: Assuming treatment increases time to event (compared with controls), the minimum detectable hazard of a study involving 500 subjects 
   ## (250 in the treatment group and 250 in the controls) is 1.33.

## FROM: https://rdrr.io/cran/gsDesign/man/nSurvival.html 
library(gsDesign)
library(ggplot2)
# consider a trial with 2 year maximum follow-up (Ts = 2) and 6 month uniform enrollment (Tr = 0.5)
# intervention/control hazard rates = 0.1/0.2 per 1 person-year
# dropout rate (eta) =  0.1 per 1 person-year
# alpha = 0.025 (1-sided)
# power = 0.9 (default beta = 0.1)

ss <- nSurvival(lambda1 = .2, lambda2 = .1, eta = .1, Ts = 2, Tr = .5, sided = 1, alpha = .025)
ss
# Results: 
# Fixed design, two-arm trial with time-to-event outcome (Lachin and Foulkes, 1986)
# Study duration (fixed):          Ts=2
# Accrual duration (fixed):        Tr=0.5
# Uniform accrual:              entry="unif"
# Control median:      log(2)/lambda1=3.5
# Experimental median: log(2)/lambda2=6.9
# Censoring median:        log(2)/eta=6.9
# Control failure rate:       lambda1=0.2
# Experimental failure rate:  lambda2=0.1
# Censoring rate:                 eta=0.1  (dropout rate = 0.1; equal in both arms)
# Power:                 100*(1-beta)=90%
# Type I error (1-sided):   100*alpha=2.5%
# Equal randomization:          ratio=1
# Sample size based on hazard ratio=0.5 (type="rr")
# Sample size (computed):           n=430
# Events required (computed): nEvents=91

## An example:
   # Slide 58 @ https://bvag.com.vn/wp-content/uploads/2013/01/k2_attachments_CT-Lecture-23.-Sample-size---part-2.pdf
   # More examples in following slides
#  • Package gsDesign with function nSurvival()
#  • Rate of death in group 1: 0.30 (lambda1)
#  • Rate of death in group 2: 0.20 (lambda2)
#  • Drop-out rate: 0.10 (eta; also known as censoring rate because dropouts are censored at the point of dropout)
#  • Duration of follow-up: 3 years (Ts - Time-study)
#  • Duration of recruitment: 1 year (Tr - Time-recruitment/accrual)

library(gsDesign)
nSurvival(lambda1 = 0.3, lambda2 = 0.20, eta = 0.1, Ts = 3, Tr = 1, alpha = 0.05, beta = 0.20) 

## Further advanced analyses (optional) 
#  group sequential translation with default bounds
#  note that delta1 is log hazard ratio; used later in gsBoundSummary summary
x <- gsDesign(k = 5, test.type = 2, n.fix = ss$nEvents, nFixSurv = ss$n, delta1 = log(ss$lambda2 / ss$lambda1))

# boundary plot
plot(x)
# effect size plot
plot(x, plottype = "hr")
# total sample size
x$nSurv
# number of events at analyses
x$n.I
# print the design
x
# overall design summary
cat(summary(x))
# tabular summary of bounds
gsBoundSummary(x, deltaname = "HR", Nname = "Events", logdelta = TRUE)

# approximate number of events required using Schoenfeld's method
# for 2 different hazard ratios
nEvents(hr = c(.5, .6), tbl = TRUE)
# vector output
nEvents(hr = c(.5, .6))

# approximate power using Schoenfeld's method
# given 2 sample sizes and hr=.6
nEvents(hr = .6, n = c(50, 100), tbl = TRUE)
# vector output
nEvents(hr = .6, n = c(50, 100))

# approximate hazard ratio corresponding to 100 events and z-statistic of 2
zn2hr(n = 100, z = 2)
# same when hr0 is 1.1
zn2hr(n = 100, z = 2, hr0 = 1.1)
# same when hr0 is .9 and hr1 is greater than hr0
zn2hr(n = 100, z = 2, hr0 = .9, hr1 = 1)

# approximate number of events corresponding to z-statistic of 2 and
# estimated hazard ratio of .5 (or 2)
hrz2n(hr = .5, z = 2)
hrz2n(hr = 2, z = 2)

# approximate z statistic corresponding to 75 events
# and estimated hazard ratio of .6 (or 1/.6)
# assuming 2-to-1 randomization of experimental to control
hrn2z(hr = .6, n = 75, ratio = 2)
hrn2z(hr = 1 / .6, n = 75, ratio = 2) 

**********
## R package ‘SampleSize4ClinicalTrials’ is another good R package to estimate sample size for four types of RCTs (equality, superiority, non-inferiority, equivalence). 
   # However, survival design is not covered. See: https://cran.r-project.org/web/packages/SampleSize4ClinicalTrials/index.html