Bioequivalence and Bioavailability Forum • Potvin – all effects fixed (PHX/WNL vs. SAS, R)

Helmut
★★★

Vienna, Austria,
2012-10-26 21:35
(4558 d 00:31 ago)

Posting: # 9462
Views: 12,933

Potvin – all effects fixed (PHX/WNL vs. SAS, R) [Two-Stage / GS Designs]

Dear all,

Jiří Hofmann brought one subtlety to my attention. I set up the pooled analysis with Sequence + Stage + Period(Stage) + Treatment [fixed] and Subject(Sequence × Stage) [random].

For Example 2 I got:

User-Specified Confidence Level for CI's and Power = 94.1200

Reference: R   LSMean= 1.133431  SE= 0.171385  GeoLSM= 3.106297

---------------------------------------------------------------

Test:      T   LSMean= 1.147870  SE= 0.171385  GeoLSM= 3.151473



    Ratio(%Ref) = 101.4544



     CI User = ( 88.4452,  116.3770)



    Average bioequivalence shown for confidence=94.12 and percent=20.0.

When we specified all effects fixed (EMA?), Phoenix/WinNonlin (tested 6.2.51, 6.3beta, 6.3) spits out:
ERROR 11050: Parsing error: The containing term must be in the model.

OK, since in PHX/WNL the BE-module sits on top of the LME-engine I set up the model directly in LME. Could reproduce the mixed model:

Least squares means

Treatment  Estimate  StdError Denom_DF T_stat P_value Conf T_crit Lower_CI Upper_CI

-----------------------------------------------------------------------------------

        R   1.13343  0.171385  18.4   6.61338  0.0000   95  2.098   0.7739    1.493

        T   1.14787  0.171385  18.4   6.69763  0.0000   95  2.098   0.7883    1.507



Estimate Statements

Title   Estimate  StdError  Denom_DF   T_stat P_value Conf T_crit Lower_CI Upper_CI

-----------------------------------------------------------------------------------

T – R  0.0144389  0.0677463    17    0.213132  0.8338 94.12 2.026  -0.1228   0.1517

… giving back-transformed 101.4544 (88.4452, 116.3770), identical to BE-wizard’s results.

I get the same values in the all-fixed model – but could not request LSMeans any more; otherwise I got the same error like in the BE-wizard. Though we don’t need these values in the BE assessment I would be interested how other software deals with it (SAS, R?).

—
Dif-tor heh smusma 🖖🏼 Довге життя Україна!
Helmut Schütz

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes

ElMaestro
★★★

Denmark,
2012-10-26 22:14
(4557 d 23:51 ago)

@ Helmut
Posting: # 9463
Views: 11,155

Potvin – all effects fixed (PHX/WNL vs. SAS, R)

Post reply

Hi Helmut,

❝ For Example 2 I got:

Where is example 2??

❝ When we specified all effects fixed (EMA?), Phoenix/WinNonlin (tested 6.2.51, 6.3beta, 6.3) spits out:

❝ ERROR 11050: Parsing error: The containing term must be in the model.

This should absolutely work, I think. Either it is the syntax that was wrong or PHX/WNL is flawed?

❝ OK, since in PHX/WNL the BE-module sits on top of the LME-engine I set up the model directly in LME. Could reproduce the mixed model:

(and blah)

❝ … giving back-transformed 101.4544 (88.4452, 116.3770), identical to BE-wizard’s results.

❝ I get the same values in the all-fixed model – but could not request LSMeans any more; otherwise I got the same error like in the BE-wizard. Though we don’t need these values in the BE assessment I would be interested how other software deals with it (SAS, R?).

You don't need a mixed model when you have only one sigma in play (this isn't a replicated study I presume?). In the simple non-replicated case fitting a mixed model with subject as random should give the same result as the normal linear model based on all fixed. If there are subtle differences at some decimal then this sounds most likely just like either a rounding phenomenon, or a convergence issue where somehow the optimiser in this case is capable of finding smaller sigma when the mixed model is applied.
The normal linear model can be solved exactly with e.g. R at least to the extent we can rely on its numerical stability for inversion of X^tX; I would take a look at sigma that way and see what the truth :-D

might look like.

—
Pass or fail!
ElMaestro

