14 replies to this topic

[krinaj@gmail.com]

Hello,

 

I am fairly new to PML coding and am trying to reproduce a pop PK model in Phoenix NLME that we had earlier run in NONMEM. I want to add assay type (categorical with values 0 or 1) as a covariate on proportional residual errors. I am able to run the model with below mentioned code; however, I am getting really high estimates for proportional error variance ~200. We are getting error of around 40% in NONMEM, so I know I am making some error here.

 

Can you please check below mentioned code and suggest if I am making any error in observe statemetn?

 

Thank you,

Krina

 

test(){

cfMicro(A1, Cl / V, Cl2 / V, Cl2 / V2)

dosepoint(A1)

C = A1 / V

error(CEps = 0.78)

observe(CObs = (C * (1 + CEps))*(dCObsdPKASY1*(PKASY==1)))

stparm(V = tvV * (WTBL/70)^dVdWTBL * exp(dVdDIAL1*(DIAL==1)) * exp(nV))

stparm(V2 = tvV2 * exp(nV2))

stparm(Cl = tvCl * (WTBL/70)^dCldWTBL * exp(dCldDIAL1*(DIAL==1)) * exp(dCldPE1*(PE==1)) * exp(nCl))

stparm(Cl2 = tvCl2 * exp(nCl2))

fcovariate(PKASY())

fcovariate(WTBL)

fcovariate(PE())

fcovariate(DIAL())

fixef(tvV = c(, 6.22151, ))

fixef(tvV2 = c(, 5.0984, ))

fixef(tvCl = c(, 0.0182788, ))

fixef(tvCl2 = c(, 0.0279412, ))

fixef(dVdWTBL(enable=c(0)) = c(, 0.816193, ))

fixef(dCldWTBL(enable=c(1)) = c(, 0.697362, ))

fixef(dVdDIAL1(enable=c(2)) = c(, 0.0156024, ))

fixef(dCldDIAL1(enable=c(3)) = c(, 0.352684, ))

fixef(dCldPE1(enable=c(4)) = c(0, 1, ))

fixef(dCObsdPKASY1(enable=c(5)) = c(, 1.27989, ))

ranef(block(nV, nCl) = c(0.0060795062, 0.035842864, 0.26366941), diag(nV2, nCl2) = c(1.2030433, 0.92190507))

}

Edited by krinaj@gmail.com, 01 November 2016 - 10:09 PM.

[serge guzy]

Dear Krina

I think your code is not correct.

What you want to say is that the cv is different for the 2 assays. This should change to me only the value of CEps.

Therefore

 

observe(CObs = (C * (1 + CEps))*(dCObsdPKASY1*(PKASY==1)))

 

is not correct and should be (my call).

 

observe(CObs = (C * (1 + CEps*dCObsdPKASY1*(PKASY==1)     )))

 

 

Start with CEps=0.1 and dCObsdPKASY1=1. This is the assumption they share the same cv.

Then for PKASY=1, you will get CEps*dCObsdPKASY1 as the second cv after fitting.

Let me know if it works better.

best REgards

Serge

[krinaj@gmail.com]

Dear Krina

I think your code is not correct.

What you want to say is that the cv is different for the 2 assays. This should change to me only the value of CEps.

Therefore

 

observe(CObs = (C * (1 + CEps))*(dCObsdPKASY1*(PKASY==1)))

 

is not correct and should be (my call).

 

observe(CObs = (C * (1 + CEps*dCObsdPKASY1*(PKASY==1)     )))

 

 

Start with CEps=0.1 and dCObsdPKASY1=1. This is the assumption they share the same cv.

Then for PKASY=1, you will get CEps*dCObsdPKASY1 as the second cv after fitting.

Let me know if it works better.

best REgards

Serge

Many thanks, Serge.

 

I have one more question to clarify, which I have been a bit confused about. When Phoenix reports stdev0 in Structure Tab or in Theta table, does it refer to variance or SD/CV?

 

Thanks,
Krina

[Ana Henry]

Many thanks, Serge.

 

I have one more question to clarify, which I have been a bit confused about. When Phoenix reports stdev0 in Structure Tab or in Theta table, does it refer to variance or SD/CV?

 

Thanks,
Krina

Phoenix reports the SD not the variance so you would need to square it to compare to your NONMEM results.

[krinaj@gmail.com]

thanks, Ana.

[krinaj@gmail.com]

Dear Krina
I think your code is not correct.
What you want to say is that the cv is different for the 2 assays. This should change to me only the value of CEps.
Therefore

observe(CObs = (C * (1 + CEps))*(dCObsdPKASY1*(PKASY==1)))

is not correct and should be (my call).

observe(CObs = (C * (1 + CEps*dCObsdPKASY1*(PKASY==1) )))


