Confidence Intervals of $f_2$ Using Bootstrap Method

The estimated $f_2$ ($\hat{f}_2$)

The estimated $f_2$ ($\hat{f}_2$, denoted by est.f2 in the function) is the one written in various regulatory guidelines. It is expressed differently, but mathematically equivalently, as $$\hat{f}_2 = 100-25\log\left(1 + \frac{1}{P}\sum_{i=1}^P\left( \bar{X}_{\mathrm{T},i} - \bar{X}_{\mathrm{R},i}\right)^2 \right)\:,$$ where $P$ is the number of time points; $\bar{X}_{\mathrm{T},i}$ and $\bar{X}_{\mathrm{R},i}$ are mean dissolution data at the $i$th time point of random samples chosen from the test and the reference population, respectively. Compared to the equation of population $f_2$ above, the only difference is that in the equation of $\hat{f}_2$ the sample means of dissolution profiles replace the population means for the approximation. In other words, a point estimate is used for the statistical inference in practice.

The Bias-corrected $f_2$ ($\hat{f}_{2,\mathrm{bc}}$)

Bias-corrected $f_2$ ($\hat{f}_{2,\mathrm{bc}}$, denoted by bc.f2 in the function) was described in the article of Shah et al⁵, as $$\hat{f}_{2,\mathrm{bc}} = 100-25\log\left(1 + \frac{1}{P} \left(\sum_{i=1}^P\left(\bar{X}_{\mathrm{T},i} - \bar{X}_{\mathrm{R},i}\right)^2 - \frac{1}{n}\sum_{i=1}^P \left(S_{\mathrm{T},i}^2 + S_{\mathrm{R},i}^2\right)\right)\right)\,,$$ where $S_{\mathrm{T},i}^2$ and $S_{\mathrm{R},i}^2$ are unbiased estimates of variance at the $i$th time points of random samples chosen from test and reference population, respectively; and $n$ is the sample size. $\bar{X}_{\mathrm{T},i}$, $\bar{X}_{\mathrm{R},i}$ and $P$ are same as described previously. As domain of the $\log(\cdot)$ function needs to be positive, when variance is sufficiently high, i.e., when inequality $$\frac{1}{n}\sum_{i=1}^P\left(S_{\mathrm{T},i}^2 + S_{\mathrm{R},i}^2\right)\ge \sum_{i=1}^P \left(\bar{X}_{\mathrm{T},i}-\bar{X}_{\mathrm{R},i}\right)^2+P$$ is valid, $\hat{f}_{2,\mathrm{bc}}$ cannot be calculated.

The variance- and bias-corrected $f_2$ ($\hat{f}_{2, \mathrm{vcbc}}$)

Bias-corrected $f_2$ described earlier assumes equal weight of variance between test and reference, which might not necessarily be true, therefore, the function also includes the following variance- and bias-corrected $f_2$ ($\hat{f}_{2, \mathrm{vcbc}}$, denoted by vc.bc.f2 in the function) $$\hat{f}_{2, \mathrm{vcbc}} = 100-25\log\left(1 + \frac{1}{P}\left(\sum_{i=1}^P \left(\bar{X}_{\mathrm{T},i} - \bar{X}_{\mathrm{R},i}\right)^2 - \frac{1}{n}\sum_{i=1}^P\left(w_{\mathrm{T},i}\cdot S_{\mathrm{T},i}^2 + w_{\mathrm{R},i}\cdot S_{\mathrm{R},i}^2\right)\right)\right)\,,$$ where $w_{\mathrm{T},i}$ and $w_{\mathrm{R},i}$ are weighting factors for variance of test and reference products, respectively, which can be calculated as follows: $$w_{\mathrm{T},i} = 0.5 + \frac{S_{\mathrm{T},i}^2} {S_{\mathrm{T},i}^2 + S_{\mathrm{R},i}^2}\,,$$ and $$w_{\mathrm{R},i} = 0.5 + \frac{S_{\mathrm{R},i}^2} {S_{\mathrm{T},i}^2 + S_{\mathrm{R},i}^2}\,.$$ Similar to $\hat{f}_{2, \mathrm{bc}}$, the $\hat{f}_{2, \mathrm{vcbc}}$ cannot be estimated when the variance is sufficiently high, i.e., when the inequality $$\frac{1}{n}\sum_{i=1}^P \left(w_{\mathrm{T},i}\cdot S_{\mathrm{T},i}^2 + w_{\mathrm{R},i}\cdot S_{\mathrm{R},i}^2\right) \ge \sum_{i=1}^P \left(\bar{X}_{\mathrm{T},i}-\bar{X}_{\mathrm{R},i}\right)^2+P$$ is valid.

The expected $f_2$ ($\hat{f}_{2, \mathrm{exp}}$)

