Bioequivalence and Bioavailability Forum • How to calculate intersubject variability in PHX WinNonlin

zan
☆

US,
2014-01-30 00:58
(3967 d 20:49 ago)

Posting: # 12289
Views: 29,078

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
★★★

Vienna, Austria,
2014-01-30 02:16
(3967 d 19:31 ago)

@ zan
Posting: # 12291
Views: 27,346

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 (variables separated by blanks, please). I have the latest beta-version of the new release of PHX on my machine.

—
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

zan ☆ US, 2014-01-30 19:11 (3967 d 02:36 ago) @ Helmut Posting: # 12296 Views: 27,184	Negative variance component Post reply
	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
(3966 d 20:31 ago)

@ Helmut
Posting: # 12299
Views: 27,209

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 (variables 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
(3966 d 12:27 ago)

@ zan
Posting: # 12300
Views: 27,131

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

:-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
★★★

Kaohsiung, Taiwan,
2014-01-31 11:26
(3966 d 10:21 ago)

@ zan
Posting: # 12302
Views: 27,293

Negative variance component

Dear all,

Sorry to cut in. Yes, non-estimable CV_inter 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(AUC_0-t) and ln(AUC_0-inf), we got

for ln(AUC_0-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(AUC_0-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.2:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here

Helmut
★★★

Vienna, Austria,
2014-02-01 17:03
(3965 d 04:45 ago)

@ yjlee168
Posting: # 12305
Views: 27,337

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

:-D

—
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

yjlee168
★★★

Kaohsiung, Taiwan,
2014-02-01 18:40
(3965 d 03:08 ago)

@ Helmut
Posting: # 12307
Views: 27,078

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.2:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here

Helmut
★★★

Vienna, Austria,
2014-02-02 03:04
(3964 d 18:43 ago)

@ yjlee168
Posting: # 12310
Views: 27,276

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 workaround, 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 CV_intra is not reported (why not?)…

—
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

yjlee168
★★★

Kaohsiung, Taiwan,
2014-02-02 08:54
(3964 d 12:54 ago)

@ Helmut
Posting: # 12313
Views: 27,169

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.2:- 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
(3965 d 04:16 ago)

@ zan
Posting: # 12306
Views: 27,107

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 SS_total=SS_within+SS_between ?

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

:-D

—
Pass or fail!
ElMaestro

yjlee168
★★★

Kaohsiung, Taiwan,
2014-02-01 18:47
(3965 d 03:00 ago)

@ ElMaestro
Posting: # 12308
Views: 27,311

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 SS_total=SS_within+SS_between ?

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.2:- 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
(3965 d 01:45 ago)

(edited on 2014-02-02 01:27)
@ yjlee168
Posting: # 12309
Views: 27,100

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:

SS_tot=SS_trt+SS_seq+SS_per+SS_subj+SS_res
where
MS_subj=SS_subj/(n-2) [n=total number of subjects counting both sequences]
and

MS_res= SS_res/(n-2).

I want to know what the color and taste of "MS_subj-MS_res" as prescribed above is:
MS_subj=SS_subj/(n-2)
or
MS_subj=(SS_tot-SS_trt-SS_seq-SS_per-SS_res)/(n-2)

Then
"MS_subj-MS_res" = (SS_tot-SS_trt-SS_seq-SS_per-SS_res)/(n-2) - SS_res/(n-2)
= (SS_tot-SS_trt-SS_seq-SS_per-2*SS_res)/(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
★★★

Vienna, Austria,
2014-02-02 03:31
(3964 d 18:16 ago)

@ ElMaestro
Posting: # 12312
Views: 27,035

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:

:confused:

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

—
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

yjlee168
★★★

Kaohsiung, Taiwan,
2014-02-02 09:07
(3964 d 12:40 ago)

(edited on 2014-02-02 16:03)
@ Helmut
Posting: # 12314
Views: 27,498

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.2:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here

Helmut
★★★

Vienna, Austria,
2014-02-02 03:19
(3964 d 18:28 ago)

@ ElMaestro
Posting: # 12311
Views: 27,203

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:

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.
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 bioequivalence studies. Int J Clin Pharmacol Ther. 1994;32(7):376–8.
Hauschke D, Steinijans V, Pigeot I. Bioequivalence Studies in Drug Development. Chichester; John Wiley: 2007.

—
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,
2014-02-02 10:56
(3964 d 10:51 ago)

@ Helmut
Posting: # 12315
Views: 26,918

References

Thanks Helmut,

once again I am sort of baffled.
I have to admit that I actually thought MS_subj 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
(3963 d 11:46 ago)

@ zan
Posting: # 12318
Views: 26,897

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 σ²_S 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 σ²):
"To avoid negative estimates a typical approach is to consider the following estimator

σ²_S=max(0,σ²_S)

σ²_e=σ²_e if MS_inter≥MS_intra

σ²_e=σ²  if MS_inter<MS_intra

where σ² = (SS_inter+SS_intra)/(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
(3963 d 10:25 ago)

@ d_labes
Posting: # 12319
Views: 27,235

Negative variance component – Chow/Liu

Hi Detleffff,

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

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

❝ σ²_S=max(0,σ²_S)

❝ σ²_e=σ²_e if MS_inter≥MS_intra

❝ σ²_e=σ² if MS_inter<MS_intra

❝ where σ² = (SS_inter+SS_intra)/(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(V_inter)=V_e+2V_S
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 V_e 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
(3963 d 08:49 ago)

(edited on 2014-02-03 13:11)
@ ElMaestro
Posting: # 12320
Views: 27,346

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   -> CV_intra= 16.6%



ln(Cmax)

   V(subject(sequence))= 0.002918  -> CV_inter=  5.4%

   V(residual)         = 0.01780   -> CV_intra= 13.4%

Variance parameters are by default in SAS restricted to ≥0 in the maximization of the likelihood.
Reasonable to me :-D

:-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  -> CV_intra= 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:

: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
(3963 d 07:50 ago)

@ d_labes
Posting: # 12321
Views: 27,018

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

:-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:

:pirate:

—
Pass or fail!
ElMaestro

d_labes ★★★ Berlin, Germany, 2014-02-03 14:16 (3963 d 07:32 ago) @ ElMaestro Posting: # 12322 Views: 26,886	Variance components – Proc mixed 90% CIs Post reply
	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 . — Regards, Detlew

Helmut
★★★

Vienna, Austria,
2014-02-03 15:14
(3963 d 06:33 ago)

@ ElMaestro
Posting: # 12323
Views: 27,036

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; 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 CV_intra. 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 variances) always conservative? I think so, but I’m lacking the intellectual horsepower to prove it.

—
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

d_labes
★★★

Berlin, Germany,
2014-02-03 16:54
(3963 d 04:53 ago)

@ Helmut
Posting: # 12324
Views: 27,063

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
★★★

Vienna, Austria,
2014-02-03 17:16
(3963 d 04:31 ago)

@ d_labes
Posting: # 12325
Views: 26,960

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 Procedures for Bioequivalence 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 🖖🏼 Довге життя Україна!
Helmut Schütz

[image]

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

yjlee168
★★★

Kaohsiung, Taiwan,
2014-02-03 21:43
(3963 d 00:04 ago)

@ Helmut
Posting: # 12327
Views: 26,985

Variance components – Proc mixed

Dear Helmut,

Yes, I noticed that too.

❝ [...] but result in values smaller than CV_intra. 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.2:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here

Helmut
★★★

Vienna, Austria,
2014-02-03 23:08
(3962 d 22:39 ago)

@ yjlee168
Posting: # 12328
Views: 27,270

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 constraint of 0, and the NOBOUND option allows their estimates to be negative.

Parameter Constraints

By default, some covariance parameters are assumed to satisfy certain boundary constraints during the Newton-Raphson algorithm. For example, variance components are constrained to be nonnegative […]. You can remove these constraints […] with the NOBOUND option in the PROC MIXED statement, but this can lead to estimates that produce an infinite likelihood. […]
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 component estimate can equal zero; in this case, you might want to drop the corresponding random effect from the model. However, be aware that changing the model in this fashion can affect degrees-of-freedom calculations.

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 🖖🏼 Довге життя Україна!
Helmut Schütz

[image]

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

yjlee168
★★★

Kaohsiung, Taiwan,
2014-02-03 21:22
(3963 d 00:25 ago)

(edited on 2014-02-03 21:55)
@ d_labes
Posting: # 12326
Views: 26,992

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 CV_inter is NaN, if(is.nan(CV_inter){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.2:- 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
(3962 d 22:36 ago)

@ yjlee168
Posting: # 12329
Views: 26,839

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
★★★

Kaohsiung, Taiwan,
2014-02-04 14:09
(3962 d 07:38 ago)

(edited on 2014-02-04 17:06)
@ ElMaestro
Posting: # 12335
Views: 26,909

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.1562085²), CV_intra(15.7%) and V(subj(seq)) (1.11*10^-06)² for AUC_0-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.2:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here

yjlee168
★★★

Kaohsiung, Taiwan,
2014-02-05 20:12
(3961 d 01:35 ago)

(edited on 2014-02-06 08:02)
@ yjlee168
Posting: # 12349
Views: 27,369

bear for 2x2x2 study with negative variance components

Dear all,

Using lme() to a 2x2x2 BE study, we can code with R (taking AUC_0-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)² ≈ 0 = V(subject(seq))



 Formula: ~1 | seq %in% subj

         (Intercept)  Residual

StdDev: 1.100023e-06 0.1562085 <- (0.156)² = 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 CV_intra 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.2:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here

Admin contact
23,336 posts in 4,902 threads, 1,669 registered users;
44 visitors (0 registered, 44 guests [including 9 identified bots]).
Forum time: 21:48 CET (Europe/Vienna)

Every man is fully satisfied that there is such a thing as truth,
or he would not ask any question. Charles Sanders Peirce

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC

Ing. Helmut Schütz