zan
☆    

US,
2014-01-30 00:58
(3711 d 10:33 ago)

Posting: # 12289
Views: 27,500
 

 How to calculate intersubject variability in PHX WinNonlin [Software]

All,
I encountered a negative Var(residual) for Cmax in PHX winnonlin when I run the 2x2 crossover test for intersubject and intrasubject CV% for Cmax and AUC.
I was able to obtain the intrasub CV for both parameters but the intersub CV for Cmax is missing. However I would still like to get the estimate for this intersub CV, I was wondering how can I calculate this value out in this case.

My limited searching tells me that this might be due to intersub CV% < intrasub CV%. What is the caution and steps to do when seeing this phenomena in the results?

Many thanks

zan


Edit: Category changed. [Helmut]
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2014-01-30 02:16
(3711 d 09:16 ago)

@ zan
Posting: # 12291
Views: 25,969
 

 Negative variance component

Hi Zan,

❝ I encountered a negative Var(residual) for Cmax in PHX winnonlin […]


❝ My limited searching tells me that this might be due to intersub CV% < intrasub CV%. What is the caution and steps to do when seeing this phenomena in the results?


Did you find this thread? A workaround for PHX at the end. If you want you can post your data set here (vari­ables separated by blanks, please). I have the latest beta-version of the new release of PHX on my machine.

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

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

US,
2014-01-30 19:11
(3710 d 16:20 ago)

@ Helmut
Posting: # 12296
Views: 25,809
 

 Negative variance component

Thank you, Helmut!
The thread really helps (although no clear answer as to how this happens). I will play with the workaround in PHX and calculate the CVinter from this.
Best regards


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! [Helmut]
zan
☆    

US,
2014-01-31 01:16
(3710 d 10:16 ago)

@ Helmut
Posting: # 12299
Views: 25,836
 

 Negative variance component

Hi Helmut,
I played with the workaround method in PHX but came up with invalid CVinter. Based on the suggested formula for PHX:
if(Hypothesis = 'Sequence*Subject', 100*sqrt(exp((MS-MSerror)/2)-1), ''), in my data MS is smaller than MSerror hence the resultant number before taking last step squareroot becomes negative.
Is it true that in the case the intersubject CV is considered nonestimable? Or is there other way to calculate this out?
Many thanks,
zan

❝ Did you find this thread? A workaround for PHX at the end. If you want you can post your data set here (vari­ables separated by blanks, please). I have the latest beta-version of the new release of PHX on my machine.

ElMaestro
★★★

Denmark,
2014-01-31 09:20
(3710 d 02:12 ago)

@ zan
Posting: # 12300
Views: 25,754
 

 Negative variance component

Hi,

❝ if(Hypothesis = 'Sequence*Subject', 100*sqrt(exp((MS-MSerror)/2)-1), ''), in my data MS is smaller than MSerror hence the resultant number before taking last step squareroot becomes negative.

❝ Is it true that in the case the intersubject CV is considered nonestimable? Or is there other way to calculate this out?


I think PHX is in mixed mode here. I'd switch to an all-fixed model. The results is the same as long as we talk standard 2,2,2-BE. If you really want to apply a mixed model for the result, then you can probably first do the all-fixed trick to get estimates of effects and between + within variances, then use the results coming from the all-fixed analysis as starter guess for the optimizer in the mixed model. That also puts less stress on the optimiser so your risk of getting trouble with the optimizer worker's union is low :-D
I am not a PHX/winnonlin user so I do not know the syntax for this approach.

Finally, without knowing anything about PHX/winnonlin, how about trying Subject everywhere rather than Sequence*Subject or Subject*Sequence?

Pass or fail!
ElMaestro
yjlee168
★★★
avatar
Homepage
Kaohsiung, Taiwan,
2014-01-31 11:26
(3710 d 00:06 ago)

@ zan
Posting: # 12302
Views: 25,907
 

 Negative variance component

Dear all,

