Bioequivalence and Bioavailability Forum

Main page Policy/Terms of Use Abbreviations Latest Posts

 Log-in |  Register |  Search

Back to the forum  Query: 2017-06-27 10:52 CEST (UTC+2h)
 
lizhao
Junior

US,
2016-01-22 00:13

Posting: # 15841
Views: 6,435
 

 R code for analyzing classical 2X2 crossover designed bioequivalence data [R for BE/BA]

Dear All,

I have a dataset that comes from a classical 2X2 crossover designed bioequivalence study.
I wonder if somebody could share with me the R code that can be used to analyze this dataset???

Should I use lm() or lme() in R (considering subject effects as random?)

Thanks a lot!!!
ElMaestro
Hero

Denmark,
2016-01-22 00:59

@ lizhao
Posting: # 15842
Views: 5,847
 

 R code for analyzing classical 2X2 crossover designed bioequivalence data

Hi Lizhao,

it is not too difficult.
You can use lm in the form of something like:

Mod = lm(lnCmax ~Trt+Seq+Subj+Per)

An easy option for the direct treatment effects is to use
Mod = lm(lnCmax ~0+Trt+Seq+Subj+Per)
and then you have the treatment effects directly out of the fixed effects or model summaries for Mod.

To generate a type III anova you can use drop1(Mod, test="F")
To generate the confidence interval, you just need the critical value at your chosen alpha and df, the fixed treatment effects (estimates thereof), and the SE which is calculated from the model's MSE and from the sequence sample sizes (which determine df too). You can find some code via this link.

Good luck. If you get stuck, send an email.

I could be wrong, but…


Best regards,
ElMaestro

- since June 2017 having an affair with the bootstrap.
lizhao
Junior

US,
2016-01-22 01:26

@ ElMaestro
Posting: # 15843
Views: 5,796
 

 R code for analyzing classical 2X2 crossover designed bioequivalence data

Hi ElMaestro,

Thank you very much for your help. I do have a question about using lm() instead of lme() function in R. Based on my understanding, subject effect should be considered to be random effect:confused:? If I use lm(), all effects are treated as fixed effects. You suggest me to use lm, is it based on FDA or EMA guidance?


Another question would be related to plot mean curves of Test and Reference product. If I want to plot a mean plot (mean and error bars) of a crossover biequivalence study, should I use within-subject confidence intervals or between-subject confidence intervals? They are apparently quite different.

Thanks again!
ElMaestro
Hero

Denmark,
2016-01-22 12:13

@ lizhao
Posting: # 15844
Views: 5,746
 

 On random effects and bogus statements

Hi Lizhao,

» Thank you very much for your help. I do have a question about using lm() instead of lme() function in R. Based on my understanding, subject effect should be considered to be random effect:confused:? If I use lm(), all effects are treated as fixed effects. You suggest me to use lm, is it based on FDA or EMA guidance?

Actually, you will see that FDA own SAS code for analysis of 222BE-studies involves PROC GLM, and this will always fit all effects as fixed. Hence lm in R. This despite the fact that Chow&Liu mentions the subject effect as random. It turns out that when you run the analysis on completers (completers only!) then PROC GLM/lm or a mixed effect model will achieve exactly the same. The beauty of PROC GLM/lm is that you get a well-defined anova along with it if you want. A well-defined anova doesn't exist for mixed models (but anovas of sorts can still be constructed).

You will also see that in FDA's code for 222BE which is based on PROC GLM there is a random statement pertaining to subject. One should think this means we are now asking PROC GLM to add a random effect (and thus to analyse as a mixed model), perhaps. It doesnt, though. PROC GLM with the random statement for subject still treats subject as a fixed effect. The bogus statement achieves another goal and compares sequence against the subject MS, and for that there is no direct equivalent in R, and this is why will see that there is a difference in the ANOVAs generated from PROC GLM and R's drop1.
Both can be argued to be correct, but since mass hysteria and brainlessness have trumped common sense in recent years everyone and his cousin just seem to want to reproduce PROC GLM's output without any questions asked.

