acfalcao
●    

2007-09-16 22:34

Posting: # 1082
Views: 15,979
 

 Three-way crossover BABE Studies [Design Issues]

Dear All

My name is Amilcar Falcão and I am professor at the University of Coimbra, Portugal. I am looking for help regarding the execution and interpretation of three-way cross-over bioequivalence trials (for academic purposes).

I would like to have an example (explicit) where a three-way crossover study is appropriately analysed. Namely, I would like to know what kind of ANOVA should be performed and how the 90%CI should be calculated for the different combinations of study formulations (i.e. A vs B and A vs C). Perhaps this is a naive question, but in fact I am trying to look for it in literature and I have seen different approaches always without the adopted statistical pathway. Therefore, I would like to be sure regarding the correct way to deal with.

I strongly appreciate if you could help me about this matter.

Thank you for your attention.

Best regards; :-)

Amílcar Falcão
Laboratory of Pharmacology
Faculty of Pharmacy
University of Coimbra
Coimbra - Portugal
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2007-09-17 20:47

@ acfalcao
Posting: # 1088
Views: 15,182
 

 Three-way crossover example data set

Dear Amilcar!

» I would like to have an example (explicit) where a three-way crossover study is appropriately analysed.

You have to use a Williams’ design (three period, six sequences); the topic was covered in previous threads.

A 6×3 design is needed in order to ‘extract’ two 2×2 tables, which are also balanced. Although the full 6×3 table will be used in the analysis of AUC and Cmax, you will need these 2×2s for the nonparametric analysis of tmax (un­fortunately there’s no confidence interval based nonparametric method available for more than two formulations/periods). The asterisks * denote pseudo-sequences and pseudo-periods, e.g. P1* means only that the treatment was given in a period prior to P2* – irrespective of the true study period:
+----+------------+  -->  +----+--------+  and  +----+--------+
|    | P1  P2  P3 |       |    | P1* P2*|       |    | P1* P2*|
+----+------------+       +----+--------+       +----+--------+
| S1 | T   R1  R2 |       | S1*| T   R1 |       | S1*| T   R2 |
| S2 | R1  R2  T  |       | S2*| R1  T  |       | S2*| R2  T  |
| S3 | R2  T   R1 |       | S1*| T   R1 |       | S2*| R2  T  |
| S4 | T   R2  R1 |       | S1*| T   R1 |       | S1*| T   R2 |
| S5 | R1  T   R2 |       | S2*| R1  T  |       | S1*| T   R2 |
| S6 | R2  R1  T  |       | S2*| R1  T  |       | S2*| R2  T  |
+----+------------+       +----+--------+       +----+--------+
                            ^   balanced          ^   balanced

A common mistake is to design the study as a set of 3×3 latin squares, which will lead (especially if the sample size is small and in the case of drop outs) to extremely imbalanced data sets:
+----+------------+  -->  +----+--------+  and  +----+--------+
|    | P1  P2  P3 |       |    | P1* P2*|       |    | P1* P2*|
+----+------------+       +----+--------+       +----+--------+
| S1 | T   R1  R2 |       | S1*| T   R1 |       | S1*| T   R2 |
| S2 | R1  R2  T  |       | S2*| R1  T  |       | S2*| R2  T  |
| S3 | R2  T   R1 |       | S1*| T   R1 |       | S2*| R2  T  |
+----+------------+       +----+--------+       +----+--------+
                            ^ imbalanced          ^ imbalanced


» Namely, I would like to know what kind of ANOVA should be performed and how the 90%CI should be calculated for the different combinations of study formulations (i.e. A vs B and A vs C).

For the design see Chapter 10 of

Chow S-C, Liu J-p.
Design and Analysis of Bioavailability and Bioequivalence Studies. New York: Marcel Dekker; 2nd ed. 2001, p. 302–32.

For a detailed discussion of variance balanced designs see Chapter 4 of

Jones B, Kenward MG.
Design and Analysis of Cross-over Trials. Boca Raton: Chapman & Hall/CRC; 2nd ed. 2003, p. 151–204.

For an example data set see Chapter 4 of

Patterson S, Jones B.
Bioequivalence and Statistics in Clinical Pharmacology. Boca Raton: Chapman & Hall/CRC; 2006, p. 79–132.

Their example 4.5 matches your design (one test formulation is compared to two reference formulations)
You may download zipped datasets and programs (SAS/S+) from CRC’s website. If you don’t have access to commercial software, S+ code will run with open-source R with minor modifications.

Patterson/Jones give results in Table 4.12 (p. 105) with:
T/R (% Test vs. Reference 1)
+----------+-------+---------------+
| Endpoint |  PE   |     90% CI    |
+----------+-------+---------------+
| AUC      | 116.2 | 109.0 , 123.9 |
| Cmax     | 130.0 | 119.1 , 141.8 |
+----------+-------+---------------+

T/S (% Test vs. Reference 2)
+----------+-------+---------------+
| Endpoint |  PE   |     90% CI    |
+----------+-------+---------------+
| AUC      |  82.8 |  77.6 ,  88.3 |
| Cmax     |  81.5 |  74.7 ,  89.0 |
+----------+-------+---------------+


