SAS and R (?) for variability comparison (FDA NTID Guidance) [General Sta­tis­tics]

posted by Shuanghe  – Spain, 2015-09-25 18:48 (2139 d 06:10 ago) – Posting: # 15471
Views: 26,520

(edited by Shuanghe on 2015-09-25 19:16)

Dear forum members,

Recent discussion on the new FDA guidance of dabigatran and rivaroxaban led me to the statistical method mentioned in the warfarin Na guidance.

It seems the method is almost the same as the one for HVDP in progesterone guidance, with difference regulatory constant and criteria of course. With the sas code example there, it's not difficult to calculate the 95% upper limit for BE evaluation.

However, there's no sas code for the variability comparison so here are my questions.

90% CI for σWTWR was expressed as {(sWT/sWR) / SQRT(F(0.05, νT, νR)), (sWT/sWR) / SQRT(F(0.95, νT, νR))}
  1. First silly question: when expressed like this I just assumed the left (blue colour) is the lower limit while the right (red colour) is the upper limit that's more important. But when I tried several data sets I always had higher value with the blue terms. Did I do anything wrong or it suppose to be so? I'm not a statistician so I don't know if F(0.05, v1, v2) is always < F(0.95, v1, v2) when v1 and v2 are kept constant.

  2. From the sas code for sWR (from Dij=R1-R2 then running with PROC MIXED, output CovParms, and calculate sWR = SQRT(estimate/2)), is it correct to assume that we need to modify the code to have something like DijT = T1 - T2, then run the same procedure as for sWR to calculate sWT?

  3. The unscaled sas code (method C in EMA's guideline) will output ANOVA residuals for both T and R, which were usually be used to calculate ISCV (ISCV=SQRT(EXP(residual)-1)). If so, it seems sWR and sWT can also be calculated use this approach. I tried this as well and the result is slight different from approach mentioned in point 2 (see below example). So which one should we use?

  4. Since FDA request that only subject who complete all 4 periods should be included in the analysis, for obtaining F-distribution value, it seems the degree of freedom vT and vR will always be equal, which is always Ntotal - 2 for full replicate. Right?

  5. I tried with EMA's sample data set I which is full replicate. There are some subjects with incomplete periods: subject 11, 20, 24, 31, 42 and 69 complete 3 periods, subject 67 and 71 complete only 2 periods. So I removed them from the data set (69 subjects left). The results are:
    • sWR = 0.45168
    • sWT = 0.34444
    • 90% CI lower = 0.62286
    • 90% CI upper = 0.93363

  6. If I use unscaled method as mentioned in point 3, then the results are:
    • residual R = 0.2097; sWR = SQRT(residual) = 0.45795
    • residual T = 0.1202; sWT = SQRT(residual) = 0.34677
    • 90% CI lower = 0.61848
    • 90% CI upper = 0.92708

  7. Why different? And, though the difference is small but in borderline case, one might be > 2.5 while the other might be < 2.5. So it's better to be clear which one we should approach. I prefer the former one.

    I put data set here. To be more precise I didn't use the last colum logpk. instead, I calculate LOG(pk) in the program.

    Could anyone check it and please let me know if there is anything wrong here? Thanks.
    What I did is add a line to calculate mean difference between Test
    /* Iij and Dij. data has previously sorted by seq and subj */
    DATA &pk._sabe;
      MERGE &pk._t1 &pk._t2 &pk._r1 &pk._r2;
      BY seq subj;
        &pk._difftr = 0.5 * (&pk._t1 + &pk._t2 - &pk._r1 - &pk._r2);
        &pk._diffr = &pk._r1 - &pk._r2;  /* <-- for s2wr as in the guidance sas example*/
        &pk._difft = &pk._t1 - &pk._t2;  /* <-- new line for s2wt later */

    and then run an additional PROC MIXED

    /* Intermediate analysis of Dij(T), the mean diff between T1 and T2. */
    PROC MIXED DATA = &pk._sabe;
        CLASS seq;
        MODEL &pk._difft = seq / DDFM = SATTERTH;
        ESTIMATE "&pk" INTERCEPT 1 seq 0.5 0.5 / E CL ALPHA = 0.1;
        ODS OUTPUT COVPARMS = &pk._difft_out1;
        ODS OUTPUT ESTIMATES = &pk._difft_out2;
        ODS OUTPUT NOBS = &pk._difft_out3;
        ODS OUTPUT TESTS3 = &pk._difft_out4;
    TITLE6 "Scaled Average BE, Intermediate Analysis of Dij(T), PROC MIXED";
    TITLE7 " "; 

    DATA &pk._difft_out1;
        SET &pk._difft_out1;
        s2wt = ESTIMATE / 2;
        KEEP s2wt;

    and then calculate 90% interval for variability comparison

    DATA &pk._final;
        MERGE &pk._difftr_out &pk._difft_out1 &pk._difft_out2 &pk._diffr_out1 &pk._diffr_out2;
        theta = (LOG(1 / 0.9) / 0.1)**2;
        y = -theta * s2wr;
        boundy = y * dfr / CINV(0.95, dfr);
        swr = SQRT(s2wr);
        swt = SQRT(s2wt);
        critbound = (x + y) + SQRT((boundx - x)**2 + (boundy - y)**2);
        varlower = (swt / swr) / SQRT(FINV(0.05, dft, dfr)); /* <-- actually upper */
        varupper = (swt / swr) / SQRT(FINV(0.95, dft, dfr)); /* <-- actually lower */


  8. Based on the description of the method in the guidance, it seems the analysis can be done in R. Has anyone tried it yet? I'll give it a try and report the result later.

  9. Helmut, this is for you. It seems LaTeX expression can not be used to write math equation here. Is it possible to implement it? Since there are so many members here and some times one has to write complicate formula it might be a good idea to be able to write equations in LaTeX. :-D

OK, that's for now. I wish you all a wonderful weekend!

All the best,

Complete thread:

 Admin contact
21,606 posts in 4,517 threads, 1,532 registered users;
online 19 (0 registered, 19 guests [including 4 identified bots]).
Forum time: Wednesday 00:58 CEST (Europe/Vienna)

Scientists cannot simply hang their subjectivities
up on a hook outside the laboratory door.    Ruth Bleier

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz