Helmut
★★★
avatar
Homepage
Vienna, Austria,
2012-10-26 19:35
(3376 d 11:02 ago)

Posting: # 9462
Views: 10,683
 

 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
[image]

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

Denmark,
2012-10-26 20:14
(3376 d 10:22 ago)

@ Helmut
Posting: # 9463
Views: 9,367
 

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

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 XtX; I would take a look at sigma that way and see what the truth :-D might look like.

Pass or fail!
ElMaestro
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2012-10-26 20:25
(3376 d 10:12 ago)

@ ElMaestro
Posting: # 9464
Views: 9,378
 

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

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 XtX; 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
[image]

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

Denmark,
2012-10-26 22:05
(3376 d 08:32 ago)

(edited by ElMaestro on 2012-10-26 22:51)
@ Helmut
Posting: # 9465
Views: 9,464
 

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

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:-D:-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 09:57
(3373 d 19:39 ago)

@ Helmut
Posting: # 9479
Views: 9,421
 

 Potvin – all effects fixed - SAS

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
★★★
avatar
Homepage
Vienna, Austria,
2012-11-06 01:43
(3366 d 03:54 ago)

@ d_labes
Posting: # 9505
Views: 9,400
 

 Potvin – all effects fixed (PHX/WNL resolved)

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
[image]

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

Denmark,
2012-11-06 11:48
(3365 d 17:48 ago)

@ Helmut
Posting: # 9507
Views: 9,228
 

 Potvin – all effects fixed (PHX/WNL resolved)

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
★★★
avatar
Homepage
Vienna, Austria,
2012-11-06 13:48
(3365 d 15:48 ago)

@ ElMaestro
Posting: # 9512
Views: 9,171
 

 Potvin – all effects fixed (PHX/WNL resolved)

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
[image]

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

Denmark,
2012-11-06 15:27
(3365 d 14:10 ago)

@ Helmut
Posting: # 9514
Views: 9,236
 

 Mysterious ways

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 :-D 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
★★★
avatar
Homepage
Vienna, Austria,
2012-11-06 17:06
(3365 d 12:31 ago)

@ ElMaestro
Posting: # 9515
Views: 9,175
 

 Mysterious ways

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
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
Activity
 Admin contact
21,844 posts in 4,571 threads, 1,555 registered users;
online 2 (0 registered, 2 guests [including 2 identified bots]).
Forum time: Monday 05:37 CET (Europe/Vienna)

Do, or do not.
There is no ‘try’.    Yoda

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