WinNonlin 5.2 comes up with:
T/R (% Test vs. Reference 1)
+----------+-------+---------------+
| Endpoint |  PE   |     90% CI    |
+----------+-------+---------------+
| AUC      | 116.2 | 109.0 , 123.8 |
| Cmax     | 129.7 | 118.4 , 141.5 |
+----------+-------+---------------+

T/S (% Test vs. Reference 2)
+----------+-------+---------------+
| Endpoint |  PE   |     90% CI    |
+----------+-------+---------------+
| AUC      |  82.6 |  77.5 ,  88.1 |
| Cmax     |  81.2 |  74.3 ,  88.6 |
+----------+-------+---------------+


Results are slightly different (although both WinNonlin and SAS use GLM – not ANOVA; implementation, rounding, etc. is different – and of course ‘proprietary information’ and not documented). I assume Patterson/Jones' results were produced by SAS; I will check the results from their S+ code the next days.

» Perhaps this is a naive question,...

Not at all; little is published - and no worked examples at all.

There’s another point which is a little bit tricky: multiplicity.
If you are testing one test (A) against two references (B, C), any impression of ‘data dredging’ must be avoided, e.g., calculation 90% CIs of A/B and A/C - and only picking out the best getting an approval.
Since in the EU the Innovator's product from any European country may be used as the reference, you may run into problems (A vs B is BE, whereas A vs C is not BE). It may be wise to use 95% CIs instead (overall Bonferroni-corrected patient’s risk: α = 1–(1–0.05/k)k, where k is the number of simultaneous comparisons).
95% CI should also be applied in testing for dose proportionality of three dose levels (or 96.67% for four levels).
IMHO the only case where 90% CIs should be used is the comparison of two test formulations against one reference, and only one of them will be further used in the approval process.

Cheers,
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. ☼
Science Quotes
acfalcao
●    

2007-09-17 23:02

@ Helmut
Posting: # 1089
Views: 14,365
 

 Three-way crossover example data set

Dear Helmut

Thank you very much for your help. The correct codification of "sequences/periods" is the key for this problem. Unfortunately the link to download the SAS files was not working. I hope it can be solve in a near future.

Best regards;

Amílcar ;-)


Edit: Full quote removed. Please delete everything from the text of the original poster which is not necessary in understanding your answer; see also this post #5[Helmut]
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2007-09-18 12:47

@ acfalcao
Posting: # 1092
Views: 14,499
 

 Three-way crossover example data set

Dear Amilcar!

» Thank you very much for your help. The correct codification of "sequences/periods" is the key for this problem. Unfortunately the link to download the SAS files was not working. I hope it can be solve in a near future.

Please try it again (it worked yesterday, and today as well): CRC’s website
Having downloaded the file C5300.zip, unzip it keeping folder's structure; in folder [chapter4] you will find:
Example4.5.sas (the SAS code)
Example45.SSC (the S+ code)
exam45.dat (the raw data for S+)

As I get from the SAS-code a mixed model (degrees of freedom according to Kenward) was used, whereas WinNonlin uses Satterthwaite's method.
The S+ code does not calculate BE, but is used for the example’s figures only.


Edit: Wrong link corrected in both posts. [Helmut]

Cheers,
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. ☼
Science Quotes
Nirali
★    

India,
2007-09-21 06:49

@ Helmut
Posting: # 1106
Views: 14,339
 

 Three-way crossover (WinNonlin)

Dear Helmut,

Which hypothesis Winnonlin will use while running BE for 2 test (T1 & T2) vs a reference (R)?
Would it be "MEAN(T1)=MEAN(T2)=MEAN(R)"?
In threeway crossover study, does we try to compare Equivalence of all three product together?

Thanks & Regards,
Nirali
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2007-09-21 13:04

@ Nirali
Posting: # 1109
Views: 14,366
 

 Three-way crossover (WinNonlin)

Dear Nirali,

» Which hypothesis Winnonlin will use while running BE for 2 test (T1 & T2) vs a reference (R)?
» Would it be "MEAN(T1)=MEAN(T2)=MEAN(R)"?
» In threeway crossover study, does we try to compare Equivalence of all three product together?

Comparisons are based on paired differences; IMHO a simultaneous test of T1=T2=R is not possible.
WinNonlin always uses one reference in a run, i.e., T1=R (one test), and T2=R (the second test).
If you are interested in T1=T2 you have to start a second run (defining T2 as the reference), and will obtain two tests (T1=T2 and R=T2) within.

BTW, why are you interested in T1=T2?

Cheers,
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. ☼
Science Quotes
Nirali
★    

India,
2007-09-25 08:05

@ Helmut
Posting: # 1127
Views: 14,310
 

 Three-way crossover (WinNonlin)

Dear Helmut,

Actully, I was under impression that WinNonlin is using the hypo:
"MEAN(Test1)=MEAN(Test2)=MEAN(Ref)".


As I understood from your reply--
(A) suppose we have two REFERENCES and single TEST [R1, R2, T]for threeway comparison: if we select R1 as REFERENCE, WNL will calculate BE with hypo:
(1) R1=T and (2) R1=R2.