» Another question would be related to plot mean curves of Test and Reference product. If I want to plot a mean plot (mean and error bars) of a crossover biequivalence study, should I use within-subject confidence intervals or between-subject confidence intervals? They are apparently quite different.

I think you can plot those curves in any way you want? The BE decision is taken on basis of the CI's.
If you plot Conc versus time and average across subjects then you would often choose to use between-subject variances for error bars.

I could be wrong, but…


Best regards,
ElMaestro

- since June 2017 having an affair with the bootstrap.
d_labes
Hero

Berlin, Germany,
2016-01-27 09:57

@ ElMaestro
Posting: # 15862
Views: 5,429
 

 FDA 'own' SAS code for 2x2

Dear ElMaestro!

Do you have a reference document for your mentioned FDA 'own' SAS code for 2x2 crossover based on Proc GLM with the "bogus" random statement?

Count up - post 999

Regards,

Detlew
nobody
Senior

2016-01-27 10:09

@ d_labes
Posting: # 15863
Views: 5,401
 

 FDA 'own' SAS code for 2x2

"If you get stuck, send an email."

totally OT: Is the email functionality of the forum back? Last time I wanted to invite myself via email for a cup of coffee in Denmark while travelling :-D I tried of no avail...

Kindest regards, nobody
Helmut
Hero
Homepage
Vienna, Austria,
2016-01-27 13:39

@ nobody
Posting: # 15864
Views: 5,340
 

 OT: PM to ElMaestro

Hi nobody,

» totally OT: Is the email functionality of the forum back?

It is.

» Last time I wanted to invite myself via email for a cup of coffee in Denmark while travelling :-D I tried of no avail...

ElMaestro’s ISP and mine have a long history of blacklisting each other. Since 2008 we exchanged 1200+ Mehls and occasionally they bounced. No idea why.
@ElMaestro: Consider changing your address to …@ymail.com.

[image]All the best,
Helmut Schütz 
[image]

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

Denmark,
2016-01-27 13:47

@ nobody
Posting: # 15865
Views: 5,357
 

 On coffee and such...

Hi nobody,

» totally OT: Is the email functionality of the forum back? Last time I wanted to invite myself via email for a cup of coffee in Denmark while travelling :-D I tried of no avail...

My coffee-maker broke down so I simply had to ask Helmut to deactivate the forum's email function temporarily. :-D:-D:-D Sorry about it. I am now using a manual drip filter (like switching from Windows 10 back to 3.1) and Helmut is trying to figure out which button to click in order to restore the forum's email capability. Have faith. Miracles will happen.

And while we are at it, here is a link to one of the better poems I read recently.

Want a coffee appointment in Shitville? Then click here for a direct email option.

I could be wrong, but…


Best regards,
ElMaestro

- since June 2017 having an affair with the bootstrap.
Helmut
Hero
Homepage
Vienna, Austria,
2016-01-27 13:51

@ ElMaestro
Posting: # 15866
Views: 5,392
 

 Geometric means ± SD

Hi ElMaestro & Lizhao,

» » If I want to plot a mean plot (mean and error bars) of a crossover biequivalence study, should I use within-subject confidence intervals or between-subject confidence intervals?
»
» I think you can plot those curves in any way you want? The BE decision is taken on basis of the CI's.
» If you plot Conc versus time and average across subjects then you would often choose to use between-subject variances for error bars.

Between – yes. (m/V)2? Ahem…
I suggest to plot geometric means* ± backtransformed standard deviations. Concentrations likely follow a log-normal distribution. Arithmetic means ± SD can lead to bizarre plots (like this one published by the FDA, page 12 of the PDF).


  • Having subjects with BQLs leads to trouble. Can be handled in an SOP.
    • Calculate the geometric mean if ≥⅔ of subjects at the given time point have concentrations ≥LLOQ and none otherwise, or
    • only (!) for the mean plots replace BQLs with LLOQ/2.
    An alternative is the median and 25/75 percentiles.

[image]All the best,
Helmut Schütz 
[image]

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

2016-01-27 15:20

@ Helmut
Posting: # 15868
Views: 5,321
 

 Geometric means ± SD

» "...like this one published by the FDA, page 12 of the PDF"