The mathematical expectation of $\hat{f}_2$ is $$\begin{align}\label{eq:mathexp} \mathbb{E}\left(\hat{f}_2\right) &= \mathbb{E}\left(100 - 25\log\left(1 + \frac{1}{P} \sum_{i=1}^P\left(\bar{X}_{\mathrm{T},i} - \bar{X}_{\mathrm{R},i}\right)^2 \right)\right)\notag\\ &\approx 100 - 25\log\left(1 + \mathbb{E}\left(\frac{1}{P} \sum_{i=1}^P \left(\bar{X}_{\mathrm{T},i} - \bar{X}_{\mathrm{R},i}\right)^2\right)\right) \notag\\ &= 100 - 25\log\left(1 + \frac{1}{P} \left(\sum_{i=1}^P\left(\mu_{\mathrm{T},i} - \mu_{\mathrm{R},i}\right)^2 + \frac{1}{n}\sum_{i=1}^P \left(\sigma_{\mathrm{T},i}^2 + \sigma_{\mathrm{R},i}^2\right)\right)\right), \end{align}$$ where $\sigma_{\mathrm{T},i}^2$ and $\sigma_{\mathrm{R},i}^2$ are population variance of the test and the reference product, respectively.

The expected $f_2$ ($\hat{f}_{2, \mathrm{exp}}$, denoted by exp.f2 in the function) is calculated based on the mathematical expectation of $\hat{f}_2$, using mean dissolution profiles and variance from samples for the approximation^6–8. $$\hat{f}_{2, \mathrm{exp}} = 100-25\log\left(1 + \frac{1}{P}\left(\sum_{i=1}^P \left(\bar{X}_{\mathrm{T},i} - \bar{X}_{\mathrm{R},i}\right)^2 + \frac{1}{n}\sum_{i=1}^P\left(S_{\mathrm{T},i}^2 + S_{\mathrm{R},i}^2\right)\right)\right)\,.$$

The variance-corrected expected $f_2$ ($\hat{f}_{2, \mathrm{vcexp}}$))

Similarly, since $\hat{f}_{2, \mathrm{exp}}$ assumes equal weight of variance between test and reference, the variance-corrected version ($\hat{f}_{2, \mathrm{vcexp}}$, denoted by vc.exp.f2 in the function) was also included in the function. $$\hat{f}_{2, \mathrm{vcexp}} = 100-25\log\left(1 + \frac{1}{P}\left(\sum_{i=1}^P \left(\bar{X}_{\mathrm{T},i} - \bar{X}_{\mathrm{R},i}\right)^2 + \frac{1}{n}\sum_{i=1}^P\left(w_{\mathrm{T},i}\cdot S_{\mathrm{T},i}^2 + w_{\mathrm{R},i}\cdot S_{\mathrm{R},i}^2\right)\right)\right)\,.$$

The Normal interval

The Normal interval (denoted by normal in the function) with bias correction was estimated according to Davison and Hinkley⁹, $$\begin{align} \hat{f}_{2, \mathrm{L}} &= \hat{f}_{2, \mathrm{S}} - E_B - \sqrt{V_B}\cdot Z_{1-\alpha}\,, \\ \hat{f}_{2, \mathrm{U}} &= \hat{f}_{2, \mathrm{S}} - E_B + \sqrt{V_B}\cdot Z_{1-\alpha}\,, \end{align}$$ where $\hat{f}_{2, \mathrm{L}}$ and $\hat{f}_{2, \mathrm{U}}$ are the lower and upper limit of the confidence interval estimated from bootstrap samples; $\hat{f}_{2, \mathrm{S}}$ denotes the estimators as described in Section Types of $f_2$ and calculated using the randomly selected samples from populations; $Z_{1-\alpha}$ is $\Phi^{-1}(1-\alpha)$, where $\Phi(\cdot)$ represents standard normal cumulative distribution function and $\Phi^{-1}(\cdot)$ denotes its inverse, the quantile function; $\alpha$ is the type I error rate and equal to 0.05 in the current study; $E_B$ and $V_B$ are the resampling estimates of bias and variance calculated as $$\begin{align} E_B &= \frac{1}{B}\sum_{b=1}^{B}\hat{f}_{2,b}^\star - \hat{f}_{2, \mathrm{S}} = \bar{f}_2^\star - \hat{f}_{2, \mathrm{S}}\label{eq:eb}\,, \\ V_B &= \frac{1}{B-1}\sum_{b=1}^{B} \left(\hat{f}_{2,b}^\star-\bar{f}_2^\star\right)^2\label{eq:vb}\,, \end{align}$$ where superscript $^\star$ denotes that the $f_2$ estimates are calculated from bootstrap samples; $B$ is the number of bootstrap samples, which is equal to 10000 in the function as default value; $\hat{f}_{2,b}^\star$ is the $f_2$ estimate with the $b$th bootstrap sample and $\bar{f}_2^\star$ is the mean value.