same way,
(B) for two TEST and single REFERENCE [T1, T2, R]: if we select R as standard, hypo will be: (1) R=T1 and (2)R=T2

pls do correnct me if i am wrong to understand.

Thanks & Regards,
Nirali
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2007-09-25 13:15

@ Nirali
Posting: # 1132
Views: 14,422
 

 Three-way crossover (WinNonlin)

Dear Nirali!

» Actully, I was under impression that WinNonlin is using the hypo:
» "MEAN(Test1)=MEAN(Test2)=MEAN(Ref)".

Again, we obtain only pairwise comparisons, not simultaneous ones.

» (A) suppose we have two REFERENCES and single TEST [R1, R2, T]for threeway comparison: if we select R1 as REFERENCE, WNL will calculate BE with hypo:
» (1) R1=T and (2) R1=R2.

Not quite, but very close. ;-)
If you run the BE wizard (R1=reference), you will get
(1a) T vs. R1, and
(1b) R2 vs. R1.
You have to run the BE wizard again (this time with R2 as reference), getting
(2a) T vs. R2, and
(2b) R1 vs. R2.
In order to get comparisons of T=R1 and T=R2, you have to run the wizard twice (results 1a and 2a), omitting 1b/2b. Of course PE and CL of 2b=1/1b.

» same way,
» (B) for two TEST and single REFERENCE [T1, T2, R]: if we select R as standard, hypo will be: (1) R=T1 and (2)R=T2

Yes this is correct; in this case only one run of the BE wizard is necessary, getting:
(1a) T1 vs. R, and
(1b) T2 vs. R.

Cheers,
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. ☼
Science Quotes
Irene_I
☆    

Indonesia,
2018-06-07 11:09

@ Helmut
Posting: # 18861
Views: 6,561
 

 Three-way crossover example data set

Dear Helmut,

I'm currently looking for statistics methods for William Design in Bioequivalence Study and I found this thread. I tried to download the SAS code from the link you've attached, but apparently the link doesn't work for now.

» Please try it again (it worked yesterday, and today as well): CRC's website
» Having downloaded the file C5300.zip, unzip it keeping folder's structure; in folder [chapter4] you will find:
» Example4.5.sas (the SAS code)
» Example45.SSC (the S+ code)
» exam45.dat (the raw data for S+)

I know that this topic was discussed years ago, but I think I really need to download these files for my further understanding about this topic. I would really appreciate if you could help me with this. Thank you. :-)

Best Regards,

Irene I
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2018-06-07 13:04

@ Irene_I
Posting: # 18862
Views: 6,567
 

 Three-way crossover example data set

Hi Irene,

» I know that this topic was discussed years ago, but I think I really need to download these files for my further understanding about this topic. I would really appreciate if you could help me with this.

The files of the 2005 edition are gone. However, here is the link for the 2016 edition.

Cheers,
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. ☼
Science Quotes
Irene_I
☆    

Indonesia,
2018-06-08 11:09

@ Helmut
Posting: # 18869
Views: 6,528
 

 Three-way crossover example data set

Dear Helmut,

Thanks a lot! It is really a great help to me.. :-)

Regards,

Irene I
Irene_I
☆    

Indonesia,
2018-06-12 09:21

@ Helmut
Posting: # 18886
Views: 6,441
 

 Three-way crossover example data set

Dear Helmut,

Could I ask you some questions about this example?

After I tried to run the SAS code (Example 4.5), I got the CI, but there is only one Residual value. I guess I suppose to calculate the intra-subject CV from this value, right? So, My question is, Do we only get one Intra-subject CV from three-comparison of the treatment (T-R, S-T, S-R) while we got three CI (T-R, S-T, S-R)?

I am looking forward for your help.
Thank you very much.

Best Regards,


Irene I
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2018-06-12 12:42

@ Irene_I
Posting: # 18888
Views: 6,491
 

 Leave-One-Out (IBD)

Hi Irene,

» After I tried to run the SAS code (Example 4.5), I got the CI, but there is only one Residual value. I guess I suppose to calculate the intra-subject CV from this value, right? So, My question is, Do we only get one Intra-subject CV from three-comparison of the treatment (T-R, S-T, S-R) while we got three CI (T-R, S-T, S-R)?

I’m not equipped with [image] … Hence, my understanding of the code is limited.
You are right, in the model you get only one (pooled) estimate of the variance. That’s not a good idea, since it will “work” only if intra-subject variances would be identical and the treatment differences at least very similar. Otherwise, the treatment estimates will be biased and the type I error is not controlled. See also this post.
At the 2nd Conference of the Global Bioequivalence Harmonisation Initiative (Rockville, September 2016) Pina D’Angelo gave a presentation “Testing for Bioequivalence in Higher‐Order Crossover Designs: Two‐at‐a‐Time Principle Versus Pooled ANOVA” showing exactly that. Here some of her slides:

Purpose

  • To determine which method of statistical analysis
    is more appropriate to conclude bioequivalence in
    higher‐order crossover studies:
    1. The two‐at‐a‐time principle using two separate
      incomplete block design ANOVAs
    2. A pooled approach using one ANOVA and a common
      error term for the two contrasts

Introduction


