gep_by_nera: Points on confidence region (found by Newton-Raphson search)

Description Usage Arguments Details Value Comparison of highly variable dissolution profiles Similarity limits in terms of MSD References See Also Examples

View source: R/utility.R

Description

The function gep_by_nera() is a function for finding points on confidence region (CR) bounds by aid of the “Method of Lagrange Multipliers” (MLM) and by “Newton-Raphson” (nera) optimisation. The multivariate confidence interval for profiles with four time points, e.g., is an “ellipse” in four dimensions.

Usage

1
gep_by_nera(n_p, K, mean_diff, S_pool, F_crit, y, max_trial, tol)

Arguments

n_p

A positive integer specifying the number of (time) points n_p.

K

A non-negative numeric value specifying the scaling factor K for the calculation of the Hotelling's T^2 statistic.

mean_diff

A vector of the mean differences between the dissolution profiles of the reference and the test batch. It must have the length specified by the parameter n_p.

S_pool

The pooled variance-covariance matrix of the dissolution profiles of the reference and the test batch. It must have the dimension n_p \times n_p.

F_crit

The critical F value (i.e. a non-negative numeric).

y

A numeric vector of y values that serve as starting points for the Newton-Raphson search, i.e. values supposed to lie on or close to the confidence interval bounds. It must have a length of n_p + 1.

max_trial

A positive integer specifying the maximum number of Newton-Raphson search rounds to be performed.

tol

A non-negative numeric specifying the accepted minimal difference between two consecutive search rounds.

Details

The function gep_by_nera() determines the points on the CR bounds for each of the n_p time points. It does so by aid of the “Method of Lagrange Multipliers” (MLM) and by “Newton-Raphson” (nera) optimisation, as proposed by Margaret Connolly (Connolly 2000).

For more information, see the sections “Comparison of highly variable dissolution profiles” and “Similarity limits in terms of MSD” below.

Value

A list with the following elements is returned:

points

A matrix with one column and n_p + 1 rows is returned, where rows 1 to n_p represent, for each time point, the points on the CR bounds. For symmetry reasons, the points on the opposite side are obtained by addition/subtraction. The last row in the matrix, with index n_p + 1, represents the λ parameter of the MLM, also known as lambda multiplier method, that is used to optimise under constraint(s). The variable λ is thus called the Lagrange multiplier.

converged

A logical stating if the NR algorithm converged or not.

n.trial

Number of trials until convergence.

max.trial

Maximal number of trials.

Comparison of highly variable dissolution profiles

When comparing the dissolution data of a post-approval change product and a reference approval product, the goal is to assess the similarity between the mean dissolution values at the observed sample time points. A widely used method is the f_2 method that was introduced by Moore & Flanner (1996). Similarity testing criteria based on f_2 can be found in several FDA guidelines and in the guideline of the European Medicines Agency (EMA) “On the investigation of bioequivalence” (EMA 2010).

In situations where within-batch variation is greater than 15%, FDA guidelines recommend use of a multivariate confidence interval as an alternative to the f_2 method. This can be done using the following stepwise procedure:

  1. Establish a similarity limit in terms of “Multivariate Statistical Distance” (MSD) based on inter-batch differences in % drug release from reference (standard approved) formulations, i.e. the so- called “Equivalence Margin” (EM).

  2. Calculate the MSD between test and reference mean dissolutions.

  3. Estimate the 90% confidence interval (CI) of the true MSD as determined in step 2.

  4. Compare the upper limit of the 90% CI with the similarity limit determined in step 1. The test formulation is declared to be similar to the reference formulation if the upper limit of the 90% CI is less than or equal to the similarity limit.

Similarity limits in terms of MSD

For the calculation of the “Multivariate Statistical Distance” (MSD), the procedure proposed by Tsong et al. (1996) can be considered as well-accepted method that is actually recommended by the FDA. According to this method, a multivariate statistical distance, called Mahalanobis distance, is used to measure the difference between two multivariate means. This distance measure is calculated as

D_M = sqrt((x_T - x_R)^{\top} S_{pooled}^{-1} (x_T - x_R)) ,

where S_{pooled} = (S_T + S_R) / 2 is the sample variance-covariance matrix pooled across the batches, x_T and x_R are the vectors of the sample means for the test (T) and reference (R) profiles, and S_T and x_R are the variance-covariance matrices of the test and reference profiles.

In order to determine the similarity limits in terms of the MSD, i.e. the Mahalanobis distance between the two multivariate means of the dissolution profiles of the formulations to be compared, Tsong et al. (1996) proposed using the equation

D_M^{max} = sqrt(d_g^{\top} S_{pooled}^{-1} d_g) ,