Start with CEps=0.1 and dCObsdPKASY1=1. This is the assumption they share the same cv.
Then for PKASY=1, you will get CEps*dCObsdPKASY1 as the second cv after fitting.
Let me know if it works better.
best REgards
Serge

Hi Serge and Ana,

I tried several different ways to run this model; I am still getting multiplicative stdev0 estimate in double digits (range 16-80).

Since Phoenix presents stdev estimate in percent format, it should be less than 1 for multiplicative model, right? It seems like I am messing up something by 2 digits. I am using the below mentioned code. Any comments/tips will be appreciated.

EDIT: I am able to resolve the issue with double digit multiplicative stdev0 estimate by using add+mult. Is there a way to give bounds to additive portion of error in combined error model?

Thanks,
Krina

test(){
cfMicro(A1, Cl / V)
dosepoint(A1)
C = A1 / V
error(CEps = 0.4)
observe(CObs = (C * (1 + CEps*dCObsdPKASY1*(PKASY==1))))
stparm(V = tvV * (WTBL/70)^dVdWTBL * (1+dVdDIAL1*(DIAL==1)) * exp(nV))
stparm(Cl = tvCl * (WTBL/70)^dCldWTBL * (1+dCldPE1*(PE==1)) * (1+dCldDIAL1*(DIAL==1)) * exp(nCl))
fcovariate(WTBL)
fcovariate(PE())
fcovariate(DIAL())
fcovariate(PKASY())
fixef(tvV = c(0, 5.02679, ))
fixef(tvCl = c(0, 0.0138065, ))
fixef(dVdWTBL(enable=c(0)) = c(0, 0.748727, ))
fixef(dCldWTBL(enable=c(1)) = c(0, 0.750154, ))
fixef(dCldPE1(enable=c(2)) = c(0, 102.508, ))
fixef(dVdDIAL1(enable=c(3)) = c(0, 0.4704, ))
fixef(dCldDIAL1(enable=c(4)) = c(0, 1.51546, ))
fixef(dCObsdPKASY1(enable=c(5)) = c(0, 0.800001, ))
ranef(diag(nV, nCl) = c(0.10007425, 0.10009075))
}

Edited by krinaj@gmail.com, 03 November 2016 - 01:59 AM.

[mittyright]

Hi Krina,

 

I'm wondering how does your model work. As I see from this statement:

observe(CObs = (C * (1 + CEps*dCObsdPKASY1*(PKASY==1))))

when PKASY is not 1 then CObs = C, that cannot be right.

I would substitute it to

observe(CObs = C * (1+ (1+dCObsdPKASY1*(PKASY==1))*CEps))

Then 

IF PKASY=0 -> CObs = C * (1+(1+0)*CEps)= C + C* CEps

IF PKASY!=0 -> CObs = C + C * CEps + C*CEps*dCObsdPKASY1

 

I kindly recommend you to review the Chapter Modeling circadian rhythm, inter-occasion variability, and two epsilons in PHX7 User's Guide if you are curious about epsilon manipulation (and many other interesting things)

 

Is there a way to give bounds to additive portion of error in combined error model?

I never heard about such a request. Sigma is randomly distributed along the data with mean 0 and fitted SD. You cannot set a bound for distribution of error.

If you can provide some NONMEM model text, it would be useful for understanding of your request. 

 

I hope it helps,

Mittyright

[krinaj@gmail.com]

Hi Krina,

 

I'm wondering how does your model work. As I see from this statement:

when PKASY is not 1 then CObs = C, that cannot be right.

I would substitute it to

observe(CObs = C * (1+ (1+dCObsdPKASY1*(PKASY==1))*CEps))

Then 

IF PKASY=0 -> CObs = C * (1+(1+0)*CEps)= C + C* CEps

IF PKASY!=0 -> CObs = C + C * CEps + C*CEps*dCObsdPKASY1

 

I kindly recommend you to review the Chapter Modeling circadian rhythm, inter-occasion variability, and two epsilons in PHX7 User's Guide if you are curious about epsilon manipulation (and many other interesting things)

 

I never heard about such a request. Sigma is randomly distributed along the data with mean 0 and fitted SD. You cannot set a bound for distribution of error.

If you can provide some NONMEM model text, it would be useful for understanding of your request. 

 

I hope it helps,

Mittyright

 

Thanks, Mittyright. 

I will surely check out user guide for the chapters you suggested.

 

The above model worked fine. I think the model suggest that when PKASY is 0, CObs = C + C*CEps, i.e. normal proportional residual error model and when PKASY is 1 there is additional coefficient, so CObs=C + C*CEps*Coefficient. 

 

observe(CObs = (C * (1 + CEps*dCObsdPKASY1*(PKASY==1))))

 