Statistical Concerns:

  1. Different means (point estimates) between formulations
  2. Different variances between formulations
  • If either situations exist, which method of analysis
    reduces bias the most: two‐at‐a‐time principle or
    pooled ANOVA?

Simulated Data: Summary

  • When all three treatments have similar means and there is
    homogeneity of variances, both methods give very similar
    results.
  • When treatment means differ but there is homogeneity of
    variances, both methods give very similar results. With higher
    variability, the power is slightly increased when using the pooled
    ANOVA method.
  • When treatment means are similar but variances are not
    homogeneous, the two‐at‐a‐time method gives higher power to
    detect BE for the treatment with lower variability
  • When treatment means differ and variances are not
    homogeneous, the two‐at‐a‐time method increases power to
    detect BE for the treatment with lower variability. Moreover,
    type I error is higher when using the pooled ANOVA method.

Closing Remarks

  • Using a two‐at‐a‐time principle for statistical analysis of a
    higher‐order pilot study will have more value for decision-
    making on which multiple tests lots will be selected for use
    in a pivotal study based on the pilot study results. The
    intra‐subject variability of a specific test‐to‐reference
    comparison can be determined using the two‐at‐a‐time
    principle, which may be an important factor in selecting a
    test product considering test products are generally
    formulated to show different characteristics for testing in a
    pilot study.
(my emphasis)

If you are interested in all pairwise comparisons, generate three data sets (values or T&R, S&R, T&S, i.e., exclude all values of the respective other treatment S, T, R) whilst keeping the codes for sequence and period. This will give you data sets which represent an IBD (incomplete block design). Run the usual model on them.
This approach is also recommended in the EMA’s BE-GL:

In studies with more than two treatment arms (e.g. a three period study including two references, one from EU and another from USA […], the analysis for each comparison should be conducted excluding the data from the treatments that are not relevant for the comparison in question.

(my emphasis)
Example 4.5 (Phoenix/WinNonlin 8, mixed effects, Satterthwaite’s degrees of freedom):
                           PE (%)    90% CI (%)       s²w    CVw (%)
pooled model T vs. R  AUC  116.15 (108.97, 123.81)  0.043954  21.20
             T vs. R  Cmax 129.65 (118.84, 141.45)  0.084569  29.71
             S vs. R  AUC  140.63 (131.93, 149.90)  0.043954  21.20
             S vs. R  Cmax 159.75 (146.42, 174.30)  0.084569  29.71
             T vs. S  AUC   82.60 ( 77.46,  88.07)  0.043954  21.20
             T vs. S  Cmax  81.16 ( 74.35,  88.56)  0.084569  29.71
IBD models   T vs. R  AUC  116.05 (108.92, 123.65)  0.042525  20.84
             T vs. R  Cmax 129.54 (119.26, 140.71)  0.074744  27.86
             S vs. R  AUC  141.09 (131.39, 151.51)  0.053706  23.49
             S vs. R  Cmax 160.48 (145.20, 177.36)  0.109492  34.02
             T vs. S  AUC   83.12 ( 78.37,  88.15)  0.035788  19.09
             T vs. S  Cmax  81.63 ( 75.33,  88.46)  0.069372  26.80

It’s clear that the variances are not identical. The ones of T/R are smaller than the ones of S/R. If we apply the pooled model the CIs of T/R will be wider than in the IBD model and the CIs of S/R narrower.

BTW, I don’t understand what the purpose of the carry variable in Byron’s code is.

Cheers,
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. ☼
Science Quotes
Irene_I
☆    

Indonesia,
2018-06-13 05:02

@ Helmut
Posting: # 18896
Views: 6,336
 

 Leave-One-Out (IBD)

Dear Helmut,

Thank you for your explanation. I am sure this will be helpful not only for me, but also others who looking for information about this topic too. :-)

Best Regards,

Irene I
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2018-06-13 15:02

@ Irene_I
Posting: # 18898
Views: 6,276
 

 Impact of pooled variance (bias, CI)

Hi Irene,

» Thank you for your explanation.

You are welcome. A picture tells more than a thousand words (AUC only). ;-)

[image]


The PE and its CI of T/R are similar. The small bias is caused by the fact that the average PE of S/R and T/S (√141.09×83.12=108.29) is close to the T/R with 116.05. The other comparisons show negative biases. The CI of S/R is narrower in the pooled model and the CI of T/S is wider.
We don’t want to go there.

Cheers,
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. ☼
Science Quotes
d_labes
★★★

Berlin, Germany,
2018-06-13 15:31

@ Helmut
Posting: # 18899
Views: 6,249
 

 carry (over?)

Dear Helmut,

» BTW, I don’t understand what the purpose of the carry variable in Byron’s code is.

As far as I see it: It serves for nothing because it is not used except in the class statement of Proc Mixed.
The class statement serves the purpose of declaring factors even if a variable is numeric. Nothing more.

The carry variable seems to me a left over from old times where carry-over was an issue.

Regards,

Detlew
ElAlumno
☆    

2019-03-14 23:40

@ Helmut
Posting: # 20036
Views: 2,920
 

 Pseudo-periods

Dear Helmut,