This plot is REALLY worth its money, quote "Study PR08210 enrolled 252 subjects and 238 completed the study". Made my day!

Kindest regards, nobody
Helmut
Hero
Homepage
Vienna, Austria,
2016-01-27 16:11

@ nobody
Posting: # 15869
Views: 5,309
 

 Money makes the world go ’round

Hi nobody,

» This plot is REALLY worth its money, quote "Study PR08210 enrolled 252 subjects and 238 completed the study". Made my day!

Did you realize that it was a 4-period full replicate study?

[image]All the best,
Helmut Schütz 
[image]

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

2016-01-27 16:31

@ Helmut
Posting: # 15870
Views: 5,253
 

 Money makes the world go ’round

» Did you realize that it was a 4-period full replicate study?

That's enough for the rest of the week! :-D ...but wait: How to get 238 completers? In prison? ooops, sorry, now I saw it...

Kindest regards, nobody
kumarnaidu
Regular

Mumbai, India,
2016-11-04 14:05

@ nobody
Posting: # 16773
Views: 3,250
 

 Fixed or random effect

Hello All,
Recently we got one query from WHO for 2x2 crossover study "Please take into account that if a mixed effect ANOVA model with the main effects of treatment, period and sequence as fixed effect and subjects nested within sequence as random effect is going to be used, as defined in the protocol, in SAS this refers to PROC MIXED as defined in the FDA guidelines. However, the protocol states at the same time that the sequence effect will be tested using the subjects nested within sequence mean square from the ANOVA as the error term. All other main effects will be tested against the residual error (mean square error/MSE) from the ANOVA as the error term. This refers to PROC GLM of SAS. This inconsistency should be corrected."

As per some experts if we want to keep sub(seq) as a error term for testing seq then it should be random rather than fixed.
But as I understand all effects in the PROC GLM are fixed including sub(seq).
Can anybody help me out?

Thanks in advance.

Kumar Naidu
mittyri
Senior

Russia,
2016-11-04 23:00

@ kumarnaidu
Posting: # 16774
Views: 3,208
 

 PROC GLM fixes Subject

Dear Kumar,

I would cite one of our Guru:
You will also see that in FDA's code for 222BE which is based on PROC GLM there is a random statement pertaining to subject. One should think this means we are now asking PROC GLM to add a random effect (and thus to analyse as a mixed model), perhaps. It doesnt, though. PROC GLM with the random statement for subject still treats subject as a fixed effect. The bogus statement achieves another goal and compares sequence against the subject MS

*--- GLM evaluation, output of T-R;
ODS output estimates=_est;
title "&file - GLM evaluation";
ods exclude  ExpectedMeanSquares;
Proc GLM data=PKmetric;
  class trt per seq subj;
  model logvar=trt per seq subj(seq)/CLparm alpha=&al2;
  random subj(seq)/test;
  estimate 'T-R' trt -1 1;
run; quit;

*--- back-transformation;
title "&file - &Ominus2al.% confidence intervals, back-transformed";
data _est;
  set _est;
  pe=exp(estimate);
  lower=exp(lowerCL);
  upper=exp(upperCL);
run;

Seems to be odd but here Subject is still fixed as ElMaestro wrote. So the experts are right (and ElMaestro too:-D)

Kind regards,
Mittyri
kumarnaidu
Regular

Mumbai, India,
2016-11-05 05:22

@ mittyri
Posting: # 16775
Views: 3,187
 

 PROC GLM fixes Subject

Hi Mittyri

Thanks for your reply.

» Seems to be odd but here Subject is still fixed as ElMaestro wrote. So the experts are right (and ElMaestro too:-D)

Ok that means sub(seq) is a random effect in proc glm. But what about above query from WHO. Also in EMA que and answer guidance they said all effects should be fixed rather than random (method A). Also as you said there is no meaning of random statement in proc glm. EMA Method also dont have such statement. My confusion is that in the protocol can we write that sub, period, sequence, treatment and sub(seq) will be a fixed effects and seq will be tested against sub(seq) if I am going to use Proc glm.

Thanks

Kumar Naidu
mittyri
Senior

Russia,
2016-11-05 10:13

