Bioequivalence and Bioavailability Forum

Dear Yung-jin & Elba,

❝ ❝ I was wondering if there is any initiative in this forum to develop "public data reference" in order to evaluate the performance of BEAR (or another package).

❝ Sounds do-able. May need more opinions and time to plan it. There are plenty BE/BA guru members in this Forum from worldwide. They should be able to figure out how-to sooner or later.

I’m a little bit short in time, but I give you some ideas on what I have done back in the dark-ages validating my 60,000-lines Pascal code. A very important point in validation is to challenge the software with extreme data. Are errors trapped, are the results reasonable, etc.
I’ll split it up into two parts: NCA and BE.

NCA

► Data-limits

• Generate a dataset containing only the headers, but do data. Does the software crash or trap the error and come up with a useful (!) error-message?
  
• Dataset with only one line (i.e, time+concentration).
  
• Dataset with empty lines or – in R-terms – ‘NAs’.
  
• Dataset with negative or non-numeric values.

► C_max/t_max

• Two ore more identical values. Is the first one identified as C_max/t_max?

► λ_z-estimation

• Less than three decreasing values.
  
• C_z ≥ C_z-1.
  
• 0 or NA within the decreasing values.
  
• Dataset with known λ_z.
  Example: t(0, 1, 2), C(1, 0.5, 0.25) should give exactly λ_z=ln(2) and t_½=1.
  
• Simple dataset (use the AUC-example below) and calculate λ_z from the last three values manually.
  1. Y_i=ln(C_i)=(4.04305, 2.94444, 1.38629)
     
  2. t=(∑t_i)/3=13.33333 and Y=(∑Y_i)/3=2.79126
     
  3. λ_z=-∑(t_i-t) (Y_i-Y)/(t_i-t)²=0.13260 and
     t_½=ln(2)/λ_z=0.69315/0.13260=5.22729.
     
  4. If you want to calculate AUC_0-∞ based on the estimated Ĉ_z (instead of the observed one) you need also the regression’s intercept:
     ln(C₀)=Y-λ_zt=2.79126-0.13260×13.33333=4.55928
     All calculations in full precision; final result rounded to the software’s precision.
     A note for cross-validators: In PHX/WNL you can request this value in
     NCA > Options > Model Settings > [×] Intermediate Output.
     Results in the Core output:
     Intermediate Output
-------------------
Value for Lambda_z: 0.1326, and intercept: 4.5593

► Trapezoidal rules

• Dataset with only one concentration >0.
  
• Dataset with two identical concentrations (lin-up/log-down should switch to linear).
• Dataset which can easily be calculated manually (KISS!). Example:
  t   C inc? linear                ∑ lin-up/log-down                      ∑
 0  0  NA                    0   0                            0         0
 1 62  y  ( 1- 0)(62+ 0)/2= 31  31 apply linear              31        31
 2 70  y  ( 2- 1)(70+62)/2= 66  97 apply linear              97        97
 4 57  n  ( 4- 2)(57+70)/2=127 224 ( 4- 2)(70-57)/ln(70/57)=126.55518 223.55518
12 19  n  (12- 4)(19+57)/2=304 528 (12- 4)(57-19)/ln(57/19)=276.71272 500.26791
24  4  n  (24-12)( 4+19)/2=138 666 (24-12)(19- 4)/ln(19/ 4)=115.52201 615.78992
  Perform the manual calculation to the same precision as the software’s output.

► AUC_0-∞

• AUC_0-∞=AUC_t+C_z,obs/λ_z (based on the observed C_z) should give
  666+4/0.13260=696.16555 (linear trapezoidal; extrapolated 4.3331%) and
  615.78992+4/0.13260=645.95547 (lin-up/log-down; extrapolated 4.6699%).
  
• For the second variant we need the predicted C_z,pred=ℯ^{ln(C₀)–λ_zt_z}=ℯ^{4.55928–0.13260×24}=3.96238. Then
  666+3.96238/0.13260=695.88220 (linear trapezoidal; extrapolated 4.2941%) and
  615.78992+3.96238/0.13260=645.67212 (lin-up/log-down; extrapolated 4.6281%).

BE (coming sooner or later…)

Dear Yung-jin,

❝ ❝ Note that NIST’s datasets do not necessarily represent the “truth”, only the “best-available” solutions. […]

❝ I see. But it has been good enough for comparing nonlinear regression software.

Sure. It’s nice to have results obtained in 128-bit precision. Important for ‘flat’ SSQ-surfaces where rounding may be an issue.