Sorry to cut in. Yes, non-estimable CVinter could happen, I guess. In bear, we had a dataset for demo run purpose (single-dose, 2x2x2 design; and the dataset can be obtained from the top menu of bear ('Generate/export all demo datatsets'; it's called 'Single2x2x2_demo.csv' or 'Single2x2x2_demo.RData'.) ln(Cmax) is OK, but for ln(AUC0-t) and ln(AUC0-inf), we got

for ln(AUC0-t)...
Intra_subj. CV = 100*sqrt(exp(MSResidual)-1) = 18.392 %
Inter_subj. CV = 100*sqrt(exp((MSSubject(seq)-MSResidual)/2)-1)
               =  NaN %
    MSResidual = 0.0332664
MSSubject(seq) = 0.02123317 ...


and

for ln(AUC0-inf)...
Intra_subj. CV = 100*sqrt(exp(MSResidual)-1) = 17.788 %
Inter_subj. CV = 100*sqrt(exp((MSSubject(seq)-MSResidual)/2)-1)
               =  NaN %
    MSResidual = 0.03114947
MSSubject(seq) = 0.01895625  ...


Happy Chinese Lunar New Year (today)!

❝ ...for PHX:

❝ if(Hypothesis = 'Sequence*Subject', 100*sqrt(exp((MS-MSerror)/2)-1), ''), in my data MS is smaller than MSerror hence the resultant number before taking last step squareroot becomes negative.


All the best,
-- Yung-jin Lee
bear v2.9.1:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2014-02-01 17:03
(3708 d 18:29 ago)

@ yjlee168
Posting: # 12305
Views: 25,940
 

 Example data set

Hi Yung-jin,

wonderful! Below the data set for non-Rusers. No way in PHX (neither subjects random or fixed + my workaround).

Subject Period Formulation Sequence Cmax   AUCt    AUCinf
   1      2         R         TR    1633 11752.19 12432.84
   1      1         T         TR    1739 13712.93 14193.31
   2      1         R         RT    1481 11933.36 12602.04
   2      2         T         RT    1837 14681.77 15592.56
   3      2         R         TR    2073 14410.17 14903.48
   3      1         T         TR    1780 14516.58 15177.15
   4      1         R         RT    1374 10564.23 11169.47
   4      2         T         RT    1629 13308.92 13977.32
   5      2         R         TR    1385 11338.69 12037.22
   5      1         T         TR    1555 13278.27 13863.78
   6      1         R         RT    1756 14503.37 15091.20
   6      2         T         RT    1522 13093.42 13598.57
   7      2         R         TR    1643 11794.74 12412.79
   7      1         T         TR    1566 12764.19 13391.09
   8      1         R         RT    1939 12749.50 13328.77
   8      2         T         RT    1615 13509.55 13849.32
   9      2         R         TR    1759 14965.49 15726.48
   9      1         T         TR    1475 11775.35 12281.42
  10      1         R         RT    1388 12675.98 13350.59
  10      2         T         RT    1483 11187.10 12019.93
  11      2         R         TR    1682 14451.41 15132.62
  11      1         T         TR    1127  8821.43  9232.76
  12      1         R         RT    1542 14207.96 14950.48
  12      2         T         RT    1247 10084.26 10561.59
  13      2         R         TR    1605 14671.99 15548.11
  13      1         T         TR    1235  9368.49  9977.93
  14      1         R         RT    1598 14262.14 15200.68
  14      2         T         RT    2018 19367.23 19510.00



My post #3000. :-D

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

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
yjlee168
★★★
avatar
Homepage
Kaohsiung, Taiwan,
2014-02-01 18:40
(3708 d 16:52 ago)

@ Helmut
Posting: # 12307
Views: 25,678
 

 Example data set

Dear Helmut,

Thanks a lot for your help. Could you present the results (part of it will be good enough) from running PHX (or its beta)? I don't have PHX.

❝ ... No way in PHX (neither subjects random or fixed + my workaround).

❝ ... My post #3000. :-D


Amazing and thanks for sharing.

All the best,
-- Yung-jin Lee
bear v2.9.1:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2014-02-02 03:04
(3708 d 08:28 ago)

@ yjlee168
Posting: # 12310
Views: 25,870
 

 PHX build 6.3.0.395 / 6.4.0.511

Hi Yung-jin,

❝ […] Could you present the results (part of it will be good enough) from running PHX (or its beta)?


There are small differences to bear, since I started from the raw data and performed NCA first. This is a complete, balanced data set; therefore, the same results for subject(sequence) random or fixed are expected.

Table Partial SS

Dependent  Hypothesis       DF        SS        MS F_stat P_value
Ln(Cmax)   Sequence          1 0.0006903 0.0006903 0.02920 0.8672
Ln(Cmax)   Sequence*Subject 12 0.2836764 0.0236397 1.32782 0.3155
Ln(Cmax)   Formulation       1 0.0181691 0.0181691 1.02054 0.3323
Ln(Cmax)   Period            1 0.0362379 0.0362379 2.03545 0.1792
Ln(Cmax)   Error            12 0.2136406 0.0178034
Ln(AUCt)   Sequence          1 0.0151077 0.0151077 0.71554 0.4142
Ln(AUCt)   Sequence*Subject 12 0.2533637 0.0211136 0.63428 0.7791
Ln(AUCt)   Formulation       1 0.0108554 0.0108554 0.32611 0.5785
Ln(AUCt)   Period            1 0.0368259 0.0368259 1.10629 0.3136
Ln(AUCt)   Error            12 0.3994534 0.0332878
Ln(AUCinf) Sequence          1 0.0143291 0.0143291 0.76043 0.4003
Ln(AUCinf) Sequence*Subject 12 0.2261202 0.0188434 0.60445 0.8022
Ln(AUCinf) Formulation       1 0.0150303 0.0150303 0.48214 0.5007
Ln(AUCinf) Period            1 0.0346915 0.0346915 1.11282 0.3122
Ln(AUCinf) Error            12 0.3740929 0.0311744

This table is missing for the fixed effects model in the current release 6.3 (build 6.3.0.395) – therefore my work­around, but available in pre-release 1.4 (build 6.4.0.511).

Of course in the mixed model for AUCt and AUCinf:

Warning 11094: Negative final variance component. Consider omitting this VC structure.


Table Final Variance parameters (mixed)

Dependent  Parameter              Estimate
Ln(Cmax)   Var(Sequence*Subject)  0.0029182
Ln(Cmax)   Var(Residual)          0.0178034
Ln(Cmax)   Intersubject CV        0.0540594
Ln(Cmax)   Intrasubject CV        0.1340254
Ln(AUCt)   Var(Sequence*Subject) -0.0060871
Ln(AUCt)   Var(Residual)          0.0332878
Ln(AUCt)   Intrasubject CV        0.1839783
Ln(AUCinf) Var(Sequence*Subject) -0.0061655
Ln(AUCinf) Var(Residual)          0.0311744
Ln(AUCinf) Intrasubject CV        0.1779478


Table Final Variance parameters (fixed)

Dependent  Parameter              Estimate
Ln(Cmax)   Var(Residual)          0.0178034
Ln(AUCt)   Var(Residual)          0.0332878
Ln(AUCinf) Var(Residual)          0.0311744

Note the CVintra is not reported (why not?)…

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

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
yjlee168
★★★
avatar
Homepage
Kaohsiung, Taiwan,
2014-02-02 08:54
(3708 d 02:38 ago)

@ Helmut
Posting: # 12313
Views: 25,797
 

 PHX build 6.3.0.395 / 6.4.0.511

Dear Helmut,

Great post and thanks for sharing again. This allows me have a chance to compare the differences between bear and PHX (and its pre-release also).

❝ [...] There are small differences to bear, since I started from the raw data and performed NCA first. [...]


All the best,
-- Yung-jin Lee
bear v2.9.1:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here
ElMaestro
★★★

Denmark,
2014-02-01 17:31
(3708 d 18:00 ago)

@ zan
Posting: # 12306
Views: 25,726
 

 Negative variance component

Hi all,

how did

❝ if(Hypothesis = 'Sequence*Subject', 100*sqrt(exp((MS-MSerror)/2)-1), '')


or

"MSSubject(seq)-MSResidual" enter the game?

Doesn't it appear funny, given SStotal=SSwithin+SSbetween ?

With MSE for Cmax being much lower than MSE for AUC's and because the MS for subject is lower than the EMS I'd be inclined to think the Bear dataset itself is also a bit ... well ... deserving an audit :-D

Pass or fail!
ElMaestro
yjlee168
★★★
avatar
Homepage
Kaohsiung, Taiwan,
2014-02-01 18:47
(3708 d 16:45 ago)

@ ElMaestro
Posting: # 12308
Views: 25,904
 

 Negative variance component

Dear Elmaestro,

❝ ...

❝ "MSSubject(seq)-MSResidual" enter the game?


Should be this (as presented values in my previous post), I guess.

❝ Doesn't it appear funny, given SStotal=SSwithin+SSbetween ?


It's true. But I don't have any explanation for this.

❝ ... I'd be inclined to think the Bear dataset itself is also a bit ... well ... deserving an audit :-D


Though it's demo dataset for bear, it is real data.

All the best,
-- Yung-jin Lee
bear v2.9.1:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here
ElMaestro
★★★

Denmark,
2014-02-01 20:02
(3708 d 15:29 ago)

(edited by ElMaestro on 2014-02-02 01:27)
@ yjlee168
Posting: # 12309
Views: 25,721
 

 Just thinking loud

Thanks Yung-jin,

❝ ❝ "MSSubject(seq)-MSResidual" enter the game?


❝ Should be this (as presented values in my previous post), I guess.



Here I am thinking loud, and usually nothing good results when I do it. But here goes:


By way of a typical anova (forget for a moment the case of imbalance and type III SS which is irrelevant anyway for the actual dataset example) we have:

SStot=SStrt+SSseq+SSper+SSsubj+SSres
where
MSsubj=SSsubj/(n-2) [n=total number of subjects counting both sequences]
and

MSres= SSres/(n-2).

I want to know what the color and taste of "MSsubj-MSres" as prescribed above is:
MSsubj=SSsubj/(n-2)
or
MSsubj=(SStot-SStrt-SSseq-SSper-SSres)/(n-2)


Then
"MSsubj-MSres" = (SStot-SStrt-SSseq-SSper-SSres)/(n-2) - SSres/(n-2)
= (SStot-SStrt-SSseq-SSper-2*SSres)/(n-2)
This is truly a mindf%cker. My brain refuses to cooperate. I am experiencing the cerebral equivalent of a
Proton M rocket launch.
Actually, I am inclined to think there is something wrong with that subtraction in the first place.

Pass or fail!
ElMaestro
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2014-02-02 03:31
(3708 d 08:00 ago)

@ ElMaestro
Posting: # 12312
Views: 25,654
 

 All models are wrong…

Hi ElMaestro,

❝ Actually, I am inclined to think there is something wrong with that subtraction in the first place.


I told you on the phone that I’m going to meet Martin. He said: “Maybe the model is wrong.” Fuck. In Chow/Liu Chapter 7.3.1 (p192 last edition) I read:

A negative estimate may indicate that model 3.1.1 is incorrect or sample size is too small. More details on negative estimates in analysis of variance components can be found in Hocking (1985).

My emphasis. 3.1.1 is the usual 2,2,2 (though with carry over).
So if this model is wrong, which one is correct? :confused:

I’m not sure whether I shall invest in a book I likely will not understand anyway.

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

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
yjlee168
★★★
avatar
Homepage
Kaohsiung, Taiwan,
2014-02-02 09:07
(3708 d 02:24 ago)

(edited by yjlee168 on 2014-02-02 16:03)
@ Helmut
Posting: # 12314
Views: 26,079
 

 another book for linear model

Dear Helmut and Elmaestro,

As I read from Faraway JJ's textbook (Extending the Linear Model with R: Generalized Linear, Mixed Effects and Nonparametric Regression Models. Julian J. Faraway, Chapman & Hall/CRC, 2006 by Taylor & Francis Group, LLC., ISBN 1-58488-424-X., pp. 170-4), it said (p.171) '...The estimates (ANOVA estimator from glm()) can take negative values... This is rather embarrassing since variances cannot be negative. Various fixes have been proposed, but these all take away from the original simplicity of the estimation method...' there are some discussions about this in the text, including using the linear mixed effect model (lme or lmer). Will lme/lmer be the answer for this if we are using a wrong model when doing a 2x2x2 BE study for this situation only?
---
[edited] The problem of negative variance components seems commonly seen with ANOVA. When googling the website, we can find a lots about this information from user forums of some other big stat apps, such as SAS, SPSS and MINITAB. Here is one of them from MINITAB.


❝ [...] I’m not sure whether I shall invest in a book I likely will not understand anyway.


All the best,
-- Yung-jin Lee
bear v2.9.1:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2014-02-02 03:19
(3708 d 08:12 ago)

@ ElMaestro
Posting: # 12311
Views: 25,789
 

 References

Hi ElMaestro,

❝ how did

❝ ❝ if(Hypothesis = 'Sequence*Subject', 100*sqrt(exp((MS-MSerror)/2)-1), '')

❝ or

❝ "MSSubject(seq)-MSResidual" enter the game?


To name some:
  1. Midha KK, Ormsby ED, Hubbard JW, McKay G, Hawes EM, Gavalas L, McGilveray IJ. Logarithmic Transformation in Bioequivalence: Application with Two Formulations of Perphenazine. J Pharm Sci. 1993;82(2):138–44.
  2. Hauschke D, Steinijans VW, Diletti E, Schall R, Luus HG, Elze M, Blume H. Presentation of the intrasubject coefficient of variation for sample size planning in bio­equi­valence studies. Int J Clin Pharmacol Ther. 1994;32(7):376–8.
  3. Hauschke D, Steinijans V, Pigeot I. Bioequivalence Studies in Drug Development. Chichester; John Wiley: 2007.

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

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

Denmark,
2014-02-02 10:56
(3708 d 00:36 ago)

@ Helmut
Posting: # 12315
Views: 25,521
 

 References

Thanks Helmut,

once again I am sort of baffled.
I have to admit that I actually thought MSsubj was a direct measure of intersubject variability. Plenty of room for improvement at my end, it seems.

RCA: I don't understand because my brain is walnut-sized.
CAPA: I shall stop trying to understand.

Can we possibly do an lm on Test only with Sequence and Period as fixed and take the MSE as some measure of total variance (between+within), then do the same for Ref only, then at the end pool (or average) the two variance estimates to obtain some kind of total variability, from which we then derive a between by subtracting the within stemming from the grand analysis with subj, seq, per, and trt?

Pass or fail!
ElMaestro
d_labes
★★★

Berlin, Germany,
2014-02-03 10:02
(3707 d 01:30 ago)

@ zan
Posting: # 12318
Views: 25,523
 

 Negative variance component – Chow/Liu

Dear "discutanten"!

Did you remember the Chow, Liu book, Chapter 7.3?
There is an in depth discussion of the possibility of σ2S becoming negative if estimated via ANOVA mean squares.

There is also a formula giving the probability for obtaining a negative estimate, which is (was?) available in bear, I think.

Ways out? Just to cite Chow and Liu (without the 'hat' above the σ2):
"To avoid negative estimates a typical approach is to consider the following estimator
σ2S=max(0,σ2S)
σ2e2e if MSinter≥MSintra
σ2e2  if MSinter<MSintra
where σ2 = (SSinter+SSintra)/(2*(n1+n2)

The above estimators are known as restricted maximum likelihood (REML) estimators". End of citation.

In SAS a Proc MIXED call with subject as random should do that, I think. Will check it for the example later.

Regards,

Detlew
ElMaestro
★★★

Denmark,
2014-02-03 11:22
(3707 d 00:10 ago)

@ d_labes
Posting: # 12319
Views: 25,796
 

 Negative variance component – Chow/Liu

Hi Detleffff,

❝ Ways out? Just to cite Chow and Liu (without the 'hat' above the σ2):

"To avoid negative estimates a typical approach is to consider the following estimator

σ2S=max(0,σ2S)

❝ σ2e2e if MSinter≥MSintra

❝ σ2e2  if MSinter<MSintra

❝ where σ2 = (SSinter+SSintra)/(2*(n1+n2)

The above estimators are known as restricted maximum likelihood (REML) estimators". End of citation.


Yes, indeed the holy scripture says so; I think much of the confusion hinges on eq.7.3.3 saying
E(Vinter)=Ve+2VS
I think I do not follow it.



When PROC MIXED or R's LME/LER fit a true MM with REML, the Al Gore Rhythm does not derive the sigmas via any of these equations but from iterartively maximising the likelihoood of the V matrix with the fixed effects. It may be that the result ends up being these same as one of the REML estimators above. But are these estimates really "maximum likelihood"-related (however restricted they are), or are they just cheap ways out of an annoying situation?? After all, the residual of an LM/GLM (traditional ANOVA) is a decomposition of the entire variability into avialable factors, so if we e.g. set Ve=V as 7.3.3 suggests on the odd occasion then I think we are solving one problem in a quick and dirty fashion and at the same time doing something that appears very dubious. I guess my confusion boils down to something like "annoying negative values aside, where's the likelihood basis behind the idea of fiddling with the model's residual which truly is some kind of maximum likelihod estimator" ?
(and if we manually tweak Ve would we then reflect that in our calc. of the CI for "likelihood" reasons?)

If or when you test with PROC MIXED, can you paste the entire co-variance matrix (not the Z or the G)? It will be the one with 14 columns if you use the dataset above; I'd expect a common sigma sq. on the diagonal and a single beween-sigma sq. elsewhere in each row.

Pass or fail!
ElMaestro
d_labes
★★★

Berlin, Germany,
2014-02-03 12:58
(3706 d 22:34 ago)

(edited by d_labes on 2014-02-03 13:11)
@ ElMaestro
Posting: # 12320
Views: 25,895
 

 Variance components – Proc mixed

Dear ElMaestro,

❝ If or when you test with PROC MIXED, can you paste the entire co-variance matrix (not the Z or the G)? It will be the one with 14 columns if you use the dataset above; I'd expect a common sigma sq. on the diagonal and a single beween-sigma sq. elsewhere in each row.


As requested: the V matrices using the following code
Proc Mixed data=BEBAC12305;
  class subject period formulation sequence;
  model &lnMetric = formulation period sequence;
  random subject(sequence);
run;


ln(AUCt)
diagonal matrix (28 rows/cols) with 0.02720 on the diagonal

ln(Cmax)
Block diagonal matrix (28 rows/cols) with blocks
 0.02072   0.002918
 0.002918  0.02072
on the diagonal


Covariance Parameter Estimates:
ln(AUCt)
   V(subject(sequence))= 0 !
   V(residual)         = 0.02720   -> CVintra= 16.6%

ln(Cmax)
   V(subject(sequence))= 0.002918  -> CVinter=  5.4%
   V(residual)         = 0.01780   -> CVintra= 13.4%


Variance parameters are by default in SAS restricted to ≥0 in the maximization of the likelihood.
Reasonable to me :-D.
I get negative variance for subject(sequence) in case of ln(AUCt) only if I use the non-standard NOBOUND option in the Proc MIXED call:
ln(AUCt)
   V(subject(sequence))= -0.00609
   V(residual)         =  0.03329  -> CVintra= 18.4%
Seems the same as PHX build 6.3.0.395 / 6.4.0.511 results. See Helmut's post above.


Edit:
Interesting! Here the 90% CI's for AUCt (PE=96.14%):
Proc GLM             85.02% ... 108.71%
Proc Mixed           86.03% ... 107.44%
Proc Mixed (nobound) 85.02% ... 108.71%


Our (at least my own :lookaround:) belief "Proc GLM and Proc MIXED give the same results for a complete, balanced data set" has to be modified!

Regards,

Detlew
ElMaestro
★★★

Denmark,
2014-02-03 13:58
(3706 d 21:34 ago)

@ d_labes
Posting: # 12321
Views: 25,629
 

 Variance components – Proc mixed

Haha, thanks Detlefffff,

❝ Variance parameters are by default in SAS restricted to ≥0 in the maximization of the likelihood.

❝ Reasonable to me :-D.

❝ I get negative variance for subject(sequence) in case of ln(AUCt) only if I use the non-standard NOBOUND option in the Proc MIXED call:

ln(AUCt)

   V(subject(sequence))= -0.00609

   V(residual)         =  0.03329

❝ Seems the same as PHX build 6.3.0.395 / 6.4.0.511 results. See Helmut's post above.


This is mindblowing. I can't say that I understand in any way, but it is clear that the unbounded (=unfiddled, native and pure) optimisation analysis results in the same residual as the all-fixed lm/anova.
The million-dollar question asked in the nastiest fashion:
Do you either

believe in negative variances between subjects

-or-

would you inflate the MSE and get wide confidence intervals?

:-D:-D:-D:-D
An agonising choice indeed.
Another wrong question: Why care about between-Vars in a 2,2,2-BE? Why not just do the anova, fetch the residual, calculate a CI on basis of it and punch any guy who asks about betweens hard in the face?

I wonder how the PK-workgroup at EMA would deal with this.

I can't say I understand any details of the stats but this thread opened my eyes to an issue that I had no idea existed. Thanks.
:pirate:

Pass or fail!
ElMaestro
d_labes
★★★

Berlin, Germany,
2014-02-03 14:16
(3706 d 21:16 ago)

@ ElMaestro
Posting: # 12322
Views: 25,495
 

 Variance components – Proc mixed 90% CIs

Dear ÖbersterGrößterMeister!

❝ Do you either


❝ believe in negative variances between subjects


-or-


❝ would you inflate the MSE and get wide confidence intervals?


See my edit above :cool:.

Regards,

Detlew
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2014-02-03 15:14
(3706 d 20:17 ago)

@ ElMaestro
Posting: # 12323
Views: 25,643
 

 Variance components – Proc mixed

Dear all,

❝ This is mindblowing. I can't say that I understand in any way, but it is clear that the unbounded (=unfiddled, native and pure) optimisation analysis results in the same residual as the all-fixed lm/anova.


I followed Yung-jin’s suggestions; [image] comes up with 1.78 mio hits. Amazing this one – especially the “ways out” mentioned on page 130.

❝ An agonising choice indeed.


Yep.

❝ Why care about between-Vars in a 2,2,2-BE?


The only possible reason I can imagine: The cross-over turned out to be not such a good idea and you want to plan the next study in a parallel design. Then you would need the total (pooled) CV for sample size estimation. The common formula

\(CV_p\% = 100\sqrt{e^{(MS_s+MS_e)/2}-1}\)

would “work”, but result in values smaller than CVintra. Forget it.
           MSe       MSs    CVintra CVinter CVpooled
Cmax    0.0178034 0.0236397  13.40%  5.41%   14.47%
AUCt    0.0332878 0.0211136  18.40%   NA     16.61%
AUCinf  0.0311744 0.0188434  17.79%   NA     15.91%


❝ Why not just do the anova, fetch the residual, calculate a CI on basis of it…


Sure.

❝ …and punch any guy who asks about betweens hard in the face?


That’s Zan’s business.

❝ I wonder how the PK-workgroup at EMA would deal with this.


Not a problem for them, I guess. Doesn’t appear in their mandatory all fixed effects model (PHX) or Proc GLM. For the rest of the world (if running Proc MIXED) I would go with the UNBOUND option (see Detlew’s post). But – hey! – that’s not the code given in FDA’s guidances. ;-)
Reading a lot of stuff it’s evident that the restriction to ≥0 results in biased estimates. With Yung-jin’s data set the resulting CI turned out to be liberal. Keeping patient’s risk in mind that’s not a good idea. Is the unrestricted method (accepting negative vari­ances) always conservative? I think so, but I’m lacking the intellectual horsepower to prove it.

Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
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,
2014-02-03 16:54
(3706 d 18:38 ago)

@ Helmut
Posting: # 12324
Views: 25,632
 

 FDA code for non-replicate crossover?

Dear Helmut,

❝ ... But – hey! – that’s not the code given in FDA’s guidances. ;-)


where does your opinion came from?
So far as I know there is no code given in the FDA's guidances for non-replicate cross-over studies.

Moreover on page 10 of the 2001 Statistical guidance they recommend Proc GLM:
"General linear model procedures available in PROC GLM in SAS or equivalent software are preferred, although linear mixed-effects model procedures can also be indicated for analysis of nonreplicated crossover studies.

For example, for a conventional two-treatment, two-period, two-sequence (2 x 2) randomized crossover design, the statistical model typically includes factors accounting for the following sources of variation: sequence, subjects nested in sequences, period, and treatment. The Estimate statement in SAS PROC GLM, or equivalent statement in other software, should be used to obtain estimates for the adjusted differences between treatment means and the standard error associated with these differences."

Regards,

Detlew
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2014-02-03 17:16
(3706 d 18:15 ago)

@ d_labes
Posting: # 12325
Views: 25,527
 

 Proc GLM rulez

Dear Detlew,

❝ ❝ ... But – hey! – that’s not the code given in FDA’s guidances. ;-)


❝ where does your opinion came from?

❝ So far as I know there is no code given in the FDA's guidances for non-replicate cross-over studies.


Correct. Fingers faster than brain. In the obsolete (July 1992) guidance “Statistical Pro­cedures for Bio­equi­va­lence Studies using a Standard Two-Treatment Crossover Design” FDA recommended that

An analysis of variance (ANOVA) should be performed […] using General Linear Models (GLM) procedures of SAS or an equivalent program.


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

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
yjlee168
★★★
avatar
Homepage
Kaohsiung, Taiwan,
2014-02-03 21:43
(3706 d 13:48 ago)

@ Helmut
Posting: # 12327
Views: 25,584
 

 Variance components – Proc mixed

Dear Helmut,

Yes, I noticed that too.

❝ [...] but result in values smaller than CVintra. Forget it.

           MSe       MSs    CVintra CVinter CVpooled

Cmax    0.0178034 0.0236397  13.40%  5.41%   14.47%

AUCt    0.0332878 0.0211136  18.40%   NA     16.61%

AUCinf  0.0311744 0.0188434  17.79%   NA     15.91%


❝ [...] Not a problem for them, I guess. Doesn’t appear in their mandatory all fixed effects model (PHX) or Proc GLM. For the rest of the world (if running Proc MIXED) I would go with the UNBOUND option (see Detlew's post).


Sorry but I think it should be NOBOUND option with Proc MIXED.

All the best,
-- Yung-jin Lee
bear v2.9.1:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2014-02-03 23:08
(3706 d 12:24 ago)

@ yjlee168
Posting: # 12328
Views: 25,884
 

 NOBOUND

Hi Yung-jin,

❝ Sorry but I think it should be NOBOUND option with Proc MIXED.


Absolutely right – I dabble in SAS. I just have the 6.3-Manual (a verbose 8,640 pages):

NOBOUND

requests the removal of boundary constraints on covariance parameters. For example, variance components have a default lower boundary con­straint of 0, and the NOBOUND option allows their estimates to be neg­a­tive.


Parameter Constraints

By default, some covariance parameters are assumed to satisfy certain bound­ary constraints during the Newton-Raphson algorithm. For example, variance components are constrained to be nonnegative […]. You can re­move these constraints […] with the NOBOUND option in the PROC MIXED statement, but this can lead to estimates that produce an infinite likeli­hood. […]
For some data sets the final estimate of a parameter might equal one of its boundary constraints. This is usually not a cause for concern, but it might lead you to consider a different model. For instance, a variance com­­po­nent estimate can equal zero; in this case, you might want to drop the cor­re­spond­ing random effect from the model. However, be aware that changing the model in this fashion can affect degrees-of-freedom cal­cu­la­tions.


Convergence Problems
  • Using the NOPROFILE and NOBOUND options in the PROC MIXED statement might help convergence, although they can produce unusual results.

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

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
yjlee168
★★★
avatar
Homepage
Kaohsiung, Taiwan,
2014-02-03 21:22
(3706 d 14:10 ago)

(edited by yjlee168 on 2014-02-03 21:55)
@ d_labes
Posting: # 12326
Views: 25,592
 

 lm() or lme() for 2x2x2 study design?

Dear Detlew and all,

Thanks for your messages. Great discussions.

❝ ...

❝ There is also a formula giving the probability for obtaining a negative estimate, which is (was?) available in bear, I think.


Do you mean VarCorr() in lme()? yes, but I didn't use it any more. The question is that in bear we use lm() first, and get the negative variance. Then we can switch to use lme(..., method="REML") to re-do it for that negative variance (i.e., when CVinter is NaN, if(is.nan(CVinter){do lme()}). Why do we need to know the probability for obtaining a negative estimate? BTW, should we consider to use lme() instead if lm() for all 2x2x2 BE study since there seems no regulatory consideration? That really can avoid negative variances, though it happens occasionally.

All the best,
-- Yung-jin Lee
bear v2.9.1:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here
ElMaestro
★★★

Denmark,
2014-02-03 23:11
(3706 d 12:21 ago)

@ yjlee168
Posting: # 12329
Views: 25,447
 

 lm() or lme() for 2x2x2 study design?

Hi Yung-jin,


EU: All fixed, no discussion.
US: All fixed, but lots of confusion because the SAS PROC GLM has the random statement, giving users the impression that subject is treated as a random effect in the model.
Other places: Same, as far as I know. Some places might accept lme, but would probably at the outset nevertheless expect a simple linear model.

So lm and not lme for 2,2,2-BE.

Pass or fail!
ElMaestro
yjlee168
★★★
avatar
Homepage
Kaohsiung, Taiwan,
2014-02-04 14:09
(3705 d 21:22 ago)

(edited by yjlee168 on 2014-02-04 17:06)
@ ElMaestro
Posting: # 12335
Views: 25,525
 

 lm() or lme() for 2x2x2 study design?

Dear ElMaestro and all,

Thanks for your messages and good to hear that. The reason I asked was that it would be related to how I would modify bear later. BTW, results obtained from R's lme() may differ from that from SAS's PROC MIXED. As I ran lme() with the dataset that I provided before, I got different MSE (0.0244=0.15620852), CVintra(15.7%) and V(subj(seq)) (1.11*10-06)2 for AUC0-tfrom that using SAS's PROC MIXED as posted by Detlew. I still check with my codes to see if I do anything wrong.

❝ [...]

❝ So lm and not lme for 2,2,2-BE.


All the best,
-- Yung-jin Lee
bear v2.9.1:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here
yjlee168
★★★
avatar
Homepage
Kaohsiung, Taiwan,
2014-02-05 20:12
(3704 d 15:20 ago)

(edited by yjlee168 on 2014-02-06 08:02)
@ yjlee168
Posting: # 12349
Views: 25,980
 

 bear for 2x2x2 study with negative variance components

Dear all,

Using lme() to a 2x2x2 BE study, we can code with R (taking AUC0-t as example of the dataset in this thread) something like
modlnAUC0t<-lme(log(AUC0t) ~ drug + seq + prd,
               random=~1|subj/seq,
               data=blablabla, method="REML")
cat("\n")
print(summary(modlnAUC0t))
cat("\n")
cat("Type I Tests of Fixed Effects\n")
print(anova(modlnAUC0t)[2:4,])
cat("\n")
cat("Type III Tests of Fixed Effects\n")
print(anova(modlnAUC0t, type="marginal")[2:4,])
cat("\n\n")
...

Then outputs will be

Linear mixed-effects model fit by REML
 Data: blablabla
       AIC      BIC   logLik
  2.163918 10.41029 5.918041

Random effects:
 Formula: ~1 | subj
         (Intercept)
StdDev: 1.106768e-06   <- (1.11*10-6)2 ≈ 0 = V(subject(seq))

 Formula: ~1 | seq %in% subj
         (Intercept)  Residual
StdDev: 1.100023e-06 0.1562085 <- (0.156)2 = 0.0244 = MSE

Fixed effects: log(AUC0t) ~ drug + seq + prd
                Value  Std.Error DF   t-value p-value
(Intercept)  9.516955 0.05904125 12 161.19163  0.0000
drug2        0.062297 0.05904125 12   1.05514  0.3121
seq2        -0.036450 0.05904125 12  -0.61736  0.5485
prd2        -0.048908 0.05904125 12  -0.82837  0.4236
...
Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max
-1.8446164 -0.7801051  0.1796820  0.6539527  1.6438227

Number of Observations: 28
Number of Groups:
         subj seq %in% subj
           14            14
            numDF denDF   F-value p-value
(Intercept)     1    12 103679.28  <.0001
drug            1    12      1.11  0.3121
seq             1    12      0.38  0.5485
prd             1    12      0.69  0.4236

Type I Tests of Fixed Effects
     numDF denDF   F-value p-value
drug     1    12 1.1133109  0.3121
seq      1    12 0.3811317  0.5485
prd      1    12 0.6862044  0.4236

Type III Tests of Fixed Effects
     numDF denDF   F-value p-value
drug     1    12 1.1133109  0.3121
seq      1    12 0.3811317  0.5485
prd      1    12 0.6862044  0.4236

Since we will still stick on lm() with 2x2x2 study, so bear will show the outputs as
...
Intra_subj. CV = 100*sqrt(exp(MSResidual)-1) = 18.038 %
Inter_subj. CV = 100*sqrt(exp((MSSubject(seq)-MSResidual)/2)-1)
               = 0 % (with a negative variance component)
*** the above CV_intra is estimated from lm() which may be different
    from than that obtained from lme().


    MSResidual = 0.03201731
MSSubject(seq) = 0.01678485

Therefore, the CVintra from lme() will not be calculated in this case. But luckily for us, Detlew has showed us that V(subject(seq)) is zero (or should be close to zero) and the 90%CI was the same as what we got from lm() using SAS PROC MIXED. Sorry for this lengthy post.

All the best,
-- Yung-jin Lee
bear v2.9.1:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here
UA Flag
Activity
 Admin contact
22,957 posts in 4,819 threads, 1,636 registered users;
112 visitors (1 registered, 111 guests [including 9 identified bots]).
Forum time: 11:32 CET (Europe/Vienna)

With four parameters I can fit an elephant,
and with five I can make him wiggle his trunk.    John von Neumann

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