Helmut
★★★

Vienna, Austria,
2012-10-26 22:25
(4557 d 23:41 ago)

@ ElMaestro
Posting: # 9464
Views: 11,175

Potvin – all effects fixed (PHX/WNL vs. SAS, R)

Post reply

Dear ElMaestro!

❝ ❝ For Example 2 I got:

❝ Where is example 2??

At the end of Madame Potvin’s paper. For convenience:

sub  <- as.factor(sort(rep(c(1:20),2)))

stg  <- as.factor(c(rep(1,12),rep(2,8)))

seq  <- as.factor(c(rep('TR',6),c(rep('RT',6)),

                    rep('TR',4),c(rep('RT',4))))

per  <- as.factor(rep(c(1,2),20))

trt  <- as.factor(c(rep(c('T','R'),6),rep(c('R','T'),6),

                    rep(c('T','R'),4),rep(c('R','T'),4)))

data <- data.frame(stg, sub, trt, per, seq)

data$resp <- c(4.8,4.6,1.4,1.3,3.9,4.4,1.5,2.4,3.8,2.9,3.3,2.5,

               1.5,1.7,1.9,1.3,2.4,2.1,13.9,11.4,2.3,1.5,2.4,2.6,

               4.0,4.2,0.7,0.5,10.3,11.9,4.7,4.8,5.5,5.9,5.4,3.8,2.1,4.3,2.4,3.6)

❝ This should absolutely work, I think. Either it is the syntax that was wrong or PHX/WNL is flawed?

Well, the flaw sits behind the LSMeans.

❝ You don't need a mixed model when you have only one sigma in play (this isn't a replicated study I presume?).

I know the Silly-O-Meter! In PHX/WNL everything goes the REML-way. The only way to mimick the behaviour of SAS Proc GLM is to specify all effects fixed, and filter incomplete data beforehand. Here this is not an issue, of course.

❝ In the simple non-replicated case fitting a mixed model with subject as random should give the same result as the normal linear model based on all fixed.

Yes.

❝ If there are subtle differences at some decimal then this sounds most likely just like either a rounding phenomenon, or a convergence issue where somehow the optimiser in this case is capable of finding smaller sigma when the mixed model is applied.

Yep – though the setup (convergence criteria, step size, :blahblah:

) was the same. Doesn’t bother me too much (yet).

❝ The normal linear model can be solved exactly with e.g. R at least to the extent we can rely on its numerical stability for inversion of X^tX; I would take a look at sigma that way and see what the truth :-D might look like.

I’m try to dig out some code from the grave.

—
Dif-tor heh smusma 🖖🏼 Довге життя Україна!
Helmut Schütz

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes

ElMaestro
★★★

Denmark,
2012-10-27 00:05
(4557 d 22:01 ago)

@ Helmut
Posting: # 9465
Views: 11,268

Potvin – all effects fixed (PHX/WNL vs. SAS, R)

Post reply

Hi Helmut,

check your coding, I thik there's an issue with seq.

I get:


Subj=as.factor(c(1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10,

         11,11,12,12,13,13,14,14,15,15,16,16,17,17,18,18,19,19,20,20))

Seq=as.factor(c("TR","TR","TR","TR","TR","TR",   

                "TR","TR","TR","TR","TR","TR",  

                "RT","RT","RT","RT","RT","RT", 

                "RT","RT","RT","RT","RT","RT",

                "TR","TR","TR","TR",

                "TR","TR","TR","TR",

                "RT","RT","RT","RT",

                "RT","RT","RT","RT"))

Per=as.factor((rep(c(1,2), 20)))

Trt=as.factor(c(1,2,1,2,1,2,1,2,1,2,1,2, 

                2,1,2,1,2,1,2,1,2,1,2,1, 

                1,2,1,2,1,2,1,2,

                2,1,2,1,2,1,2,1))

Stg=as.factor(c(rep(1,24), rep(2,16)))



PK=c(4.8, 4.6, 1.4, 1.3, 3.9, 4.4, 1.5, 2.4, 3.8, 2.9, 3.3, 2.5, 1.7, 1.5, 1.3, 1.9, 2.1, 2.4, 11.4, 13.9, 1.5, 2.3, 2.6, 2.4,

4.0,4.2,0.7, 0.5, 10.3, 11.9,4.7, 4.8, 5.9,5.5, 3.8, 5.4, 4.3, 2.1, 3.6, 2.4)