However, this model gave very high estimate for CEps (even without PKASY part) and some trend in CWRES vs PRED plot, so I ended up using combined additive+multiplicative error model with PKASY effect on multiplicative portion only. It worked well. Then, I wanted to try to have additive portion as small as possible and have the rest of error distribute towards proportional error, if possible. I could fix additive error to small number, which increased proportional portion of combined error, but I wanted to see if there is an option to give bound to additive error rather than fixing.

[mittyright]

Hi Krina,

 

The above model worked fine. I think the model suggest that when PKASY is 0, CObs = C + C*CEps, i.e. normal proportional residual error model and when PKASY is 1 there is additional coefficient, so CObs=C + C*CEps*Coefficient. 

 

observe(CObs = (C * (1 + CEps*dCObsdPKASY1*(PKASY==1))))

 

A little algebra if you don't mind using statement above::

IF PKASY=1 -> CObs = (C * (1 + CEps*dCObsdPKASY1*(PKASY==1))) = (C * (1+CEps*dCObsdPKASY1*1)) = C + C*CEps * dCObsdPKASY1

OK, proportional residual error

IF PKASY=0 -> CObs = (C * (1 + CEps*dCObsdPKASY1*(PKASY==1))) = (C * (1+CEps*dCObsdPKASY1*0)) = C + 0 = C

No residual error at all

 

BTW I'm happy that you came up with solution using combined residual model.

May be I  can answer to your question regarding error model modification if you post it here.

Thanks

 

BR,

Mittyright

[krinaj@gmail.com]

Hi Krina,

 

 

A little algebra if you don't mind using statement above::

IF PKASY=1 -> CObs = (C * (1 + CEps*dCObsdPKASY1*(PKASY==1))) = (C * (1+CEps*dCObsdPKASY1*1)) = C + C*CEps * dCObsdPKASY1

OK, proportional residual error

IF PKASY=0 -> CObs = (C * (1 + CEps*dCObsdPKASY1*(PKASY==1))) = (C * (1+CEps*dCObsdPKASY1*0)) = C + 0 = C

No residual error at all

 

BTW I'm happy that you came up with solution using combined residual model.

May be I  can answer to your question regarding error model modification if you post it here.

Thanks

 

BR,

Mittyright

 

Hi Mittyright,

 

I do not know if (PKASY==0) statement literally multiplies with 0. I was thinking it only means that when PKASY is 0.

 

Let's say if we have a covariate with more categories, i.e. Race represented by values 0, 1, 2, 3, 4 and we use above type proportional covariate model (+1 option in built in Phoenix). Are the estimates for Race == 2, 3, or 4 really derived after multiplying with 2, 3, or 4? That cant'  be right in terms of interpretation.

 

Hope someone from Certara can please clarify this.

 

Thanks,
Krina

[mittyright]

Dear Krina,

 

Hope someone from Certara can please clarify this.

 

Don't you trust my experience? 

 

BTW it works as intended (and I'm sure that NONMEM doesn't have an alternative way)

You can try it by yourself, add this line to your code:

secondary(flag_PKASY=(PKASY==0))

and switch to individual mode (do not forget to sort by id)

And then go to the secondary tab

You'll see that

if PKASY=0 -> flag_PKASY=1

if PKASY=1 -> flag_PKASY=0

if PKASY=2 -> flag_PKASY=0

if PKASY=100000 -> flag_PKASY=0

and so on

 

That's because 

(PKASY==0)

is equal to

IF (PKASY.EQ.0) THEN

FLAGPKASY=1

ELSE

FLAGPKASY=0

ENDIF

NM speaking.

 

Let's say if we have a covariate with more categories, i.e. Race represented by values 0, 1, 2, 3, 4 and we use above type proportional covariate model (+1 option in built in Phoenix). Are the estimates for Race == 2, 3, or 4 really derived after multiplying with 2, 3, or 4?

 

It would be like this:

observe(CObs = C * (1+ (1+dCObsdPKASY1*(PKASY==1)+dCObsdPKASY2*(PKASY==2)+dCObsdPKASY3*(PKASY==3))*CEps))

So you'd need to describe coeffs for all categories.

 

If you want more difficult conditions, you need to use ternary operator. For example:

flag_PKASY= (PKASY ==1 ? 10 : PKASY ==2? 100: 0 )

 

Please see more about ternary operator here

 

All the best,

Mittyright

Edited by mittyright, 06 November 2016 - 06:42 PM.

[krinaj@gmail.com]

Dear Krina,

Don't you trust my experience?

BTW it works as intended (and I'm sure that NONMEM doesn't have an alternative way)
You can try it by yourself, add this line to your code:
secondary(flag_PKASY=(PKASY==0))
and switch to individual mode (do not forget to sort by id)
And then go to the secondary tab
You'll see that

if PKASY=0 -> flag_PKASY=1

if PKASY=1 -> flag_PKASY=0