In your old but still-useful post (excerpt below), you describe generating pseudo-sequences and pseudo-periods when assessing ABE.

» A 6×3 design is needed in order to ‘extract’ two 2×2 tables, which are also balanced. Although the full 6×3 table will be used in the analysis of AUC and Cmax, you will need these 2×2s for the nonparametric analysis of tmax (unfortunately there’s no confidence interval based nonparametric method available for more than two formulations/periods). The asterisks * denote pseudo-sequences and pseudo-periods, e.g. P1* means only that the treatment was given in a period prior to P2* – irrespective of the true study period:

Your post preceded the EMA's Guideline on the Investigation of Bioequivalence by a few years, but it appears to agree with their stipulation:

In studies with more than two treatment arms (e.g. a three period study including two references, one from EU and another from USA, or a four period study including test and reference in fed and fasted states), the analysis for each comparison should be conducted excluding the data from the treatments that are not relevant for the comparison in question.

I completely understand the use of pseudo-sequences in this situation. But what about pseudo-periods? The EMA guidelines aren't too specific.

On the one hand, if you wish to display your data with pseudo-sequences, you would want to use pseudo-periods too (so that the timing of treatments is consistent within pseudo-sequence). On the other hand, the ANOVA doesn't care whether your pseudo-sequences are logically consistent with your period.

If there is no period effect, the difference between using period and pseudo-period will be minimal (use a few more df in the ANOVA table with period than pseudo-period). But what if there is a serious period effect? You would miss it (wrong sum of squares for period) by using pseudo-periods.

Opinions?
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2019-03-15 00:47

@ ElAlumno
Posting: # 20037
Views: 2,900
 

 Pseudo-periods

Hi ElAlumno (nice nick!),

» In your old but still-useful post (excerpt below), you describe generating pseudo-sequences and pseudo-periods when assessing ABE.
»
» » A 6×3 design is needed in order to ‘extract’ two 2×2 tables, [:blahblah:]
»
» Your post preceded the EMA's Guideline on the Investigation of Bioequivalence by a few years, but it appears to agree with their stipulation: […]

Not quite.

» I completely understand the use of pseudo-sequences in this situation. But what about pseudo-periods? The EMA guidelines aren't too specific.

The EMA means the ‘Two‐at‐a‐Time Principle’ (see above) and not my recoding (which I abandoned a good while ago as well).
David Brown (MHRA, member the EMA’s Biostatistics Working Party) gave a clarification at the 3rd EGA Symposium on Bioequivalence (London, June 2010). The Q&A document states:

However, the treatment, groups, sequences and periods should have their original values maintained in the analysis, and not have the values modified. For example an observation made in period 3 should still be coded as period 3, not have the period changed to “2” because the results for that subject in one of the earlier periods has now be removed.

(my emphasis)

» On the one hand, if you wish to display your data with pseudo-sequences, you would want to use pseudo-periods too …

Nope. Relevant slides:

[image] GL
[image] Q
[image] Q
[image] A
[image] A


Cheers,
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. ☼
Science Quotes
ElAlumno
☆    

2019-03-22 21:59

@ Helmut
Posting: # 20069
Views: 2,494
 

 Two‐at‐a‐Time analysis in R

Thank you Helmut for your quick and detailed response!

I'm glad to hear that pseudo-periods have gone the way of the dodo. I agree with the EMA Biostatistics Working Party that it makes more sense (statistically) to maintain the original codings when leaving out data from irrelevant treatments (not to mention being less prone to human error than recoding). Also, I love their name; I bet the EMA Biostatistics Working Party could form a coalition government with the Slightly Silly Party, the Surprise Party, and the Rent Is Too Damn High Party. But I digress.

I don't have WinNonlin, but I tried analyzing Patterson & Jones' Example 4.5 with the Two‐at‐a‐Time Principle using R 1.1.456. My results were similar to the results you posted above, but not as close as I would have expected. I will post the code & table below. Do you see anything I am doing wrong or any obvious explanation of the differences? I don't know how WinNonlin handles missing values, so I tried a few variations (subjects with complete cases in the full dataset, and subjects with complete cases in each of the 3 incomplete-block subsets), but the differences persisted.

# 3-treatment design from Patterson & Jones (2017)
dta = read.table("exam45.dat", header=T)
dta[dta==99999] = NA
dta$Subj = factor(dta$subject)
dta$Per  = factor(dta$period)
dta$Seq  = factor(dta$sequence)
dta$Trt  = factor(dta$formula)
options(contrasts=c("contr.treatment","contr.poly"), digits=4)

# Data frame to store results
IBD = data.frame(Test = paste(rep(c("T vs. R", "S vs. R", "T vs. S"), each=2),
                              rep(c(" AUC", "Cmax"), times=3)),
                 PE = NA, LCI = NA, UCI = NA, CV = NA)

# AUC T vs R
muddle = lm(log(AUC)~Trt+Per+Seq+Subj, data=dta[dta$Trt!="S",])
IBD[1, 2:5] = c(100*exp(coef(muddle)["TrtT"]),
                100*exp(confint(muddle,c("TrtT"), level=.9)),
                100*sqrt(exp(summary(muddle)$sigma^2)-1))