@ kumarnaidu
Posting: # 16776
Views: 3,156
 

 PROC GLM fixes Subject

Hi Kumar,

» Ok that means sub(seq) is a random effect in proc glm.

Nope! Subject is still fixed, because you are using PROC GLM.
Please see the guide:
Note:PROC GLM uses only the information pertaining to expected mean squares when you specify the TEST option in the RANDOM statement and, even then, only in the extra tests produced by the RANDOM statement. Other features in the GLM procedure—including the results of the LSMEANS and ESTIMATE statements—assume that all effects are fixed, so that all tests and estimability checks for these statements are based on a fixed-effects model, even when you use a RANDOM statement.

» But what about above query from WHO.

They just want to see unambiguity in your statement.

» Also in EMA que and answer guidance they said all effects should be fixed rather than random (method A). Also as you said there is no meaning of random statement in proc glm. EMA Method also dont have such statement.

Are you working with Replicate study? The things are different there

» My confusion is that in the protocol can we write that sub, period, sequence, treatment and sub(seq) will be a fixed effects and seq will be tested against sub(seq) if I am going to use Proc glm.

I think you can provide the SAS code above, so all experts will understand that all effects are still fixed, and the reason of random part there is just an F test.

Kind regards,
Mittyri
kumarnaidu
Regular

Mumbai, India,
2016-11-07 06:24

@ mittyri
Posting: # 16783
Views: 3,084
 

 PROC GLM fixes Subject

Dear Mittyri and all,

» Are you working with Replicate study? The things are different there

Its a 2x2 crossover study.

» I think you can provide the SAS code above, so all experts will understand that all effects are still fixed, and the reason of random part there is just an F test.

proc glm data=replicate;
class formulation subject period sequence;
model logDATA= sequence subject (sequence) period formulation;
estimate "test-ref" formulation -1+1;
test h=sequence e=subject(sequence);
lsmeans formulation / adjust=t pdiff=control("R") CL alpha=0.10;
run;

This is Method A from EMA Q&A guidance. Ideally in GLM all effects are fixed.
Can we call here sub(seq) as a fixed effect?

Kumar Naidu
mittyri
Senior

Russia,
2016-11-07 11:57

@ kumarnaidu
Posting: # 16784
Views: 3,031
 

 PROC GLM fixes Subject

Dear Kumar,

» This is Method A from EMA Q&A guidance. Ideally in GLM all effects are fixed.
» Can we call here sub(seq) as a fixed effect?

yes, you can.
I kindly suggest to review the statistical part of the protocol with statistician experienced in BEQ studies evaluation. It would be useful for you, for this study and for the product.
EOD from my side.

Kind regards,
Mittyri
kumarnaidu
Regular

Mumbai, India,
2016-11-07 12:42

@ mittyri
Posting: # 16785
Views: 3,051
 

 PROC GLM fixes Subject

Thanks Mittyri,
» yes, you can.
» I kindly suggest to review the statistical part of the protocol with statistician experienced in BEQ studies evaluation. It would be useful for you, for this study and for the product.
The confusion arises due to the difference in the opinion :-D. Some CRO stat experts says all effects in the above program are fixed and some says that sub(seq) cannot be a fixed because we are test seq against this. Also the EMA guidance says that all effects are fixed in the Method A (PROC GLM).

Kumar Naidu
Ben
Regular

2016-11-08 21:20

@ Helmut
Posting: # 16789
Views: 2,970
 

 Geometric means ± SD

Dear Helmut,

» I suggest to plot geometric means* ± backtransformed standard deviations. Concentrations likely follow a log-normal distribution.

We should bear in mind that when we back transform to the original scale we end up with asymmetry with respect to +/-. Therefore, it is IMHO not appropriate to do +/- back transformed standard deviation. Symmetry is however obtained with respect to * and /. So the error bars should be such that lower end = gMean / gSD and upper end = gMean * gSD, where gSD = exp(standard deviation from log-transformed values).
[Mathematical rationale: exp(m - s) = exp(m)/exp(s) = gMean / gSD and exp(m + s) = exp(m) * exp(s) = gMean * gSD), where m and s are mean and sd on log scale, respectively]