if PKASY=2 -> flag_PKASY=0

if PKASY=100000 -> flag_PKASY=0

and so on

That's because

(PKASY==0)

is equal to

IF (PKASY.EQ.0) THEN

FLAGPKASY=1

ELSE

FLAGPKASY=0

ENDIF

NM speaking.

It would be like this:

observe(CObs = C * (1+ (1+dCObsdPKASY1*(PKASY==1)+dCObsdPKASY2*(PKASY==2)+dCObsdPKASY3*(PKASY==3))*CEps))

So you'd need to describe coeffs for all categories.

If you want more difficult conditions, you need to use ternary operator. For example:

flag_PKASY= (PKASY ==1 ? 10 : PKASY ==2? 100: 0 )

Please see more about ternary operator here

All the best,
Mittyright

Many thanks, Mittyright.

No, I do not mean to not trust your answer. Sorry about that. It's just that the answer was a little more out of what I thought and shocked me.

One more clarification that I need:
With Race example below, if I use built in 1+ proportional covariate model, and lets say if we consider that Race==2, do I need to divide the coefficient of Race 2 by 2 to get correct coefficient?

CLi = tvCl * (1+ dCldRACE1 * (RACE ==1)) * (1+ dCldRACE2 * (RACE ==2)) * (1+ dCldRACE3 * (RACE ==3)) * exp (nCl)

When Race == 2,
CLi = tvCl * (1+ dCldRACE2 * 2) * exp (nCl)

Thus,
actual dCldRACE2 == dCldRACE2/2

I have never divided coefficient by flag number or seen anyone doing that to interpret coefficients. And that is why I was confused and shocked.

I will really appreciate if you can clarify.

Thanks,
Krina

CLi = tvCl * (1+ dCldRACE1 * (RACE ==1)) * (1+ dCldRACE1 * (RACE ==1)) * (1+ dCldRACE1 * (RACE ==1)) * exp (nCl)

Edited by krinaj@gmail.com, 06 November 2016 - 07:29 PM.

[krinaj@gmail.com]

Many thanks, Mittyright.

 

No, I do not mean to not trust your answer. Sorry about that. It's just that the answer was a little more out of what I thought and shocked me.

 

One more clarification that I need:

With Race example below, if I use built in 1+ proportional covariate model, and lets say if we consider that Race==2, do I need to divide the coefficient of Race 2 by 2 to get correct coefficient?

 

CLi = tvCl * (1+ dCldRACE1 * (RACE ==1)) * (1+ dCldRACE2 * (RACE ==2)) * (1+ dCldRACE3 * (RACE ==3)) * exp (nCl)

 

When Race == 2,

CLi = tvCl * (1+ dCldRACE2 * 2) *  exp (nCl)

 

Thus,

actual dCldRACE2 == dCldRACE2/2

 

I have never done above or seen anyone doing that to interpret coefficients. And that is why I was confused and shocked.

 

I will really appreciate if you can clarify.

 

Thanks,
Krina

Edited by krinaj@gmail.com, 06 November 2016 - 07:29 PM.

[mittyright]

Dear Krina,

With Race example below, if I use built in 1+ proportional covariate model, and lets say if we consider that Race==2, do I need to divide the coefficient of Race 2 by 2 to get correct coefficient?

 

CLi = tvCl * (1+ dCldRACE1 * (RACE ==1)) * (1+ dCldRACE2 * (RACE ==2)) * (1+ dCldRACE3 * (RACE ==3)) * exp (nCl)

 

When Race == 2,

CLi = tvCl * (1+ dCldRACE2 * 2) *  exp (nCl)

 

 

When Race = 2 ->

(RACE==1) -> 0

(RACE==2) -> 1

(RACE==3) -> 0

So '==' is just a condition to check. For simplicity:

when you are typing (RACE==2) you are saying to the engine: 'please, if Race is equal to 2, treat it as 1, otherwise - all other values - treat it as 0'

when you are typing (RACE==999) you are saying to the engine: 'please, if Race is equal to 999, treat it as 1, otherwise - all other values - treat it as 0'

 

So you don't need any additional coefficients.

 

BR,

Mittyright

[krinaj@gmail.com]

Dear Krina,

 

When Race = 2 ->

(RACE==1) -> 0

(RACE==2) -> 1

(RACE==3) -> 0

So '==' is just a condition to check. For simplicity:

when you are typing (RACE==2) you are saying to the engine: 'please, if Race is equal to 2, treat it as 1, otherwise - all other values - treat it as 0'

when you are typing (RACE==999) you are saying to the engine: 'please, if Race is equal to 999, treat it as 1, otherwise - all other values - treat it as 0'

 

So you don't need any additional coefficients.

 

BR,

Mittyright

Hi Mittyright,

 

Okay, that is very helpful and clear now.

Many thanks,

Krina