where d_g is a 1 x p vector with all p elements equal to an empirically defined limit d_g, e.g., 15%, for the maximum tolerable difference at all time points, and p is the number of sampling points. By assuming that the data follow a multivariate normal distribution, the 90% confidence region (CR) bounds for the true difference between the mean vectors, μ_T - μ_R, can be computed for the resultant vector μ to satisfy the following condition:

CR = sqrt((μ - (x_T - x_R))^{\top} S_{pooled}^{-1} (μ - (x_T - x_R))) ≤q F_{p, n_T + n_R - p - 1, 0.9} ,

where K is the scaling factor that is calculated as

(n_T n_R) / (n_T + n_R) * (n_T + n_R - p - 1) / ((n_T + n_R - 2) p) ,

and F_{p, n_T + n_R - p - 1, 0.9} is the 90^{th} percentile of the F distribution with degrees of freedom p and n_T + n_R - p - 1. It is obvious that (n_T + n_R) must be greater than (p + 1). The formula for CR gives a p-variate 90% confidence region for the possible true differences.

References

United States Food and Drug Administration (FDA). Guidance for industry: dissolution testing of immediate release solid oral dosage forms. 1997.
https://www.fda.gov/media/70936/download

United States Food and Drug Administration (FDA). Guidance for industry: immediate release solid oral dosage form: scale-up and post-approval changes, chemistry, manufacturing and controls, in vitro dissolution testing, and in vivo bioequivalence documentation (SUPAC-IR). 1995.
https://www.fda.gov/media/70949/download

European Medicines Agency (EMA), Committee for Medicinal Products for Human Use (CHMP). Guideline on the Investigation of Bioequivalence. 2010; CPMP/EWP/QWP/1401/98 Rev. 1.
https://www.ema.europa.eu/en/documents/scientific-guideline/guideline-investigation-bioequivalence-rev1_en.pdf

Moore, J.W., and Flanner, H.H. Mathematical comparison of curves with an emphasis on in-vitro dissolution profiles. Pharm Tech. 1996; 20(6): 64-74.

Tsong, Y., Hammerstrom, T., Sathe, P.M., and Shah, V.P. Statistical assessment of mean differences between two dissolution data sets. Drug Inf J. 1996; 30: 1105-1112.
doi: 10.1177/009286159603000427

Connolly, M. SAS(R) IML Code to calculate an upper confidence limit for multivariate statistical distance; 2000; Wyeth Lederle Vaccines, Pearl River, NY.
https://analytics.ncsu.edu/sesug/2000/p-902.pdf

See Also

mimcr, bootstrap_f2.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
# Collecting the required information
time_points <- suppressWarnings(as.numeric(gsub("([^0-9])", "",
                                                colnames(dip1))))
tico <- which(!is.na(time_points))
b1 <- dip1$type == "R"

n_tp <- length(tico)
n_b1 <- length(dip1[b1, "type"])
n_b2 <- length(dip1[!b1, "type"])
df_b1 <- n_tp
df_b2 <- n_b1 + n_b2 - n_tp - 1
K_limit <- (n_b2 * n_b1) / (n_b2 + n_b1)
K <- K_limit * df_b2 / (df_b1 * (n_b2 + n_b1 - 2))

mean_b1 <- apply(X = dip1[b1, tico], MARGIN = 2, FUN = mean)
mean_b2 <- apply(X = dip1[!b1, tico], MARGIN = 2, FUN = mean)
mean_diff <- mean_b2 - mean_b1

S_b1 <- cov(dip1[b1, tico])
S_b2 <- cov(dip1[!b1, tico])
S <- ((n_b1 - 1) * S_b1 + (n_b2 - 1) * S_b2) / (n_b1 + n_b2 - 2)

F_crit <- qf(p = (1 - 0.05), df1 = df_b1, df2 = df_b2)
y_b1 <- rep(1, times = (n_tp + 1))

# Calling gep_by_nera()
res <- gep_by_nera(n_p = n_tp, K = K, mean_diff = mean_diff, S_pool = S,
                   F_crit = F_crit, y = y_b1, max_trial = 100, tol = 1e-9)

# Expected result in res[["points"]]
#              [,1]
# t.5   -15.7600077
# t.10  -13.6501734
# t.15  -11.6689469
# t.20   -9.8429369
# t.30   -6.6632182
# t.60   -0.4634318
# t.90    2.2528551
# t.120   3.3249569
#       -17.6619995

# Rows t.5 to t.120 represent the points on the CR bounds.The unnamed last row
# represents the Lagrange multiplier lambda.

# If 'max_trial' is too small, the Newton-Raphson search may not converge.
tryCatch(
  gep_by_nera(n_p = n_tp, K = K, mean_diff = mean_diff, S_pool = S,
              F_crit = F_crit, y = y_b1, max_trial = 5, tol = 1e-9),
  warning = function(w) message(w),
  finally = message("\nMaybe increasing the number of max_trial could help."))

disprofas documentation built on Dec. 8, 2021, 5:10 p.m.