# Cmax T vs R
muddle = lm(log(CMAX)~Trt+Per+Seq+Subj, data=dta[dta$Trt!="S",])
IBD[2, 2:5] = c(100*exp(coef(muddle)["TrtT"]),
                100*exp(confint(muddle,c("TrtT"), level=.9)),
                100*sqrt(exp(summary(muddle)$sigma^2)-1))

# AUC S vs R
muddle = lm(log(AUC)~Trt+Per+Seq+Subj, data=dta[dta$Trt!="T",])
IBD[3, 2:5] = c(100*exp(coef(muddle)["TrtS"]),
                100*exp(confint(muddle,c("TrtS"), level=.9)),
                100*sqrt(exp(summary(muddle)$sigma^2)-1))

# Cmax S vs R
muddle = lm(log(CMAX)~Trt+Per+Seq+Subj, data=dta[dta$Trt!="T",])
IBD[4, 2:5] = c(100*exp(coef(muddle)["TrtS"]),
                100*exp(confint(muddle,c("TrtS"), level=.9)),
                100*sqrt(exp(summary(muddle)$sigma^2)-1))

# AUC T vs S
muddle = lm(log(AUC)~Trt+Per+Seq+Subj, data=dta[dta$Trt!="R",])
IBD[5, 2:5] = c(100*exp(coef(muddle)["TrtT"]),
                100*exp(confint(muddle,c("TrtT"), level=.9)),
                100*sqrt(exp(summary(muddle)$sigma^2)-1))

# Cmax T vs S
muddle = lm(log(CMAX)~Trt+Per+Seq+Subj, data=dta[dta$Trt!="R",])
IBD[6, 2:5] = c(100*exp(coef(muddle)["TrtT"]),
                100*exp(confint(muddle,c("TrtT"), level=.9)),
                100*sqrt(exp(summary(muddle)$sigma^2)-1))

# Print table
print(IBD)


This produced the following table:

          Test     PE    LCI    UCI    CV
1 T vs. R  AUC 116.14 108.98 123.78 20.88
2 T vs. R Cmax 129.82 119.50 141.04 27.87
3 S vs. R  AUC 140.73 131.02 151.15 23.50
4 S vs. R Cmax 160.04 144.76 176.94 34.06
5 T vs. S  AUC  83.51  78.73  88.58 19.07
6 T vs. S Cmax  82.30  75.96  89.18 26.67
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2019-03-23 01:03

@ ElAlumno
Posting: # 20070
Views: 2,520
 

 fixed & mixed (dammit!) and a request to SASians

Hi ElAlumno,

» I agree with the EMA Biostatistics Working Party […]. Also, I love their name; I bet the EMA Biostatistics Working Party could form a coalition government with the Slightly Silly Party, the Surprise Party, and the Rent Is Too Damn High Party.

:lol3: It is really existing though David Brown will be brexited (is this already a word?) unceremoniously.

» But I digress.

That’s fine. Adds some spice to that dry stuff.

» […] I tried analyzing Patterson & Jones' Example 4.5 with the Two‐at‐a‐Time Principle using R 1.1.456.

What? 1.1.456? 1.1.4 was the first one I used in 2002. Current is 3.5.3!

» My results were similar to the results you posted above, but not as close as I would have expected. I will post the code & table below. Do you see anything I am doing wrong or any obvious explanation of the differences? I don't know how WinNonlin handles missing values,…

That’s not the point. The devil is in the details. I wrote above:

» » » » Phoenix/WinNonlin 8, mixed effects

Why? Well, Patterson & Jones analyzed the data with SAS Proc MIXED and I wanted to reproduce it in PHX/WNL (restricted maximum likelihood). OK, they used in SAS ddfm=KENWARDROGER1 and yours truly Satterthwaite’s df.
When I switched in PHX/WNL to a fixed effects model (i.e., like lm() in R and SAS Proc GLM) I got exactly your results. Hence, your code is fine though you can streamline it (specify NAs and classes):

dta <- read.table("exam45.dat", header=TRUE, na.strings="99999",
                  colClasses=c(rep("factor", 4), rep("numeric", 2))
)

Checks:

sum(is.na(dta$AUC))
[1] 6
sum(is.na(dta$CMAX))
[1] 2
str(dta)
'data.frame':   186 obs. of  6 variables:
 $ subject : Factor w/ 62 levels "001","002","004",..: […]
 $ sequence: Factor w/ 6 levels "1","2","3","4",..: […]
 $ period  : Factor w/ 3 levels "1","2","3": […]
 $ formula : Factor w/ 3 levels "R","S","T": […]
 $ AUC     : num  4089 7411 5513 2077 3684 ...
 $ CMAX    : num  906 1711 1510 504 845 ...


I would not use option(..., digits=4) cause IMHO, this may affect also the precision of internal calculations. Better to round yourself at the end, f.i.

IBD[1, 2:5] <- round(c(foo, bar, baz), 2).


So what do we have (T/R of AUC)?

                                            PE        90% CI      CV%