Best regards,
Ben.
Helmut
Hero
Homepage
Vienna, Austria,
2016-11-09 15:29

@ Ben
Posting: # 16793
Views: 2,901
 

 Misunderstanding

Dear Ben,

» We should bear in mind that when we back transform to the original scale we end up with asymmetry with respect to +/-.

This is exactly the idea – reflecting the distributional properties of the data:

[image]



length of upper / lower whiskers:
0.826 0.950 0.859 0.730 0.528 0.411 0.329 0.273 0.228 0.179 0.111  0.0566
0.607 0.735 0.693 0.607 0.454 0.358 0.286 0.232 0.185 0.137 0.0758 0.0337


» Therefore, it is IMHO not appropriate to do +/- back transformed standard deviation. Symmetry is however obtained with respect to * and /. So the error bars should be such that lower end = gMean / gSD and upper end = gMean * gSD, where gSD = exp(standard deviation from log-transformed values).
» [Mathematical rationale: exp(m - s) = exp(m)/exp(s) = gMean / gSD and exp(m + s) = exp(m) * exp(s) = gMean * gSD), where m and s are mean and sd on log scale, respectively]

Can’t agree more. I did not express myself clearly. Blue your method and red the one* I meant. Of course they are equivalent.

x        <- c(1, 2, 4)
log.x    <- log(x)
mean.log <- mean(log.x)
SD.log   <- sd(log.x)
gMean    <- exp(mean.log)
gSD      <- exp(SD.log)
BL.lohi  <- c(gMean / gSD, gMean * gSD)
HS.lohi  <- exp(mean.log + c(-1, +1) * SD.log)
unique(BL.lohi == HS.lohi)
# [1] TRUE



  • Also used in Phoenix/WinNonlin.
    ──────────────
    Recipe for novices setting up the plot. I assume a table with at least the following columns: subject, scheduled, time, concentration. It is important to have a column with scheduled times. You must not use the actual times.
    Descriptive Stats:
      Mappings
        concentration ⦿ Summary
        scheduled ⦿ Sort
      Options
         Number of SD     1
    Execute. Results > Output Data > Results [right click] Send To > Plotting > XY Plot
    XY Plot:
    Mappings
        scheduled ⦿ X
        GeometricMean ⦿ Y
        GEOLowerSD Error Bars ⦿ Lower
        GEOUpperSD Error Bars ⦿ Upper
    Options
        Graphs > GeometricMean vs scheduled > Error Bars > Content
        Change the User Calculation Type Relative to Absolute

[image]All the best,
Helmut Schütz 
[image]

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

2016-11-14 19:40

@ Helmut
Posting: # 16802
Views: 2,705
 

 Misunderstanding

I see, then everything is fine :cool: Thanks for the clarification!

» Also used in Phoenix/WinNonlin.
» ──────────────
» Recipe ...

Ha! Great, thanks. The PK guys always told me that it is not possible to create such error bars in WinNonlin. From now on I can refer to this recipe :-D

Best regards,
Ben.
Helmut
Hero
Homepage
Vienna, Austria,
2016-11-14 21:50

@ Ben
Posting: # 16803
Views: 2,708
 

 RTFM

Hi Ben,

» The PK guys always told me that it is not possible to create such error bars in WinNonlin.

Possible for years! These guys should either RTFM or consult Certara/Pharsight’s forum. Last time I posted the recipe was two years ago.

[image]All the best,
Helmut Schütz 
[image]

The quality of responses received is directly proportional to the quality of the question asked. ☼
Science Quotes
Back to the forum Activity
 Thread view
Bioequivalence and Bioavailability Forum | Admin contact
17,007 Posts in 3,644 Threads, 1,043 registered users;
15 users online (0 registered, 15 guests).

We absolutely must leave room for doubt
or there is no progress and no learning.
There is no learning without having to pose a question.
And a question requires doubt.
People search for certainty.
But there is no certainty.    Richard Feynman

The BIOEQUIVALENCE / BIOAVAILABILITY FORUM is hosted by
BEBAC Ing. Helmut Schütz
XHTML/CSS RSS Feed