The basic interval

The basic interval (denoted by basic in the function) was calculated according to Davison and Hinkley⁹, $$\begin{align} \hat{f}_{2, \mathrm{L}} &= 2\hat{f}_{2, \mathrm{S}} - \hat{f}_{2,(B+1)(1-\alpha)}^\star\,, \\ \hat{f}_{2, \mathrm{U}} &= 2\hat{f}_{2, \mathrm{S}} - \hat{f}_{2,(B+1)\alpha}^\star\,, \end{align}$$ where $\hat{f}_{2,(B+1)\alpha}^\star$ and $\hat{f}_{2,(B+1)(1-\alpha)}^\star$ are the $(B+1)\alpha$th and the $(B+1)(1-\alpha)$th ordered resampling estimates of $f_2$, respectively. When $(B+1)\alpha$ is not an integer, the following equation is used for interpolation, $$\hat{f}_{2,(B+1)\alpha}^\star = \hat{f}_{2,k}^\star + \frac{\Phi^{-1}\left(\alpha\right)-\Phi^{-1}\left(\frac{k}{B+1}\right)} {\Phi^{-1}\left(\frac{k+1}{B+1}\right)-\Phi^{-1}\left(\frac{k}{B+1}\right)} \left(\hat{f}_{2,k+1}^\star-\hat{f}_{2,k}^\star\right),$$ where $k = \left[(B+1)\alpha\right]$, the integer part of $(B+1)\alpha$, $\hat{f}_{2,k+1}^\star$ and $\hat{f}_{2,k}^\star$ are the $(k+1)$th and the $k$th ordered resampling estimates of $f_2$, respectively.

The percentile interval

The percentile intervals (denoted by percentile in the function) were estimated using nine different types of quantiles, percentile Type 1 to Type 9, as summarized in Hyndman and Fan’s article¹⁰ and implemented in R’s quantile function. Using R’s boot package, program bootf2BCA outputs a percentile interval using the equation above for interpolation. To be able to compare the results among different programs, the same interval, denoted by Percentile (boot) in the function, is also calculated in this study.

The BCa intervals

The bias-corrected and accelerated intervals were estimated as follows according to literature¹¹, $$\begin{align} \hat{f}_{2, \mathrm{L}} &= \hat{f}_{2, \alpha_1}^\star\,,\\ \hat{f}_{2, \mathrm{U}} &= \hat{f}_{2, \alpha_2}^\star\,, \end{align}$$ where $\hat{f}_{2, \alpha_1}^\star$ and $\hat{f}_{2, \alpha_2}^\star$ are the $100\alpha_1$th and the $100\alpha_2$th percentile of the resampling estimates of $f_2$, respectively. Type I errors $\alpha_1$ and $\alpha_2$ were obtained as $$\begin{align} \alpha_1 &= \Phi\left(\hat{z}_0 + \frac{\hat{z}_0 + \hat{z}_\alpha} {1-\hat{a}\left(\hat{z}_0 + \hat{z}_\alpha\right)}\right), \\ \alpha_2 &= \Phi\left(\hat{z}_0 + \frac{\hat{z}_0 + \hat{z}_{1-\alpha}} {1-\hat{a}\left(\hat{z}_0 + \hat{z}_{1-\alpha}\right)}\right), \end{align}$$ where $\hat{z}_0$ and $\hat{a}$ are called bias-correction and acceleration, respectively.

There are different methods to estimate the $\hat{z}_0$ and $\hat{a}$. Shah et al.⁵ used jackknife technique (denoted by bca.jackknife) as $$\begin{equation}\label{eq:z0} \hat{z}_0 = \Phi^{-1}\left(\frac{\#\left\{\hat{f}_{2,b}^\star < \hat{f}_{2,\mathrm{S}} \right\}}{B}\right), \end{equation}$$ and $$\begin{equation}\label{eq:a0} \hat{a} = \frac{\sum_{i=1}^{n}\left(\hat{f}_{2,\mathrm{m}} - \hat{f}_{2, i}\right)^3}{6\left(\sum_{i=1}^{n} \left(\hat{f}_{2,\mathrm{m}} - \hat{f}_{2, i}\right)^2\right)^{3/2}}\,, \end{equation}$$ where $\#\left\{\cdot\right\}$ denotes the number of element in the set, $\hat{f}_{2, i}$ is the $i$th jackknife statistic, $\hat{f}_{2,\mathrm{m}}$ is the mean of the jackknife statistics, and $n$ is the sample size.

Program bootf2BCA gives a slightly different BCa interval with R’s boot package⁸. This approach, denoted by bca.boot in the function, was also implemented in the function bootf2 for estimating the interval.