pooled SAS Proc MIXED (ddfm=KENWARDROGER)  116.2   109.0  123.9   21.3 
2
       SAS Proc MIXED (ddfm=SATTERTHWAITE) 116.15  108.97 123.81  21.20 ╮
       PHX/WNL REML (DF Satterthwaite)     116.15  108.97 123.81  21.20 |
       R lmer() ddf="Satterthwaite"        116.15  108.97 123.81  21.20 |
       R lmer() ddf="Kenward-Roger"        116.15  108.97 123.81  21.20 ╯
       PHX/WNL fixed                       116.24  109.04 123.92  21.22 ╮
       R lm()                              116.24  109.04 123.92  21.22 ╯
IBD    SAS Proc MIXED (ddfm=SATTERTHWAITE) 116.05  108.92 123.65  20.84 ╮
       PHX/WNL REML (DF Satterthwaite)     116.05  108.92 123.65  20.84 ╯
       R lmer() ddf="Satterthwaite"        116.05  108.97 123.58  20.84 ╮
       R lmer() ddf="Kenward-Roger"        116.05  108.97 123.58  20.84 ╯
       PHX/WNL fixed                       116.14  108.98 123.78  20.88 ╮
       R lm()                              116.14  108.98 123.78  20.88 ╯


If Satterthwaite’s df are used, the mixed effects model in SAS and PHX/WNL agree. R lm() agrees with PHX/WNL’s fixed effects model. A little bit surprising to me that Proc MIXED with ddfm=KENWARDROGER3 gives the same result as a fixed effects model. OK, but the FDA recommends Satterthwaite’s df for ages.4 A mixed effects model may give a lower residual variance – like in this case – because information of incomplete data is recovered. R lmer() is set up according to the EMA’s ‘Method B’ which treats subject as a random effect but is not a true mixed effects model. Given only for completeness (for the code see below).
A request to the SASians of the forum: Please run Proc MIXED with DDFM=SATTERTHWAITE of the pooled data and both DDFM-options on the data where S is excluded. THX to mittyri!

Mixed effects in R are a swamp. The closest I could get was this:

library(nlme)
muddle1 <- lme(log(AUC) ~ Trt + Per + Seq, data = dta[dta$Trt != "S", ],
                          random = ~Trt-1|Subj, weights = varIdent(form = ~1|Trt),
                          method = "REML", na.action = na.exclude,
                          control = lmeControl(opt = "optim", msMaxIter = 200))
anova(muddle1, type="marginal")
summary(muddle1)


Extracting the required stuff is beyond me. We had many discussions already in the forum and the verdict was “mixed effects models for BE are not possible in R”.
So if you are happy with a fixed effect model (many jurisdictions like the EMA), fine. If you have to deal with the FDA or Health Canada, maybe not.


  1. Jones B, Kenward MG. Design and Analysis of Cross-Over Trials. Boca Raton: CRC Press; 3rd ed. 2015.
  2. Limited precision of results in the log-domain of the reference, CV not given. Estimated:
    library(PowerTOST)
    round(100*CVfromCI(lower=1.090, upper=1.239, design="3x6x3",
                       n=c(8, 11, 11, 10, 10, 10)), 1)

  3. SAS tells me: This method changes output in the following tables: Contrast, CorrB, CovB, Diffs, Estimates, InvCovB, LSMeans, Slices, SolutionF, SolutionR, Tests1–Tests3. The OUTP= and OUTPM= data sets are also affected.
  4. FDA/CDER. Guidance for Industry. Statistical Approaches to Establishing Bioequivalence. January 2001. online.

Cheers,
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. ☼
Science Quotes
mittyri
★★  

Russia,
2019-03-23 23:13

@ Helmut
Posting: # 20073
Views: 2,387
 

 Pooled vs IBD T-R in SAS from non-SASian

Dear Helmut,

» A request to the SASians of the forum: Please run Proc MIXED with DDFM=SATTERTHWAITE of the pooled data and both DDFM-options on the data where S is excluded.

Here's my (naive) attempt

<...>
ODS output estimates = _pooled;

proc mixed data=ex631 method=reml ITDETAILS maxiter=200;
class sequence subject period formula carry;
model ln_auc= sequence period formula/ddfm=SATTERTHWAITE outp=out;
random subject(sequence);
estimate 't T-R' formula -1 0 1/cl alpha=0.10;;
run;

Data _pooled;
set _pooled;
pe=round(100*exp(estimate), 0.00001);
lowerCL = lower;
upperCL = upper;
lower=round(100*exp(lowerCL), 0.00001);
upper=round(100*exp(upperCL), 0.00001);
run;

Data noS;
set ex631;
if formula = 'R' | formula = 'T';
run;

ODS output estimates = _TRonly;

proc mixed data=noS method=reml ITDETAILS maxiter=200;
class sequence subject period formula carry;
model ln_auc= sequence period formula/ddfm=SATTERTHWAITE outp=out;
random subject(sequence);
estimate 'T-R' formula -1 1/cl alpha=0.10;;
run;

Data _TRonly;
set _TRonly;
pe=round(100*exp(estimate), 0.00001);
lowerCL = lower;
upperCL = upper;
lower=round(100*exp(lowerCL), 0.00001);
upper=round(100*exp(upperCL), 0.00001);
run;