❝ ❝ Why? []

❝ […] The textbook or published papers may have final results in Tables for their dataset, but do not have whole detailed output results.

Agree. The most detailed one is #1 from this post.

❝ ❝ BTW: WinNonlin Pro 4.01, but which versions of bear and SAS?

❝ bear was v2.4.0 (p.92) and SAS should be v9.1.3 I guessed when I was at KMU (university license).

THX!

Hi ElMaestro,

❝ Looks to me like their examples also include a Potvin B run. I am not a veterinarian but ... I think they got it wrong as they plugged in the observed GMR from stage 1 and not 0.95 for stage 2 dimensioning.

Yes, that’s weird. Subsection II.D.3 of the guidance smells of copypasting from EMA’s GL (contrary to Method C which is generally suggested by the FDA). Furthermore,

The plan to use a two-stage approach should be pre-specified in the protocol along with the number of animals to be included in each stage and the adjusted signi­fi­cance levels to be used for each of the analyses.

I beg your pardon?

Figure 1 of the supp. is Potvin’s B:

[…] sample size based on variance stage 1…

But in the paragraph below:

[…] sample size based on the information derived at Stage 1.

But let’s continue:

SCENARIO: For the sake of this example, we will use the following Stage 1 assump­tions:

• We estimate a residual error of 15% CV and a ratio of the test/reference means of 0.90.
  
• We will conduct Stage 1 with 20 subjects (10 per sequence).

Nice to call animals subjects, but these guys are adventurous, aren’t they? Let’s assume they don’t know that with a T/R ratio of 0.90 an α_adj of 0.0294 is not sufficient any more (covering the grid of n₁ 12…60, CV 10…100%). My magic code tells me…

For adj. α of 0.0294 a maximum inflation
of 0.053753 is detected at CV=20% and n1=12.

… and suggests to use 0.0272 instead.*



However, let’s believe in Mr Pocock.

library(PowerTOST)
sampleN.TOST(alpha=0.0294, CV=0.15, theta0=0.9, design="2x2x2")
+++++++++++ Equivalence test - TOST +++++++++++
            Sample size estimation
-----------------------------------------------
Study design:  2x2 crossover
log-transformed data (multiplicative model)

alpha = 0.0294, target power = 0.8
BE margins        = 0.8 ... 1.25
Null (true) ratio = 0.9,  CV = 0.15

Sample size (total)
 n     power
26   0.802288

power.TOST(alpha=0.0294, CV=0.15, theta0=0.9, n=20)
[1] 0.6850327

Brilliant! They must love a second stage. Some more stuff about inflation & power:

library(Power2Stage)
power.2stage(method="B", alpha=c(0.0294, 0.0294), n1=20,
   GMR=0.9, CV=0.15, targetpower=0.8, pmethod="nct",
   usePE=FALSE, Nmax=Inf, theta0=1.25, nsims=1e6)

Method B: alpha (s1/s2)= 0.0294 0.0294
Futility criterion Nmax= Inf
CV= 0.15; n(stage 1)= 20; GMR= 0.9
BE margins = 0.8 ... 1.25
GMR= 0.9 and mse of stage 1 in sample size est. used

1e+06 sims at theta0= 1.25 (p(BE)='alpha').
p(BE)   = 0.038733
p(BE) s1= 0.029144
pct studies in stage 2= 74.27%

Distribution of n(total)
- mean (range)= 27.4 (20 ... 80)
- percentiles
 5% 50% 95%
 20  26  42

With this CV, T/R, n₁, and 0.0294 there will be no inflation. But that’s not true for other combinations. Overall power ~85% (~27 of studies will proceed to the second stage).

What if they want to adapt for the T/R of stage 1 (in the code above switch to usePE=TRUE)?

Method B: alpha (s1/s2)= 0.0294 0.0294
Futility criterion Nmax= Inf
CV= 0.15; n(stage 1)= 20; GMR= 0.9
BE margins = 0.8 ... 1.25
PE and mse of stage 1 in sample size est. used

1e+06 sims at theta0= 1.25 (p(BE)='alpha').
p(BE)   = 0.047555
p(BE) s1= 0.029144
pct studies in stage 2= 36.39%

Distribution of n(total)
- mean (range)= 203753.4 (20 ... 4375926130)
- percentiles
  5%  50%  95%
  20   20 5932

Still no inflation, but the sample sizes ring the alarm bell. Old story. Full adaption rarely ‘works’ in BE. This time the entire output for power:

