--- title: "Confidence Intervals of $f_2$ Using Bootstrap Method" author: "Zhengguo XU" date: "`r format(Sys.Date())`" output: rmarkdown::html_vignette: keep_md: true toc: false toc_depth: 3 fig_width: 7 fig_height: 4.5 highlight: "tango" bibliography: ref.bib notes-after-punctuation: false link-citations: yes csl: ref.csl vignette: > %\VignetteIndexEntry{Confidence Intervals of $f_2$ Using Bootstrap Method} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup0, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#" ) ``` ```{r setup, echo=FALSE} library(bootf2) ``` ## Introduction To compare the dissolution profiles, the most widely used method is the similarity factor $f_2$, which is usually defined in the regulatory guideline as $$f_2 = 50 \log\frac{100}{\sqrt{1 + \frac{\sum_{t=1}^{t=n}\left(\bar{R}_t-\bar{T}_t\right)^2}{n}}},$$ where $\bar{R}_t$ and $\bar{T}_t$ are mean dissolution profile of reference ad test product at time $t$; $n$ is the number of time points. Nevertheless, there are several prerequisites for the use of f2 method according to the regulatory guidelines, such as no more than one time point above 85% dissolved should be used and the variability, expressed as coefficient of variation (CV), should be no more than 20% and 10% for early and later time points, respectively. See vignette *Introduction to bootf2* for details. Recently, several guidelines recommended the use of confidence interval of $f_2$ approach using bootstrap when the traditional $f_2$ method cannot be applied [@EMA-2018-09-QA.MSD.DISSO; @Davit-2013-03-BA; @Lum-2019-05-WS; @Mandula-2019-05-WS]. However, none of the guidelines specified in details which estimators or types of confidence intervals should be used. Therefore, the function `bootf2()` will output various confidence intervals by default for several $f_2$ estimators. ### Types of $f_2$ According to Shah et al, the population $f_2$ for the inference is $$ f_2 = 100-25\log\left(1 + \frac{1}{P}\sum_{i=1}^P \left(\mu_{\mathrm{T},i} - \mu_{\mathrm{R},i}\right)^2 \right)\, , $$ where $P$ is the number of time points; $\mu_{\mathrm{T},i}$ and $\mu_{\mathrm{R},i}$, typically unknown, are population mean of test and reference product at time point $i$, respectively [@Shah-1998-06-PR]. Several estimators for the population $f_2$ has been described in the literature [@Shah-1998-06-PR; @Ma-2000-05-JBS; @Mendyk-2013-02-DT; @Mendyk-2019]. The function is able to calculate the following five estimators for $f_2$, along with their corresponding 90% confidence intervals using bootstrap. #### 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 [@Shah-1998-06-PR], 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 [@Ma-2000-05-JBS; @Mendyk-2013-02-DT; @Mendyk-2019]. $$ \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)\,. $$ ### Types of confidence intervals The following 90% confidence intervals are included in the function: #### The Normal interval The Normal interval (denoted by `normal` in the function) with bias correction was estimated according to Davison and Hinkley [@Davison-1997-Bootstrap], \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 [@Davison-1997-Bootstrap], \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 [@Hyndman-1996-AS] 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 [@Efron1993AIttB], \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. [@Shah-1998-06-PR] 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 [@Mendyk-2019]. This approach, denoted by `bca.boot` in the function, was also implemented in the function `bootf2` for estimating the interval. ## Usage The complete list of arguments of the function is as follows: ```{r bootf2-code0, eval=FALSE} bootf2(test, ref, path.in, file.in, path.out, file.out, n.boots = 10000L, seed = 306, digits = 2L, alpha = 0.05, regulation = c("EMA", "FDA", "WHO"), min.points = 1L, both.TR.85 = FALSE, print.report = TRUE, report.style = c("concise", "intermediate", "detailed"), f2.type = c("all", "est.f2", "exp.f2", "bc.f2", "vc.exp.f2", "vc.bc.f2"), ci.type = c("all", "normal", "basic", "percentile", "bca.jackknife", "bca.boot"), quantile.type = c("all", 1:9, "boot"), jackknife.type = c("nt+nr", "nt*nr", "nt=nr"), time.unit = c("min", "h"), output.to.screen = FALSE, sim.data.out = FALSE) ``` ### Notes on function arguments 1. Data input: `test`, `ref`, `path.in`, `file.in` - In the typical interactive use, `test` and `ref` are *data frames* with the time as the first column, and individual units for the rest of columns. In such cases, arguments `path.in` and `file.in` should not be used. - Data can be read directly from an Excel file with extension `.xlsx` or `.xls`. In this case, data of test and reference should be stored in separate worksheets. The first column should be time, the rest columns are individual units. *The first row* is the column head indicating the names, such as 'time', 'unit 01', unit 02', .... It doesn't matter what the names are as columns will be renamed internally by the function. The important point is that *the first row will not be read, so do not put dissolution data on the first row*. - When `path.in` and `file.in` are provided, the argument `test` and `ref` should be the *worksheet names inside quotation mark*, e.g., `"lot ABCD1234 pH 6.8"`. - `path.in` can be an absolute path such as`"/home/myname/my.project/dat/"`, or a relative path such as `"../dat/"` if the working directory is in the folder `"/home/myname/my.project/analysis/"` and the data file is in the folder `"/home/myname/my.project/dat/"`. - One more note for Windows user: As Windows use "\" instead of "/" to separate path, you have to either escape it by an additional "\", e.g., `"C:\user\myname\my.project\dat\"` *cannot* be the `path.in`, you have to changed it to `"C:\\user\\myname\\my.project\\dat\\"`, or to `"C:/user/myname/my.project/dat/"`, the same format as used in Linux system. 1. Output: `path.out` and `file.out`, `output.to.screen`: - If argument `path.out` and `file.out` are not provided, but argument `print.report` is `TRUE`, then a plain text report will be generated automatically in the current working directory with file name `test_vs_ref_TZ_YYYY-MM-DD_HHMMSS.txt`, where `test` and `ref` are data set names of test and reference, `TZ` is the time zone such as CEST, `YYYY-MM-DD` is the numeric date format and `HHMMSS` is the numeric time format for hour, minute, and second. - For a quick check, set argument `output.to.screen = TRUE`, a summary report very similar to concise style report will be printed on screen. 1. Argument `jackknife.type` - For any sample with size $n$, the jackknife estimator is obtained by estimating the parameter for each subsample omitting the $i$th observation. However, when two samples (e.g., test and reference) are involved, there are several possible ways to do it. Assuming sample size of test and reference are $n_\mathrm{T}$ and $n_\mathrm{R}$, the following three possibility are considered: - Estimated by removing one observation from both test and reference samples. In this case, the prerequisite is $n_\mathrm{T} = n_\mathrm{R}$, denoted by `nt=nr` in the function. So if there are 12 units in test and reference data sets, there will be 12 jackknife estimators. - Estimate the jackknife for test sample while keeping the reference data unchanged; and then estimate jackknife for reference sample while keeping the test sample unchanged. This is denoted by `nt+nr` in the function. This is the default method. So if there are 12 units in test and reference data sets, there will be $12 + 12 = 24$ jackknife estimators. - For each observation deleted from test sample, estimate jackknife for reference sample. This is denoted by `nt*nr` in the function. So if there are 12 units in test and reference data sets, there will be $12 \times 12 = 144$ jackknife estimators. 1. Arguments `f2.type` and `ci.type` are explained in the previous section. 1. By default, the individual bootstrapped data set are not included in the output as those data sets are only useful for validation purpose. To include the individual data sets, set the argument `sim.data.out = TRUE`. 1. Read the function manual `help("bootf2")` for more details of each argument. ## Examples The minimum required arguments are `test` and `ref`. First, let's simulate some dissolution profiles to play with. ```{r bootf2-sim01} # time points tp <- c(5, 10, 15, 20, 30, 45, 60) # model.par for reference with low variability par.r <- list(fmax = 100, fmax.cv = 3, mdt = 15, mdt.cv = 14, tlag = 0, tlag.cv = 0, beta = 1.5, beta.cv = 8) # simulate reference data dr <- sim.dp(tp, model.par = par.r, seed = 100, plot = FALSE) # model.par for test par.t <- list(fmax = 100, fmax.cv = 3, mdt = 12.29, mdt.cv = 12, tlag = 0, tlag.cv = 0, beta = 1.727, beta.cv = 9) # simulate test data with low variability dt <- sim.dp(tp, model.par = par.t, seed = 100, plot = FALSE) ``` Use default setting for most arguments. To reduce the compilation time, bootstrap number is set to 100 only. ```{r bootf2-ex01, eval=FALSE} # output file will be generated automatically bootf2(dt$sim.disso, dr$sim.disso, n.boots = 100, print.report = FALSE, output.to.screen = TRUE) ```