Results:
Option Label  Estimate StdErr DF tValue Probt Alpha Lower Upper    pe        lowerCL    upperCL
Pooled "t T-R" 0.1497 0.03848 115 3.89 0.0002 0.1  108.97 123.81 116.153  0.0859308889 0.2135523619
IBD     T-R    0.1488 0.03793 56.8 3.92 0.0002 0.1 108.92 123.65 116.047  0.0854035502 0.2122465728

Kind regards,
Mittyri
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2019-03-23 23:36

@ mittyri
Posting: # 20074
Views: 2,376
 

 Pooled vs IBD T-R in SAS from non-SASian

Hi mittyri,

» » A request to the SASians of the forum: Please run Proc MIXED with DDFM=SATTERTHWAITE of the pooled data and both DDFM-options on the data where S is excluded.
» Here's my (naive) attempt

THX a lot. I added your result to the comparison above.

Cheers,
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. ☼
Science Quotes
mittyri
★★  

Russia,
2019-03-23 23:50

@ Helmut
Posting: # 20076
Views: 2,370
 

 mixed in R

Dear Helmut,

since you are very quick with answers even Saturday night, I'm not editing my previous post but commenting it here :-D

» Extracting the required stuff is beyond me. We had many discussions already in the forum and the verdict was “mixed effects models for BE are not possible in R”.

Hmmm...
If I remember correctly we were struggling with FDA model where sophisticated 'repeated' statement exists.
We successfully crosschecked EMA method B (simple mixed effects model with a Subject as random effect). So linear mixed effects models are possible until FDA-style is required :cool:.

Kind regards,
Mittyri
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2019-03-24 01:23

@ mittyri
Posting: # 20077
Views: 2,384
 

 mixed in R (EMA B ≠ FDA)

Hi mittyri,

» If I remember correctly we were struggling with FDA model where sophisticated 'repeated' statement exists.
» We successfully crosschecked EMA method B (simple mixed effects model with a Subject as random effect).

Are you reminding Detlew and me obout our TODO-list? ;-)

» So linear mixed effects models are possible until FDA-style is required :cool:.

Yes, but that’s the point. However, recycled our code:

library(lmerTest) # (requires lme4, Matrix)
dta           <- read.table("exam45.dat", header=TRUE, na.strings="99999",
                            colClasses=c(rep("factor", 4), rep("numeric", 2)))
names(dta)[4] <- "treatment"
TR.only       <- dta[dta$treatment != "S", ]
ci            <- data.frame(rep(NA, 2), rep(NA, 4))
res           <- data.frame(method=c(rep("pooled", 2), rep("IBD", 2)),
                            PE=NA, CL.lo=NA, CL.hi=NA, CV=NA,
                            DFM=rep(c("Satterthwaite", "Kenward-Roger"), 2),
                            DF=NA, stringsAsFactors=FALSE)
for (j in 1:4) {
  if (j == 1) { # pooled (all at once)
    muddle <- lmer(log(AUC) ~ sequence + period + treatment + (1|subject),
                              data=dta)
  }
  if (j == 3) { # IBD (S excluded)
    muddle <- lmer(log(AUC) ~ sequence + period + treatment + (1|subject),
                              data=TR.only)
  }
  sum.muddle  <- summary(muddle, ddf=res$DFM[j])
  log.pe      <- sum.muddle$coefficients["treatmentT", "Estimate"]
  ci[j, 1:2]  <- round(100*exp(log.pe + c(-1, +1) *
                           qt(1-0.05, sum.muddle$coef["treatmentT", "df"]) *
                           sum.muddle$coef["treatmentT", "Std. Error"]), 2)
  res$PE[j]   <- round(100*exp(log.pe), 2)
  res$CL.lo[j] <- ci[j, 1]; res$CL.hi[j] <- ci[j, 2]
  res$CV[j]   <- round(100*sqrt(exp(sum.muddle$devcomp$cmp[["sigmaREML"]]^2)-1), 2)
  res$DF[j]   <- signif(sum.muddle$coefficients["treatmentT", "df"], 5)
}
print(res, row.names=FALSE)


Gives:

Method     PE  CL.lo  CL.hi    CV           DFM      DF
pooled 116.15 108.97 123.81 21.20 Satterthwaite 115.040
pooled 116.15 108.97 123.81 21.20 Kenward-Roger 114.630
   IBD 116.05 108.92 123.65 20.84 Satterthwaite  56.823
   IBD 116.05 108.91 123.65 20.84 Kenward-Roger  56.468


Similar same.

Though the DFs are slightly different, the CIs look only identical due to rounding.

Cheers,
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. ☼
Science Quotes
Activity
 Thread view
Bioequivalence and Bioavailability Forum |  Admin contact
19,567 posts in 4,150 threads, 1,340 registered users;
online 7 (0 registered, 7 guests [including 5 identified bots]).
Forum time (Europe/Vienna): 22:34 CEST

Mediocrity knows nothing higher than itself,
but talent instantly recognizes genius.    Arthur Conan Doyle

The BIOEQUIVALENCE / BIOAVAILABILITY FORUM is hosted by
BEBAC Ing. Helmut Schütz
HTML5