frailtyDesign: Sample Size calculation and Power Analysis using...
In frailtypack: Shared, Joint (Generalized) Frailty Models; Surrogate Endpoints
frailtyDesign R Documentation
Sample Size calculation and Power Analysis using Gamma-Frailty Models
Description
A collection of functions to calculate statistical power and required
sample sizes for survival analysis using frailty models, specifically
the Shared Frailty Model (SFM), Nested Frailty Model (NFM), Joint
Frailty Model (JFM), and General Joint Frailty Model (GJFM).
For each frailty model type (denoted by *, where * corresponds to
SFM, NFM, JFM or GJFM), the package provides two distinct functions:
-
*.power: Computes the statistical power under given study settings.
-
*.ssize: Determines the required sample size needed to achieve
a specified target power under given study settings.
Usage
######################################################
## 1. SHARED FRAILTY MODEL (SFM)
######################################################
# Compute power for a given sample size in a SFM
# --------------------------------------------------
SFM.power(
Groups = 80, ni = 8, ni.type = "max", Acc.Dur = 0,FUP = 12,
FUP.type = "UpToEnd", median.H0 = 1, beta.H0 = 0, beta.HA = log(0.75),
shape.W = 1, theta = 0.25, ratio = 1, samples.mc = 1e4, seed = 42,
timescale = "gap", data.type = "grouped",
cens.par = 5, cens.type = "Expo", statistic = "Wald",
typeIerror = 0.05, test.type = "2-sided"
)
# Compute sample size for a given power in a SFM
# --------------------------------------------------
SFM.ssize(
power = 0.8, ni = 8, ni.type = "max", Acc.Dur = 0,FUP = 12,
FUP.type = "UpToEnd", median.H0 = 1, beta.H0 = 0, beta.HA = log(0.75),
shape.W = 1, theta = 0.25, ratio = 1, samples.mc = 1e4, seed = 42,
timescale = "gap", data.type = "grouped",
cens.par = 5, cens.type = "Expo", statistic = "Wald",
typeIerror = 0.05, test.type = "2-sided"
)
######################################################
## 2. NESTED FRAILTY MODEL (NFM)
######################################################
# Compute power for a given sample size in a NFM
# --------------------------------------------------
NFM.power(
Groups = 80, ni = 8, ni.type = "max", kij = 15, kij.type = "max",
Acc.Dur = 0, FUP = 12, FUP.type = "UpToEnd", median.H0 = 1,
beta.H0 = 0, beta.HA = log(0.75), shape.W = 1, theta = 0.25, eta = 0.5,
ratio = 1, samples.mc = 1e4, seed = 42,
timescale = "gap", data.type = "grouped", cens.par = 5, cens.type = "Expo",
statistic = "Wald", typeIerror = 0.05, test.type = "2-sided"
)
# Compute sample size for a given power in a NFM
# --------------------------------------------------
NFM.ssize(
power = 0.8, ni = 8, ni.type = "max", kij = 15, kij.type = "max",
Acc.Dur = 0, FUP = 12, FUP.type = "UpToEnd", median.H0 = 1,
beta.H0 = 0, beta.HA = log(0.75), shape.W = 1, theta = 0.25, eta = 0.5,
ratio = 1, samples.mc = 1e4, seed = 42,
timescale = "gap", data.type = "grouped", cens.par = 5, cens.type = "Expo",
statistic = "Wald", typeIerror = 0.05, test.type = "2-sided"
)
######################################################
## 3. JOINT FRAILTY MODEL (JFM)
######################################################
# Compute power for a given sample size in a JFM
# --------------------------------------------------
JFM.power(
Npts = 400, ni = 8, ni.type = "max", Acc.Dur = 0, FUP = 12,
FUP.type = "UpToEnd", medianR.H0 = 3, medianD.H0 = 10, betaTest.type = "joint",
betaR.H0 = 0, betaR.HA = log(0.75), betaD.H0 = 0, betaD.HA = log(0.85),
shapeR.W = 1, shapeD.W = 1, theta = 0.25, alpha = 1, ratio = 1,
samples.mc = 1e4, seed = 42, timescale = "gap",
statistic = "Wald", typeIerror = 0.05, test.type = "2-sided"
)
# Compute sample size for a given power in a JFM
# --------------------------------------------------
JFM.ssize(
power = 0.8, ni = 8, ni.type = "max", Acc.Dur = 0, FUP = 12,
FUP.type = "UpToEnd", medianR.H0 = 3, medianD.H0 = 10, betaTest.type = "joint",
betaR.H0 = 0, betaR.HA = log(0.75), betaD.H0 = 0, betaD.HA = log(0.85),
shapeR.W = 1, shapeD.W = 1, theta = 0.25, alpha = 1, ratio = 1,
samples.mc = 1e4, seed = 42, timescale = "gap",
statistic = "Wald", typeIerror = 0.05, test.type = "2-sided"
)
######################################################
## 4. GENERAL JOINT FRAILTY MODEL (GJFM)
######################################################
# Compute power for a given sample size in a GJFM
# --------------------------------------------------
GJFM.power(
Npts = 400, ni = 8, ni.type = "max", Acc.Dur = 0, FUP = 12,
FUP.type = "UpToEnd", medianR.H0 = 3, medianD.H0 = 10,
betaTest.type = "joint", betaR.H0 = 0, betaR.HA = log(0.75),
betaD.H0 = 0, betaD.HA = log(0.85), shapeR.W = 1, shapeD.W = 1,
theta = 0.25, eta = 0.5, ratio = 1, samples.mc = 1e4,
seed = 42, timescale = "gap",
statistic = "Wald", typeIerror = 0.05, test.type = "2-sided"
)
# Compute sample size for a given power in a GJFM
# --------------------------------------------------
GJFM.ssize(
power = 0.8, ni = 8, ni.type = "max", Acc.Dur = 0, FUP = 12,
FUP.type = "UpToEnd", medianR.H0 = 3, medianD.H0 = 10,
betaTest.type = "joint", betaR.H0 = 0, betaR.HA = log(0.75),
betaD.H0 = 0, betaD.HA = log(0.85), shapeR.W = 1, shapeD.W = 1,
theta = 0.25, eta = 0.5, ratio = 1, samples.mc = 1e4,
seed = 42, timescale = "gap",
statistic = "Wald", typeIerror = 0.05, test.type = "2-sided"
)
Arguments
Groups
Only in SFM and NFM: A numeric value, where interpretation
depends on the data.type parameter and on the model:
For SFM, it corresponds to either the number of groups (grouped data)
or the number of subjects (recurrent events data).
For NFM, it corresponds to either the number of groups (grouped and
recurrent event data) or the number of subjects (multi-type recurrent
events data).
Default is 80.
ni
A numeric value or an array (dim = 2), representing
expected values or distribution parameters. Interpretation depends on
the data.type parameter and on the model:
For SFM, it corresponds to either the expected number of subjects per group (grouped data)
or the expected number of recurrent events per subject (recurrent events data).
For NFM, it corresponds to the expected number of subgroups within each group (grouped data),
the expected number of recurrent events per group (recurrent event data), or the number of
distinct recurrent event type (multi-type recurrent event data).
For JFM/GJFM, it corresponds to the expected number of recurrent events per subject.
The default value is 8.
ni.type
Character value, specifying ni. Valid options:
-
"max": ni is a fixed number.
-
"pois": ni is a mean (parameter of a Poisson distribution).
-
"unif": ni is the lower and upper bound parameters of a uniform distribution.
Options "pois" and "unif" can only be selected for recurrent event data; see Note for details.
Default is "max".
Acc.Dur
Non-negative numeric value. Parameter for a uniform accrual
from time 0 to time Acc.Dur. Default is 0.
FUP
A positive numeric value of follow-up duration as defined in the
study protocol (i.e. administrative censoring). Default is 12.
FUP.type
Character value, indicating the type of follow-up. Valid options:
-
"Fixed": each subject is followed exactly for FUP time units after enrollment.
-
"UptoEnd": global study cutoff at time FUP; individual follow-up for at most FUP.
Default is "Fixed".
median.H0
Only in SFM and NFM: A positive numeric value, used for the scale parameter (Weibull) calculation.
If recurrent event data, it is the median gap time to event under the null (excluding censoring times).
If grouped data, it is the median time to an event under the null (excluding censoring times).
Default is 1.
beta.H0
Only in SFM and NFM: log-hazard ratio parameter under the null
hypothesis (H0). Default is 0.
beta.HA
Only in SFM and NFM: log-hazard ratio parameter under the
alternative hypothesis (HA). Default is log(0.75).
shape.W
Only in SFM and NFM: A positive numeric value, corresponding to
the shape parameter (Weibull) of the baseline hazard of the recurrent event. Default is 1.
theta
A positive numeric value, corresponding to the Gamma-frailty
variance for the main random effect. Default is 0.25.
ratio
A positive numeric value, corresponding to the allocation ratio
(experimental : control). Default is 1.
samples.mc
A positive numeric value, corresponding to the number of
Monte Carlo samples used to approximate the Fisher information matrix.
Default is 1e4.
seed
Integer number for random-number generation seed. Ensures reproducibility
of the Monte-Carlo simulations. Default is 42.
timescale
Character value indicating the timescale when recurrent event data type
is considered. Can be either 'gap' or 'calendar'. See note for more detail. Default is 'gap'.
data.type
Only in SFM and NFM: Character value indicating what kind of data we
want to consider for the current frailty model. Valid options differ depending on the model:
For SFM, can be either "grouped" (corresponding to subjects included in a group)
or "rec_event" (corresponding to subjects experiencing recurrent events).
For NFM, the hierarchical structure of the data can be either "grouped"
(where subjects are included into subgroup and subgroups into groups),
"rec_event1" (where the group level corresponds to a group (e.g., hospitals)
and subgroup level to a subject) or "rec_event2" (where the group
level corresponds to a subject and subgroup level to a type of recurrent event).
Default is "grouped".
cens.par
Only in SFM and NFM: A numeric value corresponding
to the parameter of the distribution for non-administrative censoring.
Default is 10000.
cens.type
Only in SFM and NFM: Character value, specifying
the distribution for non-administrative censoring. Valid options:
-
"Expo": in this case, cens.par is the median
from an exponential distribution.
-
"Unif": in this case, cens.par is the lower and upper
bound parameters of a uniform distribution.
Default is "Expo".
statistic
Type of test statistic used. Currently, only "Wald" is available.
typeIerror
A numeric value corresponding to the type I error level. Default is 0.05.
test.type
Character value indicating whether It is a one-tailed or two-tailed test.
Valid options are either "1-sided" or "2-sided". Default is "2-sided".
power
Numeric in (0,0.99]. The target power 1 - \beta. Default is 0.8.
kij
Only in NFM: A numeric value or an array (dim = 2), representing
expected values or distribution parameters. Interpretation depends on the data.type parameter:
For grouped data: It is the number of observations per subgroup.
For recurrent events data: It is the number of observation per subjects.
For multi-type recurrent events data: It is the number of recurrences for
each distinct type of event.
Default is 15.
kij.type
Character value, specifying kij. Valid options:
-
"max": kij is a fixed number.
-
"pois": kij is a mean (parameter of a Poisson distribution).
-
"unif": kij is the lower and upper bound parameters of a uniform distribution.
Options "pois" and "unif" can only be selected for recurrent event data; see Note for details.
Default is "max".
eta
Only in NFM and GJFM: positive numeric value, corresponding to an
additional Gamma-frailty variance parameter for second-level nesting (NFM)
or inter-recurrence dependence (GJFM). Default is 0.5.
Npts
Only in JFM and GJFM: positive numeric value, corresponding to the
total number of subjects. Default is 400.
medianR.H0
Only in JFM and GJFM: positive numeric value, corresponding
to the expected median time between two recurrent events under the null (H0),
for the scale parameter (Weibull) calculation. Default is 3.
medianD.H0
Only in JFM and GJFM: positive numeric value, corresponding
to the expected median time to the terminal event under the null (H0),
for the scale parameter (Weibull) calculation. Default is 10.
betaTest.type
Only in JFM and GJFM: character value indicating which
hypothesis is tested when computing power. Our implementation allows either
power calculation or sample-size estimation, testing recurrent events alone,
terminal event alone or both. Valid options:
"joint" (for testing both \beta_R and \beta_D), "betaRtest"
(for testing only \beta_R) or "betaDtest"
(for testing only \beta_D). Default is "joint".
betaR.H0
Only in JFM and GJFM: numeric value, corresponding to the
log-hazard ratios for recurrent events under the null hypothesis (H0). Default is 0.
betaR.HA
Only in JFM and GJFM: numeric value, corresponding to the
log-hazard ratios for recurrent events under the alternative hypothesis (HA). Default is log(0.75).
betaD.H0
Only in JFM and GJFM: numeric value, corresponding to the
log-hazard ratios for terminal events under the null hypothesis (H0). Default is 0.
betaD.HA
Only in JFM and GJFM: numeric value, corresponding to the
log-hazard ratios for terminal events under the alternative hypothesis (HA). Default is log(0.85).
shapeR.W
Only in JFM and GJFM: positive numeric value, corresponding to
the shape parameter (Weibull) of the recurrent-event baseline hazard function. Default is 1.
shapeD.W
Only in JFM and GJFM: positive numeric value, corresponding to
the shape parameter (Weibull) of the terminal-event baseline hazard function. Default is 1.
alpha
Only in JFM: numeric value, corresponding to the parameter \alpha
that modulates the association between recurrent and terminal events. Default is 1.
Details
See Dinart et al. (2024) for the original article. We present here the case
where we want to assess the treatment effect only. Our null hypothesis is that
there is no treatment effect (i.e. zero log-hazard ratio).
This approach relies on the squared Wald test to assess the presence of a
treatment effect using an estimator \hat{\Theta}. Specifically, under
our null hypothesis, the test statistic
X_w = Z^2 = (\hat{\Theta} - \Theta)^2 / \mathcal{I}^{-1}(\hat{\Theta})
follows a central \chi^2_1 distribution (i.e., non-centrality parameter
\mu = 0), whereas under the alternative hypothesis, it follows a
non-central \chi^2_1(\mu) distribution with \mu > 0.
The parameter \mu is estimated algorithmically, and the Fisher
information \mathcal{I}_1(\hat{\Theta}) is obtained by simulation,
leveraging the law of large numbers. Concretely, for an M-sample
generated by simulation, the matrix \mathcal{I}(\hat{\Theta}) is
approximated via the empirical mean of the products
\partial_{\Theta_k} l(\Theta(i)) \times \partial_{\Theta_l} l(\Theta(i))^\top,
i \in [\![1, \dots, M]\!]. The algorithmic estimation for \mu
follows the three-step procedure described by Dinart et al. (2024):
Compute the \alpha-quantile (denoted q_{1,\alpha}) of a central
chi-square distribution with 1 degree of freedom (\chi^2_1),
given a specified type I error rate \alpha.
Determine a non-centrality parameter \vartheta such that
1 - P\bigl(\chi^2_1(\vartheta) < q_{1,\alpha}\bigr) > 1 - \beta,
where 1 - \beta represents the desired statistical power and
1 - P\bigl(\chi^2_1(\vartheta) < q_{1,\alpha}\bigr) the computed power.
Optimize \mu to find the smallest value satisfying that condition for
all
x \in [0,\vartheta], i.e.,
\mu
= \min_{x \in [0,\vartheta]} \Bigl\{ 1 - P\bigl(\chi^2_1(x) < q_{1,\alpha}\bigr)
- (1 - \beta)\Bigr\}.
Once \mu is estimated, the sample size n is derived from
n \,\geq\, \mu
\times \Bigl(\Theta_A^2 \times \mathcal{I}_1(\hat{\Theta})\Bigr)^{-1},
where \Theta_A denotes the parameter value under H_A.
If we are interested in the evaluation of the power, we estimate the
non-centrality parameter under H_A for a given sample size N, then compute
the power as P\bigl(\chi^2_1(\vartheta) > q_{1,\alpha}\bigr\vert H_A).
For both the joint frailty model and general joint frailty model, by following
the same methodology as in the univariate case, we can derive an expression
for the sample size from the generalized Wald statistic. Let H_0:
(\beta_R = 0) \text{ and } (\beta_D = 0) vs. H_A: (\beta_R = \beta_R^A)
\text{ or } (\beta_D = \beta_D^A), be our null and alternative hypotheses respectively.
This multivariate test then follows a \chi^2_Q distribution, where
Q is the rank of the matrix C, corresponding to the number of
constraints applied on the parameters under the null hypothesis. The test
statistic is:
X_W \;=\; n\,\bigl(C\,\Omega\bigr)^\top
\Bigl(C\,\mathcal{I}_1^{-1}(\Omega)\,C^\top\Bigr)^{-1}\bigl(C\,\Omega\bigr)
\;\sim\;\chi^2_q(\mu),
where \Omega^\top is the vector parameter from the corresponding model.
From this, we derive a sample size formula:
n \,\geq\, \mu \,\Bigl(\bigl(C\,\Omega\bigr)^\top
\bigl(C\,\mathcal{I}_1^{-1}(\Omega)\,C^\top\bigr)^{-1}\bigl(C\,\Omega\bigr)\Bigr)^{-1}.
For instance, for the JFM, we have
\Omega^\top = (\beta_R, \beta_D, r_0(.), h_0(.), \theta, \alpha)^\top).
If we want to test a treatment effect on both the recurrent and the terminal event
(i.e. H_0 : \beta_R=0 \text{ and } \beta_D=0), hence:
C\times \Omega = \begin{pmatrix}1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \end{pmatrix}
\times (\beta_R, \beta_D, r_0(.), h_0(.), \theta, \alpha)^\top
Value
For SFM and NFM, the *.power() function returns a list with:
-
estimated.power: the estimated power.
-
number.events: an array (dim = 2) containing the number of recurrent events
under the null and alternative hypotheses, respectively.
The *.ssize() function returns a list with:
-
Groups: the total number of subjects (recurrent event data)
or number of groups (clustered data).
-
number.events: an array (dim = 2) containing the number of events
under the null and alternative hypotheses, respectively.
For JFM and GJFM, the *.power() function returns a list with:
-
estimated.power: the estimated power.
-
events.rec: an array (dim = 2) containing the number of recurrent events
under the null and alternative hypotheses, respectively.
-
events.D: an array (dim = 2) containing the number of terminal events
under the null and alternative hypotheses, respectively.
The *.ssize() function returns a list with:
-
Npts: the computed total number of subjects.
-
events.rec: an array (dim = 2) containing the number of recurrent events
under the null and alternative hypotheses, respectively.
-
events.D: an array (dim = 2) containing the number of terminal events
under the null and alternative hypotheses, respectively.
All returned lists additionally include several input parameters: target.power or
Groups/Npts (depending on the called function), ni, FUP,
FUP.type, Acc.Dur, ratio, data.type, the test type testType,
alpha, theta, eta (for NFM and GJFM) and samplesMC. For SFM and NFM,
we have the corresponding hazard ratios (HR.H0, HR.HA) from the given betas,
median.H0. For JGM and GJFM, we have the corresponding hazard ratios (HR.R0,
HR.RA, HR.D0, HR.DA) from the given betas, medianR.H0, medianD.H0,
and the testing structure tested.structure. Along with that, model
("SFM", "NFM", "JFM" or "GJFM"), method ("power" or "ssize") and timescale
("gap" or "calendar") are also included.
All these parameters are utilized by the S3 methods print.frailtyDesign and
summary.frailtyDesign for further detailed output.
Note
Internally, these functions rely on extensive numerical integration using Gaussian-Laguerre
quadrature to approximate the Fisher information. As such, computations may become resource-intensive.
You may need to adjust the parameter samples.mc or other integration parameters to enhance
computational performance or precision.
In this implementation, users must provide both the median time and the shape parameter explicitly;
the scale parameter is then computed automatically. Under the Weibull distribution, the median time
t_{1/2} relates to the scale and shape parameters via:
t_{1/2} = \text{scale} \times \log(2)^{1/\text{shape}}. Consequently, the scale parameter is
calculated as: \text{scale} = \frac{t_{1/2}}{\log(2)^{1/\text{shape}}}.
For both SFM and NFM, when analyzing grouped data, the arguments ni.type
and kij.type are restricted to the value "max" to define an
exact sample size. For instance, specifying ni.type = "pois" would
represent a mean number of subgroups per group, thereby precluding
the determination of a precise sample size.
In NFM, the parameter "rec_event2" might initially appear difficult to interpret.
As clarified by Derek et al. (2024), multitype recurrent events include situations such as
transient ischemic attacks classified by anatomical location in cardiovascular studies, or
migraines differentiated according to severity in neurological research.
In survival analysis involving recurrent events, the interpretation of
regression coefficients (\beta) is contingent upon the chosen timescale.
This distinction is crucial, as the timescale directly influences the risk
assessment and the corresponding interpretation of model parameters.
When employing a gap timescale, the timescale resets after each event,
measuring the duration until the next occurrence. Consequently, the regression
coefficients represent the modification of the inter-event risk, reflecting
how the treatment influence the hazard of experiencing a subsequent event
after the previous one. This approach focuses on the conditional risk between
events.
In contrast, utilizing a calendar timescale measures the time from a fixed
origin, such as study entry, without resetting after each event. Here, the
regression coefficients pertain to the modification of the risk since the
initiation of the study, indicating how the treatment affect the hazard of
experiencing events over the entire follow-up period. This approach
focuses on the cumulative risk from the study entry.
Author(s)
Original code by Dinart Derek. Implementation by Adrien Orué.
References
Derek Dinart, Carine Bellera & Virginie Rondeau (09 Feb 2024).
Sample size estimation for recurrent event data using multifrailty
and multilevel survival models,
Journal of Biopharmaceutical Statistics,
DOI: 10.1080/10543406.2024.2310306.
See Also
frailtyPenal, print.frailtyDesign,
summary.frailtyDesign
Examples
# Example 1 (SFM): a total of 400 patients (1:1 randomization scheme),
# with a fixed number of 3 recurrent events per patient. Gamma-frailty
# variance of 0.5. Expected hazard ratio of 0.7, time-to-death are uniformly
# distributed, with a mean time to death of (3+10)/2=6.5 years. Each subject is
# followed-up for a maximum of 6 years, with a median time-to-event of 1.5 years.
# Patients are recruited over a 0.5-year period.
SFM.power(
Groups = 400, ni = 3, ni.type = "max",
FUP = 6, Acc.Dur = 0.5, median.H0 = 1.5,
beta.HA = log(0.7), theta = 0.5,
cens.par = c(3, 10), cens.type = "Unif",
data.type = "rec_event"
) # power ~ 90%
# Example 2 (NFM): same parameters as above, but we now assume that we have
# 40 hospitals, 10 subjects per hospital (10 × 40 = 400 subjects in total)
# and 3 recurrent events per subject.
NFM.power(
Groups = 40, ni = 10, ni.type = "max", kij = 3, kij.type = "max",
FUP = 6, Acc.Dur = 0.5, median.H0 = 1.5,
beta.HA = log(0.7), theta = 0.5,
cens.par = c(3, 10), cens.type = "Unif",
data.type = "rec_event1"
) # power ~ 83%
# Example 3 (NFM): we aim to compute the required sample size to achieve
# 80% power for detecting a hazard ratio of 0.75 in a neurological study,
# where migraine episodes experienced by subjects are classified into three
# severity subtypes (mild, moderate, severe). For each subject, we anticipate
# a mean number of 2 migraine episodes per severity subtype, with a median
# time-to-event of 6 months. The study duration includes a 1-year accrual period
# followed by a 5-year total follow-up. All subjects will be followed until
# the end of the study.
NFM.ssize(
power = 0.80, ni = 3, ni.type = "max", kij = 2, kij.type = "pois",
FUP = 5, Acc.Dur = 1, FUP.type = "uptoend", median.H0 = 0.5,
beta.HA = log(0.75), data.type = "rec_event2"
) # sample size ~ 363 patients
# Example 4 (JFM): power estimation, testing a treatment effect on recurrent
# events only. We assume a uniformly distributed number of recurrent events,
# ranging from 1 to 6 recurrent events per subject. The allocation ratio
# experimental:control is 2:1, and the follow-up is 10 weeks.
# The expected hazard ratio is 0.70 for recurrent event and 0.90 for the
# terminal event. We have chosen 0.5 as the variance of the frailties.
JFM.power(
Npts = 400, ni = c(1, 6), ni.type = "unif",
FUP = 10, FUP.type = "fixed", ratio = 2,
betaTest.type = "betaRtest", betaR.HA = log(.70), betaD.HA = log(.90),
theta = .5
) # power ~ 76%
# Example 5 (JFM): sample size calculation, to assess the treatment effect on
# both recurrent and terminal events. We want to achieve an 80% power.
# We anticipate a maximum of 5 recurrent events, over a 6-year period and a
# 0.5-year accrual period. We assume that the gamma-frailty variance is 0.5.
# For the control group, we expect a 2-year and a 5-year median time-to-event
# for recurrent events and terminal events, respectively. We consider a 30%
# and 20% risk reduction for recurrent events and terminal event, respectively.
JFM.ssize(
power = 0.80, ni = 9,
FUP = 6, Acc.Dur = 1.5, medianR.H0 = 2, medianD.H0 = 5,
betaTest.type = "joint", betaR.HA = log(.70), betaD.HA = log(.80), theta = .5
) # sample size ~ 445 patients / ~ approx
# Example 6: Sample size calculation for GJFM
# Same as above, but with two random effects (with two Gamma-frailty variances
# theta and eta). To ensure sample size estimation stability, we use 10000
# Monte-Carlo samples.
GJFM.ssize(
power = 0.80, ni = 5,
FUP = 6, Acc.Dur = 0.5, medianR.H0 = 2, medianD.H0 = 5,
betaR.HA = log(0.70), betaD.HA = log(0.80), theta = 0.5, eta = 0.75,
samples.mc = 1e5
) # sample size ~ 705 patients / ~ approx 4 min.
# Example 7:
# --------------------------------------------------
# Post-hoc power analysis for a Joint Frailty Model
# --------------------------------------------------
# See original article by Gonzalez et al. (2005)
data(readmission)
modJFM <- frailtyPenal(
Surv(time, event) ~ cluster(id) + as.factor(chemo) + terminal(death),
formula.terminalEvent = ~ as.factor(chemo), data = readmission,
hazard = "Weibull"
)
# Test both recurrent and death events
# # - Let us assume an underlying Poisson distribution for ni.type. The
# # empirical mean of the number of recurrent events per patients is: ni = 1.136476.
# # - For the null hypothesis, let us consider betaR.H0 = betaD.H0 = 0. For the
# # alternative hypothesis, we use the estimated parameters for betaD.HA and
# # betaR.HA.
# # - "Patients were actively followed up until June 2002" -> the follow-up
# # type is "UpToEnd".
# # - "The study took place in the Hospital de Bellvitge, [...] between
# # January 1996 and December 1998" -> the accrual time is approximately 3 years
# # - We can assume that the study duration is approximately 6 years
ni <- 1.136476
ni.type <- "Pois"
Acc.Dur <- 3 * 365.25 # time unit = days
FUP <- 6 * 365.25 # same as above
betaR.HA <- as.numeric(modJFM$coef[1]) # else "Named numeric"
betaD.HA <- as.numeric(modJFM$coef[2]) # same as above
med <- modJFM$scale.weib * log(2)^(1 / modJFM$shape.weib)
medianR.H0 <- med[1]
medianD.H0 <- med[2]
shapeR.W <- modJFM$shape[1]
shapeD.W <- modJFM$shape[2]
theta <- modJFM$theta
alpha <- modJFM$alpha
Npts <- length(unique(readmission[, "id"])) # 403 patients
nTreated <- length(unique(readmission[readmission$chemo == "Treated", "id"])) #217 treated patients
ratio <- nTreated / (Npts - nTreated)
JFM.power(
Npts = Npts, ni = ni, ni.type = ni.type,
Acc.Dur = Acc.Dur, FUP = FUP, medianR.H0 = medianR.H0, medianD.H0 = medianD.H0,
betaTest.type = "joint", betaR.HA = betaR.HA, betaD.HA = betaD.HA,
shapeR.W = shapeR.W, shapeD.W = shapeD.W, theta = theta, alpha = alpha,
ratio = ratio
) # power ~ 92%
# --------------------------------------------------
# Required sample size under the same setting
# --------------------------------------------------
# Here, let us consider that readmission is a “pilot study” with 403 patients,
# from which we estimate parameters. Under this scenario, let us compute the
# needed sample size, but to achieve an 80% power.
JFM.ssize(
power = 0.80, ni = ni, ni.type = ni.type,
Acc.Dur = Acc.Dur, FUP = FUP, medianR.H0 = medianR.H0, medianD.H0 = medianD.H0,
betaTest.type = "joint", betaR.HA = betaR.HA, betaD.HA = betaD.HA,
shapeR.W = shapeR.W, shapeD.W = shapeD.W, theta = theta, alpha = alpha,
ratio = ratio
) # 289 patients needed under the same settings vs. 403
frailtypack documentation built on April 4, 2025, 3:08 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| frailtyDesign | R Documentation |
Sample Size calculation and Power Analysis using Gamma-Frailty Models
Description
A collection of functions to calculate statistical power and required sample sizes for survival analysis using frailty models, specifically the Shared Frailty Model (SFM), Nested Frailty Model (NFM), Joint Frailty Model (JFM), and General Joint Frailty Model (GJFM).
For each frailty model type (denoted by *, where * corresponds to
SFM, NFM, JFM or GJFM), the package provides two distinct functions:
-
*.power: Computes the statistical power under given study settings. -
*.ssize: Determines the required sample size needed to achieve a specified target power under given study settings.
Usage
######################################################
## 1. SHARED FRAILTY MODEL (SFM)
######################################################
# Compute power for a given sample size in a SFM
# --------------------------------------------------
SFM.power(
Groups = 80, ni = 8, ni.type = "max", Acc.Dur = 0,FUP = 12,
FUP.type = "UpToEnd", median.H0 = 1, beta.H0 = 0, beta.HA = log(0.75),
shape.W = 1, theta = 0.25, ratio = 1, samples.mc = 1e4, seed = 42,
timescale = "gap", data.type = "grouped",
cens.par = 5, cens.type = "Expo", statistic = "Wald",
typeIerror = 0.05, test.type = "2-sided"
)
# Compute sample size for a given power in a SFM
# --------------------------------------------------
SFM.ssize(
power = 0.8, ni = 8, ni.type = "max", Acc.Dur = 0,FUP = 12,
FUP.type = "UpToEnd", median.H0 = 1, beta.H0 = 0, beta.HA = log(0.75),
shape.W = 1, theta = 0.25, ratio = 1, samples.mc = 1e4, seed = 42,
timescale = "gap", data.type = "grouped",
cens.par = 5, cens.type = "Expo", statistic = "Wald",
typeIerror = 0.05, test.type = "2-sided"
)
######################################################
## 2. NESTED FRAILTY MODEL (NFM)
######################################################
# Compute power for a given sample size in a NFM
# --------------------------------------------------
NFM.power(
Groups = 80, ni = 8, ni.type = "max", kij = 15, kij.type = "max",
Acc.Dur = 0, FUP = 12, FUP.type = "UpToEnd", median.H0 = 1,
beta.H0 = 0, beta.HA = log(0.75), shape.W = 1, theta = 0.25, eta = 0.5,
ratio = 1, samples.mc = 1e4, seed = 42,
timescale = "gap", data.type = "grouped", cens.par = 5, cens.type = "Expo",
statistic = "Wald", typeIerror = 0.05, test.type = "2-sided"
)
# Compute sample size for a given power in a NFM
# --------------------------------------------------
NFM.ssize(
power = 0.8, ni = 8, ni.type = "max", kij = 15, kij.type = "max",
Acc.Dur = 0, FUP = 12, FUP.type = "UpToEnd", median.H0 = 1,
beta.H0 = 0, beta.HA = log(0.75), shape.W = 1, theta = 0.25, eta = 0.5,
ratio = 1, samples.mc = 1e4, seed = 42,
timescale = "gap", data.type = "grouped", cens.par = 5, cens.type = "Expo",
statistic = "Wald", typeIerror = 0.05, test.type = "2-sided"
)
######################################################
## 3. JOINT FRAILTY MODEL (JFM)
######################################################
# Compute power for a given sample size in a JFM
# --------------------------------------------------
JFM.power(
Npts = 400, ni = 8, ni.type = "max", Acc.Dur = 0, FUP = 12,
FUP.type = "UpToEnd", medianR.H0 = 3, medianD.H0 = 10, betaTest.type = "joint",
betaR.H0 = 0, betaR.HA = log(0.75), betaD.H0 = 0, betaD.HA = log(0.85),
shapeR.W = 1, shapeD.W = 1, theta = 0.25, alpha = 1, ratio = 1,
samples.mc = 1e4, seed = 42, timescale = "gap",
statistic = "Wald", typeIerror = 0.05, test.type = "2-sided"
)
# Compute sample size for a given power in a JFM
# --------------------------------------------------
JFM.ssize(
power = 0.8, ni = 8, ni.type = "max", Acc.Dur = 0, FUP = 12,
FUP.type = "UpToEnd", medianR.H0 = 3, medianD.H0 = 10, betaTest.type = "joint",
betaR.H0 = 0, betaR.HA = log(0.75), betaD.H0 = 0, betaD.HA = log(0.85),
shapeR.W = 1, shapeD.W = 1, theta = 0.25, alpha = 1, ratio = 1,
samples.mc = 1e4, seed = 42, timescale = "gap",
statistic = "Wald", typeIerror = 0.05, test.type = "2-sided"
)
######################################################
## 4. GENERAL JOINT FRAILTY MODEL (GJFM)
######################################################
# Compute power for a given sample size in a GJFM
# --------------------------------------------------
GJFM.power(
Npts = 400, ni = 8, ni.type = "max", Acc.Dur = 0, FUP = 12,
FUP.type = "UpToEnd", medianR.H0 = 3, medianD.H0 = 10,
betaTest.type = "joint", betaR.H0 = 0, betaR.HA = log(0.75),
betaD.H0 = 0, betaD.HA = log(0.85), shapeR.W = 1, shapeD.W = 1,
theta = 0.25, eta = 0.5, ratio = 1, samples.mc = 1e4,
seed = 42, timescale = "gap",
statistic = "Wald", typeIerror = 0.05, test.type = "2-sided"
)
# Compute sample size for a given power in a GJFM
# --------------------------------------------------
GJFM.ssize(
power = 0.8, ni = 8, ni.type = "max", Acc.Dur = 0, FUP = 12,
FUP.type = "UpToEnd", medianR.H0 = 3, medianD.H0 = 10,
betaTest.type = "joint", betaR.H0 = 0, betaR.HA = log(0.75),
betaD.H0 = 0, betaD.HA = log(0.85), shapeR.W = 1, shapeD.W = 1,
theta = 0.25, eta = 0.5, ratio = 1, samples.mc = 1e4,
seed = 42, timescale = "gap",
statistic = "Wald", typeIerror = 0.05, test.type = "2-sided"
)
Arguments
Groups |
Only in SFM and NFM: A numeric value, where interpretation
depends on the
Default is 80. |
ni |
A numeric value or an array (dim = 2), representing
expected values or distribution parameters. Interpretation depends on
the
The default value is 8. |
ni.type |
Character value, specifying
Options |
Acc.Dur |
Non-negative numeric value. Parameter for a uniform accrual
from time 0 to time |
FUP |
A positive numeric value of follow-up duration as defined in the study protocol (i.e. administrative censoring). Default is 12. |
FUP.type |
Character value, indicating the type of follow-up. Valid options:
Default is |
median.H0 |
Only in SFM and NFM: A positive numeric value, used for the scale parameter (Weibull) calculation. If recurrent event data, it is the median gap time to event under the null (excluding censoring times). If grouped data, it is the median time to an event under the null (excluding censoring times). Default is 1. |
beta.H0 |
Only in SFM and NFM: log-hazard ratio parameter under the null hypothesis (H0). Default is 0. |
beta.HA |
Only in SFM and NFM: log-hazard ratio parameter under the alternative hypothesis (HA). Default is log(0.75). |
shape.W |
Only in SFM and NFM: A positive numeric value, corresponding to the shape parameter (Weibull) of the baseline hazard of the recurrent event. Default is 1. |
theta |
A positive numeric value, corresponding to the Gamma-frailty variance for the main random effect. Default is 0.25. |
ratio |
A positive numeric value, corresponding to the allocation ratio (experimental : control). Default is 1. |
samples.mc |
A positive numeric value, corresponding to the number of Monte Carlo samples used to approximate the Fisher information matrix. Default is 1e4. |
seed |
Integer number for random-number generation seed. Ensures reproducibility of the Monte-Carlo simulations. Default is 42. |
timescale |
Character value indicating the timescale when recurrent event data type is considered. Can be either 'gap' or 'calendar'. See note for more detail. Default is 'gap'. |
data.type |
Only in SFM and NFM: Character value indicating what kind of data we want to consider for the current frailty model. Valid options differ depending on the model:
Default is "grouped". |
cens.par |
Only in SFM and NFM: A numeric value corresponding to the parameter of the distribution for non-administrative censoring. Default is 10000. |
cens.type |
Only in SFM and NFM: Character value, specifying the distribution for non-administrative censoring. Valid options:
Default is "Expo". |
statistic |
Type of test statistic used. Currently, only |
typeIerror |
A numeric value corresponding to the type I error level. Default is 0.05. |
test.type |
Character value indicating whether It is a one-tailed or two-tailed test.
Valid options are either |
power |
Numeric in (0,0.99]. The target power |
kij |
Only in NFM: A numeric value or an array (dim = 2), representing
expected values or distribution parameters. Interpretation depends on the
Default is 15. |
kij.type |
Character value, specifying
Options |
eta |
Only in NFM and GJFM: positive numeric value, corresponding to an additional Gamma-frailty variance parameter for second-level nesting (NFM) or inter-recurrence dependence (GJFM). Default is 0.5. |
Npts |
Only in JFM and GJFM: positive numeric value, corresponding to the total number of subjects. Default is 400. |
medianR.H0 |
Only in JFM and GJFM: positive numeric value, corresponding to the expected median time between two recurrent events under the null (H0), for the scale parameter (Weibull) calculation. Default is 3. |
medianD.H0 |
Only in JFM and GJFM: positive numeric value, corresponding to the expected median time to the terminal event under the null (H0), for the scale parameter (Weibull) calculation. Default is 10. |
betaTest.type |
Only in JFM and GJFM: character value indicating which
hypothesis is tested when computing power. Our implementation allows either
power calculation or sample-size estimation, testing recurrent events alone,
terminal event alone or both. Valid options:
|
betaR.H0 |
Only in JFM and GJFM: numeric value, corresponding to the log-hazard ratios for recurrent events under the null hypothesis (H0). Default is 0. |
betaR.HA |
Only in JFM and GJFM: numeric value, corresponding to the log-hazard ratios for recurrent events under the alternative hypothesis (HA). Default is log(0.75). |
betaD.H0 |
Only in JFM and GJFM: numeric value, corresponding to the log-hazard ratios for terminal events under the null hypothesis (H0). Default is 0. |
betaD.HA |
Only in JFM and GJFM: numeric value, corresponding to the log-hazard ratios for terminal events under the alternative hypothesis (HA). Default is log(0.85). |
shapeR.W |
Only in JFM and GJFM: positive numeric value, corresponding to the shape parameter (Weibull) of the recurrent-event baseline hazard function. Default is 1. |
shapeD.W |
Only in JFM and GJFM: positive numeric value, corresponding to the shape parameter (Weibull) of the terminal-event baseline hazard function. Default is 1. |
alpha |
Only in JFM: numeric value, corresponding to the parameter |
Details
See Dinart et al. (2024) for the original article. We present here the case where we want to assess the treatment effect only. Our null hypothesis is that there is no treatment effect (i.e. zero log-hazard ratio).
This approach relies on the squared Wald test to assess the presence of a
treatment effect using an estimator \hat{\Theta}. Specifically, under
our null hypothesis, the test statistic
X_w = Z^2 = (\hat{\Theta} - \Theta)^2 / \mathcal{I}^{-1}(\hat{\Theta})
follows a central \chi^2_1 distribution (i.e., non-centrality parameter
\mu = 0), whereas under the alternative hypothesis, it follows a
non-central \chi^2_1(\mu) distribution with \mu > 0.
The parameter \mu is estimated algorithmically, and the Fisher
information \mathcal{I}_1(\hat{\Theta}) is obtained by simulation,
leveraging the law of large numbers. Concretely, for an M-sample
generated by simulation, the matrix \mathcal{I}(\hat{\Theta}) is
approximated via the empirical mean of the products
\partial_{\Theta_k} l(\Theta(i)) \times \partial_{\Theta_l} l(\Theta(i))^\top,
i \in [\![1, \dots, M]\!]. The algorithmic estimation for \mu
follows the three-step procedure described by Dinart et al. (2024):
Compute the
\alpha-quantile (denotedq_{1,\alpha}) of a central chi-square distribution with 1 degree of freedom (\chi^2_1), given a specified type I error rate\alpha.Determine a non-centrality parameter
\varthetasuch that1 - P\bigl(\chi^2_1(\vartheta) < q_{1,\alpha}\bigr) > 1 - \beta, where1 - \betarepresents the desired statistical power and1 - P\bigl(\chi^2_1(\vartheta) < q_{1,\alpha}\bigr)the computed power.Optimize
\muto find the smallest value satisfying that condition for allx \in [0,\vartheta], i.e.,\mu = \min_{x \in [0,\vartheta]} \Bigl\{ 1 - P\bigl(\chi^2_1(x) < q_{1,\alpha}\bigr) - (1 - \beta)\Bigr\}.
Once \mu is estimated, the sample size n is derived from
n \,\geq\, \mu
\times \Bigl(\Theta_A^2 \times \mathcal{I}_1(\hat{\Theta})\Bigr)^{-1},
where \Theta_A denotes the parameter value under H_A.
If we are interested in the evaluation of the power, we estimate the
non-centrality parameter under H_A for a given sample size N, then compute
the power as P\bigl(\chi^2_1(\vartheta) > q_{1,\alpha}\bigr\vert H_A).
For both the joint frailty model and general joint frailty model, by following
the same methodology as in the univariate case, we can derive an expression
for the sample size from the generalized Wald statistic. Let H_0:
(\beta_R = 0) \text{ and } (\beta_D = 0) vs. H_A: (\beta_R = \beta_R^A)
\text{ or } (\beta_D = \beta_D^A), be our null and alternative hypotheses respectively.
This multivariate test then follows a \chi^2_Q distribution, where
Q is the rank of the matrix C, corresponding to the number of
constraints applied on the parameters under the null hypothesis. The test
statistic is:
X_W \;=\; n\,\bigl(C\,\Omega\bigr)^\top
\Bigl(C\,\mathcal{I}_1^{-1}(\Omega)\,C^\top\Bigr)^{-1}\bigl(C\,\Omega\bigr)
\;\sim\;\chi^2_q(\mu),
where \Omega^\top is the vector parameter from the corresponding model.
From this, we derive a sample size formula:
n \,\geq\, \mu \,\Bigl(\bigl(C\,\Omega\bigr)^\top
\bigl(C\,\mathcal{I}_1^{-1}(\Omega)\,C^\top\bigr)^{-1}\bigl(C\,\Omega\bigr)\Bigr)^{-1}.
For instance, for the JFM, we have
\Omega^\top = (\beta_R, \beta_D, r_0(.), h_0(.), \theta, \alpha)^\top).
If we want to test a treatment effect on both the recurrent and the terminal event
(i.e. H_0 : \beta_R=0 \text{ and } \beta_D=0), hence:
C\times \Omega = \begin{pmatrix}1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \end{pmatrix}
\times (\beta_R, \beta_D, r_0(.), h_0(.), \theta, \alpha)^\top
Value
For SFM and NFM, the *.power() function returns a list with:
-
estimated.power: the estimated power. -
number.events: an array (dim = 2) containing the number of recurrent events under the null and alternative hypotheses, respectively.
The *.ssize() function returns a list with:
-
Groups: the total number of subjects (recurrent event data) or number of groups (clustered data). -
number.events: an array (dim = 2) containing the number of events under the null and alternative hypotheses, respectively.
For JFM and GJFM, the *.power() function returns a list with:
-
estimated.power: the estimated power. -
events.rec: an array (dim = 2) containing the number of recurrent events under the null and alternative hypotheses, respectively. -
events.D: an array (dim = 2) containing the number of terminal events under the null and alternative hypotheses, respectively.
The *.ssize() function returns a list with:
-
Npts: the computed total number of subjects. -
events.rec: an array (dim = 2) containing the number of recurrent events under the null and alternative hypotheses, respectively. -
events.D: an array (dim = 2) containing the number of terminal events under the null and alternative hypotheses, respectively.
All returned lists additionally include several input parameters: target.power or
Groups/Npts (depending on the called function), ni, FUP,
FUP.type, Acc.Dur, ratio, data.type, the test type testType,
alpha, theta, eta (for NFM and GJFM) and samplesMC. For SFM and NFM,
we have the corresponding hazard ratios (HR.H0, HR.HA) from the given betas,
median.H0. For JGM and GJFM, we have the corresponding hazard ratios (HR.R0,
HR.RA, HR.D0, HR.DA) from the given betas, medianR.H0, medianD.H0,
and the testing structure tested.structure. Along with that, model
("SFM", "NFM", "JFM" or "GJFM"), method ("power" or "ssize") and timescale
("gap" or "calendar") are also included.
All these parameters are utilized by the S3 methods print.frailtyDesign and
summary.frailtyDesign for further detailed output.
Note
Internally, these functions rely on extensive numerical integration using Gaussian-Laguerre
quadrature to approximate the Fisher information. As such, computations may become resource-intensive.
You may need to adjust the parameter samples.mc or other integration parameters to enhance
computational performance or precision.
In this implementation, users must provide both the median time and the shape parameter explicitly;
the scale parameter is then computed automatically. Under the Weibull distribution, the median time
t_{1/2} relates to the scale and shape parameters via:
t_{1/2} = \text{scale} \times \log(2)^{1/\text{shape}}. Consequently, the scale parameter is
calculated as: \text{scale} = \frac{t_{1/2}}{\log(2)^{1/\text{shape}}}.
For both SFM and NFM, when analyzing grouped data, the arguments ni.type
and kij.type are restricted to the value "max" to define an
exact sample size. For instance, specifying ni.type = "pois" would
represent a mean number of subgroups per group, thereby precluding
the determination of a precise sample size.
In NFM, the parameter "rec_event2" might initially appear difficult to interpret.
As clarified by Derek et al. (2024), multitype recurrent events include situations such as
transient ischemic attacks classified by anatomical location in cardiovascular studies, or
migraines differentiated according to severity in neurological research.
In survival analysis involving recurrent events, the interpretation of
regression coefficients (\beta) is contingent upon the chosen timescale.
This distinction is crucial, as the timescale directly influences the risk
assessment and the corresponding interpretation of model parameters.
When employing a gap timescale, the timescale resets after each event, measuring the duration until the next occurrence. Consequently, the regression coefficients represent the modification of the inter-event risk, reflecting how the treatment influence the hazard of experiencing a subsequent event after the previous one. This approach focuses on the conditional risk between events.
In contrast, utilizing a calendar timescale measures the time from a fixed origin, such as study entry, without resetting after each event. Here, the regression coefficients pertain to the modification of the risk since the initiation of the study, indicating how the treatment affect the hazard of experiencing events over the entire follow-up period. This approach focuses on the cumulative risk from the study entry.
Author(s)
Original code by Dinart Derek. Implementation by Adrien Orué.
References
Derek Dinart, Carine Bellera & Virginie Rondeau (09 Feb 2024). Sample size estimation for recurrent event data using multifrailty and multilevel survival models, Journal of Biopharmaceutical Statistics, DOI: 10.1080/10543406.2024.2310306.
See Also
frailtyPenal, print.frailtyDesign,
summary.frailtyDesign
Examples
# Example 1 (SFM): a total of 400 patients (1:1 randomization scheme),
# with a fixed number of 3 recurrent events per patient. Gamma-frailty
# variance of 0.5. Expected hazard ratio of 0.7, time-to-death are uniformly
# distributed, with a mean time to death of (3+10)/2=6.5 years. Each subject is
# followed-up for a maximum of 6 years, with a median time-to-event of 1.5 years.
# Patients are recruited over a 0.5-year period.
SFM.power(
Groups = 400, ni = 3, ni.type = "max",
FUP = 6, Acc.Dur = 0.5, median.H0 = 1.5,
beta.HA = log(0.7), theta = 0.5,
cens.par = c(3, 10), cens.type = "Unif",
data.type = "rec_event"
) # power ~ 90%
# Example 2 (NFM): same parameters as above, but we now assume that we have
# 40 hospitals, 10 subjects per hospital (10 × 40 = 400 subjects in total)
# and 3 recurrent events per subject.
NFM.power(
Groups = 40, ni = 10, ni.type = "max", kij = 3, kij.type = "max",
FUP = 6, Acc.Dur = 0.5, median.H0 = 1.5,
beta.HA = log(0.7), theta = 0.5,
cens.par = c(3, 10), cens.type = "Unif",
data.type = "rec_event1"
) # power ~ 83%
# Example 3 (NFM): we aim to compute the required sample size to achieve
# 80% power for detecting a hazard ratio of 0.75 in a neurological study,
# where migraine episodes experienced by subjects are classified into three
# severity subtypes (mild, moderate, severe). For each subject, we anticipate
# a mean number of 2 migraine episodes per severity subtype, with a median
# time-to-event of 6 months. The study duration includes a 1-year accrual period
# followed by a 5-year total follow-up. All subjects will be followed until
# the end of the study.
NFM.ssize(
power = 0.80, ni = 3, ni.type = "max", kij = 2, kij.type = "pois",
FUP = 5, Acc.Dur = 1, FUP.type = "uptoend", median.H0 = 0.5,
beta.HA = log(0.75), data.type = "rec_event2"
) # sample size ~ 363 patients
# Example 4 (JFM): power estimation, testing a treatment effect on recurrent
# events only. We assume a uniformly distributed number of recurrent events,
# ranging from 1 to 6 recurrent events per subject. The allocation ratio
# experimental:control is 2:1, and the follow-up is 10 weeks.
# The expected hazard ratio is 0.70 for recurrent event and 0.90 for the
# terminal event. We have chosen 0.5 as the variance of the frailties.
JFM.power(
Npts = 400, ni = c(1, 6), ni.type = "unif",
FUP = 10, FUP.type = "fixed", ratio = 2,
betaTest.type = "betaRtest", betaR.HA = log(.70), betaD.HA = log(.90),
theta = .5
) # power ~ 76%
# Example 5 (JFM): sample size calculation, to assess the treatment effect on
# both recurrent and terminal events. We want to achieve an 80% power.
# We anticipate a maximum of 5 recurrent events, over a 6-year period and a
# 0.5-year accrual period. We assume that the gamma-frailty variance is 0.5.
# For the control group, we expect a 2-year and a 5-year median time-to-event
# for recurrent events and terminal events, respectively. We consider a 30%
# and 20% risk reduction for recurrent events and terminal event, respectively.
JFM.ssize(
power = 0.80, ni = 9,
FUP = 6, Acc.Dur = 1.5, medianR.H0 = 2, medianD.H0 = 5,
betaTest.type = "joint", betaR.HA = log(.70), betaD.HA = log(.80), theta = .5
) # sample size ~ 445 patients / ~ approx
# Example 6: Sample size calculation for GJFM
# Same as above, but with two random effects (with two Gamma-frailty variances
# theta and eta). To ensure sample size estimation stability, we use 10000
# Monte-Carlo samples.
GJFM.ssize(
power = 0.80, ni = 5,
FUP = 6, Acc.Dur = 0.5, medianR.H0 = 2, medianD.H0 = 5,
betaR.HA = log(0.70), betaD.HA = log(0.80), theta = 0.5, eta = 0.75,
samples.mc = 1e5
) # sample size ~ 705 patients / ~ approx 4 min.
# Example 7:
# --------------------------------------------------
# Post-hoc power analysis for a Joint Frailty Model
# --------------------------------------------------
# See original article by Gonzalez et al. (2005)
data(readmission)
modJFM <- frailtyPenal(
Surv(time, event) ~ cluster(id) + as.factor(chemo) + terminal(death),
formula.terminalEvent = ~ as.factor(chemo), data = readmission,
hazard = "Weibull"
)
# Test both recurrent and death events
# # - Let us assume an underlying Poisson distribution for ni.type. The
# # empirical mean of the number of recurrent events per patients is: ni = 1.136476.
# # - For the null hypothesis, let us consider betaR.H0 = betaD.H0 = 0. For the
# # alternative hypothesis, we use the estimated parameters for betaD.HA and
# # betaR.HA.
# # - "Patients were actively followed up until June 2002" -> the follow-up
# # type is "UpToEnd".
# # - "The study took place in the Hospital de Bellvitge, [...] between
# # January 1996 and December 1998" -> the accrual time is approximately 3 years
# # - We can assume that the study duration is approximately 6 years
ni <- 1.136476
ni.type <- "Pois"
Acc.Dur <- 3 * 365.25 # time unit = days
FUP <- 6 * 365.25 # same as above
betaR.HA <- as.numeric(modJFM$coef[1]) # else "Named numeric"
betaD.HA <- as.numeric(modJFM$coef[2]) # same as above
med <- modJFM$scale.weib * log(2)^(1 / modJFM$shape.weib)
medianR.H0 <- med[1]
medianD.H0 <- med[2]
shapeR.W <- modJFM$shape[1]
shapeD.W <- modJFM$shape[2]
theta <- modJFM$theta
alpha <- modJFM$alpha
Npts <- length(unique(readmission[, "id"])) # 403 patients
nTreated <- length(unique(readmission[readmission$chemo == "Treated", "id"])) #217 treated patients
ratio <- nTreated / (Npts - nTreated)
JFM.power(
Npts = Npts, ni = ni, ni.type = ni.type,
Acc.Dur = Acc.Dur, FUP = FUP, medianR.H0 = medianR.H0, medianD.H0 = medianD.H0,
betaTest.type = "joint", betaR.HA = betaR.HA, betaD.HA = betaD.HA,
shapeR.W = shapeR.W, shapeD.W = shapeD.W, theta = theta, alpha = alpha,
ratio = ratio
) # power ~ 92%
# --------------------------------------------------
# Required sample size under the same setting
# --------------------------------------------------
# Here, let us consider that readmission is a “pilot study” with 403 patients,
# from which we estimate parameters. Under this scenario, let us compute the
# needed sample size, but to achieve an 80% power.
JFM.ssize(
power = 0.80, ni = ni, ni.type = ni.type,
Acc.Dur = Acc.Dur, FUP = FUP, medianR.H0 = medianR.H0, medianD.H0 = medianD.H0,
betaTest.type = "joint", betaR.HA = betaR.HA, betaD.HA = betaD.HA,
shapeR.W = shapeR.W, shapeD.W = shapeD.W, theta = theta, alpha = alpha,
ratio = ratio
) # 289 patients needed under the same settings vs. 403
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.