## 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