,
Helmut Schütz [aut] ,
Benjamin Lang [aut]| Title: | Power and Sample Size for (Bio)Equivalence Studies |
|---|---|
| Description: | Contains functions to calculate power and sample size for various study designs used in bioequivalence studies. Use known.designs() to see the designs supported. Power and sample size can be obtained based on different methods, amongst them prominently the TOST procedure (two one-sided t-tests). See README and NEWS for further information. |
| Authors: | Detlew Labes [aut, cre] ,
Helmut Schütz [aut] ,
Benjamin Lang [aut] |
| Maintainer: | Detlew Labes <[email protected]> |
| License: | GPL (>=2) |
| Version: | 1.5-6 |
| Built: | 2024-10-15 04:00:59 UTC |
| Source: | https://github.com/detlew/powertost |
This function returns the ‘design’ matrix of incomplete
block designs described by Chow & Liu. The design matrices were
recoded 1=R, 2=T1, 3=T2, ...
bib.CL(trt, p)bib.CL(trt, p)
trt |
Number of treatments ( |
p |
Number of periods ( |
Matrix containing the sequences in rows and periods in columns.
The entry (i, j) of the matrix corresponds to the treatment or
dose (index) a subject within i-th sequence gets in the
j-th period.
D. Labes
Chow SC, Liu JP. Design and Analysis of Bioavailability and Bioequivalence Studies. Boca Raton: CRC Press; 3rd edition 2009. Chapter 2.6.
# 4 treatments/doses, 3 periods bib.CL(4, 3) # gives 4 sequences # to see this in Chow & Liu's coding tmt <- c("R", "T1", "T2", "T3") matrix(tmt[bib.CL(4, 3)], ncol=3)# 4 treatments/doses, 3 periods bib.CL(4, 3) # gives 4 sequences # to see this in Chow & Liu's coding tmt <- c("R", "T1", "T2", "T3") matrix(tmt[bib.CL(4, 3)], ncol=3)
Utility function to calculate the 1–2α CI given point estimate, CV, and n for the various designs covered in this package.
CI.BE(alpha = 0.05, pe, CV, n, design = "2x2", robust = FALSE)CI.BE(alpha = 0.05, pe, CV, n, design = "2x2", robust = FALSE)
alpha |
Type I error probability, significance level. Defaults to 0.05. |
pe |
Point estimate (GMR). |
CV |
Coefficient of variation as ratio (not percent). |
n |
Total number of subjects if a scalar is given. |
design |
Character string describing the study’s design. |
robust |
Defaults to |
Returns the 1–2α
confidence interval.
Returns a vector with named elements lower, upper if
arguments pe and CV are scalars, else a matrix with
columns lower, upper is returned.
The function assumes an evaluation using log-transformed data.
The function assumes equal variances in case of design="parallel"
and the higher order crossover designs.
The implemented formula covers balanced and unbalanced designs.
Whether the function vectorizes properly is not thoroughly tested.
D. Labes
# 90% confidence interval for the 2x2 crossover # n(total) = 24 CI.BE(pe = 0.95, CV = 0.3, n = 24) # should give # lower upper # 0.8213465 1.0988055 # same total number but unequal sequences CI.BE(pe = 0.95, CV = 0.3, n = c(13, 11)) # lower upper # 0.8209294 1.0993637# 90% confidence interval for the 2x2 crossover # n(total) = 24 CI.BE(pe = 0.95, CV = 0.3, n = 24) # should give # lower upper # 0.8213465 1.0988055 # same total number but unequal sequences CI.BE(pe = 0.95, CV = 0.3, n = c(13, 11)) # lower upper # 0.8209294 1.0993637
Utility function to calculate the 1–2α Fieller confidence interval given the point estimate, CV (, CVb), and n for the parallel group and 2×2 crossover.
CI.RatioF(alpha = 0.025, pe, CV, CVb, n, design = c("2x2", "parallel"))CI.RatioF(alpha = 0.025, pe, CV, CVb, n, design = c("2x2", "parallel"))
alpha |
Type I error probability, aka significance level. |
pe |
Point estimate of T/R ratio. |
CV |
Coefficient of variation as ratio (not percent). In case of |
CVb |
CV of the between-subject variability. Only necessary for
|
n |
Total number of subjects if a scalar is given. |
design |
A character string describing the study design. |
The CV(within) and CVb(etween) in case of design="2x2" are
obtainedvia an appropriate ANOVA from the error term and from the
difference (MS(subject within sequence)-MS(error))/2.
Returns the 1–2α confidence interval.
The function assumes an evaluation using un-transformed data.
The function assumes equal variances in case of design="parallel".
The formula implemented covers balanced and unbalanced designs.
Note that when the mean of the denominator of the ratio is close to zero,
confidence intervals might be degenerated and are returned as NA. In such
a case a warning is issued.
Whether the function vectorizes properly is not thoroughly tested.
This function is intended for studies with clinical endpoints. In such
studies the 95% confidence intervals are usually used for equivalence
testing. Therefore, alpha defaults here to 0.025 (see EMA 2000).
D. Labes
Locke CS. An exact confidence interval from untransformed data for the ratio of two formulation means. J Pharmacokin Biopharm. 1984;12(6):649–55. doi:10.1007/BF01059558
Hauschke D, Steinijans VW, Pigeot I. Bioequivalence Studies in Drug Development. Chichester: John Wiley; 2007. Chapter 10. doi:10.1002/9780470094778.fmatter
European Medicines Agency, Committee for Proprietary Medicinal Products. Points to consider on switching between superiority and non-inferiority. London, 27 July 2000. CPMP/EWP/482/99
# 95% Fieller CI for the 2x2 crossover CI.RatioF(pe = 0.95, CV = 0.3, CVb = 0.6, n = 32)# 95% Fieller CI for the 2x2 crossover CI.RatioF(pe = 0.95, CV = 0.3, CVb = 0.6, n = 32)
These data.frames give sample size tables calculated with
sampleN.TOST() for the 2×2 design.
The data.frames can be accessed by their names.
| data.frame | Description |
| ct5.1 | Multiplicative model, theta1=0.8, theta2=1.25 (1/theta1), exact |
| ct5.2 | Multiplicative model, theta1=0.75, theta2=1.3333 (1/theta1), exact |
| ct5.3 | Multiplicative model, theta1=0.9, theta2=1.1111 (1/theta1), exact |
| ct5.4.1 | Additive model, theta1=--0.2, theta2=+0.2 (BE limits 0.80 -- 1.20), exact |
Scripts for creation of these data.frames can be found in the /tests
sub-directory of the package.
Comparing the results of these scripts to the corresponding data.frames can
be used for validation purposes.
PowerTOST
| data.frame | Origin | Details |
| ct5.1 | Hauschke et al. | Table 5.1 (p 113--114) |
| ct5.2 | Hauschke et al. | Table 5.2 (p 115--116) |
| ct5.3 | Hauschke et al. | Table 5.3 (p 118) |
| ct5.4.1 | Chow & Liu | Table 5.4.1 (p 158) |
Hauschke D, Steinijans VW, Pigeot I. Bioequivalence Studies in Drug Development. Chichester: John Wiley; 2007.
Chow SC, Liu JP. Design and Analysis of Bioavailability and Bioequivalence Studies. Boca Raton: CRC Press; 3rd edition 2009.
ct5.1 ct5.2 ct5.3 ct5.4.1ct5.1 ct5.2 ct5.3 ct5.4.1
These data.frames give sample size tables calculated with
sampleN.TOST() for the 2×2×3 replicate crossover design
(2-treatment 2-sequence 3-period design.
The data.frames can be accessed by their names.
| data.frame | Description |
| ct9.6.2 | Additive model, theta1=--0.2, theta2=+0.2 (BE limits 0.80 -- 1.20) |
| approximate power via shifted non-central t-distribution | |
| ct9.6.6 | Multiplicative model, theta1=0.8, theta2=1.25 (1/theta1) |
| approximate power via shifted non-central t-distribution |
Attention! CV is se (standard error) of residuals.
Scripts for creation of these data.frames can be found in the /tests
sub-directory of the package.
Comparing the results of these scripts to the corresponding data.frames can
be used for validation purposes.
PowerTOST
| data.frame | Origin | Details |
| ct9.6.2 | Chow & Liu | Table 9.6.2 (p 292) |
| ct9.6.6 | Chow & Liu | Table 9.6.6 (p 293) |
Chow SC, Liu JP. Design and Analysis of Bioavailability and Bioequivalence Studies. Boca Raton: CRC Press; 3rd edition 2009.
ct9.6.2 ct9.6.6ct9.6.2 ct9.6.6
These data.frames give sample size tables calculated with
sampleN.TOST() for the 2×4×4 replicate crossover design
(2-treatment 4-sequence 4-period design).
The data.frames can be accessed by their names.
| data.frame | Description |
| ct9.6.4 | Additive model, theta1=--0.2, theta2=+0.2 (BE limits 0.80 -- 1.20) |
| approximate power via shifted non-central t-distribution | |
| ct9.6.8 | Multiplicative model, theta1=0.8, theta2=1.25 (1/theta1) |
| approximate power via shifted non-central t-distribution |
Attention! CV is se (standard error) of residuals.
Scripts for creation of these data.frames can be found in the /tests
sub-directory of the package.
Comparing the results of these scripts to the corresponding data.frames can
be used for validation purposes.
PowerTOST
| data.frame | Origin | Details |
| ct9.6.4 | Chow & Liu | Table 9.6.4 (p 294) |
| ct9.6.8 | Chow & Liu | Table 9.6.8 (p 298) |
Chow SC, Liu JP. Design and Analysis of Bioavailability and Bioequivalence Studies. Boca Raton: CRC Press; 3rd edition 2009.
ct9.6.4 ct9.6.8ct9.6.4 ct9.6.8
These data.frames give sample size tables calculated with
sampleN.TOST() for the parallel group design (2 groups).
The data.frames can be accessed by their names.
| data.frame | Description |
| ctSJ.VIII.10 | Multiplicative model, theta1=0.9, theta2=1.1111 (1/theta1), target power=90% |
| approximate power via non-central t-distribution | |
| ctSJ.VIII.20 | Multiplicative model, theta1=0.8, theta2=1.25 (1/theta1), target power=90% |
| approximate power via non-central t-distribution | |
| ctCW.III | Additive model, theta1=--0.2, theta2=+0.2 (BE limits 0.80 -- 1.20), exact |
Attention! Julious gives sample size per group.
Scripts for creation of these data.frames can be found in the /tests
sub-directory of the package.
Comparing the results of these scripts to the corresponding data.frames can
be used for validation purposes.
PowerTOST
| data.frame | Origin | Details |
| ctSJ.VIII.10 | Julious | Table VIII (p. 1972), column ‘Level of bioequivalence 10%’ |
| ctSJ.VIII.20 | Julious | Table VIII (p. 1972), column ‘Level of bioequivalence 20%’ |
| ctCW.III | Chow & Wang | Table III (p. 164) |
Seems the last reference is not very reliable (compare to the table in the paper).
Julious SA. Tutorial in Biostatistics. Sample sizes for clinical trials with Normal data. Stat Med. 2004;23(12):1921–86. doi:10.1002/sim.1783
Chow SC, Wang H. On Sample Size Calculation in Bioequivalence Trials. J Pharmacokinet Pharmacodyn. 2001;28(2):155–69. doi:10.1023/A:1011503032353
ctSJ.VIII.10 ctSJ.VIII.20 ctCW.IIIctSJ.VIII.10 ctSJ.VIII.20 ctCW.III
Calculates the standard error or the mean squared error from a given CV and vice versa for log-normal data.
CV2se(CV) se2CV(se) CV2mse(CV) mse2CV(mse)CV2se(CV) se2CV(se) CV2mse(CV) mse2CV(mse)
CV |
coefficient of variatio as ratio (not percent) |
se |
standard error |
mse |
mean squared error (aka residual variance) |
Returns se = sqrt(log(CV^2+1)) CV = sqrt(exp(se^2)-1) mse = log(CV^2+1) CV = sqrt(exp(mse)-1)
These functions were originally intended for internal use only but may be useful for others.
D. Labes
# these functions are one liners: CV2se <- function(CV) return(sqrt(log(1.0 + CV^2))) se2CV <- function(se) return(sqrt(exp(se^2)-1)) CV2se(0.3) # should give: # [1] 0.2935604 se2CV(0.2935604) # [1] 0.3# these functions are one liners: CV2se <- function(CV) return(sqrt(log(1.0 + CV^2))) se2CV <- function(se) return(sqrt(exp(se^2)-1)) CV2se(0.3) # should give: # [1] 0.2935604 se2CV(0.2935604) # [1] 0.3
The function calculates the 1–α confidence limits (either 1-sided or 2-sided) via the χ2 distribution of the error variance the CV is based on.
CVCL(CV, df, side = c("upper", "lower", "2-sided"), alpha = 0.05)CVCL(CV, df, side = c("upper", "lower", "2-sided"), alpha = 0.05)
CV |
Coefficient of variation as ratio (not percent) |
df |
degrees of freedom of the CV (error variance) |
side |
Side(s) to calculate the confidence limits for, defaults to |
alpha |
Type I error probability, aka significance level |
Numeric vector of the confidence limits named as lower CL and upper CL.
In case of the one-sided upper confidence limit the lower CL is = 0.
In case of the one-sided lower confidence limit the upper CL is = Inf.
D. Labes
# upper one-sided 95% CL of a CV=0.3 # from a study with df=22 (f.i. a 2x2 crossover with n=24) # default side="upper" since not explicitly given CVCL(0.3, df = 22) # should give: # lower CL upper CL # 0.0000000 0.4075525# upper one-sided 95% CL of a CV=0.3 # from a study with df=22 (f.i. a 2x2 crossover with n=24) # default side="upper" since not explicitly given CVCL(0.3, df = 22) # should give: # lower CL upper CL # 0.0000000 0.4075525
Calculates the CV (coefficient of variation) from a known confidence interval
of a BE study.
Useful if no CV but the 90% CI was given in literature.
CVfromCI(pe, lower, upper, n, design = "2x2", alpha = 0.05, robust = FALSE) CI2CV(pe, lower, upper, n, design = "2x2", alpha = 0.05, robust = FALSE)CVfromCI(pe, lower, upper, n, design = "2x2", alpha = 0.05, robust = FALSE) CI2CV(pe, lower, upper, n, design = "2x2", alpha = 0.05, robust = FALSE)
pe |
Point estimate of the T/R ratio. |
lower |
Lower confidence limit of the BE ratio. |
upper |
Upper confidence limit of the BE ratio. |
n |
Total number of subjects under study if given as scalar. |
design |
Character string describing the study design. |
alpha |
Error probability. Set it to |
robust |
With |
See Helmut Schütz’ presentation for the algebra underlying this function.
Numeric value of the CV as ratio.
The calculations are based on the assumption of evaluation via log-transformed values.
The calculations are further based on a common variance of Test and Reference
treatments in replicate crossover studies or parallel group study, respectively.
In case of argument n given as n(total) and is not divisible by the number
of (sequence) groups the total sample size is partitioned to the (sequence) groups
to have small imbalance only. A message is given in such cases.
The estimated CV is conservative (i.e., higher than actually observed) in case of
unbalancedness.CI2CV() is simply an alias to CVfromCI().
Original by D. Labes with suggestions by H. Schütz.
Reworked and adapted to unbalanced studies by B. Lang.
Yuan J, Tong T, Tang M-L. Sample Size Calculation for Bioequivalence Studies Assessing Drug Effect and Food Effect at the Same Time With a 3-Treatment Williams Design. Regul Sci. 2013;47(2):242–7. doi:10.1177/2168479012474273
# Given a 90% confidence interval (without point estimate) # from a classical 2x2 crossover with 22 subjects CVfromCI(lower=0.91, upper=1.15, n=22, design="2x2") # will give [1] 0.2279405, i.e a CV ~ 23% # # unbalanced 2x2 crossover study, but not reported as such CI2CV(lower=0.89, upper=1.15, n=24) # will give a CV ~ 26.3% # unbalancedness accounted for CI2CV(lower=0.89, upper=1.15, n=c(16,8)) # should give CV ~ 24.7%# Given a 90% confidence interval (without point estimate) # from a classical 2x2 crossover with 22 subjects CVfromCI(lower=0.91, upper=1.15, n=22, design="2x2") # will give [1] 0.2279405, i.e a CV ~ 23% # # unbalanced 2x2 crossover study, but not reported as such CI2CV(lower=0.89, upper=1.15, n=24) # will give a CV ~ 26.3% # unbalancedness accounted for CI2CV(lower=0.89, upper=1.15, n=c(16,8)) # should give CV ~ 24.7%
Helper function to calculate CV(T) and CV(R) from a pooled CV(T/R) assuming a ratio of the intra-subject variances.
CVp2CV(CV, ratio = 1.5)CVp2CV(CV, ratio = 1.5)
CV |
‘pooled’ CV of T and R (as ratio, not percent). |
ratio |
Ratio of the intra-subject variances |
In case of knowing only the CV(T/R) f.i. from an ordinary cross-over you can
calculate the components CV(T) and CV(R) assuming a ratio of the
intra-subject variances.
The formula the function is based on:log(1.0 + CV^2) = (sWT^2 + sWR^2)/2
Insert sWT^2 = ratio * sWR^2 and solve for sWR^2.
Returns a numeric vector of the CV values for Test and Reference if only
one ratio is given.
Returns a matrix with named columns CVwT and CVwR if ratio is given as vector.
D. Labes
CVp2CV(0.4, ratio=2) # gives # [1] 0.4677952 0.3225018CVp2CV(0.4, ratio=2) # gives # [1] 0.4677952 0.3225018
This function pools CVs of several studies.
CVpooled(CVdata, alpha = 0.2, logscale = TRUE, robust = FALSE) ## S3 method for class 'CVp' print(x, digits = 4, verbose = FALSE, ...)CVpooled(CVdata, alpha = 0.2, logscale = TRUE, robust = FALSE) ## S3 method for class 'CVp' print(x, digits = 4, verbose = FALSE, ...)
CVdata |
A data.frame that must contain the columns |
alpha |
Error probability for calculating an upper confidence limit of the pooled CV. |
logscale |
Should the calculations be done for log-transformed data? Defaults to |
robust |
Defaults to |
x |
An object of class |
digits |
Number of significant digits for the |
verbose |
Defaults to |
... |
More args to print(). None used. |
The pooled CV is obtained from the weighted average of the error variances obtained from
the CVs of the single studies, weights are the degrees of freedom df.
If only n is given in the input CVdata, the dfs are calculated via
the formulas given in known.designs(). If both n and df are given
the df column precedes.
If logscale=TRUE the error variances are obtained via function CV2se().
Otherwise the pooled CV is obtained via pooling the CV^2.
A list of class "CVp" with components
CV |
value of the pooled CV |
df |
pooled degrees of freedom |
CVupper |
upper confidence interval of the pooled CV |
alpha |
input value |
The class "CVp" has a S3 methods print.CVp.
Pooling of CVs from parallel group and crossover designs does not make any sense.
Also the function does not throw an error if you do so.
The calculations for logscale=FALSE are not described in the references.
They are implemented by analogy to the case via log-transformed data.
The calculations are based on a common variance of Test and Reference formulations
in replicate crossover studies or a parallel group study, respectively.
D. Labes
H. Schütz’ presentations about sample size challenges.
Patterson S, Jones B. Bioequivalence and Statistics in Clinical Pharmacology. Boca Raton: Chapman & Hall / CRC Press; 2nd edition 2017. Chapter 5.7 “Determining Trial Size”.
# some data: # the values for AUC, study 1 and study 2 are Example 3 of H. Schuetz' presentation CVs <- (" PKmetric | CV | n |design| source AUC | 0.20 | 24 | 2x2 | study 1 Cmax | 0.25 | 24 | 2x2 | study 1 AUC | 0.30 | 12 | 2x2 | study 2 Cmax | 0.31 | 12 | 2x2 | study 2 AUC | 0.25 | 12 | 2x2x4| study 3 (full replicate) ") txtcon <- textConnection(CVs) CVdata <- read.table(txtcon, header = TRUE, sep = "|", strip.white = TRUE, as.is = TRUE) close(txtcon) # evaluation of the AUC CVs CVsAUC <- subset(CVdata, PKmetric == "AUC") CVpooled(CVsAUC, alpha = 0.2, logscale = TRUE) # df of the 'robust' evaluation CVpooled(CVsAUC, alpha = 0.2, logscale = TRUE, robust = TRUE) # print also the upper CL, data example 3 CVsAUC3 <- subset(CVsAUC,design != "2x2x4") print(CVpooled(CVsAUC3, alpha = 0.2, robust = TRUE), digits = 3, verbose = TRUE) # will give the output: # Pooled CV = 0.235 with 32 degrees of freedom (robust dfs) # Upper 80% confidence limit of CV = 0.266 # # Combining CVs from studies evaluated by ANOVA (robust=FALSE) and # by a mixed effects model (robust=TRUE). dfs have to be provided! CVs <- (" CV | n |design| source | model | df 0.212 | 24 | 2x2 | study 1 | fixed | 22 0.157 | 27 | 3x3 | study 2 | fixed | 50 0.148 | 27 | 3x3 | study 3 | mixed | 24 ") txtcon <- textConnection(CVs) CVdata <- read.table(txtcon, header = TRUE, sep = "|", strip.white = TRUE, as.is = TRUE) close(txtcon) print(CVpooled(CVdata, alpha = 0.2), digits = 3, verbose = TRUE) # will give the output: # Pooled CV = 0.169 with 96 degrees of freedom # Upper 80% confidence limit of CV = 0.181# some data: # the values for AUC, study 1 and study 2 are Example 3 of H. Schuetz' presentation CVs <- (" PKmetric | CV | n |design| source AUC | 0.20 | 24 | 2x2 | study 1 Cmax | 0.25 | 24 | 2x2 | study 1 AUC | 0.30 | 12 | 2x2 | study 2 Cmax | 0.31 | 12 | 2x2 | study 2 AUC | 0.25 | 12 | 2x2x4| study 3 (full replicate) ") txtcon <- textConnection(CVs) CVdata <- read.table(txtcon, header = TRUE, sep = "|", strip.white = TRUE, as.is = TRUE) close(txtcon) # evaluation of the AUC CVs CVsAUC <- subset(CVdata, PKmetric == "AUC") CVpooled(CVsAUC, alpha = 0.2, logscale = TRUE) # df of the 'robust' evaluation CVpooled(CVsAUC, alpha = 0.2, logscale = TRUE, robust = TRUE) # print also the upper CL, data example 3 CVsAUC3 <- subset(CVsAUC,design != "2x2x4") print(CVpooled(CVsAUC3, alpha = 0.2, robust = TRUE), digits = 3, verbose = TRUE) # will give the output: # Pooled CV = 0.235 with 32 degrees of freedom (robust dfs) # Upper 80% confidence limit of CV = 0.266 # # Combining CVs from studies evaluated by ANOVA (robust=FALSE) and # by a mixed effects model (robust=TRUE). dfs have to be provided! CVs <- (" CV | n |design| source | model | df 0.212 | 24 | 2x2 | study 1 | fixed | 22 0.157 | 27 | 3x3 | study 2 | fixed | 50 0.148 | 27 | 3x3 | study 3 | mixed | 24 ") txtcon <- textConnection(CVs) CVdata <- read.table(txtcon, header = TRUE, sep = "|", strip.white = TRUE, as.is = TRUE) close(txtcon) print(CVpooled(CVdata, alpha = 0.2), digits = 3, verbose = TRUE) # will give the output: # Pooled CV = 0.169 with 96 degrees of freedom # Upper 80% confidence limit of CV = 0.181
Calculates the intra-subject CV (coefficient of variation) of the reference from the upper expanded limit of a BE study (replicate design for ABEL). Useful if no CVwR but the expanded limits were given.
CVwRfromU(U, regulator = "EMA") U2CVwR(U, regulator = "EMA")CVwRfromU(U, regulator = "EMA") U2CVwR(U, regulator = "EMA")
U |
Upper expanded limit. |
regulator |
Regulatory body’s settings for expanding the BE acceptance limits,
given as a string from the choices |
Only the upper expanded limit is supported since it offers one more significant digit than the lower expanded limit.
Numeric value of the CVwR as ratio, where
CVwR = sqrt(exp((log(U)/r_const)^2)-1).
U2CVwR() is simply an alias to CVwRfromU().
H. Schütz
# Given the upper expanded limit and using the defaults CVwRfromU(U = 1.38) # should give [1] 0.44355, i.e., a CVwR ~ 44% # Upper limit from a study according the Health Canada’s rules CVwRfromU(U = 1.48, regulator = "HC") # should give [1] 0.55214# Given the upper expanded limit and using the defaults CVwRfromU(U = 1.38) # should give [1] 0.44355, i.e., a CVwR ~ 44% # Upper limit from a study according the Health Canada’s rules CVwRfromU(U = 1.48, regulator = "HC") # should give [1] 0.55214
Details about removed functions in PowerTOST.
type1error.2TOST was designed to calculate the type I error (TIE) rate of two simultaneous TOST procedures with some specified correlation of the parameters. It suffered from poor precision to obtain the TIE via simulations.
Due to the intersection-union principle the TIE is always upper bounded to α by theory and hence, type1error.2TOST was removed in version 1.4-7.
power.scABEL2 and sampleN.scABEL2 were deprecated in version 1.4-1. In power.scABEL and sampleN.scABEL the regulator component "est_method" is used for switching between simulations based on the EMA’s ANOVA evaluation or ISC evaluation, respectively.power.scABEL2 and sampleN.scABEL2 were removed in version 1.4-3.
power2.TOST was deprecated in version 1.2-6 since power.TOST was modified in order to handle unbalanced sequences. power2.TOST was removed in version in version 1.2-7.
The functions power.NTIDFDA, sampleN.NTIDFDA and pa.NTIDFDA were deprecated in
version 1.5-4 and were removed in version 1.5-6.
Berger RL, Hsu JC. Bioequivalence Trials, Intersection-Union Tests and Equivalence Confidence Sets. Stat Sci. 1996;11(4):283–302. JSTOR:2246021
power.TOST, power.scABEL, sampleN.scABEL, power.NTID, sampleN.NTID, and pa.NTID for the new functions.
Details about deprecated functions in PowerTOST.
The functions power.NTIDFDA, sampleN.NTIDFDA, and pa.NTIDFDA were deprecated in version 1.5-4.
power.NTID, sampleN.NTID, pa.NTID for the new functions.
Calculates the so-called expected, i.e., unconditional, power for a variety of study designs used in bioequivalence studies.
exppower.noninf(alpha = 0.025, logscale = TRUE, theta0, margin, CV, n, design = "2x2", robust = FALSE, prior.type = c("CV", "theta0", "both"), prior.parm = list(), method = c("exact", "approx"))exppower.noninf(alpha = 0.025, logscale = TRUE, theta0, margin, CV, n, design = "2x2", robust = FALSE, prior.type = c("CV", "theta0", "both"), prior.parm = list(), method = c("exact", "approx"))
alpha |
Significance level (one-sided). Defaults here to 0.025. |
logscale |
Should the data be used on log-transformed or on original scale? |
theta0 |
Assumed ‘true’ (or ‘observed’ in case of |
margin |
Non-inferiority margin. |
CV |
In case of In case of cross-over studies this is the within-subject CV, in case of a parallel-group design the CV of the total variability. |
n |
Number of subjects under study. |
design |
Character string describing the study design.
See |
robust |
Defaults to FALSE. Set to |
prior.type |
Specifies which parameter uncertainty should be accounted for. In case of
|
prior.parm |
A list of parameters expressing the prior information about the
variability and/or treatment effect. Possible components are |
method |
Defaults to |
This function calculates the so-called expected power taking into account that
usually the parameters (CV and/or theta0) are not known but estimated
from a prior study with some uncertainty. The expected power is an unconditional
power and can therefore be seen as probability for success. See references for
further details.
The prior.parm argument is a list that can supply any of the following
components:
dfError degrees of freedom from the prior trial (>4, maybe non-integer).
df = Inf is allowed and for method = "exact" the result will then
coincide with power.noninf(...).
Note: This corresponds to the df of both the CV and the difference of means.
SEMStandard error of the difference of means from the prior trial; must always be on additive scale (i.e., usually log-scale).
mNumber of subjects from prior trial. Specification is analogous to
the main argument n.
designStudy design of prior trial. Specification is analogous to the
main argument design.
For prior.parm, the combination consisting of df and SEM
requires a somewhat advanced knowledge of the prior trial (provided in the raw
output from for example the software SAS, or may be obtained via
emmeans of package emmeans. However, it has the
advantage that if there were missing data the exact degrees of freedom and
standard error of the difference can be used, the former possibly being
non-integer valued (e.g. if the Kenward-Roger method was used).
Details on argument prior.type:
CVThe expectation is calculated with respect to the Inverse-gamma distribution.
theta0The expectation is calculated with respect to the
conditional distribution theta0 | = s^2
of the posteriori distribution of (theta0, ) from the prior
trial.
bothThe expectation is calculated with respect to the posteriori
distribution of (theta0, ) from the prior trial. Numerical calculation
of the two-dimensional integral is performed via cubature::hcubature.
Notes on the underlying hypotheses
If the supplied margin is < 0 (logscale=FALSE) or < 1 (logscale=TRUE),
then it is assumed higher response values are better. The hypotheses are H0: theta0 <= margin vs. H1: theta0 > margin
where theta0 = mean(test)-mean(reference) if logscale=FALSE
or H0: log(theta0) <= log(margin) vs. H1: log(theta0) > log(margin)
where theta0 = mean(test)/mean(reference) if logscale=TRUE.
If the supplied margin is > 0 (logscale=FALSE) or > 1 (logscale=TRUE),
then it is assumed lower response values are better. The hypotheses are H0: theta0 >= margin vs. H1: theta0 < margin
where theta0 = mean(test)-mean(reference) if logscale=FALSE
or H0: log(theta0) >= log(margin) vs. H1: log(theta0) < log(margin)
where theta0 = mean(test)/mean(reference) if logscale=TRUE.
This latter case may also be considered as ‘non-superiority’.
Value of expected power according to the input.
B. Lang, D. Labes
Grieve AP. Confidence Intervals and Sample Sizes. Biometrics. 1991;47:1597–603. doi:10.2307/2532411
O’Hagan, Stevens, JW, Campell MJ. Assurance in Clinical Trial Design. Pharm Stat. 2005;4:187–201. doi:10.1002/pst.175
Julious SA, Owen RJ. Sample size calculations for clinical studies allowing for uncertainty in variance. Pharm Stat. 2006;5:29–37. doi:10.1002/pst.197
Julious SA. Sample sizes for Clinical Trials. Boca Raton: CRC Press / Chapman & Hall; 2010.
Bertsche A, Nehmitz G, Beyersmann J, Grieve AP. The predictive distribution of the residual variability in the linear-fixed effects model for clinical cross-over trials. Biom J. 2016;58(4):797–809. doi:10.1002/bimj.201500245
Box GEP, Tiao GC. Bayesian Inference in Statistical Analysis. Boston: Addison-Wesley; 1992.
Held L, Sabanes Bove D. Applied Statistical Inference. Likelihood and Bayes. Berlin, Heidelberg: Springer; 2014. doi:10.1007/978-3-642-37887-4
Senn S. Cross-over Trials in Clinical Research. Chichester: John Wiley & Sons; 2nd edition 2002.
Zierhut ML, Bycott P, Gibbs MA, Smith BP, Vicini P. Ignorance is not bliss: Statistical power is not probability of trial success. Clin Pharmacol Ther. 2015;99:356–9. doi:10.1002/cpt.257
expsampleN.noninf, power.noninf
# Expected power for non-inferiority test for a 2x2 crossover # with 40 subjects. CV 30% known from a pilot 2x2 study with # 12 subjects # using all the defaults for other parameters (theta0 carved in stone) # should give: [1] 0.6761068 exppower.noninf(CV = 0.3, n = 40, prior.parm = list(df = 12-2)) # or equivalently exppower.noninf(CV = 0.3, n = 40, prior.parm = list(m = 12, design = "2x2")) # May be also calculated via exppower.TOST() after setting upper acceptance limit # to Inf and alpha=0.025 exppower.TOST(CV = 0.3, n = 40, prior.parm = list(df = 10), theta2 = Inf, alpha=0.025) # In contrast: Julious approximation exppower.noninf(CV = 0.3, n = 40, prior.parm = list(df = 10), method = "approx") # should give: [1] 0.6751358 # Compare this to the usual (conditional) power (CV known, "carved in stone") power.noninf(CV = 0.3, n = 40) # should give: [1] 0.7228685 # same as if setting df = Inf in function exppower.noninf() exppower.noninf(CV = 0.3, n = 40, prior.parm = list(df = Inf)) # Expected power for a 2x2 crossover with 40 subjects # CV 30% and theta0 = 0.95 known from a pilot 2x2 study with 12 subjects # using uncertainty with respect to both CV and theta0 exppower.noninf(CV = 0.3, theta0 = 0.95, n = 40, prior.parm = list(m = 12, design = "2x2"), prior.type = "both") # should give a decrease of expected power to 0.5982852# Expected power for non-inferiority test for a 2x2 crossover # with 40 subjects. CV 30% known from a pilot 2x2 study with # 12 subjects # using all the defaults for other parameters (theta0 carved in stone) # should give: [1] 0.6761068 exppower.noninf(CV = 0.3, n = 40, prior.parm = list(df = 12-2)) # or equivalently exppower.noninf(CV = 0.3, n = 40, prior.parm = list(m = 12, design = "2x2")) # May be also calculated via exppower.TOST() after setting upper acceptance limit # to Inf and alpha=0.025 exppower.TOST(CV = 0.3, n = 40, prior.parm = list(df = 10), theta2 = Inf, alpha=0.025) # In contrast: Julious approximation exppower.noninf(CV = 0.3, n = 40, prior.parm = list(df = 10), method = "approx") # should give: [1] 0.6751358 # Compare this to the usual (conditional) power (CV known, "carved in stone") power.noninf(CV = 0.3, n = 40) # should give: [1] 0.7228685 # same as if setting df = Inf in function exppower.noninf() exppower.noninf(CV = 0.3, n = 40, prior.parm = list(df = Inf)) # Expected power for a 2x2 crossover with 40 subjects # CV 30% and theta0 = 0.95 known from a pilot 2x2 study with 12 subjects # using uncertainty with respect to both CV and theta0 exppower.noninf(CV = 0.3, theta0 = 0.95, n = 40, prior.parm = list(m = 12, design = "2x2"), prior.type = "both") # should give a decrease of expected power to 0.5982852
Calculates the so-called expected, i.e., unconditional, power for a variety of study designs used in bioequivalence studies.
exppower.TOST(alpha = 0.05, logscale = TRUE, theta0, theta1, theta2, CV, n, design = "2x2", robust = FALSE, prior.type = c("CV", "theta0", "both"), prior.parm = list(), method = c("exact", "approx"))exppower.TOST(alpha = 0.05, logscale = TRUE, theta0, theta1, theta2, CV, n, design = "2x2", robust = FALSE, prior.type = c("CV", "theta0", "both"), prior.parm = list(), method = c("exact", "approx"))
alpha |
Significance level (one-sided). Commonly set to 0.05. |
logscale |
Should the data be used on log-transformed or on original scale? |
theta0 |
Assumed ‘true’ (or ‘observed’ in case of |
theta1 |
Lower (bio-)equivalence limit. |
theta2 |
Upper (bio-)equivalence limit. |
CV |
In case of In case of cross-over studies this is the within-subject CV, in case of a parallel-group design the CV of the total variability. |
n |
Number of subjects under study. |
design |
Character string describing the study design. See known.designs for designs covered in this package. |
robust |
Defaults to |
prior.type |
Specifies which parameter uncertainty should be accounted for. In case of
|
prior.parm |
A list of parameters expressing the prior information about the
variability and/or treatment effect. Possible components are |
method |
Defaults to |
This function calculates the so-called expected power taking into account that usually the parameters (CV and/or theta0) are not known but estimated from a prior study with some uncertainty. The expected power is an unconditional power and can therefore be seen as probability for success. See references for further details.
The prior.parm argument is a list that can supply any of the following
components:
dfError degrees of freedom from the prior trial (>4, maybe non-integer).
df = Inf is allowed and for method = "exact" the result will then
coincide with power.TOST(...).
Note: This corresponds to the df of both the CV and the difference of means.
SEMStandard error of the difference of means from the prior trial; must always be on additive scale (i.e., usually log-scale).
mNumber of subjects from prior trial. Specification is analogous to
the main argument n.
designStudy design of prior trial. Specification is analogous to the
main argument design.
For prior.parm, the combination consisting of df and SEM
requires a somewhat advanced knowledge of the prior trial (provided in the raw
output from for example the software SAS, or may be obtained via
emmeans of package emmeans.
However, it has the advantage that if there were
missing data the exact degrees of freedom and standard error of the difference
can be used, the former possibly being non-integer valued (e.g., if the
Kenward-Roger method was used).
Details on argument prior.type:
CVThe expectation is calculated with respect to the Inverse-gamma distribution.
theta0The expectation is calculated with respect to the
conditional distribution theta0 | = s^2
of the posteriori distribution of (theta0, ) from the prior
trial.
bothThe expectation is calculated with respect to the posteriori
distribution of (theta0, ) from the prior trial. Numerical calculation
of the two-dimensional integral is performed via hcubature.
Value of expected power according to the input.
B. Lang (thanks to G. Nehmiz for the helpful discussions), D. Labes
Grieve AP. Confidence Intervals and Sample Sizes. Biometrics. 1991;47:1597–603. doi:10.2307/2532411
O’Hagan, Stevens, JW, Campell MJ. Assurance in Clinical Trial Design. Pharm Stat. 2005;4:187–201. doi:10.1002/pst.175
Julious SA, Owen RJ. Sample size calculations for clinical studies allowing for uncertainty in variance. Pharm Stat. 2006;5:29–37. doi:10.1002/pst.197
Julious SA. Sample sizes for Clinical Trials. Boca Raton: CRC Press / Chapman & Hall; 2010.
Bertsche A, Nehmitz G, Beyersmann J, Grieve AP. The predictive distribution of the residual variability in the linear-fixed effects model for clinical cross-over trials. Biom J. 2016;58(4):797–809. doi:10.1002/bimj.201500245
Box GEP, Tiao GC. Bayesian Inference in Statistical Analysis. Boston: Addison-Wesley; 1992.
Held L, Sabanes Bove D. Applied Statistical Inference. Likelihood and Bayes. Berlin, Heidelberg: Springer; 2014. doi:10.1007/978-3-642-37887-4
Senn S. Cross-over Trials in Clinical Research. Chichester: John Wiley & Sons; 2nd edition 2002.
Zierhut ML, Bycott P, Gibbs MA, Smith BP, Vicini P. Ignorance is not bliss: Statistical power is not probability of trial success. Clin Pharmacol Ther. 2015;99:356–9. doi:10.1002/cpt.257
# Expected power for a 2x2 crossover with 40 subjects # CV 30% known from a pilot 2x2 study with 12 subjects # using all the defaults for other parameters (theta0 carved in stone) exppower.TOST(CV = 0.3, n = 40, prior.parm = list(df = 12-2)) # should give: [1] 0.7365519 # or equivalently exppower.TOST(CV = 0.3, n = 40, prior.parm = list(m = 12, design = "2x2")) # In contrast: Julious approximation exppower.TOST(CV = 0.3, n = 40, prior.parm = list(df = 10), method = "approx") # should give: [1] 0.7359771 # Compare this to the usual (conditional) power (CV known, "carved in stone") power.TOST(CV = 0.3, n = 40) # should give: [1] 0.8158453 # same as if setting df = Inf in function exppower.TOST() exppower.TOST(CV = 0.3, n = 40, prior.parm = list(df = Inf)) # Expected power for a 2x2 crossover with 40 subjects # CV 30% and theta0 = 0.95 known from a pilot 2x2 study with 12 subjects # using uncertainty with respect to both CV and theta0 exppower.TOST(CV = 0.3, theta0 = 0.95, n = 40, prior.parm = list(m = 12, design = "2x2"), prior.type = "both") # should give [1] 0.5114685# Expected power for a 2x2 crossover with 40 subjects # CV 30% known from a pilot 2x2 study with 12 subjects # using all the defaults for other parameters (theta0 carved in stone) exppower.TOST(CV = 0.3, n = 40, prior.parm = list(df = 12-2)) # should give: [1] 0.7365519 # or equivalently exppower.TOST(CV = 0.3, n = 40, prior.parm = list(m = 12, design = "2x2")) # In contrast: Julious approximation exppower.TOST(CV = 0.3, n = 40, prior.parm = list(df = 10), method = "approx") # should give: [1] 0.7359771 # Compare this to the usual (conditional) power (CV known, "carved in stone") power.TOST(CV = 0.3, n = 40) # should give: [1] 0.8158453 # same as if setting df = Inf in function exppower.TOST() exppower.TOST(CV = 0.3, n = 40, prior.parm = list(df = Inf)) # Expected power for a 2x2 crossover with 40 subjects # CV 30% and theta0 = 0.95 known from a pilot 2x2 study with 12 subjects # using uncertainty with respect to both CV and theta0 exppower.TOST(CV = 0.3, theta0 = 0.95, n = 40, prior.parm = list(m = 12, design = "2x2"), prior.type = "both") # should give [1] 0.5114685
Estimates the sample size based on the expected power for a variety of designs used in bioequivalence studies. See known.designs for the study designs covered.
expsampleN.noninf(alpha = 0.025, targetpower = 0.8, logscale = TRUE, theta0, margin, CV, design = "2x2", robust = FALSE, prior.type = c("CV", "theta0", "both"), prior.parm = list(), method = c("exact", "approx"), print = TRUE, details)expsampleN.noninf(alpha = 0.025, targetpower = 0.8, logscale = TRUE, theta0, margin, CV, design = "2x2", robust = FALSE, prior.type = c("CV", "theta0", "both"), prior.parm = list(), method = c("exact", "approx"), print = TRUE, details)
alpha |
Significance level (one-sided). Defaults here to 0.025. |
targetpower |
Power to achieve at least. Must be |
logscale |
Should the data used on log-transformed or on original scale?
|
theta0 |
Assumed ‘true’ (or ‘observed’ in case of |
margin |
Non-inferiority margin. |
CV |
In case of If In case of cross-over studies this is the within-subject CV, in case of a parallel-group design the CV of the total variability. |
design |
Character string describing the study design. |
robust |
Defaults to |
prior.type |
Specifies which parameter uncertainty should be accounted for.
In case of |
prior.parm |
A list of parameters expressing the prior information about the
variability and/or treatment effect. Possible components are |
method |
Defaults to |
print |
If |
details |
If |
The sample size is estimated based on iterative evaluation of
expected power. The starting value of the sample size search is
taken from a large sample approximation if prior.type="CV".
Else an empirical start value is obtained. Note that in case of
prior.type="both" the calculation may still take several seconds.
Note also that the expected power is always bounded above by the
so-called probability of technical success (PTS) which
may be a value less than 1.Therefore, it may be possible that it
is either not possible to calculate the required sample size at
all or that the sample size gets very large if the given targetpower
is less but close to the PTS.
Notes on the underlying hypotheses
If the supplied margin is < 0 (logscale=FALSE) or
< 1 (logscale=TRUE), then it is assumed higher response
values are better. The hypotheses are H0: theta0 <= margin H1: theta0 > margin
where theta0 = mean(test)-mean(reference) if logscale=FALSE
or H0: log(theta0) <= log(margin) H1: log(theta0) > log(margin)
where theta0 = mean(test)/mean(reference) if logscale=TRUE.
If the supplied margin is > 0 (logscale=FALSE) or
> 1 (logscale=TRUE), then it is assumed lower response
values are better. The hypotheses are H0: theta0 >= margin H1: theta0 < margin
where theta0 = mean(test)-mean(reference) if logscale=FALSE
or H0: log(theta0) >= log(margin) H1: log(theta0) < log(margin)
where theta0 = mean(test)/mean(reference) if logscale=TRUE.
This latter case may also be considered as ‘non-superiority’.
A data.frame with the input values and the result of the sample
size estimation.
The Sample size column contains the total sample
size in case of all designs implemented.
B. Lang, D. Labes
Grieve AP. Confidence Intervals and Sample Sizes. Biometrics. 1991;47:1597–603. doi:10.2307/2532411
O’Hagan, Stevens, JW, Campell MJ. Assurance in Clinical Trial Design. Pharm Stat. 2005;4:187–201. doi:10.1002/pst.175
Julious SA, Owen RJ. Sample size calculations for clinical studies allowing for uncertainty in variance. Pharm Stat. 2006;5:29–37. doi:10.1002/pst.197
Julious SA. Sample sizes for Clinical Trials. Boca Raton: CRC Press; 2010.
Bertsche A, Nehmitz G, Beyersmann J, Grieve AP. The predictive distribution of the residual variability in the linear-fixed effects model for clinical cross-over trials. Biom J. 2016;58(4):797–809. doi:10.1002/bimj.201500245
Box GEP, Tiao GC. Bayesian Inference in Statistical Analysis. Boston: Addison-Wesley; 1992.
Held L, Sabanes Bove D. Applied Statistical Inference. Likelihood and Bayes. Berlin, Heidelberg: Springer; 2014. doi:10.1007/978-3-642-37887-4
Senn S. Cross-over Trials in Clinical Research. Chichester: John Wiley & Sons; 2nd edition 2002.
Zierhut ML, Bycott P, Gibbs MA, Smith BP, Vicini P. Ignorance is not bliss: Statistical power is not probability of trial success. Clin Pharmacol Ther. 2015;99:356–9. doi:10.1002/cpt.257
exppower.noninf, known.designs, sampleN.noninf
# Classical 2x2 cross-over, target power = 80%, # assumed true ratio = 95%, margin = 0.8, # intra-subject CV=30% estimated from prior 2x2 trial # with m = 12 subjects expsampleN.noninf(theta0 = 0.95, margin = 0.8, CV = 0.3, design = "2x2", prior.parm = list(m = 12, design = "2x2")) # gives n = 58 with achieved expected power 0.809148 # Compare this to the usual sample size with CV assumed # as 'carved in stone' sampleN.noninf(theta0 = 0.95, margin = 0.8, CV = 0.3) # Perform 'non-superiority' (lower is better) with assumed # true ratio = 105% and margin 125% expsampleN.noninf(theta0 = 1.05, margin = 1.25, CV = 0.3, design = "2x2", prior.parm = list(m = 12, design = "2x2")) # should give n = 56 with achieved expected power 0.806862 # More than one CV with corresponding degrees of freedom # other settings as above in first example CVs <- c(0.25, 0.3) dfs <- c(22, 10) expsampleN.noninf(theta0 = 0.95, margin = 0.8, CV = CVs, prior.parm = list(df = dfs)) # should give a pooled CV=0.2664927 with 32 df and a sample # size n=42 with achieved expected power 0.814073 exact # achieved expected power 0.816163 approximate acc. to Julious # Uncertainty is accounted for CV and theta0 expsampleN.noninf(CV = 0.3, prior.type = "both", prior.parm = list(m = 12, design = "2x2")) # gives a dramatic increase in sample size (n = 194) # due to small pilot trial# Classical 2x2 cross-over, target power = 80%, # assumed true ratio = 95%, margin = 0.8, # intra-subject CV=30% estimated from prior 2x2 trial # with m = 12 subjects expsampleN.noninf(theta0 = 0.95, margin = 0.8, CV = 0.3, design = "2x2", prior.parm = list(m = 12, design = "2x2")) # gives n = 58 with achieved expected power 0.809148 # Compare this to the usual sample size with CV assumed # as 'carved in stone' sampleN.noninf(theta0 = 0.95, margin = 0.8, CV = 0.3) # Perform 'non-superiority' (lower is better) with assumed # true ratio = 105% and margin 125% expsampleN.noninf(theta0 = 1.05, margin = 1.25, CV = 0.3, design = "2x2", prior.parm = list(m = 12, design = "2x2")) # should give n = 56 with achieved expected power 0.806862 # More than one CV with corresponding degrees of freedom # other settings as above in first example CVs <- c(0.25, 0.3) dfs <- c(22, 10) expsampleN.noninf(theta0 = 0.95, margin = 0.8, CV = CVs, prior.parm = list(df = dfs)) # should give a pooled CV=0.2664927 with 32 df and a sample # size n=42 with achieved expected power 0.814073 exact # achieved expected power 0.816163 approximate acc. to Julious # Uncertainty is accounted for CV and theta0 expsampleN.noninf(CV = 0.3, prior.type = "both", prior.parm = list(m = 12, design = "2x2")) # gives a dramatic increase in sample size (n = 194) # due to small pilot trial
Estimates the sample size based on the expected power for a variety of study designs used in bioequivalence studies. See known.designs for the study designs covered.
expsampleN.TOST(alpha = 0.05, targetpower = 0.8, logscale=TRUE, theta0, theta1, theta2, CV, design = "2x2", robust = FALSE, prior.type = c("CV", "theta0", "both"), prior.parm = list(), method = c("exact", "approx"), print = TRUE, details)expsampleN.TOST(alpha = 0.05, targetpower = 0.8, logscale=TRUE, theta0, theta1, theta2, CV, design = "2x2", robust = FALSE, prior.type = c("CV", "theta0", "both"), prior.parm = list(), method = c("exact", "approx"), print = TRUE, details)
alpha |
Significance level (one-sided). Commonly set to 0.05. |
targetpower |
Power to achieve at least. Must be >0 and <1. Typical values are 0.8 or 0.9. |
logscale |
Should the data used on log-transformed or on original scale? |
theta0 |
Assumed ‘true’ (or ‘observed’ in case of |
theta1 |
Lower (bio-)equivalence limit. |
theta2 |
Upper (bio-)equivalence limit. |
CV |
In case of If In case of cross-over studies this is the within-subject CV, in case of a parallel-group design the CV of the total variability. |
design |
Character string describing the study design. |
robust |
Defaults to FALSE. With that value the usual degrees of freedom will be used. |
prior.type |
Specifies which parameter uncertainty should be accounted for. In case of
|
prior.parm |
A list of parameters expressing the prior information about the
variability and/or treatment effect. Possible components are |
method |
Defaults to |
print |
If |
details |
If |
The sample size is calculated based on iterative evaluation of expected power.
The starting value of the sample size search is taken from a large sample
approximation if prior.type = "CV". Otherwise, an empirical start value is
obtained. Note that in case of prior.type = "both" the calculation may
still take several seconds.
Note also that the expected power is always bounded above by the so-called
probability of technical success (PTS) which may be a value less than 1.
Therefore, it may be possible that it is either not possible to calculate the
required sample size at all or that the sample size gets very large
if the given targetpower is less but close to the PTS.
The estimated sample size gives always the total number of subjects (not subject/sequence in crossovers or subjects/group in parallel designs – like in some other software packages).
A data.frame with the input values and the result of the sample size estimation.
The Sample size column contains the total sample size in case of all
designs implemented.
B. Lang, D. Labes
Grieve AP. Confidence Intervals and Sample Sizes. Biometrics. 1991;47:1597–603. doi:10.2307/2532411
O’Hagan, Stevens, JW, Campell MJ. Assurance in Clinical Trial Design. Pharm Stat. 2005;4:187–201. doi:10.1002/pst.175
Julious SA, Owen RJ. Sample size calculations for clinical studies allowing for uncertainty in variance. Pharm Stat. 2006;5:29–37. doi:10.1002/pst.197
Julious SA. Sample sizes for Clinical Trials. Boca Raton: CRC Press / Chapman & Hall; 2010.
Bertsche A, Nehmitz G, Beyersmann J, Grieve AP. The predictive distribution of the residual variability in the linear-fixed effects model for clinical cross-over trials. Biom J. 2016;58(4):797–809. doi:10.1002/bimj.201500245
Box GEP, Tiao GC. Bayesian Inference in Statistical Analysis. Boston: Addison-Wesley; 1992.
Held L, Sabanes Bove D. Applied Statistical Inference. Likelihood and Bayes. Berlin, Heidelberg: Springer; 2014. doi:10.1007/978-3-642-37887-4
Senn S. Cross-over Trials in Clinical Research. Chichester: John Wiley & Sons; 2nd edition 2002.
Zierhut ML, Bycott P, Gibbs MA, Smith BP, Vicini P. Ignorance is not bliss: Statistical power is not probability of trial success. Clin Pharmacol Ther. 2015;99:356–9. doi:10.1002/cpt.257
exppower.TOST, known.designs, sampleN.TOST
# Classical 2x2 cross-over, target power = 80%, # BE limits 80 ... 125%, assumed true BE ratio = 95%, # intra-subject CV=30% estimated from prior 2x2 trial # with m = 30 subjects expsampleN.TOST(CV=0.3, prior.parm = list(m = 30, design = "2x2")) # -> gives n = 42 with achieved expected power 0.806262 # Compare this to the usual sample size with CV assumed known ('carved in stone') sampleN.TOST(CV=0.3) # -> gives n = 40 subjects # Compare this to the case where uncertainty is accounted for CV and theta0 # Not run due to timing policy of CRAN - may run several seconds expsampleN.TOST(CV=0.3, prior.parm = list(m = 30, design = "2x2"), prior.type = "both") # -> gives n = 72 subjects # More than one CV with corresponding degrees of freedom # other settings as above in first example CVs <- c(0.25, 0.3) dfs <- c(22, 10) expsampleN.TOST(CV=CVs, prior.parm = list(df = dfs)) # -> gives a pooled CV=0.2664927 with df=32 # and a sample size n=34 with achieved expected power 0.812653 exact # achieved expected power 0.815019 approximate acc. Julious# Classical 2x2 cross-over, target power = 80%, # BE limits 80 ... 125%, assumed true BE ratio = 95%, # intra-subject CV=30% estimated from prior 2x2 trial # with m = 30 subjects expsampleN.TOST(CV=0.3, prior.parm = list(m = 30, design = "2x2")) # -> gives n = 42 with achieved expected power 0.806262 # Compare this to the usual sample size with CV assumed known ('carved in stone') sampleN.TOST(CV=0.3) # -> gives n = 40 subjects # Compare this to the case where uncertainty is accounted for CV and theta0 # Not run due to timing policy of CRAN - may run several seconds expsampleN.TOST(CV=0.3, prior.parm = list(m = 30, design = "2x2"), prior.type = "both") # -> gives n = 72 subjects # More than one CV with corresponding degrees of freedom # other settings as above in first example CVs <- c(0.25, 0.3) dfs <- c(22, 10) expsampleN.TOST(CV=CVs, prior.parm = list(df = dfs)) # -> gives a pooled CV=0.2664927 with df=32 # and a sample size n=34 with achieved expected power 0.812653 exact # achieved expected power 0.815019 approximate acc. Julious
Returns the known study designs for which power and sample size can be calculated within this package.
known.designs()known.designs()
This function is for informal purposes and used internally for obtaining characteristics of the designs used in calculation formulas.
Returns a data.frame with
no |
number of the design |
design |
character string for identifying the design |
df |
degrees of freedom of the design |
df2 |
‘robust’ degrees of freedom of the design |
steps |
step width in the iterative sample size estimation |
bk |
so-called design constant in terms of total n |
bkni |
design constant in terms of number of subjects in (sequence) groups |
The design character string has to be used in the functions calls for power and sample size.
The design string for higher order crossover designs is named as:treatments x sequences x periods in case of replicate designs andtreatments x periods in case of crossover designs for more then 2 treatments
with number of sequences equal number of treatments.
The df for the replicate crossover designs are those without carry-over in the model.
Chen et al. used models with carry-over, i.e., one df lower than here.
The design constant bk in case of design 2x2x4 is here bk=1.
Chen et al. used bk=1.1 due to carry-over in the model.
n is the total number of subjects for all designs implemented.
df2 = degrees of freedom for the so-called ‘robust’ analysis (aka Senn’s basic estimator).
These degrees of freedom are often also more appropriate in case of evaluation via a ‘true’
mixed model (e.g. the FDA’ for replicate designs).
The design 2x2x2r is the 2-treatment-2-sequence-2-period design with
2 repeated targets determined in each period (sequences TT|RR or RR|TT) described
by Liu. Implemented are the characteristics of this design for the evaluation
via assuming no S×F interaction and equal variabilities of Test and Reference.
D. Labes
Chen KW, Chow SC, Liu G. A Note on Sample Size Determination for Bioequivalence Studies with Higher-order Crossover Designs. J Pharmacokin Biopharm. 1997;25(6):753–65.
Senn S. Cross-over Trials in Clinical Research. Chichester: John Wiley & Sons; 2nd edition 2002.
U.S. Department of Health and Human Services, Food and Drug Administration, Center for Drug Evaluation and Research (CDER). Guidance for Industry. Statistical Approaches to Establishing Bioequivalence. January 2001. download
Liu J-p. Use of the Repeated Crossover design in Assessing Bioequivalence. Stat Med. 1995;14(9-10):1067–78.
known.designs()known.designs()
Calculates Owen’s Q function.
OwensQ(nu, t, delta, a=0, b)OwensQ(nu, t, delta, a=0, b)
nu |
degree of Owen’s Q |
t |
parameter t |
delta |
parameter delta |
a |
lower integration limit, only a=0 implemented |
b |
upper integration limit |
Uses the relationship to non-central t-distribution (see Chou) OwensQ = pt(t, df=nu, ncp=delta) - Integal_b_Inf(Q_integrand)
The definite integral is numerically evaluated using integrate
after a variables transformation resulting in the integration range from 0 to 1
instead of the semi-infinite original range.
This may result in higher precision and better numerical stability.
The arguments to the function must be scalars. No vectors allowed.
Numeric value of Owen’s Q-function at given input arguments.
This function is intended for internal use in the power calculations.
But may be useful for others.
D. Labes
Owen DB. A special case of a bivariate non-central t-distribution. Biometrika. 1965;52(3/4):437–46. doi:10.2307/2333696
Chou YM. A bivariate noncentral T-distibution with applications. Commun Stat Theory Methods. 1992;21(12):3427–62. doi:10.1080/03610929208830988
# This function is mainly intended for internal use. OwensQ(10, 2.5, 5, 0, 2) #should give [1] 9.388137e-06 OwensQ(10, -2.5, -5, 0, 2) #should give [1] 0.05264363# This function is mainly intended for internal use. OwensQ(10, 2.5, 5, 0, 2) #should give [1] 9.388137e-06 OwensQ(10, -2.5, -5, 0, 2) #should give [1] 0.05264363
This is an implementation of the algorithm given by Owen via repeated integration by parts.
OwensQOwen(nu, t, delta, a=0, b)OwensQOwen(nu, t, delta, a=0, b)
nu |
degree of Owen’s Q |
t |
parameter t |
delta |
parameter delta |
a |
lower integration limit. |
b |
upper integration limit |
Numeric value of Owen’s Q function.
The argument a=0 could be dropped but is retained for sake of completeness.
This function is mainly for comparative / validation purposes.
The function requireds OwensT function.
D. Labes
Owen DB. A special case of a bivariate non-central t-distribution. Biometrika. 1965;52(3/4):437–46. doi:10.2307/2333696
# comparison of the results of both implementations # both should give [1] 0.0731726 OwensQ(2, 2.92, 4.2135, 0, 2.0407) OwensQOwen(2, 2.92, 4.2135, 0, 2.0407)# comparison of the results of both implementations # both should give [1] 0.0731726 OwensQ(2, 2.92, 4.2135, 0, 2.0407) OwensQOwen(2, 2.92, 4.2135, 0, 2.0407)
Calculates the definite integral from 0 to a of exp(-0.5*h^2*(1+x^2))/(1+x^2)/(2*pi).
OwensT(h, a)OwensT(h, a)
h |
parameter h |
a |
upper limit of integration |
The function is an R port of FORTRAN code given in the references and MATLAB
code given by John Burkardt under the GNU LGPL license.
The arguments of OwensT() have to be scalars because the implementation
doesn’t vectorize.
Numerical value of the definite integral.
This function is only needed as auxiliary in OwensQOwen.
But may be useful for others.
MATLAB code by J. Burkardt, R port by D. Labes
Goedhart PW, Jansen MJW. Remark AS R89: A Remark on Algorithm AS 76: An Integral Useful in Calculating Central t and Bivariate Normal Probabilities. J Royal Stat Soc C. 1992;41(2):496–7. doi:10.2307/2347586
Boys R. Algorithm AS R80: A Remark on Algorithm AS 76: An Integral Useful in Calculating Noncentral t and Bivariate Normal Probabilities. J Royal Stat Soc C. 1989;38(3):580–2. doi:10.2307/2347755
Thomas GE. Remark ASR 65: A Remark on Algorithm AS76: An Integral Useful in Calculating Non-Central t and Bivariate Normal Probabilities. J Royal Stat Soc C. 1986;35(3):310–2. doi:10.2307/2348031
Chou Y-M. Remark AS R55: A Remark on Algorithm AS 76: An Integral Useful in Calculating Noncentral T and Bivariate Normal Probabilities. J Royal Stat Soc C. 1985;34(1):100–1. doi:10.2307/2347894
Thomas GE. Remark AS R30: A Remark on Algorithm AS 76: An Integral Useful in Calculating Non-Central t and Bivariate Normal Probabilities. J Royal Stat Soc C. 1979;28(1):113. doi:10.2307/2346833
Young JC, Minder C. Algorithm AS 76: An Integral Useful in Calculating Non-Central t and Bivariate Normal Probabilities. J Royal Stat Soc C. 1974;23(3):455–7. doi:10.2307/2347148
Burkardt J. ASA076. Owen's T Function. https://people.math.sc.edu/Burkardt/f_src/asa076/asa076.html
Owen DB. Tables for Computing Bivariate Normal Probabilities. Ann Math Stat. 1956;27(4):1075–90. doi:10.1214/aoms/1177728074
OwensT(2.5, 0.75) # should give [1] 0.002986697 # value from Owen's tables is 0.002987 OwensT(2.5, -0.75) # should give [1] -0.002986697OwensT(2.5, 0.75) # should give [1] 0.002986697 # value from Owen's tables is 0.002987 OwensT(2.5, -0.75) # should give [1] -0.002986697
An analysis tool for exploration/visualization of the impact of expected values (CV, theta0, reduced sample size due to drop-outs) on power of BE decision via ABE if these values deviate from the ones assumed in planning the sample size of the study.
pa.ABE(CV, theta0 = 0.95, targetpower = 0.8, minpower = 0.7, design = "2x2", ...) ## S3 method for class 'pwrA' print(x, digits = 4, plotit = TRUE, ...) ## S3 method for class 'pwrA' plot(x, pct = TRUE, ratiolabel = "theta0", cols = c("blue", "red"), ...)pa.ABE(CV, theta0 = 0.95, targetpower = 0.8, minpower = 0.7, design = "2x2", ...) ## S3 method for class 'pwrA' print(x, digits = 4, plotit = TRUE, ...) ## S3 method for class 'pwrA' plot(x, pct = TRUE, ratiolabel = "theta0", cols = c("blue", "red"), ...)
CV |
Coefficient of variation as ratio (not percent). |
theta0 |
‘True’ or assumed T/R ratio. Often named GMR. |
targetpower |
Power to achieve at least in sample size estimation. Must be >0 and <1. |
minpower |
Minimum acceptable power to have if deviating from assumptions for sample size plan. |
design |
Character string describing the study design. |
... |
More arguments to pass to |
Additional arguments of the S3 methods:
x |
Object of class |
digits |
Digits for rounding power in printing. The '...' argument is currently ignored
in |
plotit |
If set to |
pct |
If set to |
ratiolabel |
Label of the T/R-ratio. Can be set to any string, e.g. to |
cols |
Colors for the plots. |
Power calculations are done via power.TOST() and calculations of CV and theta0
which gave a power=minpower are derived via R base uniroot.
While one of the parameters (CV, theta0, n) is varied, the respective two others are
kept constant. The tool shows the relative impact of single parameters on power.
The tool takes a minimum of 12 subjects as required in most BE guidances into account.
It should be kept in mind that this is not a substitute for the ‘Sensitivity Analysis’
recommended in ICH-E9. In a real study a combination of all effects occurs simultaneously.
It is up to you to decide on reasonable combinations and analyze their respective power.
Returns a list with class "pwrA" with the components
plan |
A data.frame with the result of the sample size estimation.
See output of |
paCV |
A data.frame with value pairs CV, pwr for impact of deviations from CV. |
paGMR |
A data.frame with value pairs theta0, pwr for impact of deviations from theta0 (GMR). |
paN |
A data.frame with value pairs N, pwr for impact of deviations from planned N (dropouts). |
method |
Method of BE decision. Here "ABE". |
minpower |
Minimum acceptable power. |
The class 'pwrA' has the S3 methods print() and plot().
See pa.scABE for usage.
The code of deviations from planned sample size tries to keep the degree of imbalance as low as possible between (sequence) groups. This results in a lesser decrease of power than more extreme dropout-patterns.
Idea and original code by H. Schütz with modifications by D. Labes to use PowerTOST infrastructure.
Schütz H. Deviating from assumptions. August 08, 2014. BEBA Forum
power.TOST, known.designs, pa.scABE, pa.NTIDFDA
# using the defaults # design="2x2", targetpower=0.8, minpower=0.7, theta0/GMR=0.95 # BE margins from defaults of sampleN.TOST() 0.8 ... 1.25 # print & plot implicitly pa.ABE(CV = 0.2) # print & plot res <- pa.ABE(CV = 0.2) print(res, plotit = FALSE) # print only plot(res, pct = FALSE, ratiolabel = "GMR") # changed from defaults# using the defaults # design="2x2", targetpower=0.8, minpower=0.7, theta0/GMR=0.95 # BE margins from defaults of sampleN.TOST() 0.8 ... 1.25 # print & plot implicitly pa.ABE(CV = 0.2) # print & plot res <- pa.ABE(CV = 0.2) print(res, plotit = FALSE) # print only plot(res, pct = FALSE, ratiolabel = "GMR") # changed from defaults
An analysis tool for exploration/visualization of the impact of expected values
(CV, theta0, reduced sample size due to drop-outs) on power of BE decision via
scABE for narrow therapeutic drugs (NTIDs) if these values deviate from the ones
assumed in planning the sample size of the study.
The only implemented design is the full replicate design "2x2x4" according to the
FDA’s warfarin guidance.
pa.NTID(CV, theta0 = 0.975, targetpower = 0.8, minpower = 0.7, ...)pa.NTID(CV, theta0 = 0.975, targetpower = 0.8, minpower = 0.7, ...)
CV |
Coefficient of variation of the intra-subject variabilities of Test and Reference
as ratio (not percent). |
theta0 |
‘True’ or assumed T/R ratio. Often named GMR. |
targetpower |
Power to achieve at least in sample size estimation. Must be >0 and <1. |
minpower |
Minimum acceptable power to have if deviating from assumptions for sample size plan. |
... |
More arguments to pass to |
Power calculations are done via power.NTID() and
calculations of CV and theta0 which result in minpower are obtained via uniroot().
While one of the parameters (CV, theta0, n) is varied, the respective two others are
kept constant. The tool shows the relative impact of single parameters on power.
The tool takes a minimum of 12 subjects into account as demanded in most BE guidances.
However, it should be kept in mind that the FDA requires at least 24 subjects to be enrolled
in studies intended for reference-scaling.
It should be kept in mind that this is not a substitute for the ‘Sensitivity Analysis’
recommended in ICH-E9. In a real study a combination of all effects occurs simultaneously.
It is up to you to decide on reasonable combinations and analyze their respective power.
Returns a list with class 'pwrA' with the components
plan |
A data.frame with the result of the sample size estimation.
See output of |
paCV |
A data.frame with value pairs CV, pwr for impact of deviations from CV. |
paGMR |
A data.frame with value pairs theta0, pwr for impact of deviations from theta0 (GMR). |
paN |
A data.frame with value pairs N, pwr for impact of deviations from planned N (dropouts). |
method |
Method of BE decision. Here "NTID". |
regulator |
Here "FDA". |
minpower |
Minimum acceptable power from the call of the function. |
The class 'pwrA' has the S3 methods print() and plot().
See pa.ABE for usage.
Be extremly carefull if your sample size plan has extremly small CV near or
below 0.05 (5%). Adapt in that case your expected true ratio (theta0)
to values nearer to 1 to not run into errors and/or long execution times.
The code for impact of deviations from planned sample size tries to keep the degree of imbalance as low as possible between (sequence) groups. This results in a lesser decrease of power than more extreme dropout-patterns.
D. Labes according to code by H. Schütz for pa.ABE() and pa.scABE().
Food and Drug Administration, Office of Generic Drugs (OGD). Draft Guidance on Warfarin Sodium. Recommended Dec 2012. download
Food and Drug Administration, Center for Drug Evaluation and Research (CDER). Draft Guidance for Industry. Bioequivalence Studies with Pharmacokinetic Endpoints for Drugs Submitted Under an ANDA. August 2021. download
Yu LX, Jiang W, Zhang X, Lionberger R, Makhlouf F, Schuirmann DJ, Muldowney L, Chen ML, Davit B, Conner D, Woodcock J. Novel bioequivalence approach for narrow therapeutic index drugs. Clin Pharmacol Ther. 2015;97(3):286–91. doi:10.1002/cpt.28
Jiang W, Makhlouf F, Schuirmann DJ, Zhang X, Zheng N, Conner D, Yu LX, Lionberger R. A Bioequivalence Approach for Generic Narrow Therapeutic Index Drugs: Evaluation of the Reference-Scaled Approach and Variability Comparison Criterion. AAPS J. 2015;17(4):891–901. doi:10.1208/s12248-015-9753-5
Endrényi L, Tóthfalusi L. Determination of Bioequivalence for Drugs with Narrow Therapeutic Index: Reduction of the Regulatory Burden. J Pharm Pharm Sci. 2013;16(5):676–82. open access
power.NTID, print.pwrA, plot.pwrA, pa.ABE, pa.scABE
# using the defaults: # targetpower=0.8, minpower=0.7, theta0/GMR=0.975 # BE margins from defaults of sampleN.NTID() 0.9002 ... 1.1108 # 1E5 sims in power.NTID() # not run due to timing policy of CRAN for examples # may run some ten seconds or more plot(pa.NTID(CV=0.1))# using the defaults: # targetpower=0.8, minpower=0.7, theta0/GMR=0.975 # BE margins from defaults of sampleN.NTID() 0.9002 ... 1.1108 # 1E5 sims in power.NTID() # not run due to timing policy of CRAN for examples # may run some ten seconds or more plot(pa.NTID(CV=0.1))
An analysis tool for exploration/visualization of the impact of expected values (CV, theta0, reduced sample size due to drop-outs) on power of BE decision via scABE (for highly variable drugs) if these values deviate from the ones assumed in planning the sample size of the study.
pa.scABE(CV, theta0 = 0.9, targetpower = 0.8, minpower = 0.7, design = c("2x3x3", "2x2x4", "2x2x3"), regulator = c("EMA", "HC", "FDA", "GCC"), ...)pa.scABE(CV, theta0 = 0.9, targetpower = 0.8, minpower = 0.7, design = c("2x3x3", "2x2x4", "2x2x3"), regulator = c("EMA", "HC", "FDA", "GCC"), ...)
CV |
Coefficient of variation of the intra-subject variability as ratio (not percent). |
theta0 |
‘True’ or assumed T/R ratio. Often named GMR. |
targetpower |
Power to achieve at least in sample size estimation. Must be >0 and <1. |
minpower |
Minimum acceptable power to have if deviating from assumptions for sample size plan. |
design |
Character string describing the study design. |
regulator |
Character string describing the scaled ABE method recommended by the regulatory
bodies |
... |
More arguments to pass to |
Power calculations are done via power.scABEL() or power.RSABE() and
calculations of CV and theta0 which result in minpower derived via R base uniroot.
While one of the parameters (CV, GMR, N) is varied, the respective two others are
kept constant. The tool shows the relative impact of single parameters on power.
The tool takes a minimum of 12 subjects as required in most BE guidances into account.
However, the sample size will be increased from the estimated one if one of the
following conditions is applicable:
The FDA requires at least 24 subjects enrolled in studies intended for reference-scaling.
The EMA requires at least 12 eligible subjects in the sequence RTR of the TRT|RTR-design (hence the minimum sample size is 24).
You should be aware that this is not a substitute for the “Sensitivity Analysis” recommended in ICH-E9. In a real study a combination of all effects occurs simultaneously. It is up to you to decide on reasonable combinations and analyze their respective power.
Returns a list with class 'pwrA' with the components
plan |
A data.frame with the result of the sample size estimation. |
paCV |
A data.frame with value pairs CV, pwr for impact of deviations from CV. |
paGMR |
A data.frame with value pairs theta0, pwr for impact of deviations from theta0 (GMR). |
paN |
A data.frame with value pairs N, pwr for impact of deviations from planned N (dropouts). |
method |
Method of BE decision. Here fix = "scABE". |
regulator |
"EMA", "HC", or "FDA". |
minpower |
Minimum acceptable power from the call of the function. |
The class 'pwrA' has the S3 methods print() and plot().
See pa.ABE for usage.
The code for impact of deviations from planned sample size tries to keep the degree of imbalance as low as possible between (sequence) groups. This results in a lesser decrease of power than more extreme dropout-patterns.
Idea and original code by H. Schütz with modifications by D. Labes to use PowerTOST infrastructure.
Schütz H. Deviating from assumptions. August 08, 2014. BEBA Forum
power.scABEL, power.RSABE, known.designs,
print.pwrA, plot.pwrA, pa.ABE, pa.NTIDFDA
# Implicitely using the defaults: # design = "2x3x3", targetpower = 0.8, minpower = 0.7, # theta0 = 0.9, GMR = 0.90, regulator = "EMA" # widened BE margins from defaults of sampleN.scABEL() 0.7462 ... 1.3402 # 1E5 sims in power.scABEL() # not run due to timing policy of CRAN, may run some ten seconds # Implicit print & plot pa.scABE(CV = 0.4)# Implicitely using the defaults: # design = "2x3x3", targetpower = 0.8, minpower = 0.7, # theta0 = 0.9, GMR = 0.90, regulator = "EMA" # widened BE margins from defaults of sampleN.scABEL() 0.7462 ... 1.3402 # 1E5 sims in power.scABEL() # not run due to timing policy of CRAN, may run some ten seconds # Implicit print & plot pa.scABE(CV = 0.4)
Calculates the exact power of two simultaneous TOST procedures (where the two parameters of the two TOSTs are correlated with some correlation) for various study designs used in BE studies
power.2TOST(alpha = c(0.05, 0.05), logscale = TRUE, theta0, theta1, theta2, CV, n, rho, design = "2x2", robust = FALSE, nsims, setseed = TRUE, details = FALSE)power.2TOST(alpha = c(0.05, 0.05), logscale = TRUE, theta0, theta1, theta2, CV, n, rho, design = "2x2", robust = FALSE, nsims, setseed = TRUE, details = FALSE)
alpha |
Vector; contains one-sided significance level for each of the two TOSTs. |
logscale |
Should the data used on log-transformed ( |
theta1 |
Vector; contains lower bioequivalence limit for each of the two TOSTs. |
theta2 |
Vector; contains upper bioequivalence limit for each of the two TOSTS. |
theta0 |
Vector; contains ‘true’ assumed bioequivalence ratio for each of the two TOSTs. |
CV |
Vector of coefficient of variations (given as as ratio, e.g., 0.2 for 20%). |
n |
Number of subjects under study. |
rho |
Correlation between the two PK metrics (e.g., AUC and Cmax) under consideration. This is defined as correlation between the estimator of the treatment difference of PK metric one and the estimator of the treatment difference of PK metric two. Has to be within {–1, +1}. |
design |
Character string describing the study design. |
robust |
Defaults to |
nsims |
Number of studies to simulate. Defaults to 1E5. |
setseed |
Logical; if |
details |
Logical; if |
Calculations are based on simulations and follow the distributional
properties as described in Phillips. This is in contrast to the calculations
via the 4-dimensional non-central t-distribution as described in Hua et al.
which was implemented in versions up to 1.4-6.
The formulas cover balanced and unbalanced studies w.r.t (sequence) groups.
In case of parallel group design and higher order crossover designs
(replicate crossover or crossover with more than two treatments) the calculations
are based on the assumption of equal variances for Test and Reference products
under consideration.
The formulas for the paired means 'design' do not take an additional correlation
parameter into account. They are solely based on the paired t-test
(TOST of differences = zero).
Value of power.
If n is given as scalar (total sample size) and this number is not
divisible by the number of (sequence) groups of the design an unbalanced design
with small imbalance is assumed. A corresponding message is thrown showing the
assumed numbers of subjects in (sequence) groups.
The function does not vectorize properly if design is a vector. Moreover,
theta0 and CV must be of length two, thus further vectorizing is not possible.
Other vector input is not tested yet.
B. Lang, D. Labes
Phillips KF. Power for Testing Multiple Instances of the Two One-Sided Tests Procedure. Int J Biostat. 2009;5(1):Article 15.
Hua SY, Xu S, D’Agostino RB Sr. Multiplicity adjustments in testing for bioequivalence. Stat Med. 2015;34(2):215–31. doi:10.1002/sim.6247
Lang B, Fleischer F. Letter to the Editor: Comments on ‘Multiplicity adjustments in testing for bioequivalence.’ Stat Med. 2016;35(14):2479–80. doi:10.1002/sim.6488
# Power for the 2x2x2 cross-over design with 24 subjects, intra-subject # standard deviation of 0.3 (CV = 30.7%) and assumed ratios of 1.05 for both # parameters, and correlation 0.75 between parameters (using all the other # default values) power.2TOST(theta0 = rep(1.05, 2), CV = rep(se2CV(0.3), 2), n = 24, rho = 0.75) # should give: 0.38849 # Setting as before but use rho 1 and high number of simulations # to reproduce result of power.TOST() p1 <- power.2TOST(theta0 = rep(1.05, 2), CV = rep(se2CV(0.3), 2), n = 24, rho = 1, nsims=1E7) p2 <- power.TOST(theta0 = 1.05, CV = se2CV(0.3), n = 24) all.equal(p1, p2, tolerance = 1e-04)# Power for the 2x2x2 cross-over design with 24 subjects, intra-subject # standard deviation of 0.3 (CV = 30.7%) and assumed ratios of 1.05 for both # parameters, and correlation 0.75 between parameters (using all the other # default values) power.2TOST(theta0 = rep(1.05, 2), CV = rep(se2CV(0.3), 2), n = 24, rho = 0.75) # should give: 0.38849 # Setting as before but use rho 1 and high number of simulations # to reproduce result of power.TOST() p1 <- power.2TOST(theta0 = rep(1.05, 2), CV = rep(se2CV(0.3), 2), n = 24, rho = 1, nsims=1E7) p2 <- power.TOST(theta0 = 1.05, CV = se2CV(0.3), n = 24) all.equal(p1, p2, tolerance = 1e-04)
Calculates the power of dose-proportionality studies using the power model for crossover (Latin square) or parallel group designs via a confidence interval equivalence criterion.
power.dp(alpha = 0.05, CV, doses, n, beta0, theta1 = 0.8, theta2 = 1/theta1, design = c("crossover", "parallel", "IBD"), dm = NULL, CVb)power.dp(alpha = 0.05, CV, doses, n, beta0, theta1 = 0.8, theta2 = 1/theta1, design = c("crossover", "parallel", "IBD"), dm = NULL, CVb)
alpha |
Type 1 error. Commonly set to 0.05. |
CV |
Coefficient of variation for intra-subject variability if |
doses |
Vector of dose levels. At least two doses have to be given. |
n |
Number of subjects. Is total number if given as scalar, else number of subjects
in the (sequence) groups. In the latter case the length of n vector has to be
the same as length of vector doses. |
beta0 |
‘True’ slope of power model. If missing defaults to |
theta1 |
Lower acceptance limit for the ratio of dose normalized means (Rdmn). |
theta2 |
Upper acceptance limit for the ratio of dose normalized means (Rdmn). |
design |
Crossover design (default), parallel group design or incomplete block design (IBD). |
dm |
'Design matrix' of the incomplete block design (IBD) if |
CVb |
Coefficient of variation of the between-subject variability. |
The power calculations are based on TOST for testing equivalence of the slope
of the power model with alternativ hypothesis slope = 1.
Power is calculated via non-central t-approximation only.
The calculations are based on mixed effects model (random intercept aka
random subject effect). For design="cossover" or design="parallel"
the results coincide with all-effects-fixed model.
Value of power according to the input arguments.
This function is ‘experimental’ only since it is not thorougly tested yet.
Especially for design="IBD" reliable test cases are missing.
D. Labes
Patterson S, Jones B. Bioequivalence and Statistics in Clinical Pharmacology. Boca Raton: Chapman & Hall/CRC: 2006. p. 239.
(contains presumably a bug)
Sethuraman VS, Leonov S, Squassante L, Mitchell TR, Hale MD. Sample size calculation for the Power Model for dose proportionality studies. Pharm Stat. 2007;6(1):35–41. doi:10.1002/pst.241
Hummel J, McKendrick S, Brindley C, French R. Exploratory assessment of dose proportionality: review of current approaches and proposal for a practical criterion. Pharm. Stat. 2009;8(1):38–49. doi:10.1002/pst.326
# using all the defaults, i.e. latin square crossover design, alpha=0.05, # beta0=1+log(0.95)/log(rd), theta1=0.8, theta2=1.25 power.dp(CV = 0.2, doses = c(1,2,8), n = 15) # # period balanced IBD with 3 doses, 2 periods and 3 sequences, ibd <- matrix(c(1, 2, 3, 2, 3, 1), nrow = 3, ncol = 2) power.dp(CV = 0.2, doses = c(1,2,8), n = 12, design = "IBD", dm = ibd) # considerably lower than 3x3 Latin square# using all the defaults, i.e. latin square crossover design, alpha=0.05, # beta0=1+log(0.95)/log(rd), theta1=0.8, theta2=1.25 power.dp(CV = 0.2, doses = c(1,2,8), n = 15) # # period balanced IBD with 3 doses, 2 periods and 3 sequences, ibd <- matrix(c(1, 2, 3, 2, 3, 1), nrow = 3, ncol = 2) power.dp(CV = 0.2, doses = c(1,2,8), n = 12, design = "IBD", dm = ibd) # considerably lower than 3x3 Latin square
This function performs the power calculation of the BE decision via the FDA’s method for highly variable narrow therapeutic index drugs (NTIDs) as described in respective guidances based on simulations. The study design could be the full replicate design 2x2x4 with 4-periods (TRTR|RTRT) or the 2x2x3 replicate design with 3-periods and sequences TRT|RTR.
power.HVNTID(alpha = 0.05, theta1, theta2, theta0, CV, n, design = c("2x2x4", "2x2x3"), nsims = 1e+05, details = FALSE, setseed = TRUE)power.HVNTID(alpha = 0.05, theta1, theta2, theta0, CV, n, design = c("2x2x4", "2x2x3"), nsims = 1e+05, details = FALSE, setseed = TRUE)
alpha |
Type I error probability, significance level. Commonly set to 0.05. |
theta1 |
Conventional lower ABE limit to be applied in the FDA procedure. |
theta2 |
Conventional upper ABE limit to be applied in the FDA procedure. |
theta0 |
‘True’ or assumed T/R ratio. |
CV |
Intra-subject coefficient(s) of variation as ratio (not percent).
|
n |
Number of subjects under study. |
design |
Design of the study to be planned. |
nsims |
Number of simulations to be performed to obtain the empirical power. Defaults to 100,000 = 1e+5. |
details |
If set to |
setseed |
Simulations are dependent on the starting point of the (pseudo) random number
generator. To avoid differences in power for different runs a |
For deciding BE the study must pass the conventional ABE test (90% CI within the
acceptance range) and additional the test that the ratio of sWT/sWR is < 2.5.
The simulations are done via the distributional properties of the statistical
quantities necessary for deciding BE based on this method.
Details can be found in a document Implementation_scaledABE_sims located in
the /doc sub-directory of the package.
Returns the value of the (empirical) power if argument details=FALSE.
Returns a named vector if argument details=TRUE.
p(BE) is the power, p(BE-ABE) is the power of the ABE test alone and p(BE-sratio)
is the power of the criterion 'ratio of sWT/sWR is <= 2.5' alone.
The FD’s guidances recommend only the full replicate design "2x2x4" (TRTR|RTRT).
The results for the design "2x2x3" (TRT|RTR) are to be considered as experimental since
at present not thorougly tested.
The method is also required by China’s Center of Drug Evaluation.
D. Labes
Food and Drug Administration, Office of Generic Drugs (OGD). Draft Guidance on Dabigatran Etexilate Mesylate. Recommended Jun 2012; Revised Sep 2015, Jul 2017. download
Food and Drug Administration, Office of Generic Drugs (OGD). Draft Guidance on Rivaroxaban. Recommended Sep 2015. download
Food and Drug Administration, Office of Generic Drugs (OGD). Draft Guidance on Edoxaban Tosylate. Recommended May 2017; Revised Mar 2020. download
sampleN.HVNTID
and power.NTIDFDA, sampleN.NTIDFDA for NTIDs with
low variability
# using the defaults: # GMR=0.95, theta1=0.8, theta2=1.25, full replicate design 2x2x4, 100,000 simulations # and a CV of 0.3 (=30%) for both Reference and Test, with 24 subjects, balanced power.HVNTID(CV = 0.3, n = 24) # should give a power of 0.86354# using the defaults: # GMR=0.95, theta1=0.8, theta2=1.25, full replicate design 2x2x4, 100,000 simulations # and a CV of 0.3 (=30%) for both Reference and Test, with 24 subjects, balanced power.HVNTID(CV = 0.3, n = 24) # should give a power of 0.86354
Function calculates of the power of the one-sided non-inferiority t-test for normal or log-normal distributed data.
power.noninf(alpha = 0.025, logscale = TRUE, margin, theta0, CV, n, design = "2x2", robust = FALSE)power.noninf(alpha = 0.025, logscale = TRUE, margin, theta0, CV, n, design = "2x2", robust = FALSE)
alpha |
Significance level (one-sided). Defaults here to 0.025. |
logscale |
Should the data used on log-transformed or on original scale? |
theta0 |
‘True’ or assumed T/R ratio or difference. |
margin |
Non-inferiority margin. |
CV |
In case of In case of cross-over studies this is the within-subject CV, in case of a parallel-group design the CV of the total variability. |
n |
Number of subjects under study. |
design |
Character string describing the study design. |
robust |
Defaults to |
The power is calculated exact via non-central t-distribution.
Notes on the underlying hypotheses
If the supplied margin is < 0 (logscale=FALSE) or < 1 (logscale=TRUE),
then it is assumed higher response values are better. The hypotheses are H0: theta0 <= margin vs. H1: theta0 > margin
where theta0 = mean(test)-mean(reference) if logscale=FALSE
or H0: log(theta0) <= log(margin) vs. H1: log(theta0) > log(margin)
where theta0 = mean(test)/mean(reference) if logscale=TRUE.
If the supplied margin is > 0 (logscale=FALSE) or > 1 (logscale=TRUE),
then it is assumed lower response values are better. The hypotheses are H0: theta0 >= margin vs. H1: theta0 < margin
where theta0 = mean(test)-mean(reference) if logscale=FALSE
or H0: log(theta0) >= log(margin) vs. H1: log(theta0) < log(margin)
where theta0 = mean(test)/mean(reference) if logscale=TRUE.
This latter case may also be considered as ‘non-superiority’.
Value of power according to the input arguments.
The function does not vectorize if design is a vector.
The function vectorize properly if CV or theta0 are vectors.
Other vector input is not tested yet.
This function does not rely on TOST but may be useful in planning BE studies
if the question is not equivalence but ‘non-superiority’.
Hint: Evaluation of Fluctuation in the EMEA’s Note for Guidance between a modified release
formulation and an immediate release product.
D. Labes
Julious SA. Sample sizes for clinical trials with Normal data. Stat Med. 2004;23(12):1921–86. doi:10.1002/sim.1783
# using all the defaults: margin=0.8, theta0=0.95, alpha=0.025 # log-transformed, design="2x2" # should give: 0.4916748 power.noninf(CV=0.3, n=24)# using all the defaults: margin=0.8, theta0=0.95, alpha=0.025 # log-transformed, design="2x2" # should give: 0.4916748 power.noninf(CV=0.3, n=24)
This function performs the power calculation of the BE decision via the FDA’s method for narrow therapeutic index drugs (NTIDs) by simulations. The study design could be the full replicate design 2x2x4 with 4-periods (TRTR|RTRT) or the 2x2x3 replicate design with sequences TRT|RTR.
power.NTID(alpha = 0.05, theta1, theta2, theta0, CV, n, design=c("2x2x4", "2x2x3"), nsims = 1e+05, details = FALSE, setseed = TRUE)power.NTID(alpha = 0.05, theta1, theta2, theta0, CV, n, design=c("2x2x4", "2x2x3"), nsims = 1e+05, details = FALSE, setseed = TRUE)
alpha |
Type I error probability, significance level. Conventionally mostly set to 0.05. |
theta1 |
Conventional lower ABE limit to be applied in the FDA procedure. |
theta2 |
Conventional upper ABE limit to be applied in the FDA procedure. |
theta0 |
‘True’ or assumed T/R ratio. |
CV |
Intra-subject coefficient(s) of variation as ratio (not percent).
|
n |
Number of subjects under study. |
design |
Design of the study to be planned. |
nsims |
Number of simulations to be performed to obtain the empirical power. Defaults to 100,000 = 1e+5. |
details |
If set to |
setseed |
Simulations are dependent on the starting point of the (pseudo) random number
generator. To avoid differences in power for different runs a |
The linearized scaled ABE criterion is calculated according to the SAS code
given in the FDA’s guidances. For deciding BE the study must pass that criterion,
the conventional ABE test, and that the upper confidence limit of
σwT/σwR ≤ 2.5.
The simulations are done via the distributional properties of the statistical
quantities necessary for deciding BE based on these method.
Details can be found in a document Implementation_scaledABE_sims located in
the /doc sub-directory of the package.
Returns the value of the (empirical) power if argument details = FALSE.
Returns a named vector if argument details = TRUE, where
p(BE) is the (overall) power, p(BE-sABEc) is the power of the BE test via scaled ABE criterion alone,
p(BE-ABE) is the power of the conventional ABE test alone, and p(BE-sratio)
is the power of the criterion ‘upper confidence limit of
σwT/σwR ≤ 2.5’ alone.
The FDA’s method is described for the ABE limits 0.8 ... 1.25 only. Setting theta1,
theta2 to other values may not be reasonable and is not tested.
The results for the design "2x2x3" are to be considered as experimental since at present not thorougly tested.
The method is also required by China’s Center of Drug Evaluation.
D. Labes
Food and Drug Administration, Office of Generic Drugs (OGD). Draft Guidance on Warfarin Sodium. Recommended Dec 2012. download
Food and Drug Administration, Center for Drug Evaluation and Research (CDER). Draft Guidance for Industry. Bioequivalence Studies with Pharmacokinetic Endpoints for Drugs Submitted Under an ANDA. August 2021. download
Yu LX, Jiang W, Zhang X, Lionberger R, Makhlouf F, Schuirmann DJ, Muldowney L, Chen ML, Davit B, Conner D, Woodcock J. Novel bioequivalence approach for narrow therapeutic index drugs. Clin Pharmacol Ther. 2015;97(3):286–91. doi:10.1002/cpt.28
Jiang W, Makhlouf F, Schuirmann DJ, Zhang X, Zheng N, Conner D, Yu LX, Lionberger R. A Bioequivalence Approach for Generic Narrow Therapeutic Index Drugs: Evaluation of the Reference-Scaled Approach and Variability Comparison Criterion. AAPS J. 2015;17(4):891–901. doi:10.1208/s12248-015-9753-5
Endrényi L, Tóthfalusi L. Determination of Bioequivalence for Drugs with Narrow Therapeutic Index: Reduction of the Regulatory Burden. J Pharm Pharm Sci. 2013;16(5):676–82. open access
sampleN.NTID
and power.HVNTID, sampleN.HVNTID for NTIDs with
high variability
# using the all defaults: # GMR=0.975, theta1=0.8, theta2=1.25, 100,000 simulations # and a CV of 0.1 (= 10%) with 12 subjects, balanced power.NTID(CV = 0.1, n = 12) # should give a power of 0.62553# using the all defaults: # GMR=0.975, theta1=0.8, theta2=1.25, 100,000 simulations # and a CV of 0.1 (= 10%) with 12 subjects, balanced power.NTID(CV = 0.1, n = 12) # should give a power of 0.62553
Calculates the power of the test of equivalence of the ratio of two means
with normality on original scale.
This test is based on Fieller’s confidence (‘fiducial’) interval and Sasabuchi’s
test (a TOST procedure as well).
power.RatioF(alpha = 0.025, theta1 = 0.8, theta2, theta0 = 0.95, CV, CVb, n, design = "2x2", setseed=TRUE)power.RatioF(alpha = 0.025, theta1 = 0.8, theta2, theta0 = 0.95, CV, CVb, n, design = "2x2", setseed=TRUE)
alpha |
Type I error probability, aka significance level. |
theta1 |
Lower bioequivalence limit. Typically 0.8 (default). |
theta2 |
Upper bioequivalence limit. Typically 1.25. |
theta0 |
‘True’ or assumed T/R ratio. Typically set to 0.95 for planning. |
CV |
Coefficient of variation as ratio. In case of |
CVb |
CV of the between-subject variability. Only necessary for |
n |
Number of subjects to be planned. |
design |
A character string describing the study design. |
setseed |
If set to |
The power is calculated exact using the bivariate non-central t-distribution
via function pmvt of the package mvtnorm.
Due to the calculation method of the used package mvtnorm – randomized
Quasi-Monte-Carlo – these probabilities are dependent from the state of the
random number generator within the precision of the power.
See argument setseed.
Value of power according to the input.
This function is intended for studies with clinical endpoints where the 95% confidence intervals are usually used for equivalence testing.
Therefore, alpha defaults here to 0.025 (see EMEA 2000).
The formulas given in the references rely on the assumption of equal variances
in the two treatment groups for the parallel group design or on assuming equal
within-subject and between-subject variabilities for the 2×2 crossover design.
D. Labes
Fieller EC. Some Problems in Interval Estimation. J Royal Stat Soc B. 1954;16(2):175–85. doi:10.1111/j.2517-6161.1954.tb00159.x
Sasabuchi S. A test of a multivariate normal mean with composite hypotheses determined by linear inequalities. Biometrika. 1980;67(2):429–39. doi:10.1093/biomet/67.2.429
Hauschke D, Kieser M, Diletti E, Burke M. Sample size determination for proving equivalence based on the ratio of two means for normally distributed data. Stat Med. 1999;18(1):93–105.
Hauschke D, Steinijans V, Pigeot I. Bioequivalence Studies in Drug Development. Chichester: Wiley; 2007. Chapter 10.
European Agency for the Evaluation of Medicinal Products, CPMP. Points to Consider on Switching between Superiority and Non-Inferiority. London, 27 July 2000. CPMP/EWP/482/99
# power for alpha=0.025, ratio0=0.95, theta1=0.8, theta2=1/theta1=1.25 # within-subject CV=0.2, between-subject CV=0.4 # 2x2 crossover study, n=24 # using all the defaults: power.RatioF(CV = 0.2, CVb = 0.4, n = 24) # gives [1] 0.7315357# power for alpha=0.025, ratio0=0.95, theta1=0.8, theta2=1/theta1=1.25 # within-subject CV=0.2, between-subject CV=0.4 # 2x2 crossover study, n=24 # using all the defaults: power.RatioF(CV = 0.2, CVb = 0.4, n = 24) # gives [1] 0.7315357
This function performs the power calculation of the BE decision via linearized scaled ABE criterion by simulations as recommended by the FDA.
power.RSABE(alpha = 0.05, theta1, theta2, theta0, CV, n, design = c("2x3x3", "2x2x4", "2x2x3"), regulator, nsims = 1e+05, details = FALSE, setseed=TRUE)power.RSABE(alpha = 0.05, theta1, theta2, theta0, CV, n, design = c("2x3x3", "2x2x4", "2x2x3"), regulator, nsims = 1e+05, details = FALSE, setseed=TRUE)
alpha |
Type I error probability, significance level. Conventionally mostly set to 0.05. |
theta1 |
Conventional lower ABE limit to be applied in the mixed procedure if
|
theta2 |
Conventional upper ABE limit to be applied in the mixed procedure if
|
theta0 |
‘True’ or assumed T/R ratio. |
CV |
Intra-subject coefficient(s) of variation as ratio (not percent).
|
n |
Number of subjects under study. |
design |
Design of the study. |
regulator |
Regulatory settings for RSABE. |
nsims |
Number of simulations to be performed to obtain the empirical power.
Defaults to 100,000 = 1e+5. |
details |
If set to |
setseed |
Simulations are dependent on the starting point of the (pseudo) random number
generator. To avoid differences in power for different runs a |
The linearized scaled ABE criterion is calculated according to the SAS code
given in the FDA’s progesterone guidance.
The simulations are done via the distributional properties of the statistical
quantities necessary for deciding BE based on scaled ABE criterion.
Details can be found in a document Implementation_scaledABE_simsVx.yy.pdf
located in the /doc sub-directory of the package.
If a CVcap is defined for the regulator, the BE decision is based on the inclusion
of the CI in the capped widened acceptance limits in case of CVwR > CVcap. This
resembles method ‘Howe-EMA’ in Muñoz et al. and is the standard behavior now if
regulator="EMA" is choosen.
Returns the value of the (empirical) power if argument details=FALSE.
Returns a named vector if argument details=TRUE.
p(BE) is the power, p(BE-sABEc) is the power of the scaled ABE criterion alone
and p(BE-pe) is the power of the criterion ‘point estimat within acceptance
range’ alone.
p(BE-ABE) is the power of the conventional ABE test given for comparative purposes.
Although some designs are more ‘popular’ than others, power calculations are valid for all of the following designs:
"2x2x4" |
TRTR | RTRT |
| TRRT | RTTR | |
| TTRR | RRTT | |
"2x2x3" |
TRT | RTR |
| TRR | RTT | |
"2x3x3" |
TRR | RTR | RRT |
In case of the design "2x2x3"" heteroscedasticity (i.e., CVwT != CVwR) may
lead to poor agreement of the power values compared to those calculated via the
‘classical’ way of subject data simulations if the design is unbalanced in respect
to the number of subjects in the sequence groups. Therefore, the function
issues a warning for that cases.
D. Labes
Food and Drug Administration, Office of Generic Drugs (OGD). Draft Guidance on Progesterone. Recommended Apr 2010. Revised Feb 2011. download
Tóthfalusi, L, Endrényi, L. Sample Sizes for Designing Bioequivalence Studies for Highly Variable Drugs. J Pharm Pharmaceut Sci. 2011;15(1):73–84. open access
Tóthfalusi L, Endrényi L, García Arieta A. Evaluation of Bioequivalence for Highly Variable Drugs with Scaled Average Bioequivalence. Clin Pharmacokin. 2009;48(11):725–43. doi:10.2165/11318040-000000000-00000
Muñoz J, Alcaide D, Ocaña J. Consumer’s risk in the EMA and FDA regulatory approaches for bioequivalence in highly variable drugs. Stat Med. 2015;35(12):1933–43. doi:10.1002/sim.6834
# using all the defaults: # design="2x3x3" = partial replicate # ABE limits, PE constraint 0.8-1.25 # true ratio = 0.90, 1E+5 simulations power.RSABE(CV = 0.4, n = 36) # should give # [1] 0.83634 # # to explore the simulation error due to the state of the # random number generator power.RSABE(CV = 0.4, n = 36, setseed = FALSE) # will give something like # [1] 0.83725 # # explore pure RSABE (without mixed method, without pe constraint) rs <- reg_const("FDA") rs$CVswitch <- 0 rs$pe_constr <- FALSE power.RSABE(CV = 0.4, n = 36, regulator = rs) # should give # [1] 0.84644# using all the defaults: # design="2x3x3" = partial replicate # ABE limits, PE constraint 0.8-1.25 # true ratio = 0.90, 1E+5 simulations power.RSABE(CV = 0.4, n = 36) # should give # [1] 0.83634 # # to explore the simulation error due to the state of the # random number generator power.RSABE(CV = 0.4, n = 36, setseed = FALSE) # will give something like # [1] 0.83725 # # explore pure RSABE (without mixed method, without pe constraint) rs <- reg_const("FDA") rs$CVswitch <- 0 rs$pe_constr <- FALSE power.RSABE(CV = 0.4, n = 36, regulator = rs) # should give # [1] 0.84644
These function performs the power calculation of the BE decision via
the reference scaled ABE based on subject data simulations.
Implemented are the methods ABEL, Hyslop and ‘exact’ as described in the
references.
The estimation method of the key statistics needed to perform the RSABE decision
is the usual ANOVA.
power.RSABE2L.sdsims(alpha = 0.05, theta1, theta2, theta0, CV, n, design = c("2x3x3", "2x2x4", "2x2x3"), design_dta = NULL, SABE_test = "exact", regulator, nsims = 1e+05, details = FALSE, setseed = TRUE, progress)power.RSABE2L.sdsims(alpha = 0.05, theta1, theta2, theta0, CV, n, design = c("2x3x3", "2x2x4", "2x2x3"), design_dta = NULL, SABE_test = "exact", regulator, nsims = 1e+05, details = FALSE, setseed = TRUE, progress)
alpha |
Type I error probability, significance level. Conventionally mostly set to 0.05. |
theta1 |
Conventional lower ABE (Average Bioequivalence) limit to be applied
in the mixed procedure if |
theta2 |
Conventional upper ABE limit to be applied in the mixed procedure if
|
theta0 |
‘True’ or assumed T/R ratio. |
CV |
Intra-subject coefficient(s) of variation as ratio (not percent).
|
n |
Number of subjects under study. |
design |
Design of the study. |
design_dta |
Alternatively to using the arguments |
SABE_test |
This argument specifies the test method to be used for the reference scaled
ABE decision. |
regulator |
Regulatory settings for the widening of the BE acceptance limits. |
nsims |
Number of simulations to be performed to obtain the empirical power.
Defaults to 100,000 = 1e+05. |
details |
If set to |
setseed |
Simulations are dependent on the starting point of the (pseudo) random number
generator. To avoid differences in power for different runs a |
progress |
Should a progressbar be shown? Defaults to |
The methods rely on the analysis of log-transformed data, i.e., assumes a
log-normal distribution on the original scale.
The data.frame with columns subject, sequence, period and tmt
necessary for evalution of simulated subject data is constructed internally from
the arguments design and n or may be given user defined via the argument
design_dta. The last option is usefull if missing data have to be considered
or if designs have to be evaluated which are not in the list of argument
design.
The estimation method for obtaining the statistics necessary to perform the
reference scaled ABE decision is the usual ANOVA with effects treatment, period,
sequence and subject within sequence for the evaluation of all data and period,
sequence and subject within sequence for the evaluation of the Reference formulation
data only.
The SABE tests implemented are:
"exact" |
‘exact’ based method of the two Laszlós (see references, called there ‘ncTOST’) |
"ABEL" |
Average bioequivalence with expanding limits |
"hyslop" |
BE decision via the linearized RSABE criterion and its upper 95% CI |
"fda" |
Hyslop with an additional bias correction term as implemented in the SAS code of the |
| FDA’s Guidance on Progesterone. |
Returns the value of the (empirical) power if argument details=FALSE.
Returns a named vector if argument details=TRUE.p(BE) is the power, p(BE-RSABE) is the power of using the reference
scaled ABE alone, and p(BE-pe) is the power of the criterion
‘point estimate within acceptance range’ alone. p(BE-ABE) is the power of
the conventional ABE test given for comparative purposes.
Although some designs are more ‘popular’ than others, power calculations are valid for all of the following designs:
"2x2x4" |
TRTR | RTRT |
| TRRT | RTTR | |
| TTRR | RRTT | |
"2x2x3" |
TRT | RTR |
| TRR | RTT | |
"2x3x3" |
TRR | RTR | RRT |
The function is relatively slow. The run-time for 1 Mio. simulations
is between ~ 1 up to 6 minutes for n=12 or n=120 and 1 Mio. sim’s
(see the call under examples) on a machine with an Intel core i7 processor.
Thus be patient and go for a cup of coffee if you use this function with higher
sample sizes and aim for estimating the type 1 error!
D. Labes
Food and Drug Administration, Office of Generic Drugs (OGD). Draft Guidance on Progesterone. Recommended Apr 2010. Revised Feb 2011. download
Tóthfalusi L, Endrényi L. An Exact Procedure for the Evaluation of Reference-Scaled Average Bioequivalence. AAPS J. 2016;18(2):476–89. doi:10.1208/s12248-016-9873-6.
Tóthfalusi L, Endrényi L. Algorithms for evaluating reference scaled average bioequivalence: power, bias, and consumer risk. Stat Med. 2017;36(27):4378–4390. doi:10.1002/sim.7440
# Not run due to timing policy of CRAN # pure EMA settings without mixed procedure, cap on widening and PE constraint # as in the reference from 2017 reg <- reg_const("EMA") reg$CVswitch <- 0 reg$CVcap <- Inf reg$pe_constr <- FALSE reg$name <- "EMA pure" power.RSABE2L.sds(CV = 0.4, n = 12, theta0 = exp(0.05), design = "2x2x4", regulator = reg, nsims = 50000) # should give: # [1] 0.46504 (compared to 47.1% in the 2017 reference)# Not run due to timing policy of CRAN # pure EMA settings without mixed procedure, cap on widening and PE constraint # as in the reference from 2017 reg <- reg_const("EMA") reg$CVswitch <- 0 reg$CVcap <- Inf reg$pe_constr <- FALSE reg$name <- "EMA pure" power.RSABE2L.sds(CV = 0.4, n = 12, theta0 = exp(0.05), design = "2x2x4", regulator = reg, nsims = 50000) # should give: # [1] 0.46504 (compared to 47.1% in the 2017 reference)
These function performs the power calculation of the BE decision via scaled (widened) BE acceptance limits by simulations.
power.scABEL(alpha = 0.05, theta1, theta2, theta0, CV, n, design = c("2x3x3", "2x2x4", "2x2x3"), regulator, nsims, details = FALSE, setseed = TRUE)power.scABEL(alpha = 0.05, theta1, theta2, theta0, CV, n, design = c("2x3x3", "2x2x4", "2x2x3"), regulator, nsims, details = FALSE, setseed = TRUE)
alpha |
Type I error probability, significance level. Conventionally mostly set to 0.05. |
theta1 |
Conventional lower ABE limit to be applied in the mixed procedure if
|
theta2 |
Conventional upper ABE limit to be applied in the mixed procedure if
|
theta0 |
‘True’ or assumed T/R ratio. |
CV |
Intra-subject coefficient(s) of variation as ratio (not percent).
|
n |
Number of subjects under study. |
design |
Design of the study. |
regulator |
Regulatory settings for the widening of the BE acceptance limits. |
nsims |
Number of simulations to be performed to obtain the empirical power.
Defaults to 100,000 = 1e+05. |
details |
If set to |
setseed |
Simulations are dependent on the starting point of the (pseudo) random number
generator. To avoid differences in power for different runs a |
The methods rely on the analysis of log-transformed data, i.e., assume a
log-normal distribution on the original scale.
The widened BE acceptance limits will be calculated by the formula [L, U] = exp(-/+ r_const * sWR)
with r_const the regulatory constant and sWR the standard deviation of the within
subjects variability of the Reference. r_const = 0.76 (~log(1.25)/0.29356) is used
in case of regulator="EMA" or regulator="HC" and in case of
regulator="FDA" r_const = 0.89257... (log(1.25)/0.25).
If the CVwR of the Reference is < CVswitch=0.3 the conventional ABE limits
apply (mixed procedure).
In case of regulator="EMA" a cap is placed on the widened limits if
CVwR>0.5, i.e., the widened limits are held at value calculated for CVwR=0.5.
In case of regulator="HC" the capping is done such that the acceptance
limits are 0.6666 ... 1.5 at maximum.
The case of regulator="GCC" is treatd as special case of ABEL with
CVswitch = CVcap = 0.3. The r_const = log(1.25)/CV2se(0.3) assures that for CV>0.3
the widened BE limits of 0.7 ... 1.3333 are used.
The simulations are done via the distributional properties of the statistical
quantities necessary for deciding BE based on widened ABEL.
For more details see the document Implementation_scaledABE_simsVx.yy.pdf in the
/doc sub-directory of the package.
Function power.scABEL() implements the simulation via distributional
characteristics of the ‘key’ statistics obtained from the EMA recommended
evaluation via ANOVA if regulator="EMA" or if the regulator component
est_method is set to "ANOVA" if regulator is an object of class 'regSet'.
Otherwise the simulations are based on the distributional characteristis of the
‘key’ statistics obtained from evaluation via intra-subject contrasts (ISC),
as recommended by the FDA.
Returns the value of the (empirical) power if argument details=FALSE.
Returns a named vector if argument details=TRUE.
p(BE) is the power, p(BE-ABEL) is the power of the widened ABEL criterion alone
and p(BE-pe) is the power of the criterion ‘point estimate within acceptance
range’ alone. p(BE-ABE) is the power of the conventional ABE test given for
comparative purposes.
Although some designs are more ‘popular’ than others, power calculations are valid for all of the following designs:
"2x2x4" |
TRTR | RTRT |
| TRRT | RTTR | |
| TTRR | RRTT | |
"2x2x3" |
TRT | RTR |
| TRR | RTT | |
"2x3x3" |
TRR | RTR | RRT |
Cross-validation of the simulations as implemented here and via the ‘classical’
subject data simulation have shown somewhat unsatisfactory results for the
2x3x3 design if the variabilities for Test and Reference are different and/or sequences exteremly unbalanced.
The function power.scABEL() therefore gives a warning if calculations
with different CVwT and CVwR are requested for the 2x3x3 partial replicate design. For "EMA" subject simulations are provided in power.scABEL.sdsims.
For more details see the above mentioned document Implementation_scaledABE_simsVy.xx.pdf.
In case of regulator="FDA" the (empirical) power is only approximate since
the BE decision method is not exactly what is expected by the FDA. But the “Two Laszlós” state that the scABEL method should be ‘operational equivalent’ to the
FDA method.
To get the power for the FDA favored method via linearized scaled ABE criterion
use function power.RSABE.
In case of regulator="HC" (based on ISC), power is also only approximative since Health Canada recommends an evaluation via mixed model approach. This could only implemented via
subject data simulations which are very time consuming. But ISC may be a good
substitute.
D. Labes
Tóthfalusi L, Endrényi L. Sample Sizes for Designing Bioequivalence Studies for Highly Variable Drugs. J Pharm Pharmaceut Sci. 2011;15(1):73–84. open source
sampleN.scABEL, power.RSABE, reg_const
# using all the defaults: # design="2x3x3", EMA regulatory settings # PE constraint 0.8-1.25, cap on widening if CV>0.5 # true ratio=0.90, 1E+6 simulations power.scABEL(CV = 0.4, n = 29) # should give: # Unbalanced design. n(i)=10/10/9 assumed. # [1] 0.66113 # # with details=TRUE to view the computational time and components power.scABEL(CV = 0.5, n = 54, theta0 = 1.15, details = TRUE) # should give (times may differ depending on your machine): # 1e+05sims. Time elapsed (sec): 0.07 # # p(BE) p(BE-wABEL) p(BE-pe) p(BE-ABE) # 0.81727 0.82078 0.85385 0.27542 # # exploring 'pure ABEL' with the EMA regulatory constant # (without mixed method, without capping, without pe constraint) rs <- reg_const("EMA") rs$CVswitch <- 0 rs$CVcap <- Inf rs$pe_constr <- FALSE power.scABEL(CV = 0.5, n = 54, theta0 = 1.15, regulator = rs) # should give # [1] 0.8519# using all the defaults: # design="2x3x3", EMA regulatory settings # PE constraint 0.8-1.25, cap on widening if CV>0.5 # true ratio=0.90, 1E+6 simulations power.scABEL(CV = 0.4, n = 29) # should give: # Unbalanced design. n(i)=10/10/9 assumed. # [1] 0.66113 # # with details=TRUE to view the computational time and components power.scABEL(CV = 0.5, n = 54, theta0 = 1.15, details = TRUE) # should give (times may differ depending on your machine): # 1e+05sims. Time elapsed (sec): 0.07 # # p(BE) p(BE-wABEL) p(BE-pe) p(BE-ABE) # 0.81727 0.82078 0.85385 0.27542 # # exploring 'pure ABEL' with the EMA regulatory constant # (without mixed method, without capping, without pe constraint) rs <- reg_const("EMA") rs$CVswitch <- 0 rs$CVcap <- Inf rs$pe_constr <- FALSE power.scABEL(CV = 0.5, n = 54, theta0 = 1.15, regulator = rs) # should give # [1] 0.8519
These function performs the power calculation of the BE decision via
scaled (widened) BE acceptance limits based on subject data simulations.
This function has an alias power.scABEL.sds().
power.scABEL.sdsims(alpha = 0.05, theta1, theta2, theta0, CV, n, design = c("2x3x3", "2x2x4", "2x2x3"), design_dta=NULL, regulator, nsims = 1e+05, details = FALSE, setseed = TRUE, progress)power.scABEL.sdsims(alpha = 0.05, theta1, theta2, theta0, CV, n, design = c("2x3x3", "2x2x4", "2x2x3"), design_dta=NULL, regulator, nsims = 1e+05, details = FALSE, setseed = TRUE, progress)
alpha |
Type I error probability, significance level. Conventionally mostly set to 0.05. |
theta1 |
Conventional lower ABE limit to be applied in the mixed procedure if
|
theta2 |
Conventional upper ABE limit to be applied in the mixed procedure if
|
theta0 |
‘True’ or assumed T/R ratio. |
CV |
Intra-subject coefficient(s) of variation as ratio (not percent).
|
n |
Number of subjects under study. |
design |
Design of the study to be planned. |
design_dta |
Alternatively to using the arguments |
regulator |
Regulatory settings for the widening of the BE acceptance limits. |
nsims |
Number of simulations to be performed to obtain the empirical power.
Defaults to 100,000 = 1e+05. |
details |
If set to |
setseed |
Simulations are dependent on the starting point of the (pseudo) random number
generator. To avoid differences in power for different runs a |
progress |
Should a progressbar be shown? Defaults to |
The methods rely on the analysis of log-transformed data, i.e., assume a
log-normal distribution on the original scale.
The widened BE acceptance limits will be calculated by the formula [L, U] = exp(± r_const * sWR)
with r_const the regulatory constant and sWR the standard deviation of the within
subjects variability of the Reference. r_const = 0.76 (~log(1.25)/0.29356) is used
in case of regulator="EMA".
If the CVwR of the Reference is < CVswitch=0.3 the conventional ABE limits
apply (mixed procedure).
In case of regulator="EMA" a cap is placed on the widened limits if
CVwr>0.5, i.e., the widened limits are held at value calculated for CVwR=0.5.
The simulations are done by simulating subject data (all effects fixed except the
residuals) and evaluating these data via ANOVA of all data to get the point estimate
of T vs. R along with its 90% CI and an ANOVA of the data under R(eference) only
to get an estimate of s2wR.
The data.frame with columns subject, sequence, period and tmt
necessary for evalution of simulated subject data is constructed internally from
the arguments design and n or may be given user defined via the argument
design_dta. The last option is usefull if missing data have to be considered
or if designs have to be evaluated which are not in the list of argument
design.
This feature is experimental in the sense that the data.frame is not checked
for complying with the assumed structure.
Returns the value of the (empirical) power if argument details=FALSE.
Returns a named vector if argument details=TRUE.
p(BE) is the power, p(BE-wABEL) is the power of the widened ABEL criterion alone
and p(BE-pe) is the power of the criterion 'point estimat within acceptance
range' alone. p(BE-ABE) is the power of the conventional ABE test given for
comparative purposes.
The function is mainly intended for crosscheck of power.scABEL() results.
But may be mandatory for cases where power.scABEL() results are inaccurate
(low sample sizes and/or heteroscedasticity).
It is relatively slow. The run-time of this function doing 1 Mio sims is between
~ 7-8 sec for n=12 and ~ 3-4 min for n=120 on a machine with an Intel core i7 processor.
Thus be patient and go for a cup of coffee if you use this function with high
sample sizes!
D. Labes, B. Lang
Tóthfalusi L, Endrényi L. Sample Sizes for Designing Bioequivalence Studies for Highly Variable Drugs. J Pharm Pharmaceut Sci. 2011;15(1):73–84. open source
# using all the defaults: # design="2x3x3", EMA regulatory settings # PE constraint 0.8-1.25, cap on widening if CV>0.5 # true ratio=0.90, 1E+5 simulations power.scABEL.sdsims(CV = 0.4, n = 36) # should give: # [1] 0.74321# using all the defaults: # design="2x3x3", EMA regulatory settings # PE constraint 0.8-1.25, cap on widening if CV>0.5 # true ratio=0.90, 1E+5 simulations power.scABEL.sdsims(CV = 0.4, n = 36) # should give: # [1] 0.74321
Calculates the exact or approximate power of the two-one-sided t-tests (TOST) procedure for various study designs used in BE studies.
power.TOST(alpha = 0.05, logscale = TRUE, theta1, theta2, theta0, CV, n, design = "2x2", method="exact", robust=FALSE)power.TOST(alpha = 0.05, logscale = TRUE, theta1, theta2, theta0, CV, n, design = "2x2", method="exact", robust=FALSE)
alpha |
Significance level (one-sided). Commonly set to 0.05. |
logscale |
Should the data used on log-transformed or on original scale? |
theta0 |
‘True’ or assumed T/R ratio or difference. |
theta1 |
Lower (bio-)equivalence limit. |
theta2 |
Upper (bio-)equivalence limit. |
CV |
In case of In case of cross-over studies this is the within-subject CV, in case of a parallel-group design the CV of the total variability. |
n |
Number of subjects under study. |
design |
Character string describing the study design. |
method |
Method for calculation of the power. |
robust |
Defaults to |
The exact calculations of power are based on Owen’s Q-function or by direct
integration of the bivariate non-central t-distribution via function
pmvt of package mvtnorm.
Approximate power is implemented via the non-central t-distribution
or the ‘shifted’ central t-distribution.
The formulas cover balanced and unbalanced studies w.r.t (sequence) groups.
In case of parallel group design and higher order crossover designs
(replicate crossover or crossover with more than two treatments) the calculations
are based on the assumption of equal variances for Test and Reference products
under consideration.
The formulas for the paired means 'design' do not take a correlation parameter
into account. They are solely based on the paired t-test (TOST of differences = zero).
Value of power according to the input arguments.
Of course it is highly recommended to use the default method="exact" :-).
There is no reason beside testing and for comparative purposes to use an
approximation if the exact method is available.
If n is given as scalar (total sample size) and this number is not
divisible by the number of (sequence) groups of the design an unbalanced design
with small imbalance is assumed. A corresponding message is thrown showing the
assumed numbers of subjects in (sequence) groups.
The function does not vectorize properly if design is a vector.
The function vectorizes properly if CV or theta0 are vectors.
Other vector input is not tested yet.
The former function power2.TOST() designd to handle unbalanced studies is
defunct since power.TOST() handles balanced as well as unbalanced designs.
D. Labes, direct integration of bivariate non-central t-distribution by B. Lang
Phillips KF. Power of the Two One-Sided Tests Procedure in Bioequivalence. J Pharmacokin Biopharm. 1990;18(2):137–44. doi:10.1007/BF01063556
Diletti D, Hauschke D, Steinijans VW. Sample Size Determination for Bioequivalence Assessment by Means of Confidence Intervals. Int J Clin Pharmacol Ther Toxicol. 1991;29(1):1–8.
# power for the 2x2 cross-over design with 24 subjects and CV 25% # using all the other default values power.TOST(CV = 0.25, n = 24) # should give: [1] 0.7391155 # nct approximation very good for this configuration power.TOST(CV = 0.25, n = 24, method = "nct") # gives also: [1] 0.7391155 # shifted-central-t approximation power.TOST(CV = 0.25, n = 24, method = "shifted") # gives: [1] 0.7328894 # power for the 2x2 cross-over study with 24 subjects, CV 25% # and 2 drop-outs in the same sequence group (unbalanced study) power.TOST(CV=0.25, n=c(10,12)) # should give: [1] 0.6912935 # not the same compared to the balanced setting power.TOST(CV=0.25, n=22) # should give: [1] 0.6953401# power for the 2x2 cross-over design with 24 subjects and CV 25% # using all the other default values power.TOST(CV = 0.25, n = 24) # should give: [1] 0.7391155 # nct approximation very good for this configuration power.TOST(CV = 0.25, n = 24, method = "nct") # gives also: [1] 0.7391155 # shifted-central-t approximation power.TOST(CV = 0.25, n = 24, method = "shifted") # gives: [1] 0.7328894 # power for the 2x2 cross-over study with 24 subjects, CV 25% # and 2 drop-outs in the same sequence group (unbalanced study) power.TOST(CV=0.25, n=c(10,12)) # should give: [1] 0.6912935 # not the same compared to the balanced setting power.TOST(CV=0.25, n=22) # should give: [1] 0.6953401
The power is obtained via subject data simulations.
Three models are implemented:
gmodel==1 is full FDA model for testing group-by-treatment interaction followed by gmodel==2 or gmodel==3 with data of the biggest group depending on the test of the treatment by group interaction
gmodel==2 is full FDA model but without group-by-treatment interaction
gmodel==3 is model with pooled groups, i.e. without any group term
power.TOST.sds(alpha = 0.05, theta1, theta2, theta0, CV, n, design = c("2x2", "2x2x2", "2x3x3", "2x2x4", "2x2x3"), design_dta = NULL, grps = 2, ngrp = NULL, gmodel = 2, p.level=0.1, nsims = 1e+05, details = FALSE, setseed = TRUE, progress)power.TOST.sds(alpha = 0.05, theta1, theta2, theta0, CV, n, design = c("2x2", "2x2x2", "2x3x3", "2x2x4", "2x2x3"), design_dta = NULL, grps = 2, ngrp = NULL, gmodel = 2, p.level=0.1, nsims = 1e+05, details = FALSE, setseed = TRUE, progress)
alpha |
Type I error probability, significance level. Conventionally mostly set to 0.05. |
theta1 |
Lower BE limit. Defaults to 0.8 if not given explicitely. |
theta2 |
Upper BE limit. Defaults to 1.25 if not given explicitely. |
theta0 |
‘True’ or assumed T/R ratio. |
CV |
Intra-subject coefficient(s) of variation as ratio (not percent).
|
n |
Number of subjects under study. |
design |
Design of the study to be planned. |
design_dta |
Alternatively to using the arguments |
grps |
Number of (logistical) groups. Defaults to 2. |
ngrp |
Vector of number of subjects in groups. |
gmodel |
Number describing the model incorporating group effects
Defaults to |
p.level |
Significance level of the test of a group-by-treatment interaction.
Defaults to |
nsims |
Number of simulations to be performed to obtain the empirical power.
Defaults to 100,000 = 1e+05. |
details |
If set to |
setseed |
Simulations are dependent on the starting point of the (pseudo) random number
generator. To avoid differences in power for different runs a |
progress |
Should a progressbar be shown? Defaults to |
The power is calculated via subject data sims.
The evaluation of BE is done via 1-2*alpha confidence interval using classical ANOVA
for the models with group effects.
The data.frame with columns subject, sequence, period and tmt
necessary for evaluation of simulated subject data is constructed internally from
the arguments design and n or may be given user defined via the argument
design_dta. The last option is usefull if missing data have to be considered
or if designs have to be evaluated which are not in the list of argument
design.
This feature is experimental in the sense that the data.frame is not checked
for complying with the assumed structure.
The p.value of the test of the group-by-treatment interaction in case of gmodel=1
defaults to p.level = 0.1, the value originally used by the FDA. Later on a value of
p.level = 0.05 was used.
If the group-by-treatment interaction is significant the subsequent BE decision is done with the data of the largest group. If there are more than one with the same size, one gets a warning that this feature – showing BE in all that groups – is not implemented yet. Only the first of the largest groups is tested for BE.
Returns the value of the (empirical) power
The run time of the function may be relatively long.
Take a cup of coffee and be patient.
D. Labes
Schütz H.
Multi-Group Studies in Bioequivalence. To pool or not to pool?
Presentation at BioBriges 2018, Prague. https://bebac.at/lectures/Prague2018.pdf
# power for gmodel=2, 2x2 crossover, grps=3 with even number of subjects power.TOST.sds(CV=0.2, n=18, grps=3) # gives [1] 0.78404 # without considering groups power.TOST.sds(CV=0.2, n=18, gmodel=3) # gives [1] 0.7887# power for gmodel=2, 2x2 crossover, grps=3 with even number of subjects power.TOST.sds(CV=0.2, n=18, grps=3) # gives [1] 0.78404 # without considering groups power.TOST.sds(CV=0.2, n=18, gmodel=3) # gives [1] 0.7887
Power is calculated by simulations of studies (PE via its normal distribution, MSE via its associated χ2 distribution) and application of the two one-sided t-tests. Power is obtained via ratio of studies found BE to the number of simulated studies.
power.TOST.sim(alpha = 0.05, logscale = TRUE, theta1, theta2, theta0, CV, n, design = "2x2", robust = FALSE, setseed = TRUE, nsims = 1e+05)power.TOST.sim(alpha = 0.05, logscale = TRUE, theta1, theta2, theta0, CV, n, design = "2x2", robust = FALSE, setseed = TRUE, nsims = 1e+05)
alpha |
Significance level (one-sided). Commonly set to 0.05. |
logscale |
Should the data used on log-transformed or on original scale? |
theta0 |
‘True’ or assumed T/R ratio or difference. |
theta1 |
Lower (bio-)equivalence limit. |
theta2 |
Upper (bio-)equivalence limit. |
CV |
In case of In case of cross-over studies this is the within-subject CV, in case of a parallel-group design the CV of the total variability. |
n |
Number of subjects under study. |
design |
Character string describing the study design. |
robust |
Defaults to |
setseed |
Simulations are dependent on the starting point of the (pseudo) random number
generator. To avoid differences in power for different runs a |
nsims |
Number of studies to simulate. Defaults to 100,000 = 1E5. |
Value of power according to the input arguments.
This function was intended for internal check of the analytical power
calculation methods. Use of the analytical power calculation methods
(power.TOST()) for real problems is recommended.
For sufficient precision nsims > 1E5 (default) may be necessary.
Be patient if using nsims=1E6. May take some seconds.
D. Labes
# using the default design 2x2, BE range 0.8 ... 1.25, logscale, theta0=0.95 power.TOST.sim(alpha = 0.05, CV = 0.3, n = 12) # should give 0.15054, with nsims=1E6 it will be 0.148533 # exact analytical is power.TOST(alpha = 0.05, CV = 0.3, n = 12) # should give 0.1484695 # very unusual alpha setting power.TOST.sim(alpha = 0.9, CV = 0.3, n = 12) # should give the same (within certain precision) as power.TOST(alpha = 0.95, CV = 0.3, n = 12) # or also within certain precision equal to power.TOST(alpha = 0.95, CV = 0.3, n = 12, method = "mvt") # SAS Proc Power gives here the incorrect value 0.60525# using the default design 2x2, BE range 0.8 ... 1.25, logscale, theta0=0.95 power.TOST.sim(alpha = 0.05, CV = 0.3, n = 12) # should give 0.15054, with nsims=1E6 it will be 0.148533 # exact analytical is power.TOST(alpha = 0.05, CV = 0.3, n = 12) # should give 0.1484695 # very unusual alpha setting power.TOST.sim(alpha = 0.9, CV = 0.3, n = 12) # should give the same (within certain precision) as power.TOST(alpha = 0.95, CV = 0.3, n = 12) # or also within certain precision equal to power.TOST(alpha = 0.95, CV = 0.3, n = 12, method = "mvt") # SAS Proc Power gives here the incorrect value 0.60525
Calculates the p-value(s) of the TOST procedure via students t-distribution given pe, CV and n.
pvalue.TOST(pe, CV, n, logscale = TRUE, theta1, theta2, design = "2x2", robust = FALSE, both = FALSE) pvalues.TOST(pe, CV, n, logscale = TRUE, theta1, theta2, design = "2x2", robust = FALSE, both = TRUE)pvalue.TOST(pe, CV, n, logscale = TRUE, theta1, theta2, design = "2x2", robust = FALSE, both = FALSE) pvalues.TOST(pe, CV, n, logscale = TRUE, theta1, theta2, design = "2x2", robust = FALSE, both = TRUE)
pe |
Observed point estimate of the T/R ratio or difference. |
CV |
In case of In case of cross-over studies this is the within-subject CV, in case of a parallel-group design the CV of the total variability. |
n |
Total number of subjects if given as scalar. |
logscale |
Should the data be used after log-transformation or on original scale? |
theta1 |
Lower (bio-)equivalence limit. |
theta2 |
Upper (bio-)equivalence limit. |
design |
Character string describing the study design. |
robust |
If set to |
both |
Indicates if both p-values (t-tests of pe >= theta1 and pe <= theta2) shall be given
back or only the maximum. |
Returns the p-value(s).
Returns a vector with named elements p.left, p.right if arguments pe and CV
are scalars, else a matrix with columns p.left, p.right. p.left gives the p-value of testing HA1: theta >= theta1
and p.right the p-value of testing HA2: theta <= theta2
against their respective Nulls.
The formulas implemented cover balanced and unbalanced designs.
In case of argument n given as n(total) and is not divisible by the number
of (sequence) groups the total sample size is partitioned to the (sequence)
groups to have small imbalance only. A message is given in such cases.
SAS procedure TTEST with the TOST option names p.left = Upper, p.right= Lower
according to the tail of the t-distribution to be used for obtaining the
p-values.
B. Lang, man page by D. Labes
Schuirmann DJ. A comparison of the two one-sided tests procedure and the power approach for assessing the equivalence of average bioavailability. J Pharmacokin Biopharm. 1987;15:657–80. doi:10.1007/BF01068419
Hauschke D, Steinijans V, Pigeot I. Bioequivalence Studies in Drug Development. Chichester: Wiley; 2007.
# Defaults: 2x2 crossover, log-transformation # BE acceptance limits 0.8 ... 1.25, usual dfs # interested in both p-values pvalues.TOST(pe = 0.95, CV = 0.3, n = 12) # gives the vector (named elements) # p.left p.right # 0.09105601 0.02250985 # i.e. 'left' hypothesis H01: theta<=theta1 can't be rejected # 'right' hypothesis H02: theta>=theta2 can be rejected # max. p-value only as 'overall' pvalue, preferred by Benjamin pvalue.TOST(pe = 0.912, CV = 0.333, n = 24) # should give 0.08777621, i.e., inequivalence can't be rejected # this is operationally identical to CI.BE(pe = 0.912, CV = .333, n = 24) # lower limit = 0.7766 outside 0.8 ... 1.25, i.e., inequivalence can't be rejected# Defaults: 2x2 crossover, log-transformation # BE acceptance limits 0.8 ... 1.25, usual dfs # interested in both p-values pvalues.TOST(pe = 0.95, CV = 0.3, n = 12) # gives the vector (named elements) # p.left p.right # 0.09105601 0.02250985 # i.e. 'left' hypothesis H01: theta<=theta1 can't be rejected # 'right' hypothesis H02: theta>=theta2 can be rejected # max. p-value only as 'overall' pvalue, preferred by Benjamin pvalue.TOST(pe = 0.912, CV = 0.333, n = 24) # should give 0.08777621, i.e., inequivalence can't be rejected # this is operationally identical to CI.BE(pe = 0.912, CV = .333, n = 24) # lower limit = 0.7766 outside 0.8 ... 1.25, i.e., inequivalence can't be rejected
This function may be used to define regulatory settings not implemented yet in PowerTOST.
reg_const(regulator, r_const, CVswitch, CVcap, pe_constr)reg_const(regulator, r_const, CVswitch, CVcap, pe_constr)
regulator |
Name of the regulatory body as a string. Implemented settings are for |
r_const |
Regulatory constant. |
CVswitch |
CV to switch to the widened acceptance limits. |
CVcap |
CV for capping the widening of the acceptance limits. |
pe_constr |
Logical. Shall pe constraint be applied? Defaults to |
Returns an object of class 'regSet', a list with components
name |
Name of the settings |
CVswitch |
see arguments |
r_const |
Regulatory constant |
CVcap |
see arguments |
pe_constr |
see arguments |
est_method |
|
Class 'regSet' has a S3 print method.
The component est_method is automatically set to "ANOVA", except for
regulator="FDA" or regulator="HC" where "ISC" is used.
The former inofficial regulatory settings for regulator="ANVISA" are covered by regulator="EMA" (BEBA Forum).
The settings for CVcap of Health Canada (regulator="HC") were chosen in such a way
that the limits of the acceptance range are capped nearly exact to 1/1.5
up to 1.5. Literally it is given rounded to 3 significant digits.
D. Labes
BEBA Forum. May 2016. online
Health Canada, Therapeutic Products Directorate. Comparative Bioavailability Standards: Formulations Used for Systemic Effects, 2.1.1.8 Highly variable drug products Ottawa, 08 June 2018. online
# to retrieve the EMA settings reg_const("EMA") # to define the old ANVISA settings reg <- reg_const("USER", r_const = 0.76, CVswitch = 0.4, CVcap = 0.5) reg$name <- "Old ANVISA" # Use reg as argument in the scaled ABEL power / sample size functions for# to retrieve the EMA settings reg_const("EMA") # to define the old ANVISA settings reg <- reg_const("USER", r_const = 0.76, CVswitch = 0.4, CVcap = 0.5) reg$name <- "Old ANVISA" # Use reg as argument in the scaled ABEL power / sample size functions for
Estimates the necessary sample size to have at least a given power when two parameters are tested simultaneously.
sampleN.2TOST(alpha = c(0.05, 0.05), targetpower = 0.8, logscale = TRUE, theta0, theta1, theta2, CV, rho, design = "2x2", setseed = TRUE, robust = FALSE, print = TRUE, details = FALSE, imax = 100, nsims = 1e+05)sampleN.2TOST(alpha = c(0.05, 0.05), targetpower = 0.8, logscale = TRUE, theta0, theta1, theta2, CV, rho, design = "2x2", setseed = TRUE, robust = FALSE, print = TRUE, details = FALSE, imax = 100, nsims = 1e+05)
alpha |
Vector; contains one-sided significance level for each of the two TOSTs. |
targetpower |
Power to achieve at least. Must be >0 and <1. |
logscale |
Should the data used on log-transformed or on original scale? |
theta0 |
Vector; contains ‘true’ assumed T/R ratio for each of the two TOSTs. |
theta1 |
Vector; contains lower bioequivalence limit for each of the two TOSTs. |
theta2 |
Vector; contains upper bioequivalence limit for each of the two TOSTs. |
CV |
Vector of coefficient of variations (given as as ratio, e.g., 0.2 for 20%). |
rho |
Correlation between the two PK metrics (e.g., AUC and Cmax) under consideration. This is defined as correlation between the estimator of the treatment difference of PK metric one and the estimator of the treatment difference of PK metric two. Has to be within {–1, +1}. |
design |
Character string describing the study design. |
setseed |
Logical; if |
robust |
Defaults to |
print |
If |
details |
If |
imax |
Maximum number of steps in sample size search. |
nsims |
Number of studies to simulate. Defaults to 100,000 = 1E5. |
The sample size is estimated via iterative evaluation of power of the two TOSTs.
Start value for the sample size search is taken from a large sample approximation
(one TOST) according to Zhang, modified.
The sample size is bound to 4 as minimum.
The estimated sample size gives always the total number of subjects (not subject/sequence in crossovers or subjects/group in parallel designs – like in some other software packages).
A list with the input and results will be returned.
The element name "Sample size" contains the total sample size.
The function does not vectorize properly.
If you need sample sizes with varying CVs, use f.i. for-loops or the apply-family.
If both theta0 are near the acceptance limits then the starting value may not
be a good approximation resulting in a lot of iteration steps; imax may need
to be increased to obtain the required sample size.
B. Lang, D. Labes
Phillips KF. Power for Testing Multiple Instances of the Two One-Sided Tests Procedure. Int J Biostat. 2009;5(1):Article 15.
Hua SY, Xu S, D’Agostino RB Sr. Multiplicity adjustments in testing for bioequivalence. Stat Med. 2015;34(2):215–31. doi:10.1002/sim.6247
Lang B, Fleischer F. Letter to the Editor: Comments on ‘Multiplicity adjustments in testing for bioequivalence’. Stat Med. 2016;35(14):2479–80. doi:10.1002/sim.6488
Zhang P. A Simple Formula for Sample Size Calculation in Equivalence Studies. J Biopharm Stat. 2003;13(3):529–538. doi:10.1081/BIP-120022772
# Sample size for 2x2x2 cross-over design, intra-subject CV = 30% and assumed # ratios of 0.95 for both parameters, and correlation 0.9 between parameters # (using all the other default values) # Should give n=44 with power=0.80879 sampleN.2TOST(theta0 = rep(0.95, 2), CV = rep(0.3, 2), rho = 0.9) # Sample size for a parallel group design, # evaluation on the original (untransformed) scale # BE limits 80 ... 120% = -20% ... +20% of reference, # assumed true BE ratio 0.95% = -5% to reference mean for both parameters, # total CV=20% for both parameters, and correlation 0.9 between parameters # should give n=52 with power=0.80094 sampleN.2TOST(logscale=FALSE, theta0 = rep(-0.05, 2), CV = c(0.2, 0.2), rho = 0.9, design = "parallel")# Sample size for 2x2x2 cross-over design, intra-subject CV = 30% and assumed # ratios of 0.95 for both parameters, and correlation 0.9 between parameters # (using all the other default values) # Should give n=44 with power=0.80879 sampleN.2TOST(theta0 = rep(0.95, 2), CV = rep(0.3, 2), rho = 0.9) # Sample size for a parallel group design, # evaluation on the original (untransformed) scale # BE limits 80 ... 120% = -20% ... +20% of reference, # assumed true BE ratio 0.95% = -5% to reference mean for both parameters, # total CV=20% for both parameters, and correlation 0.9 between parameters # should give n=52 with power=0.80094 sampleN.2TOST(logscale=FALSE, theta0 = rep(-0.05, 2), CV = c(0.2, 0.2), rho = 0.9, design = "parallel")
Performes a sample size estimation for dose-proportionality studies using the power model for cossover (Latin square), parallel group designs or incomplete block designs via a confidence interval equivalence criterion.
sampleN.dp(alpha = 0.05, CV, doses, targetpower = 0.8, beta0, theta1 = 0.8, theta2 = 1/theta1, design = c("crossover", "parallel", "IBD"), dm=NULL, CVb, print = TRUE, details = FALSE, imax = 100)sampleN.dp(alpha = 0.05, CV, doses, targetpower = 0.8, beta0, theta1 = 0.8, theta2 = 1/theta1, design = c("crossover", "parallel", "IBD"), dm=NULL, CVb, print = TRUE, details = FALSE, imax = 100)
alpha |
Type 1 error. Usually set to 0.05. |
CV |
Coefficient of variation. Is intra-subject CV for |
doses |
Vector of dose values under study. At least two doses have to be given. |
targetpower |
Power to achieve at least. Must be >0 and <1. |
beta0 |
‘True’ or assumed slope of the power model. If missing defaults to |
theta1 |
Lower acceptance limit for the ratio of dose normalized means (Rdmn). |
theta2 |
Upper acceptance limit for the ratio of dose normalized means (Rdmn). |
design |
Crossover design (default), parallel group design or incomplete block design (IBD). |
dm |
'Design matrix' of the incomplete block design (IBD) if |
CVb |
Coefficient of variation of the between-subject variability. |
print |
If |
details |
If |
imax |
Maximum number of steps in sample size search. |
The sample size is estimated via iterative evaluation of power.dp().
Start value for the sample size search is taken from a large sample approximation.
The sample size is bound to number of dose or sequence groups as minimum.
Balanced designs are used although this is not absolutely necessary.
The estimated sample size gives always the total number of subjects (not subject/sequence in crossovers or subjects/group in parallel designs – like in some other software packages).
A data.frame with the input and results will be returned.
The Sample size column contains the total sample size.
This function is ‘experimental’ only, since it is not thorougly tested yet.
Especially for design="IBD" reliable test cases are missing.
D. Labes
Chow SC, Liu JP. Design and Analysis of Bioavailability and Bioequivalence Studies. Boca Raton: CRC Press; 3rd edition 2009.
Patterson S, Jones B. Bioequivalence and Statistics in Clinical Pharmacology. Boca Raton: Chapman & Hall/CRC: 2006. p. 239.
(contains presumably a bug)
Sethuraman VS, Leonov S, Squassante L, Mitchell TR, Hale MD. Sample size calculation for the Power Model for dose proportionality studies. Pharm Stat. 2007;6(1):35–41. doi:10.1002/pst.241
Hummel J, McKendrick S, Brindley C, French R. Exploratory assessment of dose proportionality: review of current approaches and proposal for a practical criterion. Pharm. Stat. 2009;8(1):38–49. doi:10.1002/pst.326
# using all the defaults, i.e. crossover design, alpha=0.05 # theta1=0.8, theta2=1.25 but true slope slightly off 1 sampleN.dp(CV = 0.2, doses = c(1, 2, 8), beta0 = 1.02) # should give n=18, power=0.854528 # incomplete block design with 5 doses, 3 periods # from library(crossdes) doses <- c(5, 25, 50, 100, 200) CVb <- mse2CV(0.8) levels <- length(doses) per <- 3 block <- levels*(levels-1)/(per-1) # IBD based on balanced minimal repeated measurements design # gives n=30 and 10 sequences ibd <- crossdes::balmin.RMD(levels, block, per) sampleN.dp(CV = 0.2, doses = doses, beta0 = 1, design = "IBD", dm = ibd, CVb = CVb, targetpower=0.9)# using all the defaults, i.e. crossover design, alpha=0.05 # theta1=0.8, theta2=1.25 but true slope slightly off 1 sampleN.dp(CV = 0.2, doses = c(1, 2, 8), beta0 = 1.02) # should give n=18, power=0.854528 # incomplete block design with 5 doses, 3 periods # from library(crossdes) doses <- c(5, 25, 50, 100, 200) CVb <- mse2CV(0.8) levels <- length(doses) per <- 3 block <- levels*(levels-1)/(per-1) # IBD based on balanced minimal repeated measurements design # gives n=30 and 10 sequences ibd <- crossdes::balmin.RMD(levels, block, per) sampleN.dp(CV = 0.2, doses = doses, beta0 = 1, design = "IBD", dm = ibd, CVb = CVb, targetpower=0.9)
This function performs the sample size estimation for the BE decision via
the FDA’s method for highly variable NTIDs as described in respective guidances based on simulations.
The study designs may be the full replicate design 2x2x4 with 4 periods (TRTR|RTRT) and
the 3-period replicate design 2x2x3 with sequences RTR|TRT.
sampleN.HVNTID(alpha = 0.05, targetpower = 0.8, theta0, theta1, theta2, CV, design = c("2x2x4", "2x2x3"), nsims = 1e+05, nstart, imax = 100, print = TRUE, details = TRUE, setseed = TRUE)sampleN.HVNTID(alpha = 0.05, targetpower = 0.8, theta0, theta1, theta2, CV, design = c("2x2x4", "2x2x3"), nsims = 1e+05, nstart, imax = 100, print = TRUE, details = TRUE, setseed = TRUE)
alpha |
Type I error probability. Per convention mostly set to 0.05. |
targetpower |
Power to achieve at least. Must be >0 and <1. |
theta0 |
‘True’ or assumed T/R ratio. |
theta1 |
Conventional lower ABE limit to be applied in the FDA procedure. |
theta2 |
Conventional upper ABE limit to be applied in the FDA procedure. |
CV |
Intra-subject coefficient(s) of variation as ratio (not percent).
|
design |
Design of the study to be planned. |
nsims |
Number of simulations to be performed to obtain the empirical power. Defaults to 100,000 = 1e+5. |
nstart |
Set this to a start value for the sample size if a previous run failed. |
imax |
Maximum number of steps in sample size search. Defaults to 100. |
print |
If |
details |
If set to |
setseed |
Simulations are dependent on the starting point of the (pseudo) random number
generator. To avoid differences in power values for different runs a
|
For deciding BE the study must pass the conventional ABE test and additionally
the test that the ratio of sWT/sWR is <= 2.5.
The simulations are done via the distributional properties of the statistical
quantities necessary for deciding BE based on these method.
Details can be found in a document Implementation_scaledABE_sims located in
the /doc sub-directory of the package.
The estimated sample size gives always the total number of subjects (not subject/sequence – like in some other software packages).
Returns a data.frame with the input and sample size results.
The Sample size column contains the total sample size.
The nlast column contains the last n value. May be useful for re-starting.
For some input constellations the sample size search may be very time
consuming and will eventually also fail since the start values chosen may
not really reasonable for them.
In case of a failed sample size search you may restart with setting the argument
nstart.
The design recommended by the FDA is the full replicate design "2x2x4".
The sample size estimation is done only for balanced studies since the
break down of the total subject number in case of unbalanced sequence groups
is not unique. Moreover the formulas used are only valid for balanced designs.
The FDA method is described for the ABE limits 0.8 ... 1.25 only. Setting theta1, theta2
to other values may not be reasonable and is not tested.
The minimum sample size is 6, even if the power is higher than the intended targetpower.
The method is also required by China’s Center of Drug Evaluation.
D. Labes
Food and Drug Administration, Office of Generic Drugs (OGD). Draft Guidance on Dabigatran Etexilate Mesylate. Recommended Jun 2012; Revised Sep 2015, Jul 2017. download
Food and Drug Administration, Office of Generic Drugs (OGD). Draft Guidance on Rivaroxaban. Recommended Sep 2015. download
Food and Drug Administration, Office of Generic Drugs (OGD). Draft Guidance on Edoxaban Tosylate. Recommended May 2017; Revised Mar 2020. download
power.HVNTID
and power.NTIDFDA, sampleN.NTIDFDA for NTIDs with
low variability
# using all defaults but CV sampleN.HVNTID(CV = 0.3) # should give # n=22 with an (empirical) power of 0.829700 # Test formulation with lower variability but same pooled CV CVs <- CVp2CV(0.3, ratio = 0.25) sampleN.HVNTID(CV = CVs) # should give (no distinct difference to example above) # n=22 with an (empirical) power of 0.837520# using all defaults but CV sampleN.HVNTID(CV = 0.3) # should give # n=22 with an (empirical) power of 0.829700 # Test formulation with lower variability but same pooled CV CVs <- CVp2CV(0.3, ratio = 0.25) sampleN.HVNTID(CV = CVs) # should give (no distinct difference to example above) # n=22 with an (empirical) power of 0.837520
Function for estimating the sample size needed to have a pre-specified power for the one-sided non-inferiority t-test for normal or log-normal distributed data.
sampleN.noninf(alpha = 0.025, targetpower = 0.8, logscale = TRUE, margin,theta0, CV, design = "2x2", robust = FALSE, details = FALSE, print = TRUE, imax=100)sampleN.noninf(alpha = 0.025, targetpower = 0.8, logscale = TRUE, margin,theta0, CV, design = "2x2", robust = FALSE, details = FALSE, print = TRUE, imax=100)
alpha |
Significance level (one-sided). Defaults here to 0.025. |
targetpower |
Power to achieve at least. Must be >0 and <1. |
logscale |
Should the data used on log-transformed or on original scale? |
theta0 |
‘True’ or assumed T/R ratio or difference. |
margin |
Non-inferiority margin. |
CV |
In case of In case of cross-over studies this is the within-subject CV, in case of a parallel-group design the CV of the total variability. |
design |
Character string describing the study design. |
robust |
Defaults to FALSE. With that value the usual degrees of freedom will be used. |
details |
If |
print |
If |
imax |
Maximum number of steps in sample size search. |
The sample size is calculated via iterative evaluation of power.noninf().
Start value for the sample size search is taken from a large sample approximation.
The sample size is bound to 4 as minimum.
The estimated sample size gives always the total number of subjects (not subject/sequence in crossovers or subjects/group in parallel designs – like in some other software packages).
Notes on the underlying hypotheses
If the supplied margin is < 0 (logscale=FALSE) or < 1 (logscale=TRUE),
then it is assumed higher response values are better. The hypotheses are H0: theta0 <= margin vs. H1: theta0 > margin,
where theta0 = mean(test)-mean(reference) if logscale=FALSE
or H0: log(theta0) <= log(margin) vs. H1: log(theta0) > log(margin),
where theta0 = mean(test)/mean(reference) if logscale=TRUE.
If the supplied margin is > 0 (logscale=FALSE) or > 1 (logscale=TRUE),
then it is assumed lower response values are better. The hypotheses are H0: theta0 >= margin vs. H1: theta0 < margin
where theta0 = mean(test)-mean(reference) if logscale=FALSE
or H0: log(theta0) >= log(margin) vs. H1: log(theta0) < log(margin)
where theta0 = mean(test)/mean(reference) if logscale=TRUE.
This latter case may also be considered as ‘non-superiority’.
A data.frame with the input settings and results will be returned.
Explore it with str(sampleN.noninf(...)
The function does not vectorize properly.
If you need sample sizes with varying CVs, use f.i. for-loops or the apply-family.
D. Labes
Julious SA. Sample sizes for clinical trials with Normal data. Stat Med. 2004;23(12):1921–86. doi:10.1002/sim.1783
# using all the defaults: margin=0.8, theta0=0.95, alpha=0.025 # log-transformed, design="2x2" sampleN.noninf(CV = 0.3) # should give n=48 # # 'non-superiority' case, log-transformed data # with assumed 'true' ratio somewhat above 1 sampleN.noninf(CV = 0.3, targetpower = 0.9, margin = 1.25, theta0 = 1.05) # should give n=62# using all the defaults: margin=0.8, theta0=0.95, alpha=0.025 # log-transformed, design="2x2" sampleN.noninf(CV = 0.3) # should give n=48 # # 'non-superiority' case, log-transformed data # with assumed 'true' ratio somewhat above 1 sampleN.noninf(CV = 0.3, targetpower = 0.9, margin = 1.25, theta0 = 1.05) # should give n=62
This function performs the sample size estimation for the BE decision for the FDA’s method for NTIDs based on simulations. The study design is the full replicate design 2x2x4 (TRTR|RTRT) or the 3-period replicate design with sequences TRT|RTR.
sampleN.NTID(alpha = 0.05, targetpower = 0.8, theta0, theta1, theta2, CV, design = c("2x2x4", "2x2x3"), nsims = 1e+05, nstart, imax = 100, print = TRUE, details = TRUE, setseed = TRUE)sampleN.NTID(alpha = 0.05, targetpower = 0.8, theta0, theta1, theta2, CV, design = c("2x2x4", "2x2x3"), nsims = 1e+05, nstart, imax = 100, print = TRUE, details = TRUE, setseed = TRUE)
alpha |
Type I error probability. Per convention mostly set to 0.05. |
targetpower |
Power to achieve at least. Must be >0 and <1. |
theta0 |
‘True’ or assumed T/R ratio. |
theta1 |
Conventional lower ABE limit to be applied in the FDA procedure. |
theta2 |
Conventional upper ABE limit to be applied in the FDA procedure. |
CV |
Intra-subject coefficient(s) of variation as ratio (not percent).
|
design |
Design of the study to be planned. |
nsims |
Number of simulations to be performed to obtain the empirical power. Defaults to 100,000 = 1e+5. |
nstart |
Set this to a start value for the sample size if a previous run failed. |
imax |
Maximum number of steps in sample size search. Defaults to 100. |
print |
If |
details |
If set to |
setseed |
Simulations are dependent on the starting point of the (pseudo) random number
generator. To avoid differences in power values for different runs a
|
The linearized scaled ABE criterion is calculated according to the SAS code
given in the FDA’s guidances. For deciding BE the study must pass that criterion,
the conventional ABE test, and that the upper confidence limit of
σwT/σwR ≤ 2.5.
The simulations are done via the distributional properties of the statistical
quantities necessary for deciding BE based on this method.
Details can be found in a document Implementation_scaledABE_sims located in
the /doc sub-directory of the package.
The estimated sample size gives always the total number of subjects (not subject / sequence – like in some other software packages).
Returns a data frame with the input settings and sample size results.
The Sample size column contains the total sample size.
The nlast column contains the last n value. May be useful for re-starting.
For some input combinations the sample size search may be very time
consuming and will eventually even fail since the start values chosen may
not really reasonable for them. This applies especially for theta0 values
close to the implied scaled limits according to
exp(∓1.053605 × swR).
In case of a failed sample size search you may restart with setting the argument
nstart.
In case of theta0 values outside the implied scaled (tightened/widened) ABE limits
no sample size estimation is possible and the function throws an error
(f.i. CV = 0.04, theta0 = 0.95).
The design recommended by the FDA is the full replicate design "2x2x4".
The sample size estimation is done only for balanced studies since the
break down of the total subject number in case of unbalanced sequence groups
is not unique. Moreover the formulas used are only valid for balanced designs.
The FDA method is described for the ABE limits 0.8 ... 1.25 only. Setting theta1,
theta2 to other values may not be reasonable and is not tested.
The results for the design "2x2x3" are to be considered as experimental since
at present not thorougly tested.
The minimum sample size is 6, even if the power is higher than the intended
targetpower.
The method is also required by China’s Center of Drug Evaluation.
D. Labes
Food and Drug Administration, Office of Generic Drugs (OGD). Draft Guidance on Warfarin Sodium. Recommended Dec 2012. download
Food and Drug Administration, Center for Drug Evaluation and Research (CDER). Draft Guidance for Industry. Bioequivalence Studies with Pharmacokinetic Endpoints for Drugs Submitted Under an ANDA. August 2021. download
Yu LX, Jiang W, Zhang X, Lionberger R, Makhlouf F, Schuirmann DJ, Muldowney L, Chen ML, Davit B, Conner D, Woodcock J. Novel bioequivalence approach for narrow therapeutic index drugs. Clin Pharmacol Ther. 2015;97(3):286–91. doi:10.1002/cpt.28
Jiang W, Makhlouf F, Schuirmann DJ, Zhang X, Zheng N, Conner D, Yu LX, Lionberger R. A Bioequivalence Approach for Generic Narrow Therapeutic Index Drugs: Evaluation of the Reference-Scaled Approach and Variability Comparison Criterion. AAPS J. 2015;17(4):891–901. doi:10.1208/s12248-015-9753-5
Endrényi L, Tóthfalusi L. Determination of Bioequivalence for Drugs with Narrow Therapeutic Index: Reduction of the Regulatory Burden. J Pharm Pharm Sci. 2013;16(5):676–82. open access
power.NTID
and power.HVNTID, sampleN.HVNTID for NTIDs with
high variability
sampleN.NTID(CV = 0.04,theta0 = 0.975) # should give # n=54 with an (empirical) power of 0.809590 # # Test formulation with lower variability sampleN.NTID(CV = c(0.04,0.06),theta0 = 0.975) # should give # n=20 with an (empirical) power of 0.0.814610 # # alternative 3-period design sampleN.NTID(CV = 0.04,theta0 = 0.975, design="2x2x3") # should give # n=86 with power = 0.80364sampleN.NTID(CV = 0.04,theta0 = 0.975) # should give # n=54 with an (empirical) power of 0.809590 # # Test formulation with lower variability sampleN.NTID(CV = c(0.04,0.06),theta0 = 0.975) # should give # n=20 with an (empirical) power of 0.0.814610 # # alternative 3-period design sampleN.NTID(CV = 0.04,theta0 = 0.975, design="2x2x3") # should give # n=86 with power = 0.80364
Estimates the necessary sample size to have at least a given power based on Fieller’s confidence (‘fiducial’) interval.
sampleN.RatioF(alpha = 0.025, targetpower = 0.8, theta1 = 0.8, theta2, theta0 = 0.95, CV, CVb, design = "2x2", print = TRUE, details = FALSE, imax=100, setseed=TRUE)sampleN.RatioF(alpha = 0.025, targetpower = 0.8, theta1 = 0.8, theta2, theta0 = 0.95, CV, CVb, design = "2x2", print = TRUE, details = FALSE, imax=100, setseed=TRUE)
alpha |
Type I error probability. |
targetpower |
Power to achieve at least. Must be >0 and <1. Typical values are 0.8 or 0.9. |
theta1 |
Lower bioequivalence limit. Typically 0.8 (default). |
theta2 |
Upper bioequivalence limit. Typically 1.25. |
theta0 |
‘True’ or assumed T/R ratio. Typically set to 0.95. |
CV |
Coefficient of variation as ratio. In case of |
CVb |
CV of the between-subject variability. Only necessary for |
design |
A character string describing the study design. |
print |
If |
details |
If |
imax |
Maximum number of steps in sample size search. |
setseed |
If set to |
The sample size is based on exact power calculated using the bivariate
non-central t-distribution via function pmvt of the package mvtnorm.
Due to the calculation method used in package mvtnorm these
probabilities are dependent from the state of the random number generator
within the precision of the power.
The CV(within) and CVb(etween) in case of design="2x2" are obtained
via an appropriate ANOVA from the error term and from the difference
(MS(subject within sequence)-MS(error))/2.
The estimated sample size gives always the total number of subjects (not subject/sequence in crossovers or subjects/group in parallel designs – like in some other software packages).
A data.frame with the input values and results will be returned.
The sample size n returned is the total sample size for both designs.
This function is intended for studies with clinical endpoints.
In such studies
the 95% confidence intervals are usually used for equivalence testing.
Therefore, alpha defaults here to 0.025 (see EMEA 2000).
D. Labes
Fieller EC. Some Problems in Interval Estimation. J Royal Stat Soc B. 1954;16(2):175–85. doi:10.1111/j.2517-6161.1954.tb00159.x
Sasabuchi S. A test of a multivariate normal mean with composite hypotheses determined by linear inequalities. Biometrika. 1980;67(2):429–39. doi:10.1093/biomet/67.2.429
Hauschke D, Kieser M, Diletti E, Burke M. Sample size determination for proving equivalence based on the ratio of two means for normally distributed data. Stat Med. 1999;18(1):93–105.
Hauschke D, Steinijans V, Pigeot I. Bioequivalence Studies in Drug Development. Chichester: Wiley; 2007. Chapter 10.
European Agency for the Evaluation of Medicinal Products, CPMP. Points to Consider on Switching between Superiority and Non-Inferiority. London, 27 July 2000. CPMP/EWP/482/99
# sample size for a 2x2 cross-over study # with CVw=0.2, CVb=0.4 # alpha=0.025 (95% CIs), target power = 80% # 'true' ratio = 95%, BE acceptance limits 80-125% # using all the defaults: sampleN.RatioF(CV = 0.2, CVb = 0.4) # gives n=28 with an achieved power of 0.807774 # see Hauschke et.al. (2007) Table 10.3a # sample size for a 2-group parallel study # with CV=0.4 (total variability) # alpha=0.025 (95% CIs), target power = 90% # 'true' ratio = 90%, BE acceptance limits 75-133.33% sampleN.RatioF(targetpower = 0.9, theta1 = 0.75, theta0 = 0.90, CV = 0.4, design = "parallel") # gives n=236 with an achieved power of 0.900685 # see Hauschke et.al. (2007) Table 10.2 # a rather strange setting of ratio0! have a look at n. # it would be better this is not the sample size but your account balance ;-). sampleN.RatioF(theta0 = 0.801, CV = 0.2, CVb = 0.4)# sample size for a 2x2 cross-over study # with CVw=0.2, CVb=0.4 # alpha=0.025 (95% CIs), target power = 80% # 'true' ratio = 95%, BE acceptance limits 80-125% # using all the defaults: sampleN.RatioF(CV = 0.2, CVb = 0.4) # gives n=28 with an achieved power of 0.807774 # see Hauschke et.al. (2007) Table 10.3a # sample size for a 2-group parallel study # with CV=0.4 (total variability) # alpha=0.025 (95% CIs), target power = 90% # 'true' ratio = 90%, BE acceptance limits 75-133.33% sampleN.RatioF(targetpower = 0.9, theta1 = 0.75, theta0 = 0.90, CV = 0.4, design = "parallel") # gives n=236 with an achieved power of 0.900685 # see Hauschke et.al. (2007) Table 10.2 # a rather strange setting of ratio0! have a look at n. # it would be better this is not the sample size but your account balance ;-). sampleN.RatioF(theta0 = 0.801, CV = 0.2, CVb = 0.4)
This function performs the sample size estimation for the BE decision via linearized scaled ABE criterion based on simulations.
sampleN.RSABE(alpha = 0.05, targetpower = 0.8, theta0, theta1, theta2, CV, design = c("2x3x3", "2x2x4", "2x2x3"), regulator = c("FDA", "EMA"), nsims = 1e+05, nstart, imax=100, print = TRUE, details = TRUE, setseed=TRUE)sampleN.RSABE(alpha = 0.05, targetpower = 0.8, theta0, theta1, theta2, CV, design = c("2x3x3", "2x2x4", "2x2x3"), regulator = c("FDA", "EMA"), nsims = 1e+05, nstart, imax=100, print = TRUE, details = TRUE, setseed=TRUE)
alpha |
Type I error probability. Per convention mostly set to 0.05. |
targetpower |
Power to achieve at least. Must be >0 and <1. |
theta0 |
‘True’ or assumed T/R ratio. |
theta1 |
Conventional lower ABE limit to be applied in the mixed procedure if
|
theta2 |
Conventional upper ABE limit to be applied in the mixed procedure if
|
CV |
Intra-subject coefficient(s) of variation as ratio (not percent).
|
design |
Design of the study to be planned. |
regulator |
Regulatory body settings for the scaled ABE criterion. |
nsims |
Number of simulations to be performed to obtain the (empirical) power. |
nstart |
Set this to a start for the sample size search if a previous run failed. |
imax |
Maximum number of steps in sample size search. Defaults to 100. |
print |
If |
details |
If set to |
setseed |
Simulations are dependent on the starting point of the (pseudo) random number
generator. To avoid differences in power for different runs a |
The linearized scaled ABE criterion is calculated according to the SAS code
given in the FDA progesterone guidance.
The simulations are done via the distributional properties of the statistical
quantities necessary for deciding BE based on scaled ABE.
For more details see a document Implementation_scaledABE_simsVx.yy.pdf in the /doc sub-directory of
the package.
If a CVcap is defined for the regulator, the BE decision is based on the inclusion
of the CI in the capped widened acceptance limits in case of CVwR > CVcap. This
resembles method ‘Howe-EMA’ in Muñoz et al. and is the standard behavior now if
regulator="EMA" is choosen.
The estimated sample size gives always the total number of subjects (not subject/sequence – like in some other software packages).
Returns a data.frame with the input and sample size results.
The Sample size column contains the total sample size.
The nlast column contains the last n value. May be useful for restarting.
Although some designs are more ‘popular’ than others, sample size estimations are valid for all of the following designs:
"2x2x4" |
TRTR | RTRT |
| TRRT | RTTR | |
| TTRR | RRTT | |
"2x2x3" |
TRT | RTR |
| TRR | RTT | |
"2x3x3" |
TRR | RTR | RRT |
The sample size estimation for theta0 >1.2 and <0.85 may be very time consuming
and will eventually also fail since the start values chosen are not really
reasonable in that ranges. This is especially true in the range about CV = 0.3
and regulatory constant according to FDA.
If you really need sample sizes in that range be prepared to restart the sample
size estimation via the argument nstart.
Since the dependence of power from n is very flat in the mentioned region you
may also consider to adapt the number of simulations not to tap in the simulation
error trap.
The sample size estimation is done only for balanced designs since the
break down of the total subject number in case of unbalanced sequence groups
is not unique. Moreover the formulas used are only for balanced designs.
The minimum sample size is n=6, even if the power is higher than the intended
targetpower.
D. Labes
Food and Drug Administration, Office of Generic Drugs (OGD). Draft Guidance on Progesterone. Recommended Apr 2010. Revised Feb 2011. download
Tóthfalusi, L, Endrényi, L. Sample Sizes for Designing Bioequivalence Studies for Highly Variable Drugs. J Pharm Pharmaceut Sci. 2011;15(1):73–84. open access
Tóthfalusi L, Endrényi L, García Arieta A. Evaluation of Bioequivalence for Highly Variable Drugs with Scaled Average Bioequivalence. Clin Pharmacokin. 2009;48(11):725–43. doi:10.2165/11318040-000000000-00000
Muñoz J, Alcaide D, Ocaña J. Consumer’s risk in the EMA and FDA regulatory approaches for bioequivalence in highly variable drugs. Stat Med. 2015;35(12):1933–43. doi:10.1002/sim.6834
# using all the defaults: # design=2x3x3 (partial replicate design), theta0=0.90, # ABE limits, PE constraint 0.8 - 1.25 # targetpower=80%, alpha=0.05, 1E5 simulations sampleN.RSABE(CV = 0.3) # should result in a sample size n=45, power=0.80344# using all the defaults: # design=2x3x3 (partial replicate design), theta0=0.90, # ABE limits, PE constraint 0.8 - 1.25 # targetpower=80%, alpha=0.05, 1E5 simulations sampleN.RSABE(CV = 0.3) # should result in a sample size n=45, power=0.80344
These functions performs the sample size estimation of the BE decision via
the reference scaled ABE based on subject data simulations.
Implemented are the methods ABEL, Hyslop and ‘exact’ (see the References in
power.RSABE2L.sdsims).
The estimation method of the key statistics needed to perform the RSABE decision
is the usual ANOVA.
This function has an alias sampleN.RSABE2L.sds().
sampleN.RSABE2L.sdsims(alpha = 0.05, targetpower = 0.8, theta0, theta1, theta2, CV, design = c("2x3x3", "2x2x4", "2x2x3"), SABE_test = "exact", regulator, nsims=1e5, nstart, imax = 100, print = TRUE, details = TRUE, setseed = TRUE, progress)sampleN.RSABE2L.sdsims(alpha = 0.05, targetpower = 0.8, theta0, theta1, theta2, CV, design = c("2x3x3", "2x2x4", "2x2x3"), SABE_test = "exact", regulator, nsims=1e5, nstart, imax = 100, print = TRUE, details = TRUE, setseed = TRUE, progress)
alpha |
Type I error probability. Per convention mostly set to 0.05. |
targetpower |
Power to achieve at least. Must be >0 and <1. |
theta0 |
‘True’ or assumed T/R ratio. |
theta1 |
Conventional lower ABE limit to be applied in the mixed procedure if
|
theta2 |
Conventional upper ABE limit to be applied in the mixed procedure if
|
CV |
Intra-subject coefficient(s) of variation as ratio (not percent).
|
design |
Design of the study to be planned. |
SABE_test |
This argument specifies the test method to be used for the reference scaled
ABE decision. |
regulator |
Regulatory settings for the widening of the BE acceptance limits. |
nsims |
Number of simulations to be performed to obtain the (empirical) power. The default value 100,000 = 1e+5 is usually sufficient. Consider to rise this value if theta0<=0.85 or >=1.25. But see the warning section. |
nstart |
Set this to a start for the sample size search if a previous run failed. |
imax |
Maximum number of steps in sample size search. Defaults to 100. |
print |
If |
details |
If set to |
setseed |
Simulations are dependent on the starting point of the (pseudo) random number
generator. To avoid differences in power for different runs a |
progress |
Should a progressbar be shown? Defaults to |
The methods rely on the analysis of log-transformed data, i.e., assume a
log-normal distribution on the original scale.
The widened BE acceptance limits will be calculated by the formula [L, U] = exp(± r_const * sWR)
with r_const the regulatory constant and sWR the standard deviation of the within
subjects variability of the Reference. r_const = 0.76 (~log(1.25)/0.29356) is used
in case of regulator="EMA".
If the CVwR of the Reference is < CVswitch=0.3 the conventional ABE limits
apply (mixed procedure).
In case of regulator="EMA" a cap is placed on the widened limits if
CVwr>0.5, i.e., the widened limits are held at value calculated for CVwR=0.5.
The simulations are done by simulating subject data (all effects fixed except the
residuals) and evaluating these data via ANOVA of all data to get the point estimate
of T vs. R along with its 90% CI and an ANOVA of the data under R(eference) only
to get an estimate of s2wR.
The estimated sample size gives always the total number of subjects (not subject/sequence – like in some other software packages).
Returns a data.frame with the input settings and sample size results.
The Sample size column contains the total sample size.
The nlast column contains the last n value. May be useful for restarting.
Although some designs are more ‘popular’ than others, sample size estimations are valid for all of the following designs:
"2x2x4" |
TRTR | RTRT |
| TRRT | RTTR | |
| TTRR | RRTT | |
"2x2x3" |
TRT | RTR |
| TRR | RTT | |
"2x3x3" |
TRR | RTR | RRT |
The sample size estimation for very extreme theta0 (<0.83 or >1.21) may be very
time consuming and will eventually also fail since the start values chosen are
not really reasonable in that ranges.
If you really need sample sizes in that range be prepared to restart the sample
size estimation via the argument nstart.
Since the dependence of power from n is very flat in the mentioned region you may
also consider to adapt the number of simulations not to tap in the simulation
error trap.
We are doing the sample size estimation only for balanced designs since the
break down of the total subject number in case of unbalanced sequence groups
is not unique. Moreover the formulas used are only for balanced designs.
The minimum sample size is 6, even if the power is higher than the intended
targetpower.
Subject simulations are relatively slow. Thus be patient and go for a cup of coffee if you use this function with high sample sizes!
H. Schütz
power.RSABE2L.sdsims, sampleN.scABEL, reg_const
# using the defaults: # partial replicate design, targetpower=80%, # true assumed ratio = 0.90, 1E+5 simulated studies # ABE limits, PE constraint 0.8 - 1.25 # EMA regulatory settings # compare results CV <- 0.4 method <- c("exact", "abel", "hyslop", "fda") res <- data.frame(SABE_test = c("ncTOST", "ABEL", "Hyslop", "FDA"), n = NA, power = NA) for (i in 1:nrow(res)) { res[i, 2:3] <- sampleN.RSABE2L.sdsims(CV = CV, SABE_test = method[i], details = FALSE, print = FALSE)[8:9] } print(res, digits = 4, row.names = FALSE) # should result in a sample size n=48 with all methods, # power=0.8197 (ncTOST), 0.8411 (ABEL), 0.8089 (Hyslop), 0.8113 (FDA)# using the defaults: # partial replicate design, targetpower=80%, # true assumed ratio = 0.90, 1E+5 simulated studies # ABE limits, PE constraint 0.8 - 1.25 # EMA regulatory settings # compare results CV <- 0.4 method <- c("exact", "abel", "hyslop", "fda") res <- data.frame(SABE_test = c("ncTOST", "ABEL", "Hyslop", "FDA"), n = NA, power = NA) for (i in 1:nrow(res)) { res[i, 2:3] <- sampleN.RSABE2L.sdsims(CV = CV, SABE_test = method[i], details = FALSE, print = FALSE)[8:9] } print(res, digits = 4, row.names = FALSE) # should result in a sample size n=48 with all methods, # power=0.8197 (ncTOST), 0.8411 (ABEL), 0.8089 (Hyslop), 0.8113 (FDA)
This function performs the sample size estimation via power calculations of the BE decision via scaled (expanded) BE acceptance limits, based on simulations.
sampleN.scABEL(alpha = 0.05, targetpower = 0.8, theta0, theta1, theta2, CV, design = c("2x3x3", "2x2x4", "2x2x3"), regulator, nsims = 1e+05, nstart, imax = 100, print = TRUE, details = TRUE, setseed = TRUE)sampleN.scABEL(alpha = 0.05, targetpower = 0.8, theta0, theta1, theta2, CV, design = c("2x3x3", "2x2x4", "2x2x3"), regulator, nsims = 1e+05, nstart, imax = 100, print = TRUE, details = TRUE, setseed = TRUE)
alpha |
Type I error probability. Per convention mostly set to 0.05. |
targetpower |
Power to achieve at least. Must be >0 and <1. |
theta0 |
‘True’ or assumed T/R ratio. |
theta1 |
Conventional lower ABE limit to be applied in the mixed procedure if
|
theta2 |
Conventional upper ABE limit to be applied in the mixed procedure if
|
CV |
Intra-subject coefficient(s) of variation as ratio (not percent).
|
design |
Design of the study to be planned. |
regulator |
Regulatory settings for the widening of the BE acceptance limits. |
nsims |
Number of simulations to be performed to obtain the (empirical) power.
The default value 100,000 = 1e+5 is usually sufficient. Consider to rise
this value if |
nstart |
Set this to a start for the sample size search if a previous run failed. |
imax |
Maximum number of steps in sample size search. Defaults to 100. |
print |
If |
details |
If set to |
setseed |
Simulations are dependent on the starting point of the (pseudo) random number
generator. To avoid differences in power for different runs a |
The simulations are done via the distributional properties of the statistical
quantities necessary for deciding BE based on ABEL (‘Average Bioequivalence with Expanded Limits’). For more details see a description in the /doc sub-directory of the package.
Function sampleN.scABEL() is based on power calculations via simulations
using the distributional characteristics of the ‘key’ statistics obtained from
the EMA recommended evaluation via ANOVA if regulator="EMA" or if the
regulator component est_method is set to "ANOVA" if regulator is an object
of class 'regSet'.
Otherwise, the simulations are based on the distributional characteristis of the
‘key’ statistics obtained from evaluation via intra-subject contrasts (ISC),
as recommended by the FDA.
The estimated sample size gives always the total number of subjects
(not subject/sequence – like in some other software packages).
Function sampleN.scABEL2() is solely based on power calculations via
simulation using the distributional characteristics of the ‘key’ statistics
obtained from evaluation via intra-subject contrasts (ISC). This function is deprecated.
Returns a data.frame with the input settings and sample size results.
The Sample size column contains the total sample size.
The nlast column contains the last n value. May be useful for restarting.
Although some designs are more ‘popular’ than others, sample size estimations are valid for all of the following designs:
"2x2x4" |
TRTR | RTRT |
| TRRT | RTTR | |
| TTRR | RRTT | |
"2x2x3" |
TRT | RTR |
| TRR | RTT | |
"2x3x3" |
TRR | RTR | RRT |
The sample size estimation for extreme theta0 (<0.83 or >1.21) may be very
time consuming and will eventually also fail since the start values chosen are
not really reasonable in that ranges. This is especially true in the range around
CV = 0.3 and regulatory constant according to FDA.
If you really need sample sizes in that range be prepared to restart the sample
size estimation via the argument nstart.
Since the dependence of power from n is very flat in the mentioned region you may
also consider to adapt the number of simulations not to get caught in the simulation
error trap.
If results of power.scABEL are expected to be inaccurate (partial
replicate design with unbalanced sequences and/or heteroscedasticity in the case of CVwT > CVwR, subject data via sampleN.scABEL.sdsims should be simulated instead. Very time consuming (easily 100times slower)! Subject data simulations are only supported for regulator="EMA" and regulator="GCC".
We are doing the sample size estimation only for balanced designs since the
break down of the total subject number in case of unbalanced sequence groups
is not unique. Moreover the formulas used are only for balanced designs.
In case of regulator="FDA" the sample size is only approximate since
the BE decision method is not exactly what is expected by the FDA. But the two Lászlós state that the scABEL method should be ‘operationally’ equivalent to the
FDA method. Thus the sample size should be comparable.
Consider in case of regulator="FDA" to use the function
sampleN.RSABE() instead.
In case of regulator="HC" the underlying power is only approximative
since the Health Canada recommends evaluation by a mixed model approach.
But this could only implemented via subject data simulations which are very
time consuming.
The minimum sample size is 6, even if the power is higher than the intended
targetpower.
D. Labes
Tóthfalusi L, Endrényi L. Sample Sizes for Designing Bioequivalence Studies for Highly Variable Drugs. J Pharm Pharmaceut Sci. 2011;15(1):73–84. open access
power.scABEL, sampleN.scABEL.sdsims, sampleN.RSABE,
reg_const
# using all the defaults: # partial replicate design, targetpower=80%, # true assumed ratio = 0.90, 1E+5 simulated studies # ABE limits, PE constraint 0.8 - 1.25 # EMA regulatory settings sampleN.scABEL(CV = 0.3) # should result in a sample size n=54, power=0.8159 # Now with former (inofficial) ANVISA settings, CVswitch=40% # (since 2016 ANVISA uses the same settings as EMA) reg <- reg_const("USER", r_const = 0.76, CVswitch = 0.4, CVcap = 0.5) reg$name <- "Old ANVISA" sampleN.scABEL(CV = 0.3, regulator = reg) # should result in a sample size n=60, power=0.8101 # For the full replicate design, target power = 90% # true assumed ratio = 0.9, FDA regulatory settings # sims based on evalaution via ISC sampleN.scABEL(CV = 0.4, targetpower = 0.9, design = "2x2x4", regulator = "FDA") # should result in a sample size n=32, power=0.9125 # Fixed wider limits (0.7500 - 1.3333) for the GCC sampleN.scABEL(CV = 0.4, targetpower = 0.9, design = "2x2x4", regulator = "GCC") # should result in a sample size n=40, power=0.9039# using all the defaults: # partial replicate design, targetpower=80%, # true assumed ratio = 0.90, 1E+5 simulated studies # ABE limits, PE constraint 0.8 - 1.25 # EMA regulatory settings sampleN.scABEL(CV = 0.3) # should result in a sample size n=54, power=0.8159 # Now with former (inofficial) ANVISA settings, CVswitch=40% # (since 2016 ANVISA uses the same settings as EMA) reg <- reg_const("USER", r_const = 0.76, CVswitch = 0.4, CVcap = 0.5) reg$name <- "Old ANVISA" sampleN.scABEL(CV = 0.3, regulator = reg) # should result in a sample size n=60, power=0.8101 # For the full replicate design, target power = 90% # true assumed ratio = 0.9, FDA regulatory settings # sims based on evalaution via ISC sampleN.scABEL(CV = 0.4, targetpower = 0.9, design = "2x2x4", regulator = "FDA") # should result in a sample size n=32, power=0.9125 # Fixed wider limits (0.7500 - 1.3333) for the GCC sampleN.scABEL(CV = 0.4, targetpower = 0.9, design = "2x2x4", regulator = "GCC") # should result in a sample size n=40, power=0.9039
This function performs a sample size estimation for the BE decision via Average Bioequivalenc with Expanding Limits (ABEL) based on simulations. Simultaneously alpha is iteratively adjusted in order to maintain the consumer risk at the nominal level.
sampleN.scABEL.ad(alpha = 0.05, targetpower = 0.8, theta0, theta1, theta2, CV, design = c("2x3x3", "2x2x4", "2x2x3"), regulator, nstart = NA, nsims = 1e+06, imax = 100, tol, print = TRUE, details = FALSE, alpha.pre = 0.05, setseed = TRUE, sdsims = FALSE, progress)sampleN.scABEL.ad(alpha = 0.05, targetpower = 0.8, theta0, theta1, theta2, CV, design = c("2x3x3", "2x2x4", "2x2x3"), regulator, nstart = NA, nsims = 1e+06, imax = 100, tol, print = TRUE, details = FALSE, alpha.pre = 0.05, setseed = TRUE, sdsims = FALSE, progress)
alpha |
Type I error (TIE) probability (nominal level of the test). Per convention commonly set to 0.05. |
targetpower |
Power to achieve at least. Must be >0 and <1. Typical values are 0.80 to 0.90 (i.e., 80% to 90%). Defaults to 0.80 if not given explicitly. |
theta0 |
‘True’ or assumed T/R ratio. Defaults to 0.90 according to the two Lászlós if not given explicitly. |
theta1 |
Conventional lower ABE limit to be applied in the mixed procedure
if |
theta2 |
Conventional upper ABE limit to be applied in the mixed procedure
if |
CV |
Intra-subject coefficient(s) of variation as ratio (not percent).
|
design |
Design of the study to be planned. |
regulator |
Regulatory settings for the widening of the BE acceptance limits.
Choose from |
nstart |
Best “guess” sample size. If not given (default), simulations
start with the sample size estimated for |
nsims |
Number of simulations to be performed to estimate the (empirical) TIE and in each iteration of adjusting alpha. The default value 1,000,000 = 1E+6 should not be lowered. |
imax |
Maximum number of steps in sample size search. Defaults to 100. |
tol |
Desired accuracy (convergence tolerance). Defaults to 1E-6. |
print |
If |
details |
If |
alpha.pre |
Pre-specified alpha (optional). Must be |
setseed |
Simulations are dependent on the starting point of the (pseudo)
random number generator. To avoid differences in power for different
runs |
sdsims |
If |
progress |
Set to |
The simulations are done via the distributional properties of the statistical
quantities necessary for assessing BE based on ABEL.
Simulations of the TIE are performed at the upper (expanded) limit U
of the acceptance range. Due to the symmetry around 1 results are valid for the lower
(expanded) limit L as well.
U at the EMA’s and Health Canada’s CVcap, the GCC’s for any CVwR > 0.30:
scABEL(CV = 0.5, reg = "EMA")[["upper"]] [1] 1.43191 scABEL(CV = 0.57382, reg = "HC")[["upper"]] [1] 1.5 scABEL(CV = 0.31, reg = "GCC")[["upper"]] [1] 1.333333
Simulated studies are evaluated by ANOVA (Method A) as recommended in the
EMA’s Q&A-document and by intra-subject contrasts if regulator="HC".
Health Canada requires a mixed effects model which cannot be implemented in R. However,
intra-subjects contrasts are a sufficiently close approximation.
If an inflation of the TIE is expected (i.e., >alpha), alpha is
iteratively adjusted until at least the target power is reached and the consumer
risk is maintained (<=alpha). For details about the algorithm see the
respective section of scABEL.ad.
The estimated sample size gives always the total number of subjects (not subject/sequence – like in some other software packages).
Returns a data.frame with the input and results for adjusted alpha,
type I error, sample size, and achieved power.
The Sample size column contains the total sample size.
If no adjustment is necessary, NA will be returned in the
alpha.adj column and other results are identical to the ones
obtained by sampleN.scABEL.
Although some designs are more ‘popular’ than others, sample size estimations are valid for all of the following designs:
"2x2x4" |
TRTR | RTRT |
| TRRT | RTTR | |
| TTRR | RRTT | |
"2x2x3" |
TRT | RTR |
| TRR | RTT | |
"2x3x3" |
TRR | RTR | RRT |
The sample size estimation for extreme theta0 (<0.83 or >1.21) may be time
consuming and will eventually also fail since the start values chosen are
not really reasonable in that ranges.
If you really need sample sizes in that range be prepared to restart the sample
size estimation with nstart above the last one before failure.
Since the dependence of power from n is very flat in the mentioned region you may
also consider to adapt the number of simulations not to tap in the simulation
error trap.
See also the Warning section of the function power.scABEL concerning
the power value agreement to those obtained from simulations via subject data.
For the GCC and CVwR ≤ 0.30 simulations will be time consuming and may result in large sample sizes.
We are doing the sample size estimation only for balanced designs since the break down of the total subject number in case of unbalanced sequences is not unique. Moreover the formulas used are only for balanced designs.
H. Schütz
Tóthfalusi L, Endrényi L. Sample Sizes for Designing Bioequivalence Studies for Highly Variable Drugs. J Pharm Pharmaceut Sci. 2011;15(1):73–84. open access
Wonnemann M, Frömke C, Koch A. Inflation of the Type I Error: Investigations on Regulatory Recommendations for Bioequivalence of Highly Variable Drugs. Pharm Res. 2015;32(1):135–43. doi:10.1007/s11095-014-1450-z
Muñoz J, Alcaide D, Ocaña J. Consumer’s risk in the EMA and FDA regulatory approaches for bioequivalence in highly variable drugs. Stat Med. 2015;35(12):1933–43. doi:10.1002/sim.6834
Labes D, Schütz H. Inflation of Type I Error in the Evaluation of Scaled Average Bioequivalence, and a Method for its Control. Pharm Res. 2016;33(11):2805–14. doi:10.1007/s11095-016-2006-1
scABEL.ad, sampleN.scABEL, power.scABEL,
scABEL
# --- Not run due to timing policy of CRAN for examples # each may run some ten seconds or more # using all the defaults: # TRR|RTR|RRT, target power 80%, assumed ratio 0.90, 1E+6 simulated studies, # EMA regulatory settings (ABE limits, PE constraint 0.8 - 1.25) sampleN.scABEL.ad(CV = 0.3) # should result in n 60, power 0.8022. # Note: Without adjustment by sampleN.scABEL(): n 54, power 0.8159 # Easier to show the details: sampleN.scABEL.ad(CV = 0.3, details = TRUE) # # TRTR|RTRT, target power 90%, pre-specified alpha 0.025 sampleN.scABEL.ad(CV = 0.3, targetpower = 0.9, design = "2x2x4", alpha.pre = 0.025) # should result in n 60, power 0.9021; pre-specified alpha justified.# --- Not run due to timing policy of CRAN for examples # each may run some ten seconds or more # using all the defaults: # TRR|RTR|RRT, target power 80%, assumed ratio 0.90, 1E+6 simulated studies, # EMA regulatory settings (ABE limits, PE constraint 0.8 - 1.25) sampleN.scABEL.ad(CV = 0.3) # should result in n 60, power 0.8022. # Note: Without adjustment by sampleN.scABEL(): n 54, power 0.8159 # Easier to show the details: sampleN.scABEL.ad(CV = 0.3, details = TRUE) # # TRTR|RTRT, target power 90%, pre-specified alpha 0.025 sampleN.scABEL.ad(CV = 0.3, targetpower = 0.9, design = "2x2x4", alpha.pre = 0.025) # should result in n 60, power 0.9021; pre-specified alpha justified.
These functions performs the sample size estimation via power calculations of the BE decision via scaled (expanded) BE acceptance limits, based on subject data simulations.
This function has an alias sampleN.scABEL.sds().
sampleN.scABEL.sdsims(alpha = 0.05, targetpower = 0.8, theta0, theta1, theta2, CV, design = c("2x3x3", "2x2x4", "2x2x3"), regulator, nsims = 1e5, nstart, imax = 100, print = TRUE, details = TRUE, setseed = TRUE, progress)sampleN.scABEL.sdsims(alpha = 0.05, targetpower = 0.8, theta0, theta1, theta2, CV, design = c("2x3x3", "2x2x4", "2x2x3"), regulator, nsims = 1e5, nstart, imax = 100, print = TRUE, details = TRUE, setseed = TRUE, progress)
alpha |
Type I error probability. Per convention mostly set to 0.05. |
targetpower |
Power to achieve at least. Must be >0 and <1. |
theta0 |
‘True’ or assumed T/R ratio. |
theta1 |
Conventional lower ABE limit to be applied in the mixed procedure if
|
theta2 |
Conventional upper ABE limit to be applied in the mixed procedure if
|
CV |
Intra-subject coefficient(s) of variation as ratio (not percent).
|
design |
Design of the study to be planned. |
regulator |
Regulatory settings for the widening of the BE acceptance limits. |
nsims |
Number of simulations to be performed to obtain the (empirical) power.
The default value 100,000 = 1e+5 is usually sufficient. Consider to rise
this value if |
nstart |
Set this to a start for the sample size search if a previous run failed. |
imax |
Maximum number of steps in sample size search. Defaults to 100. |
print |
If |
details |
If set to |
setseed |
Simulations are dependent on the starting point of the (pseudo) random number generator. To avoid differences in power for different runs a |
progress |
Should a progressbar be shown? Defaults to |
The methods rely on the analysis of log-transformed data, i.e., assume a
log-normal distribution on the original scale.
The expanded BE acceptance limits will be calculated by the formula [L, U] = exp(± r_const * sWR)
with r_const the regulatory constant and sWR the standard deviation of the within
subjects variability of the Reference. r_const = 0.76 (~log(1.25)/0.29356) is used
in case of regulator = "EMA".
If the CVwR is < CVswitch=0.3 the conventional ABE limits apply (mixed procedure).
In case of regulator="EMA" a cap is placed on the widened limits if
CVwR > 0.50, i.e., the widened limits are held at value calculated for CVwR = 0.50.
In case of regulator="GCC" fixed wider limits of 0.7500 – 1.3333 for CVwR > 0.30 are applied and the conventional limits otherwise.
The simulations are done by simulating subject data (all effects fixed except the
residuals) and evaluating these data via ANOVA of all data to get the point estimate
of T vs. R along with its 90% CI and an ANOVA of the data under R(eference) only
to get an estimate of s²wR.
The estimated sample size gives always the total number of subjects (not subject/sequence – like in some other software packages).
Returns a data.frame with the input settings and sample size results.
The Sample size column contains the total sample size.
The nlast column contains the last n value. May be useful for restarting.
Although some designs are more ‘popular’ than others, sample size estimations are valid for all of the following designs:
"2x2x4" |
TRTR | RTRT |
| TRRT | RTTR | |
| TTRR | RRTT | |
"2x2x3" |
TRT | RTR |
| TRR | RTT | |
"2x3x3" |
TRR | RTR | RRT |
The sample size estimation for very extreme theta0 (<0.83 or >1.21) may be very
time consuming and will eventually also fail since the start values chosen are
not really reasonable in that ranges.
If you really need sample sizes in that range be prepared to restart the sample
size estimation via the argument nstart.
Since the dependence of power from n is very flat in the mentioned region you may
also consider to adapt the number of simulations not to get caught in the simulation
error trap.
We are doing the sample size estimation only for balanced designs since the
break down of the total subject number in case of unbalanced sequence groups
is not unique. Moreover the formulas used are only for balanced designs.
The minimum sample size is 6, even if the power is higher than the intended
targetpower.
Subject simulations are easily more than 100times slower than simulations based
on the ‘key’ statistics. We recommend this function only for the partial
replicate design (TRR|RTR|RRT) assuming heteroscedasticity in the case of CVwT > CVwR.
Thus be patient and go for a cup of coffee if you use this function with high
sample sizes!
H. Schütz
Tóthfalusi L, Endrényi L. Sample Sizes for Designing Bioequivalence Studies for Highly Variable Drugs. J Pharm Pharmaceut Sci. 2011;15(1):73–84. open access
power.scABEL.sdsims, sampleN.scABEL, reg_const
# using the defaults: # partial replicate design, targetpower=80%, # true assumed ratio = 0.90, 1E+5 simulated studies # ABE limits, PE constraint 0.8 - 1.25 # EMA regulatory settings # Heteroscedasticity (CVwT 0.4, CVwR 0.3) # compare results and run times CV <- c(0.4, 0.3) expl <- data.frame(method = c("subject simulations", "\'key\' statistics"), n = NA, power = NA, seconds = NA) start <- proc.time()[[3]] expl[1, 2:3] <- sampleN.scABEL.sdsims(CV = CV, print = FALSE, details = FALSE)[8:9] expl[1, 4] <- proc.time()[[3]] - start start <- proc.time()[[3]] expl[2, 2:3] <- sampleN.scABEL(CV = CV, print = FALSE, details = FALSE)[8:9] expl[2, 4] <- proc.time()[[3]] - start print(expl, row.names = FALSE) # should result in a sample size n=69, power=0.80198 for # the subject simulations and n=66, power=0.80775 for the # 'key' statistics# using the defaults: # partial replicate design, targetpower=80%, # true assumed ratio = 0.90, 1E+5 simulated studies # ABE limits, PE constraint 0.8 - 1.25 # EMA regulatory settings # Heteroscedasticity (CVwT 0.4, CVwR 0.3) # compare results and run times CV <- c(0.4, 0.3) expl <- data.frame(method = c("subject simulations", "\'key\' statistics"), n = NA, power = NA, seconds = NA) start <- proc.time()[[3]] expl[1, 2:3] <- sampleN.scABEL.sdsims(CV = CV, print = FALSE, details = FALSE)[8:9] expl[1, 4] <- proc.time()[[3]] - start start <- proc.time()[[3]] expl[2, 2:3] <- sampleN.scABEL(CV = CV, print = FALSE, details = FALSE)[8:9] expl[2, 4] <- proc.time()[[3]] - start print(expl, row.names = FALSE) # should result in a sample size n=69, power=0.80198 for # the subject simulations and n=66, power=0.80775 for the # 'key' statistics
Estimates the necessary sample size to obtain at least a target (desired) power.
sampleN.TOST(alpha = 0.05, targetpower = 0.8, logscale = TRUE, theta0, theta1, theta2, CV, design = "2x2", method = "exact", robust = FALSE, print = TRUE, details = FALSE, imax=100)sampleN.TOST(alpha = 0.05, targetpower = 0.8, logscale = TRUE, theta0, theta1, theta2, CV, design = "2x2", method = "exact", robust = FALSE, print = TRUE, details = FALSE, imax=100)
alpha |
Significance level (one-sided). Commonly set to 0.05. |
targetpower |
Power to achieve at least. Must be >0 and <1. |
logscale |
Should the data used on log-transformed ( |
theta0 |
‘True’ or assumed T/R ratio or difference. |
theta1 |
Lower (bio-)equivalence limit. |
theta2 |
Upper (bio-)equivalence limit. |
CV |
In case of In case of cross-over studies this is the within-subject CV, in case of a parallel-group design the CV of the total variability. |
design |
Character string describing the study design. |
method |
Method for calculation of the power. |
robust |
Defaults to |
print |
If |
details |
If |
imax |
Maximum number of steps in sample size search. |
The sample size is estimated via iterative evaluation of power of the TOST procedure.
Start value for the sample size search is taken from a large sample approximation
according to Zhang, modified.
The sample size is bound to 4 as minimum.
The estimated sample size gives always the total number of subjects (not subject/sequence in crossovers or subjects/group in parallel designs – like in some other software packages).
A data frame with the input and results will be returned.
The Sample size column contains the total sample size.
The function does not vectorize properly.
If you need sample sizes with varying CVs, use f.i. for-loops or the apply-family.
Of course it is highly recommended to use the default method = "exact". :-)
There is no reason besides testing and for comparative purposes to use an
approximation if the exact method is available at no extra costs.
D. Labes
Phillips KF. Power of the Two One-Sided Tests Procedure in Bioequivalence. J Pharmacokin Biopharm. 1990;18:137–44. doi:10.1007/BF01063556
Diletti D, Hauschke D, Steinijans VW. Sample Size Determination for Bioequivalence Assessment by Means of Confidence Intervals. Int J Clin Pharmacol Ther Toxicol. 1991;29(1):1–8.
Diletti D, Hauschke D, Steinijans VW. Sample size determination: Extended tables for the multiplicative model and bioequivalence ranges of 0.9 to 1.11 and 0.7 to 1.43. Int J Clin Pharmacol Ther Toxicol. 1992;30(Suppl 1):S59–62.
Zhang P. A Simple Formula for Sample Size Calculation in Equivalence Studies. J Biopharm Stat. 2003;13(3):529–38. doi:10.1081/BIP-120022772
# Exact calculation for a classical 2x2 cross-over (TR/RT), # BE limits 80 ... 125%, assumed true BE ratio 0.95, intra-subject CV=30%, # using all the default values # should give n=40 power=0.815845 sampleN.TOST(CV = 0.3) # Exact calculation for a parallel group design # evaluation on the original (untransformed) scale # BE limits 80 ... 120% = -20% ... +20% of reference, # assumed true BE ratio 0.95% = -5% to reference mean, # total CV=20% # should give n=48 (total) power=0.815435 sampleN.TOST(logscale = FALSE, theta1 = -0.2, theta0 = -0.05, CV = 0.2, design = "parallel") # A rather strange setting of theta0! Have a look at n. # It would be better this is not the sample size but the running total # of my bank account. But the first million is the hardest. ;-) sampleN.TOST(CV = 0.2, theta0 = 0.8005)# Exact calculation for a classical 2x2 cross-over (TR/RT), # BE limits 80 ... 125%, assumed true BE ratio 0.95, intra-subject CV=30%, # using all the default values # should give n=40 power=0.815845 sampleN.TOST(CV = 0.3) # Exact calculation for a parallel group design # evaluation on the original (untransformed) scale # BE limits 80 ... 120% = -20% ... +20% of reference, # assumed true BE ratio 0.95% = -5% to reference mean, # total CV=20% # should give n=48 (total) power=0.815435 sampleN.TOST(logscale = FALSE, theta1 = -0.2, theta0 = -0.05, CV = 0.2, design = "parallel") # A rather strange setting of theta0! Have a look at n. # It would be better this is not the sample size but the running total # of my bank account. But the first million is the hardest. ;-) sampleN.TOST(CV = 0.2, theta0 = 0.8005)
The (widened) scaled BE acceptance limits are calculated according to the regulatory settings of EMA, HC, FDA or via user defined regulatory settings.
scABEL(CV, regulator)scABEL(CV, regulator)
CV |
Coefficient of variation (of the Reference) as ratio. |
regulator |
Regulatory body settings for the widening of the BE acceptance limits. |
The widened BE acceptance limits are calculated by the formula [L, U] = exp(-/+ r_const * sWR)
with r_const the regulatory constant and sWR the standard deviation of the within
subjects variability of the Reference.
regulator="EMA" or regulator="HC"r_const = 0.76 (~ log(1.25)/sqrt(log(0.30^2+1)))
regulator="GCC"r_const = 0.97997... (= log(1/0.75)/sqrt(log(0.30^2+1)))
regulator="FDA"r_const = 0.89257... (= log(1.25)/0.25)
If the CVwR of the Reference is < CVswitch=0.3 the conventional ABE limits
apply (mixed procedure).
In case of regulator="EMA" a cap is placed on the widened limits if
CVwR>0.5, i.e., the widened limits are held at the value calculated for CVwR=0.5.
In case of regulator="HC" the capping is done such that the acceptance
limits are 0.6666 ... 1.5 at maximum, i.e., CVcap=0.57382.
Literally it is given by Health Canada rounded to three significant digits as 57.4%.
Returns a vector of lenghth 2 if one CV is given or a matrix if CV is given as vector
with named components lower and upper of the scaled acceptance limits.
The scaled acceptance limits (coined ‘implied limits’ by Davit et al.) are not directly used in the BE evaluation for HVDP(s) recommended by the FDA. They are included here for comparative purposes. Moreover, there are controversies where to locate the ‘implied limits’ and whether the so-called ‘desired consumer-risk model’ should be used.
D. Labes
Davit BM, Chen ML, Conner DP, Haidar SH, Kim S, Lee CH, Lionberger RA, Makhlouf FT, Nwakama PE, Patel DT, Schuirmann DJ, Yu LX. Implementation of a Reference-Scaled Average Bioequivalence Approach for Highly Variable Generic Drug Products by the US Food and Drug Administration. AAPS J. 2012;14(4):915–24. doi:10.1208/s12248-012-9406-x
Health Canada, Therapeutic Products Directorate. Guidance Document. Comparative Bioavailability Standards: Formulations Used for Systemic Effects. 2018/06/08. ISBN: 978-0-660-25514-9
power.scABEL, sampleN.scABEL, reg_const
scABEL(CV = 0.3, regulator = "EMA") # should give the conventional (unscaled) BE limits: # lower upper # 0.80 1.25 scABEL(CV = 0.5, regulator = "EMA") # should give the (maximum) expanded limits: # lower upper # 0.6983678 1.4319102 # define old ANVISA settings via reg_const() rc <- reg_const("USER", r_const = 0.76, CVswitch = 0.4, CVcap = 0.5) rc$name <- "ANVISAold" scABEL(CV = 0.4, regulator = rc) # should give the conventional (not expanded) limits: # lower upper # 0.80 1.25 scABEL(CV = 0.55, regulator = "HC") # should give the widened limits: # lower upper # 0.6765789 1.4780241 scABEL(CV = 0.55, regulator = "GCC") # should give the widened limits: # lower upper # 0.750000 1.333333 scABEL(CV = 0.55, regulator = "FDA") # should give the 'implied' limits: # lower upper # 0.6320032 1.5822705scABEL(CV = 0.3, regulator = "EMA") # should give the conventional (unscaled) BE limits: # lower upper # 0.80 1.25 scABEL(CV = 0.5, regulator = "EMA") # should give the (maximum) expanded limits: # lower upper # 0.6983678 1.4319102 # define old ANVISA settings via reg_const() rc <- reg_const("USER", r_const = 0.76, CVswitch = 0.4, CVcap = 0.5) rc$name <- "ANVISAold" scABEL(CV = 0.4, regulator = rc) # should give the conventional (not expanded) limits: # lower upper # 0.80 1.25 scABEL(CV = 0.55, regulator = "HC") # should give the widened limits: # lower upper # 0.6765789 1.4780241 scABEL(CV = 0.55, regulator = "GCC") # should give the widened limits: # lower upper # 0.750000 1.333333 scABEL(CV = 0.55, regulator = "FDA") # should give the 'implied' limits: # lower upper # 0.6320032 1.5822705
This function iteratively adjusts alpha for the BE decision via Average Bioequivalence with Expanding Limits (ABEL) based on simulations in order to maintain the consumer risk at the nominal level.
scABEL.ad(alpha = 0.05, theta0, theta1, theta2, CV, design = c("2x3x3", "2x2x4", "2x2x3"), regulator, n, alpha.pre = 0.05, imax = 100, tol, print = TRUE, details = FALSE, setseed = TRUE, nsims = 1e+06, sdsims = FALSE, progress)scABEL.ad(alpha = 0.05, theta0, theta1, theta2, CV, design = c("2x3x3", "2x2x4", "2x2x3"), regulator, n, alpha.pre = 0.05, imax = 100, tol, print = TRUE, details = FALSE, setseed = TRUE, nsims = 1e+06, sdsims = FALSE, progress)
alpha |
Type I Error (TIE) probability (nominal level of the test). Per convention commonly set to 0.05. |
theta0 |
‘True’ or assumed T/R ratio. Defaults to 0.90 according to the two Lászlós if not given explicitly. |
theta1 |
Conventional lower ABE limit to be applied in the mixed procedure
if |
theta2 |
Conventional upper ABE limit to be applied in the mixed procedure
if |
CV |
Intra-subject coefficient(s) of variation as ratio (not percent).
|
design |
Design of the study. |
regulator |
Regulatory settings for the expanding of the BE acceptance limits.
Choose from |
n |
Total sample size of the study or a vector of sample size / sequences.
If |
alpha.pre |
Pre-specified alpha (optional). Must be |
imax |
Maximum number of steps in sample size search. Defaults to 100. |
tol |
Desired accuracy (convergence tolerance). Defaults to 1E-6. |
print |
If |
details |
If |
setseed |
Simulations are dependent on the starting point of the (pseudo)
random number generator. To avoid differences in power for different
runs |
nsims |
Number of simulations to be performed to estimate the (empirical) TIE error and in each iteration of adjusting alpha. The default value 1,000,000 = 1E+6 should not be lowered. |
sdsims |
If |
progress |
Set to |
The simulations are done via the distributional properties of the statistical
quantities necessary for assessing BE based on ABEL.
Simulations for the TIE are performed at the upper (expanded) limit U
of the acceptance range. Due to the symmetry around 1 results are valid for the lower
(expanded) limit L as well.
U at the EMA’s and Health Canada’s CVcap, the GCC’s for any CVwR > 0.3:
scABEL(CV = 0.5, reg = "EMA")[["upper"]] [1] 1.43191 scABEL(CV = 0.57382, reg = "HC")[["upper"]] [1] 1.5 scABEL(CV = 0.5, reg = "GCC")[["upper"]] [1] 1.333333
Simulated studies are evaluated by ANOVA (Method A) as recommended in the
EMA’ Q&A-document and by intra-subject contrasts if regulator = "HC".
Health Canada requires a mixed-effects model which cannot be implemented in R.
However, intra-subjects contrasts are a sufficiently close approximation.
The Type I Error in ABEL depends only on CVwR and – to a
minor degree – the sample size. Algorithm:
The TIE is assessed based on alpha (or alpha.pre)
and compared to the nominal level of the test alpha.
If no inflation of the TIE is found, the algorithm stops.
Otherwise, alpha is iteratively adjusted (i.e., alpha.adj <alpha)
until no more relevant inflation of the TIE is detected (i.e.,
abs(TIE - alpha) <= tol).
Sends results to the console if argument print=TRUE (default).
Returns a list with the input, adjusted alpha, and Type I Error (for nominal
and adjusted alpha) if argument print=FALSE.
If no adjustment is necessary, NAs will be returned for the respective
variables (alpha.adj, TIE.adj, rel.change, pwr.adj, rel.loss).
Although some designs are more ‘popular’ than others, power calculations are valid for all of the following designs:
"2x2x4" |
TRTR | RTRT |
| TRRT | RTTR | |
| TTRR | RRTT | |
"2x2x3" |
TRT | RTR |
| TRR | RTT | |
"2x3x3" |
TRR | RTR | RRT |
See the Warning section of the function power.scABEL concerning
the power value agreement to the one obtained by simulations via subject data.
Specifying theta0 is not necessary.
If theta0 is not given, achievable power for the common target
of 0.80 (both for alpha and adjusted alpha) will be estimated. If
theta0 is specified, its value will be used; again for target power 0.80.
If you are interested in other levels of power, use sampleN.scABEL.ad.
The EMA’s method is currently recommended in other jurisdictions as well (e.g., by the WHO;
in ASEAN States, Australia, Brazil, Egypt, the Eurasian Economic Union, New Zealand, and the East African Community).
If CVwR > 30%, fixed wider limits of 0.7500–1.3333 are recommended by the Gulf Cooperation Council (Bahrain, Kuwait, Oman, Qatar, Saudi Arabia, United Arab Emirates).
H. Schütz
Wonnemann M, Frömke C, Koch A. Inflation of the Type I Error: Investigations on Regulatory Recommendations for Bioequivalence of Highly Variable Drugs. Pharm Res. 2015;32(1):135–43. doi:10.1007/s11095-014-1450-z
Muñoz J, Alcaide D, Ocaña J. Consumer’s risk in the EMA and FDA regulatory approaches for bioequivalence in highly variable drugs. Stat Med. 2015;35(12):1933–43. doi:10.1002/sim.6834
Labes D, Schütz H. Inflation of Type I Error in the Evaluation of Scaled Average Bioequivalence, and a Method for its Control. Pharm Res. 2016;33(11):2805–14. doi:10.1007/s11095-016-2006-1
Tóthfalusi L, Endrényi L. Algorithms for Evaluating Reference Scaled Average Bioequivalence: Power, Bias, and Consumer Risk. Stat Med. 2017;36(27):4378–90. doi:10.1002/sim.7440
Molins E, Cobo E, Ocaña J. Two-Stage Designs Versus European Scaled Average Designs in Bioequivalence Studies for Highly Variable Drugs: Which to Choose? Stat Med. 2017;36(30):4777–88. doi:10.1002/sim.7452
European Medicines Agency, Committee for Medicinal Products for Human Use. Guideline on the Investigation of Bioequivalence. London, 20 January 2010. CPMP/EWP/QWP/1401/98 Rev. 1/ Corr **
European Medicines Agency, Committee for Medicinal Products for Human Use. Questions & Answers: positions on specific questions addressed to the Pharmacokinetics Working Party (PKWP). London, 19 November 2015. EMA/618604/2008 Rev. 13
Health Canada, Therapeutic Products Directorate. Comparative Bioavailability Standards: Formulations Used for Systemic Effects, 2.1.1.8 Highly variable drug products Ottawa, 08 June 2018. online
Executive Board of the Health Ministers’ Council for GCC States. The GCC Guidelines for Bioequivalence. May 2021. Version 3.0 Saudi Food & Drug Authority The GCC Guidelines for Bioequivalence. Version 3.1
sampleN.scABEL.ad, power.scABEL, power.scABEL.sdsims, scABEL
# Using all defaults: # TRR|RTR|RRT, target power 80% for assumed ratio 0.90 (estimated sample size 54), # EMA regulatory settings (ABE limits and PE constraint 0.80 - 1.25), # 1E+6 simulated studies. # Not run: due to timing policy of CRAN for examples scABEL.ad(CV = 0.3) # Should result in adjusted alpha 0.03389 (TIE 0.5000, TIE for nominal alpha 0.07189). # # As above but subject data simulations. scABEL.ad(CV = 0.3, sdsims = TRUE) # Should result in adjusted alpha 0.03336 (TIE 0.5000, TIE for nominal alpha 0.07237). # # TRT|RTR, heteroscedasticity, sample size 48 (unbalanced), subject data simulations. scABEL.ad(CV = c(0.25, 0.3), design = "2x2x3", n = c(23, 25), sdsims = TRUE) # Should result in adjusted alpha 0.02465 (TIE 0.5000, TIE for nominal alpha 0.09050). # # TRTR|RTRT, CV 0.35, sample size 33 (unbalanced). scABEL.ad(CV = 0.35, design = "2x2x4", n = c(16, 17)) # Should result in adjusted alpha 0.03632 (TIE 0.5000, TIE for nominal alpha 0.06544).# Using all defaults: # TRR|RTR|RRT, target power 80% for assumed ratio 0.90 (estimated sample size 54), # EMA regulatory settings (ABE limits and PE constraint 0.80 - 1.25), # 1E+6 simulated studies. # Not run: due to timing policy of CRAN for examples scABEL.ad(CV = 0.3) # Should result in adjusted alpha 0.03389 (TIE 0.5000, TIE for nominal alpha 0.07189). # # As above but subject data simulations. scABEL.ad(CV = 0.3, sdsims = TRUE) # Should result in adjusted alpha 0.03336 (TIE 0.5000, TIE for nominal alpha 0.07237). # # TRT|RTR, heteroscedasticity, sample size 48 (unbalanced), subject data simulations. scABEL.ad(CV = c(0.25, 0.3), design = "2x2x3", n = c(23, 25), sdsims = TRUE) # Should result in adjusted alpha 0.02465 (TIE 0.5000, TIE for nominal alpha 0.09050). # # TRTR|RTRT, CV 0.35, sample size 33 (unbalanced). scABEL.ad(CV = 0.35, design = "2x2x4", n = c(16, 17)) # Should result in adjusted alpha 0.03632 (TIE 0.5000, TIE for nominal alpha 0.06544).
Was designed to calculate the type I error rate of two simultaneous TOST procedures (where the two parameters of the two TOSTs are correlated with some correlation) for various study designs used in BE studies.
Is defunct now since it suffers from insufficient precision to obtain the
type 1 error (TIE) via simulations.
Due to the intersection-union principle the TIE is always upper bounded
to alpha by theory.