Method B: alpha (s1/s2)= 0.0294 0.0294
Futility criterion Nmax= Inf
CV= 0.15; n(stage 1)= 20; GMR= 0.9
BE margins = 0.8 ... 1.25
PE and mse of stage 1 in sample size est. used

1e+05 sims at theta0= 0.9 (p(BE)='power').
p(BE)   = 0.948
p(BE) s1= 0.68365
pct studies in stage 2= 26.73%

Distribution of n(total)
- mean (range)= 17250.8 (20 ... 1117221800)
- percentiles
 5% 50% 95%
 20  20 160

Oops! In the worst case almost the entire population of India? In the spirit of Karalis / Macheras let’s add a futility criterion limiting the total sample size. OK, why not 100? No inflation; power:

Method B: alpha (s1/s2)= 0.0294 0.0294
Futility criterion Nmax= 100
CV= 0.15; n(stage 1)= 20; GMR= 0.9
BE margins = 0.8 ... 1.25
PE and mse of stage 1 in sample size est. used

1e+05 sims at theta0= 0.9 (p(BE)='power').
p(BE)   = 0.86378
p(BE) s1= 0.68365
pct studies in stage 2= 18.29%

Distribution of n(total)
- mean (range)= 27.3 (20 ... 100)
- percentiles
 5% 50% 95%
 20  20  68

Might work. Still not what the guidance wants – pre-specified sample size in both stages. Will they accept a maximum total sample size as a futility criterion? Does one have to perform the second stage in the pre-specified sample size even if the calculated one is lower?

---

• Validation of chosen adjusted α (0.0272):
“Method B”, assumed T/R 90%, target power 80%
           12       20       24       36       48       60
0.1  0.031667 0.027057 0.026990 0.027151 0.027106 0.027074
0.15 0.046101 0.036577 0.032393 0.027345 0.027106 0.027074
0.2  0.049838 0.046210 0.043916 0.036062 0.029193 0.027194
0.3  0.039896 0.048036 0.049104 0.047415 0.045185 0.042275
0.4  0.030938 0.035728 0.040070 0.048512 0.048743 0.047437
0.5  0.028471 0.029026 0.030577 0.038588 0.046841 0.049009
0.6  0.027897 0.027882 0.028254 0.030001 0.036409 0.044488
0.7  0.027576 0.027528 0.027789 0.027686 0.029420 0.034456
0.8  0.027430 0.027329 0.027578 0.027432 0.027694 0.028885
0.9  0.027352 0.027308 0.027633 0.027173 0.027348 0.027694
1    0.027275 0.027276 0.027640 0.027302 0.027268 0.027485
With an adj. α of 0.0272 a maximum inflation of 0.049838 is seen at CV=20% and n1=12.
Significance limit for 1e+06 simulations: 0.05036; no significant inflation seen.
The chosen adjusted α is justified within the assessed range of stage 1 sample sizes (12-60) and CVs (10-100%).

          12      20      24      36      48      60
0.1  0.88303 0.94968 0.97642 0.99834 0.99992 1.00000
0.15 0.83337 0.84816 0.85456 0.91422 0.96897 0.99007
0.2  0.81034 0.82734 0.83347 0.83663 0.85491 0.90205
0.3  0.77531 0.80634 0.80897 0.82105 0.82504 0.82862
0.4  0.75923 0.78497 0.79476 0.81014 0.81775 0.82105
0.5  0.75205 0.77630 0.78163 0.79811 0.80779 0.81514
0.6  0.74976 0.77267 0.77753 0.78933 0.79713 0.80734
0.7  0.74873 0.77240 0.77691 0.78657 0.79120 0.79915
0.8  0.74768 0.77115 0.77602 0.78450 0.78980 0.79423
0.9  0.74564 0.77041 0.77521 0.78571 0.78838 0.79215
1    0.74505 0.77060 0.77473 0.78502 0.78649 0.79281
With an adj. α of 0.0272 a minimum power of 0.74505 is seen at CV=100% and n1=12.

...additional info: to make the R script in the .rtf file work, release one of the two options for mse from the comments

Hi Yung-jin,

❝ […] Some posts about differences of replicate BE data analysis using WNL or SAS can be found in this forum.

❝ […] Still remains unsolved with R.

So far I do not know how to set up reference-scaling (FDA+EMA) in any software except SAS and Phoenix/WinNonlin. I’m not overly optimistic for R (see here; still current as of July 2014).