options(digits=12)

M=lm(log(PK)~0+Trt+Seq+Subj%in%(Stg*Seq)+Per%in%Stg+Stg)

lnPE=coef(M)[1]-coef(M)[2]


df=anova(M)[6,1]


ct=qt(1-0.0294,df)

nps=length(Trt)/4

foo=sqrt(0.5*(anova(M)[6,3])* (1/nps+1/nps)   )  

U=lnPE+ct*foo

L=lnPE-ct*foo  

cat("Difference in treatment means is ", coef(M)[1]-coef(M)[2], "\n")

Difference in treatment means is  0.01443888498 

cat("PE: ", exp(lnPE), "\n")

PE:  1.01454362920 

> cat("The residual is ", anova(M)[6,3], "\n")

The residual is  0.0458955516532 

cat("Upper ", exp(U), "\n")

Upper  1.16377016091 

cat("Lower ", exp(L), "\n")

Lower  0.88445193916

We have a winner?

PS: With that coding for subject we will of course not talk about subject nested in sequence or stage :-D

. There is no subject 1 in sequence "RT" and no subject 1 in stage 2 and so forth. So we might just as well do:

M=lm(log(PK)~0+Trt+Stg+Seq+Subj+Per%in%Stg)

—
Pass or fail!
ElMaestro

d_labes
★★★

Berlin, Germany,
2012-10-29 10:57
(4555 d 10:08 ago)

@ Helmut
Posting: # 9479
Views: 11,211

Potvin – all effects fixed - SAS

Post reply

Dear Helmut!

Here the sasophylistic answer (to what question ever):

Used code:

title "Long form analysis after stage 2";

Proc GLM data=Potvin2_long;

  class treatment period subject sequence stage;



  model lnPK = treatment period(stage) sequence stage subject(stage*sequence)

               /ss3 CLParm alpha=0.0588; * =2*0.0294;

  * random serves only for correct tests for sequence stage;

  random stage sequence subject(sequence*stage) /test;

  estimate 'T-R' treatment -1 1;

  lsmeans treatment / stderr cl pdiff alpha=0.0588;

run; quit;

As we know meanwhile the random statement doesn't change the "all fixed effects" solution of Proc GLM.

Output (only interested parts):

 ...

                    Least Squares Means



                                                             H0:LSMean1=

                                  Standard    H0:LSMEAN=0      LSMean2

 treatment     lnPK LSMEAN           Error       Pr > |t|       Pr > |t|



 R              1.13343123      0.04840026         <.0001         0.8338

 T              1.14787011      0.04840026         <.0001





 treatment     lnPK LSMEAN     94.12% Confidence Limits



 R                1.133431        1.035392     1.231470

 T                1.147870        1.049831     1.245909





         Least Squares Means for Effect treatment



               Difference

                  Between     94.12% Confidence Limits

   i    j           Means     for LSMean(i)-LSMean(j)



   1    2       -0.014439       -0.151665     0.122787





         94.12% (1-2*alpha) CIs - backtransformed



                   point         lower         upper



              101.454363     88.445194    116.377016

Attention! Must had change the sign (swap upper, lower) in the log domain since R-T is calculated due to lexical order of treatments coded as R and T.

Remark for the nitpickers and silly-o-meter inventors:
Within SAS you must use the subject effect nested within sequence and stage, even if the subjects are coded unique. Otherwise you don't obtain a type III ANOVA and no LSMeans.

—
Regards,

Detlew

Helmut
★★★

Vienna, Austria,
2012-11-06 02:43
(4547 d 18:23 ago)

@ d_labes
Posting: # 9505
Views: 11,204

Potvin – all effects fixed (PHX/WNL resolved)

Post reply

Dear Detlew, ElMaestro, Jiří, and all,

Linda Hughes of Pharsight clarified the setup. Actually

ERROR 11050: Parsing error: The containing term must be in the model.

gives the initiated (not me) a hint of what’s going on. Phoenix/WinNonlin requires the inner nested term in order to use it. Strange. However, the all fixed effects coding is:
Sequence+Stage+Period(Stage)+Treatment+Sequence*Stage+Subject(Sequence*Stage)

Now I get the same results (PE, CI, Means, SEs, :blahblah:

) as in PHX’ mixed model, SAS, and R. Don’t run the silly-o-meter on this setup.

—
Dif-tor heh smusma 🖖🏼 Довге життя Україна!
Helmut Schütz

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes

ElMaestro
★★★

Denmark,
2012-11-06 12:48
(4547 d 08:17 ago)

@ Helmut
Posting: # 9507
Views: 11,039

Potvin – all effects fixed (PHX/WNL resolved)

Post reply

Dear Helmut,

❝ Sequence+Stage+Period(Stage)+Treatment+Sequence*Stage+Subject(Sequence*Stage)

❝

❝ Now I get the same results (PE, CI, Means, SEs, :blahblah: ) as in PHX’ mixed model, SAS, and R. Don’t run the silly-o-meter on this setup.

Looks like this would kick the the silly-o-meter directly into a cosmic overdrive. I would intuitively say you should get the same results if you remove the Sequence (first) term. Can check?

By the way: If models need to be written anything like that then I would guess period would need to be specified as period(sequence) in the standard 2,2,2-BE model. Just trying to follow the logic that I don't understand. Nicht?

—
Pass or fail!
ElMaestro

Helmut
★★★

Vienna, Austria,
2012-11-06 14:48
(4547 d 06:17 ago)

@ ElMaestro
Posting: # 9512
Views: 10,965

Potvin – all effects fixed (PHX/WNL resolved)

Post reply

Hi ElMaestro!

❝ ❝ Sequence+Stage+Period(Stage)+Treatment+Sequence*Stage+Subject(Sequence*Stage)

❝ ❝ […] Don’t run the silly-o-meter on this setup.

❝

❝ Looks like this would kick the the silly-o-meter directly into a cosmic overdrive.

Definitely.

❝ I would intuitively say you should get the same results if you remove the Sequence (first) term. Can check?

Correct.

❝ By the way: If models need to be written anything like that then I would guess period would need to be specified as period(sequence) in the standard 2,2,2-BE model. Just trying to follow the logic that I don't understand. Nicht?

Nö (God moves in mysterious ways). Sequence+Treatment+Period+Subject(Sequence) works, whereas Sequence+Treatment+Period(Sequence)+Subject(Sequence) does not. In the magic BE-wizard LSMeans “not estimable” – end of story.

—
Dif-tor heh smusma 🖖🏼 Довге життя Україна!
Helmut Schütz

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes

ElMaestro
★★★

Denmark,
2012-11-06 16:27
(4547 d 04:39 ago)

@ Helmut
Posting: # 9514
Views: 11,012

Mysterious ways

Post reply

Hi Helmut,

❝ Nö (God moves in mysterious ways). Sequence+Treatment+Period+Subject(Sequence) works, whereas Sequence+Treatment+Period(Sequence)+Subject(Sequence) does not. In the magic BE-wizard LSMeans “not estimable” – end of story.

Excellent

Now the confusion should be total: it doesn't work, but it isn't clear why the logic doesn't or shouldn't apply to this case. Next time I talk to someone from Pharsight I will ask if this is a bug or a feature, or if I can somehow progress to the next level of cosmic insight into the logic behind. Period%in%Sequence syntactically works in R (like Sequence%in%Period, haha), although it makes no difference. Thanks for the test.

—
Pass or fail!
ElMaestro

Helmut
★★★

Vienna, Austria,
2012-11-06 18:06
(4547 d 03:00 ago)

@ ElMaestro
Posting: # 9515
Views: 10,983

Mysterious ways

Post reply

Hi ElMaestro!

❝ […] it doesn't work, but it isn't clear why the logic doesn't or shouldn't apply to this case.

Yep.

❝ Next time I talk to someone from Pharsight I will ask if this is a bug or a feature, or if I can somehow progress to the next level of cosmic insight into the logic behind.

Software developers will always tell you that’s a feature, although:

If debugging is the process of removing bugs,
then programming must be the process of putting them in. Edsger W. Dijkstra

—
Dif-tor heh smusma 🖖🏼 Довге життя Україна!
Helmut Schütz

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes