- On 9 Jul 2013 at 21:46:48, Ayyappa Chaturvedula sent the message

Back to the Top

Dear group,

I am trying to use combination error model in Phoenix. I did not find an easier way than using in

text editor mode. In the manual it mentions Cobs=C*(1+eps1)+eps1*scale1 as correct format where we

declare scale1 as a fixed effect. But this gives an error saying a repeat of eps1 in observ(). I

just tried Cobs=C*(1+eps1)+scale1 instead and worked fine. Now my questions are:

1. Is there an easier way using custom option for combination residual error?

2. That way I coded additive portion, is it like coding epsilon as theta in NONMEM? Can I use the

value directly as SD?

Thanks for help.

Regards,

Ayyappa - On 9 Jul 2013 at 23:12:59, Mukul Minocha sent the message

Back to the Top

Hi Ayyappa

In the Graphic user interface of phoenix NLME, under STRUCTURE tab you will

find 'residual error' structure drop down menu, there you will have several

options to select like additive, multiplicative or 'mixed' error model. I

guess'mixed' (combined additive and proportional) is the one that you are

looking for.

The code is:

observe(CObs = C + CEps * (1 + C * CMixRatio))

where CMixratio is estimated.

Hope this helps

Best

Mukul Minocha

Research Scientist

CTM

University of Maryland Baltimore - On 10 Jul 2013 at 08:51:04, Nick Holford sent the message

Back to the Top

Awappa.

You wrote:

"I am trying to use combination error model in Phoenix. I did not find an easier way than using in

text editor mode. In the manual it mentions Cobs=C*(1+eps1)+eps1*scale1 as correct format where we

declare scale1 as a fixed effect. But this gives an error saying a repeat of eps1 in observ(). I

just tried Cobs=C*(1+eps1)+scale1 instead and worked fine."

The model you coded may have run successfully but it is not a model for a combined residual error

model that you are trying to achieve. In NM-TRAN terms what you coded was:

C=F ; concentration prediction from PREDPP

Y=C*(1+EPS(1)) + THETA(x)

i.e. you predict concentration (C) with a PK model that has a proportional residual error and also

estimate an additional baseline concentration that is independent of dose (THETA(x)).

Phoenix NLME offers a way of coding a combined residual error model but it is not easily comparable

to what you may be familiar with using NM-TRAN eg.

observe(CObs = C + CEps * (1 + C * CMixRatio)

In NM-TRAN this might be written like this:

$SIGMA 1 FIX ; EPS1

C=F ; concentration prediction from PREDPP

SD=THETA(1)

SDCV=THETA(2)

Y=C + SD*EPS(1)*(1 + C*CV2SD)

The SD parameter is the std deviation of the additive residual error.

The CV2SD parameter is the ratio of CV/SD where CV is the co-efficient of variation of the

proportional residual error.

That is why Phoenix uses the name CMixRatio because what is estimated is the ratio of the CV to the

SD. For example if the SD was 2 mg/L and the CV was 0.1 then the CV2SD parameter would have a value

of 0.05. This means that when you estimate CMixRatio you need to do a bit of extra calculation to

get a CV value that would be understood in the mixed effects modelling literature (i.e.

CV=SD*CV2SD).

Best wishes,

Nick

--

Nick Holford, Professor Clinical Pharmacology

Dept Pharmacology & Clinical Pharmacology, Bldg 503 Room 302A

University of Auckland,85 Park Rd,Private Bag 92019,Auckland,New Zealand

email: n.holford.-a-.auckland.ac.nz

http://holford.fmhs.auckland.ac.nz/ - On 10 Jul 2013 at 08:59:59, Simon Davis sent the message

Back to the Top

Dear Ayyappa,

Mukul is indeed correct, if you want to discuss further the Phoenix discussion forum has the

advantage that you can post data files, screen shots and Phoenix projects to help debug a problem

There is a related discussion here;

http://www.pharsight.com/extranet/index.php?option=com_kunena&Itemid=55&func=view&catid=51&id=1001&

limit=6&limitstart=12#1358

Note when you choose the drop down, you will see the code updated in the statements. The term

CMixRatio can be your own choice of variable name.

[by Simon - db]

Best regards,

Simon.

PS for others who are interested in the forum it is a free but separate sign-up currently to our

support system;

http://www.pharsight.com/extranet/index.php?option=com_comprofiler&task=registers

If you have registered in the past but not been granted access, please check your spam folder and

please click the authorization /activation emails you have received to confirm your address. IF

they are in your junk folder then please mark as Not Junk to ensure future notifications come

through.

Password and user name reminders can be accessed via;

http://www.pharsight.com/extranet/index.php?option=com_comprofiler&task=lostpassword

If you still have problem then you may contact me directly

--

Simon.Davis.-at-.certara.com

Senior Scientific Consultant

Pharsight- A Certara™ Company - On 10 Jul 2013 at 09:01:23, Ayyappa Chaturvedula sent the message

Back to the Top

Mukul,

Thank you, how simple that is. I see now that is a simple rearrangement that avoids repeat of eps.

Regards,

Ayyappa - On 10 Jul 2013 at 12:49:32, Brendan Johnson sent the message

Back to the Top

So I have always wondered about this, using a single 'eps' for a combined proportional and additive

residual error model (sometimes seen when using SIGMA like THETAs in NONMEM, and using one FIXED

SIGMA=1 and two THETAs to code a combined error model).

These approaches force the sign of the eps value to be the same for both the proportional and

additive portions of the model. Should this (or would this) always be the case? Would it not make

more sense to estimate a different SIGMA and use two independent eps values when using these

combined models?

Brendan - On 10 Jul 2013 at 13:30:57, Devin Pastoor sent the message

Back to the Top

Hi Brendan,

From the standpoint of phoenix, you cannot estimate two residual error terms due to a constraint of

the software. You can only have one error (eps) term per observe statement, therefore the

work-around is to use a scaling factor such that you can derive one component as a proportion of the

other component (to allow for 'different' eps values). To estimate as such, the scaling factor is

introduced as a fixed effect to allow for estimation.

Ie: C*(1+eps1) + Eps2

With Eps2 written as Eps1*scalingfactor.

Again, as mentioned by Nick, the resulting values for eps1 and eps2 are NOT comparable to nonmem due

to phoenix being in terms of standard deviation while nonmem handles in terms of CV.

Devin Pastoor

Clinical Research Scientist

University of Maryland, Baltimore

Center for Translational Medicine - On 11 Jul 2013 at 09:03:04, Simon Davis sent the message

Back to the Top

Following user feedback, the next release of Phoenix will have two separate built in choices related

to mixed additive proportional residual models, namely the existing ‘mixratio’ and also ‘add+mult’

‘Mixratio’ results in

observe(EObs(C) = E + EEps * (1 + E * EMixRatio))

where EEps is the ‘additive’ component of the mixture std deviation, and

EEps*Emixratio (where Emixratio is a fixed effect) is the ’multiplicative’ std deviation component.

Note that PHX reports Emixratio as a fixed effect, so you have to do the extra multiplication

EEps*Emixratio yourself to get the multiplicative std deviation. The is the ‘extra’ computation

that Nick Holford refer to that you have to do and could simply be added as a Secondary parameter so

would be calculated at runtime.

Note that in NM if want this particular model, as Nick points out, you could formulate it as

$SIGMA 1 FIX ; EPS1

SD=THETA(1)

CV2SD = THETA(2)

Y=E + SD*EPS(1)*(1 + E*CV2SD)

Here SD is now the standard deviation of the additive component, and SD*CV2SD is the std deviation

of the multiplicative component.

But you would still need to do the extra multiplication SD*CV2SD to get the std deviation of the

multiplicative component.

In the next release (Phoenix 1.4), if you choose add+mix, you now get variance mixing equivalent to

the usual NM construction Y= E+eps(1) + E*eps(2)

The PML that is generated is

observe(EObs(C) = E + EEps * sqrt(1 + E^2 * (MultStdev/sigma())^2))

stparm(MultStdev = tvMultStdev)

fixef(tvMultStdev = c(, 1, ))

Note here sigma() is a special reserved word representing the std deviation of EEps, and was added

to avoid the necessity of doing the extra multiplication - now the value reported for EEps is the

additive standard deviation of the additive component (equivalent to the sqrt(sigma(1)) where

sigma(1)=var(eps(1) in the NM formulation), and the value reported for the fixed effect;

tvMultStdev is the multiplicative standard deviation (equivalent to sqrt(sigma(2) in the NM

formulation).

Best regards,

Simon. - On 11 Jul 2013 at 09:27:33, Brendan Johnson sent the message

Back to the Top

Devin, understood, but my question was related to the sign of this single Eps per observed value

In this case

C*(1+Eps1) + Eps2

Eps1 and Eps2 can have different signs, positive/negative/any combination

If Eps2 is written as Eps1*scalingfactor. Then you force Eps2 to have the same sign as Eps1

Question is, would that or should that always be the case? Intuitively, I would say no, but since

this method of writing a combined error model is fairly common, wanted to know what others thought

Brendan Johnson - On 11 Jul 2013 at 11:56:49, Devin Pastoor sent the message

Back to the Top

Hi Brendan,

I think you may be misunderstanding what the epsilon value points to. The epsilon component of the

residual error (in phoenix) corresponds to the standard deviation of the distribution of residual

error with a mean of zero and a variance --> stdev^2.

As for the reason for having multiple eps, this can be very situational. For example, let us say you

know your assay only has a certain level of precision for measurements. This residual variability

could be handled by a additive component as it this error is homoscedastic (constant). On the other

hand the other residual error component that covers things like within subject variability, model

misspecification, etc may be heteroscedastic (for example often seen in PK models where you see

concentration dependence). In this situation, two independent error distributions better describe

your data, thus why you would want two Eps.

Does this answer your question as to why Eps1, etc…. Will always have a positive sign – they don’t

correspond to the difference between Obs and pred, instead they describe the distribution of

potential differences which take on positive and negative values).

Devin Pastoor

Clinical Research Scientist

University of Maryland, Baltimore

Center for Translational Medicine - On 11 Jul 2013 at 14:50:02, Nick Holford sent the message

Back to the Top

Brendan,

You ask an excellent question.

Depending on the details of the estimation algorithm it may make no difference when estimating

parameters (just a guess).

But for simulation then I would think that NONMEM would sample independent Eps1 and Eps2 values

which could have different signs and sizes (~N(0,Eps1) and ~N(0,Eps2)).

On the other hand it seems that Phoenix NLME would have to use the same random sample of epsilon and

therefore the sign and size of the random effect would be the same for both the additive and

proportional components of the residual error.

The same thing would happen when simulating using NONMEM if only one EPS was used and the random

effect model was implemented with THETAs.

Best wishes,

Nick

--

Nick Holford, Professor Clinical Pharmacology

Dept Pharmacology & Clinical Pharmacology, Bldg 503 Room 302A

University of Auckland,85 Park Rd,Private Bag 92019,Auckland,New Zealand

email: n.holford.-at-.auckland.ac.nz

http://holford.fmhs.auckland.ac.nz/ - On 11 Jul 2013 at 15:38:58, Devin Pastoor sent the message

Back to the Top

Dear Nick,

Do you mind clarifying what you mean by different signs for Eps1/Eps2.

My understanding is that with all MLE methods you cannot get a negative variance given it is equal

to sigma^2, and even when using techniques like ANOVA that may give a negative variance as output it

would indicative of some issues such as sample size, outliers, not the right model for the

covariance structure, or that the true variance is 0.

Thanks,

Devin Pastoor - On 11 Jul 2013 at 23:20:12, Brendan Johnson sent the message

Back to the Top

Devon, sorry, I am using NONMEM-speak.

So as Nick refers to - I am using Eps1 and Eps2 as samples from N(0,sigma1) and N(0,sigma2)

distributions. So yes, sigma1 and sigma2 must be positive. I see the nomenclature is different in

PHX (sigma=Eps in PHX speak).

Brendan - On 11 Jul 2013 at 23:20:55, Devin Pastoor sent the message

Back to the Top

Dear Nick, Brendan and Simon,

I believe I may have misread Brendan's initial question. To make sure we are on the same page:

In nonmem given a situation such that

Yobs = Ypred*(1+ Err1) + Err2

Err1 and Err2 correspond to the residual random effect distributions (Eps) such that

Eps1 --> N(0, Err1)

Eps2 --> N(0, Err2)

As Eps1 and Eps2 are independent, when a sample is taken from said distributions, the sample from

Eps1 could be positive or negative, as would the sample from Eps2. These two values would not be

forced to share the same sign due to complete independence.

On the other hand, in phoenix, given the same equation parameterized in terms of one Error component

Yobs = Ypred*(1 + Err1) + Err1*scalingfactor

Given Eps --> N(0, Err1^2)

The independence of the proportional and additive components would not hold in that if only 1 sample

was taken from the eps distribution at each observation phoenix would force each component to

maintain the same sign (one could never be positive while the other negative).

This is an interesting question. To me, I would hypothesize that the results of using this approach

would result in different estimates of the variance (with phoenix underpredicting in comparison to

nonmem) due to the independence of eps in nonmem allowing the positive/negative combos whereas

phoenix is only positive/positive or negative/negative.

On the flip side, for simulation given the same distributions to sample from, phoenix would result

in more pronounced residual variability given that when samples from either extreme of the

distribution were picked, the resulting variability would always compound and you'd never see the

'dampening' of one portion being positive and one negative.

Please feel free to critique my thoughts,

Devin Pastoor

Clinical Research Scientist

University of Maryland, Baltimore

Center for Translational Medicine - On 12 Jul 2013 at 08:50:42, Leonid Gibiansky sent the message

Back to the Top

Dear All,

Concerning Nonmem/Phoenix differences, the difference is in notations only. Using Nonmem notations,

you can use two identical models

Y = IPRED*(1 + EPS1) + EPS2

(EPS1 and EPS2 are independent, with variances SIGMA1, SIGMA2) or

Y = IPRED*(1 + EPS3) + THETA(1)*EPS3

(with variance SIGMA3; I think it make sense to use different notations in these models so we can

compare them).

For these models, SD of the observations can be expressed as

SD = SQRT(IPRED**2*SIGMA1+SIGMA2)

or

SD = SQRT(IPRED**2*SIGMA3+THETA(1)**2*SIGMA3)

From here it is immediately follows that

SIGMA1=SIGMA3 and SIGMA2=THETA(1)**2*SIGMA3 (1)

If you estimate the same data with two models, you will get SIGMA1-SIGMA3 and THETA(1) estimates

that satisfy (1). If you simulate from these models, you will get exactly the same results if you

use SIGMA1-SIGMA3 and THETA(1) that satisfy (1).

THETA(1) is the same as scaling factor in Phoenix.

Thanks

Leonid

--

Leonid Gibiansky, Ph.D.

President, QuantPharm LLC

web: www.quantpharm.com

e-mail: LGibiansky at quantpharm.com - On 12 Jul 2013 at 09:27:49, Bob Leary (bob.leary.-at-.CERTARA.COM) sent the message

Back to the Top

Dear Devin -

I think its a little misleading read anything into the fact that the usual NONMEM style

additive/proportional residual error model

Yobs = Ypred*(1+ Err1) + Err2

uses two independent residual errors wheres the roughly similar (at least in intent, but actually

mathematically

distinct) Phoenix residual error model

Yobs = Ypred*(1 + Err1) + Err1*scalingfactor

uses only one.

In fact, it is easily shown that the NM (Err1,Err2) formulation is mathematically equivalent to an

implementation

that, like Phoenix, uses just one Err1 value along with a scaling factor :

Yobs=Ypred +sqrt(Ypred^2 + scalingfactor^2)*Err1

In fact, the real difference between the NM and PHX models is that the NM model in effect

expands the residual variance in a first order taylor series in Ypred^2 where the zero order

coefficient

is the 'addtive' compoment and the first order coefficient is the proportional component.

Similarly, the PHX model expands the residual standard deviation (not the variance) in a first order

Taylor series in

Ypred where the zero order coefficient can be interpreted as the 'additive' component' and the first

order coefficient

as the proportional component.

As Nick Holford has noted, the Phoenix std deviation expansion model is easily implementable in

NMTRAN, and similarly the NM variance expansion model is easily implementable in Phoenix PML. Its

just that currently the NM variance expansion version using (Err1, Err2) is the most compact and

convenient version to write in NMTRAN, while the PHX standard deviation expansion version using

(Err1, scalingfactor) is the one that is currently built into the Phoenix UI. In the next release

of PHX, both versions of the additive/proportional model will be available. - On 12 Jul 2013 at 13:12:32, Devin Pastoor sent the message

Back to the Top

Dear Bob and Leonid,

Thank you for your inputs. My main question is regarding how the software samples from the

distribution at the observation level. To use a simulation case

Given Y = IPRED*(1 + EPS1) + EPS2

It seems that for each Y a sample will be taken independently from the distributions N(0, EPS1),

N(0, EPS2). (please correct me if I am incorrect in this assumption)

However I am unsure what happens when

Y = IPRED*(1 + EPS3) + THETA(1)*EPS3

Is EPS resampled at the second time it is called for the additive component?

If not, are the additive and proportional components coerced into the same sign/place on the

curvature of the distribution? This is where my statistical understanding breaks down.

Thanks for the input!

Devin Pastoor

University of Maryland, Baltimore

Center for Translational Medicine - On 12 Jul 2013 at 13:48:40, Leonid Gibiansky sent the message

Back to the Top

> However I am unsure what happens when

> >

> Y = IPRED*(1 + EPS3) + THETA(1)*EPS3

> >

> Is EPS resampled at the second time it is called for the additive component?

> >

> If not, are the additive and proportional components coerced into the same sign/place on the

> curvature of the distribution? This is where my statistical understanding breaks down.

>

Devin,

You may re-write the model as

Y = IPRED+(IPRED + THETA(1))*EPS3

so it is sampled only once.

Leonid

--

Leonid Gibiansky, Ph.D.

President, QuantPharm LLC

web: www.quantpharm.com

e-mail: LGibiansky at quantpharm.com - On 12 Jul 2013 at 13:51:00, Nick Holford sent the message

Back to the Top

Devin,

I think we are now on the same page but let me just clarify the notation (I was a bit sloppy in my

earlier posting).

Given:

Yobs = Ypred*(1+ eps1) + eps2

I understand this to mean that eps1 and eps2 are random variables sampled from a normal distribution

such that:

eps1 ~N(0,variance1)

eps2~N(0,variance2)

The value of eps1 and eps2 can have be positive or negative anywhere in the range -infinity to

+infinity. When NONMEM performs simulation the values of eps1 and eps2 are sampled independently

from their respective normal distributions which means they have independent signs and sizes.

If in NONMEM you used this kind of model

Yobs = Ypred + eps1*sqrt((SD**2 + (CV*Ypred)**2))

then there is only one eps1 and simulating with this model would mean that the additive component

defined by SD and the proportional component defined by CV must have the same sign and size

determined by eps1.

As I noted previously I do not know if this will be true when using NONMEM for estimation but it

seems that it should be the same because only one variance is estimated.

Nick - On 12 Jul 2013 at 13:57:10, Bob Leary sent the message

Back to the Top

Devin:

If you run

Y = IPRED*(1 + EPS3) + THETA(1)*EPS3

In simulation mode, you will get only 1 sample EPS3 for each observation

(coded either in PHOENIX PML or NMTRAN)

I think you are making a fundamental mistake in assuming that there really are separate additive

and proportional components in this case -

There is only one random residual error with a normal distribution that has a standard deviation

equal to abs((IPRED + THETA1)) *stddev(EPS3).

The interpretation that there is a separate additive and proportional components is a conceptual

fiction - these separate components play no real

part in the model - it is only the combined residual which has a normal distribution with a certain

modeled std deviation value that matters.

In the NM case where you use

Y = IPRED*(1 + EPS1) + EPS2

In simulation mode in NM you will actually get separate independent samples for EPS1 and EPS2.

But if EPS1 and EPS2 appear nowhere else in the model, again the idea of a separate additive and

proportional parts is simply a conceptual

fiction - in this case the overall residual error is comes from a Normal distribution that has a

certain modeled variance . The model

does not care whether you sample EPS1 ~N(0,SIGMA1) SIGMA1 and EPS2~N(0,SIGMA2) separately and

combine them as a residual

error IPRED*EPS1 + EPS2, or whether you sample from a single normal distribution EPS3

~N(0,SIGMA1*IPRED^2 + SIGMA2)

In this case the use of EPS1 and EPS2 separately is just a convenient way to generate the final

residual error with the right variance -

Whether or not we know the separate additive and proportional components makes no difference and

this knowledge, if indeed we do know it,

is never used anywhere else in the model. - On 12 Jul 2013 at 23:25:27, Roger Jelliffe sent the message

Back to the Top

Dear All:

About the error model. We at the USC LAPK have been using another approach which we have thought

to be useful. We break the error up into that which is due to the assay and that which is due to the

uncertainties in the clinical environment in which the study has been done. So we first separately

estimate the assay error as a polynomial in which the assay SD = AoCo + A1C1 + A2C2 = A3CC3 if need

be. A0, A1, A2, and A3 are coefficients, and C0, C1, C2, and C3 are the measured serum

concentration, for example, raised to the power of 0, 1, 2, or 3. (I can't do superscripts on this

email). In addition, we used to multiply this polynomial by another parameter, Gamma, which we

estimated separately as a point estimate [1]. We found, though, that gamma, when used this way,

multiplied up the curvature in the assay error polynomial,, and was probably not the corr4ect thing

to do. Instead, we now estimate another additive term which we call lambda, which we then add to the

assay error polynomial to provide an estimate of the overall noise with which the clinical part of

the study was done.

In this way we first get the assay error polynomial totally apart from the study itself [1], and

then estimate the noise of the clinical part as lambda. Good or poor assay, that is self-evident.

Small lambda, a good clinical study, with small environmental noise due to errors is the preparation

and administration of doses, the errors of recording the timing of doses and the samples, the model

misspecification, and any possible changes in model parameter values during the period of data

analysis. Large lambda, maybe one ought to see what can be done to make the clinical part of the

study more precise. In this way we then have a means to estimate both components of this overall

error. We think this is a more informative approach rather than simply estimating the overall error

and never knowing how much is due to what. One might consider using this approach not only in the

Pmetrics software for nonparametric PK modeling, but also in your various approaches to error

estimation.

Another benefit of this approach is that we weight each assay measurement by the reciprocal of

its variance at that value. In doing this, we avoid the errors which are inevitable as the

measurements approach zero. While the CV% gets bigger, and becomes infinite at zero, the assay

variance becomes the machine noise of the assay is still always finite and can provide the best

estimate of the measurement at any value. This then avoids the illusion of the apparent need to

censor low values which have CV% values are felt to be unacceptable. In this way one has the correct

weight of each assay measurement plus the environmental noise, no need to censor low assay

measurements, and in our view, a more informative study. Some relevant references are below.

Lastly, you might consider using nonparametric pop modeling, as you can either with Pmetrics or

Phoenix. It avoids the constraints of assuming parametric shapes of the parameter distributions and

thus gets results that actually are more likely [3], and in contrast to what is possible with

parametric methods, one can develop dosage regimens which actually hit therapeutic targets with

maximal precision (minimum expected weighted squared error), using the method of Multiple Model

dosage design [4,5].

References

1. Jelliffe RW, Schumitzky A, Van Guilder M, Liu M, Hu L, Maire P, Gomis P, Barbaut X, and Tahani

B: Individualizing Drug Dosage Regimens: Roles of Population Pharmacokinetic and Dynamic Models,

Bayesian Fitting, and Adaptive Control. Therapeutic Drug Monitoring, 15: 380-393, 1993.

2. Neely M, van Guilder M, Yamada W, Schumitzky A, and Jelliffe R: Accurate Detection of Outliers

and Subpopulations with Pmetrics, a Nonparametric and Parametric Pharmacometric Modeling and

Simulation Package for R. Therap. Drug Monit. 34: 467-476, 2012.

3. Bustad A, Terziivanov D, Leary R, Port R, Schumitzky A, and Jelliffe R: Parametric and

Nonparametric Population Methods: Their Comparative Performance in Analysing a Clinical Data Set and

Two Monte Carlo Simulation Studies. Clin. Pharmacokinet., 45: 365-383, 2006.

4. Jelliffe R, Bayard D, Milman M, Van Guilder M, and Schumitzky A: Achieving Target Goals most

Precisely using Nonparametric Compartmental Models and "Multiple Model" Design of Dosage Regimens.

Therap. Drug Monit. 22: 346-353, 2000.

5. Jelliffe R, Schumitzky A, Bayard D, Milman M, Van Guilder M, Wang X, Jiang F, Barbaut X, and

Maire P: Model-Based, Goal-Oriented, Individualized Drug Therapy: Linkage of Population Modeling,

New "Multiple Model" Dosage Design, Bayesian Feedback, and Individualized Target Goals. Clin.

Pharmacokinet. 34: 57-77, 1998.

Any comments?

All the best,

Roger Jelliffe - On 12 Jul 2013 at 23:28:13, Vijay Ivaturi sent the message

Back to the Top

Dear All,

Can the estimate of the scaling factor be "negative" to account for the residual between the

observed and predicted thus driving the scaled value (THETA*EPS) for the additive error component

to be negative?

Vijay - On 13 Jul 2013 at 08:37:51, Bob Leary sent the message

Back to the Top

Vijay -

yes - the scaling factor is just estimated as a fixed effect theta. If you want to exclude the

possibility of the additive component of residual error standard deviation being negative, just put

a lower bound of 0 on it. If you want to allow this, leave the bound off.

Note that nominally here the standard deviation is being modeled, not the variance. So for small

concentrations, a negative theta on the additive component will result in an overall negative

residual error standard deviation. On the surface this appears nonsensical, but it is actually OK

because the square of the standard deviation enters the likelihood. In effect, this captures a type

of behavior that is sometimes seen when assay error magnitudes are fit to polynomials in

concentration (see Roger Jelliffe's recent message on this) - at small concentrations, assay error

magnitude can initially decrease with increasing concentration, but at even larger concentrations

this reverses itself and error magnitude starts to increase. Note this type of behavior cannot be

captured in the usual Nonmem style of proportional + additive residual error Ypred*Eps1 + Eps2,

since the residual error variance here is necessarily an increasing function of predicted

concentration Ypred. - On 13 Jul 2013 at 08:38:51, Peiming Ma sent the message

Back to the Top

Dear all,

The scaling factor (=SD(Err2)/SD(Err1)) could be large or small; so the two error models

Yobs = Ypred*(1+ Err1) + Err2 = Ypred + sqrt(Ypred^2 + scalingfactor^2)*Err

Yobs = Ypred*(1 + Err1) + Err1*scalingfactor = Ypred + (Ypred + scalingfactor)*Err

could be quite different. This is the case, for example, when Ypred has similar magnitude as the

scaling factor.

Leonid, in your email when you change both Eps1 (=Err1) and Eps2 to Eps3, I think you made the

errors 100% correlated.

Regards,

Peiming - On 13 Jul 2013 at 10:54:24, Leonid Gibiansky (lgibiansky.aaa.QUANTPHARM.COM) sent the message

Back to the Top

Hi Peiming,

Nice to hear from you. I have not assumed anything; those two error

models were considered separately, and then I showed how they relate to

each other: these are two different parametrization of the same model,

and my earlier e-mail derived the conditions when they are equivalent.

Leonid

--

Leonid Gibiansky, Ph.D.

President, QuantPharm LLC

web: www.quantpharm.com

e-mail: LGibiansky at quantpharm.com - On 14 Jul 2013 at 08:24:23, Peiming Ma sent the message

Back to the Top

Hello, Leonid, Good to hear from you too. Sorry if I missed your point. J

For your 2nd model Y = IPRED*(1 + EPS3) + THETA(1)*EPS3, the SD should be

abs(IPRED+THETA(1))*SIGMA3, not SD = SQRT(IPRED**2*SIGMA3+THETA(1)**2*SIGMA3).

What I tried to say was that the two models are not just different in notations. Their assumed error

SDs could differ by 41% (=sqrt(2)) if the concentrations happen to be around the scaling factor

(with other things being equal, and no good approximation here by expansion of the terms because of

no clear dominant terms in the square root).

Thank you,

Peiming - On 14 Jul 2013 at 11:31:43, Leonid Gibiansky sent the message

Back to the Top

Hi Peiming,

In fact, you are absolutely correct, and my original post was a mistake,

sorry for the confusion. Indeed, the derivation assumed EPS3 to be

uncorrelated with EPS3, that is not the case :(

For the model

Y = IPRED*(1 + EPS1) + EPS2

(EPS1 and EPS2 are independent, with variances SIGMA1, SIGMA2)

SD = SQRT(IPRED**2*SIGMA1+SIGMA2)

and for the model

Y = IPRED*(1 + EPS3) + THETA(1)*EPS3

(with variance SIGMA3).

SD is NOT SD = SQRT(IPRED**2*SIGMA3+THETA(1)**2*SIGMA3)

but indeed given by

abs(IPRED+THETA(1))*sqrt(SIGMA3)

Thanks

Leonid

--

Leonid Gibiansky, Ph.D.

President, QuantPharm LLC

web: www.quantpharm.com

e-mail: LGibiansky at quantpharm.com - On 25 Jul 2013 at 09:03:11, Hans Proost sent the message

Back to the Top

Dear Roger,

Thank you for your explanation of the USC LAPK view on the error model. I have a few comments and

questions:

"In addition, we used to multiply this polynomial by another parameter, Gamma, which we estimated

separately as a point estimate [1]. We found, though, that gamma, when used this way, multiplied up

the curvature in the assay error polynomial, and was probably not the correct thing to do."

I'm not sure why this would not the correct (or at least, appropriate) thing to do. Actually, the

curvature of the assay error polynomial is not changed, it is only multiplied by a parameter. Since

assay precision is the reciprocal of the variance of the assay measurement, this factor gamma will

imply that the assay precision is multiplied by 1/gamma^2, so a constant factor independent of

concentration. Why would this not be the correct thing to do?

"Instead, we now estimate another additive term which we call lambda, which we then add to the assay

error polynomial to provide an estimate of the overall noise with which the clinical part of the

study was done. In this way we first get the assay error polynomial totally apart from the study

itself [1], and then estimate the noise of the clinical part as lambda. Good or poor assay, that is

self-evident. Small lambda, a good clinical study, with small environmental noise due to errors is

the preparation and administration of doses, the errors of recording the timing of doses and the

samples, the model misspecification, and any possible changes in model parameter values during the

period of data analysis. Large lambda, maybe one ought to see what can be done to make the clinical

part of the study more precise. In this way we then have a means to estimate both components of this

overall error. We think this is a more informative approach rather than simply estimating the

overall error and never knowing how much is due to what."

The basic idea is excellent. However, there might be good arguments to assume that this type of

errors is also dependent on the concentration. For example, a dosing error of +10% will increase all

concentrations by 10%, and so the contribution of this error in de SD of peak concentrations is

higher than in the SD of trough concentrations or during elimination phase. For timing errors the

same may apply. This implies that one needs again an additional and a proportional term for this

error model (as was the origin of this thread).

best regards,

Hans Proost

Johannes H. Proost

Dept. of Pharmacokinetics, Toxicology and Targeting

University Centre for Pharmacy

Antonius Deusinglaan 1

9713 AV Groningen, The Netherlands - On 25 Jul 2013 at 09:04:30, Hans Proost sent the message

Back to the Top

Dear Bob,

To Vijay's question

"Can the estimate of the scaling factor be "negative" to account for the residual between the

observed and predicted thus driving the scaled value (THETA*EPS) for the additive error component

to be negative?"

you wrote:

"yes - the scaling factor is just estimated as a fixed effect theta. If you want to exclude the

possibility of the additive component of residual error standard deviation being negative, just put

a lower bound of 0 on it. If you want to allow this, leave the bound off.

Note that nominally here the standard deviation is being modeled, not the variance. So for small

concentrations, a negative theta on the additive component will result in an overall negative

residual error standard deviation. On the surface this appears nonsensical, but it is actually OK

because the square of the standard deviation enters the likelihood."

A negative SD sounds weird, and in my opinion it is not also weird, but also invalid. Using the

square in the calculation of the likelihood solves the numerical problem (indeed, no runtime error,

it works nicely), but not the fundamental point that an SD cannot be negative by definition because

it is derived from the variance. In my view, a lower bound of 0 is mandatory.

"In effect, this captures a type of behavior that is sometimes seen when assay error magnitudes are

fit to polynomials in concentration (see Roger Jelliffe's recent message on this) - at small

concentrations, assay error magnitude can initially decrease with increasing concentration, but at

even larger concentrations this reverses itself and error magnitude starts to increase."

If one uses a polynomial that results in negative values over a certain range of concentrations,

there is at least one concentration associated with SD = 0, and a range of concentrations around

this with very low SD values, resulting in observations with very high weight. This cannot be the

desired characteristic of the assay error polynomial. One should check (manually, e.g using excel)

that the SD has the desired characteristic over the entire concentration range of the data.

best regards,

Hans Proost

Johannes H. Proost

Dept. of Pharmacokinetics, Toxicology and Targeting

University Centre for Pharmacy

Antonius Deusinglaan 1

9713 AV Groningen, The Netherlands - On 25 Jul 2013 at 11:23:54, Bob Leary sent the message

Back to the Top

Dear Hans,

There are two ways to write the 'mixratio' type of additive -proportional model, depending on

whether you put

the theta coefficient on the additive term

a) sddev(residual) = (Pred+theta)*stddev(eps)

or on the multiplicative term

b) stddev(residual) = (1 + theta*Pred)*stddev(eps)

Normally we expect theta to be positive, and both a) and b) will give the same maximum log

likelihood , but the interpretation of stddev(eps) and theta will

be different. For example, in a), stddev(eps) is interpreted as the standard deviation of the

multiplicative component, and in b) of the additive component.

If theta is negative, then we immediately have the negative stddev(residual) problem at low

concentrations in formulation a), which was I think was the version being discussed in

the original thread . In this case, we simply need to interpret this as stddev(residual) =

abs(Pred+theta)*stddev(eps)

In formulation b) with a negative theta we start with a positive std deviation at low

concentrations and then move to a negative std deviation at high concentrations. Again, we should

just use the absolute value

Interpretation at high concentrations.

Phoenix in fact by default uses case b), so low concentrations are not a problem with negative

theta. Normally in this case, when we encounter a negative theta, it usually has a small magnitude

so

the actual predicted standard deviation is positive over the entire range of concentrations

actually encountered.

But you are correct that in both cases there is a problem in that with a linear mixing between the

additive and proportional terms with a negative theta,

there is a critical concentration where the standard deviation of the residual falls to zero. If

this critical concentration indeed is within the range of concentrations actually observed, then

this is indeed a serious

problem, and indeed one would question the validity of such a mixed model, and possibly restrict

theta to be positive (or even take a look at entirely different

types of residual error models). But if the critical concentration is well outside this range,

then in fact the model captures a type of not impossible behavior

where the std deviation of the residual error falls with increasing concentration. This type of

behavior is not possible to capture with the usual NM Ypred*eps1 + eps2

type of residual model.

My allusion to the LAPK work was probably somewhat misplaced and inappropriate. There the error

polynomials are usually at least quadratic in concentration and the predicted standard deviation

Is indeed positive (not possible with a std deviation that is only first order in the concentration

if the mixing coefficient theta is negative). What I was really trying to get at was the LAPK work

shows that a decreasing residual error magnitude with increasing concentration at low concentrations

is indeed possible . On the other hand, my point that the error model in the negative theta

linear case first falls but eventually starts growing in magnitude was poorly judged and not really

relevant, since it ignores the fact that then the zero std deviation problem at the critical

concentration then comes into play if we indeed enter the concentration region where the magnitude

starts growing again.

In any event, a negative theta result should not be categorically rejected, but also should be

examined for reasonableness. - On 26 Jul 2013 at 09:08:17, Hans Proost sent the message

Back to the Top

Dear Bob,

Thank you for your extensive reply. Your equations for cases a) and b) are

quite helpful for the discussion (moving away from the NONMEM- and Phoenix

dialect). However, I'm still worried about the negative theta. Using

absolute values is a nice trick to avoid runtime errors, but it also

obscures a fundamental mistake in the interpretation. A few comments:

1) In your case b), a negative theta indeed would imply an sd decreasing

with concentration.

"the model captures a type of not impossible behavior where the std

deviation of the residual error falls with increasing concentration"

I don't see any situation where this can be true, except for the situation

described in several papers from LAPK, with a second- or third-order

polynomial where the sd at C=0 is larger than at low concentration, but (of

course) increasing with concentration at higer levels; no problem here for a

negative coefficient in the polynomial, if one is sure that the sd is always

positive. But in the case of additional and proportional error (equivalend

to first-order polynomial) this does not make much sense, and a lower bound

of 0 seems the natural way to avoid this.

2) Even if the sd is positive over the entire range of concentrations

actually observed, this does not guarantee that there is no problem. The

predicted concentration (Pred) may be outside this range, and it probably

will do so because of the bonus in the objective function as a result of

small sd's within a certain predicted concentration range. This may cause

problems that are noticed only on careful visual inspection. So, why making

this troubles where the true solution is obvious, i.e. avoiding unacceptable

values for the parameter theta?

3) "the model captures a type of not impossible behavior where the std

deviation of the residual error falls with increasing concentration. This

type of behavior is not possible to capture with the usual NM Ypred*eps1 +

eps2"

Indeed. But you may interpret this statement in different ways. You suggest:

'so, case a) or case b) is the better way to model additional and

proportional error'. Indeed it is a more flexible way to do so. However,

this flexibility is towards the parameter space where I don't want to end

up.

best regards,

Hans

Johannes H. Proost

Dept. of Pharmacokinetics, Toxicology and Targeting

University Centre for Pharmacy

Antonius Deusinglaan 1

9713 AV Groningen, The Netherlands

Email: j.h.proost.-at-.rug.nl - On 29 Jul 2013 at 20:34:50, Roger Jelliffe sent the message

Back to the Top

Dear Hans:

Thanks for your comments on Gamma and lambda. You are correct, that with gamma the curvature of

the assay polynomial is not changed. However, it is multiplied up by gamma, the parameter reflecting

the other components of the noise in the system. Because of this, the shape of the overall noise

polynomial does get altered when it is multiplied up by Gamma. We think this is not correct any more

to do this.

The assay error polynomial gives you no information about the noise in the clinical part of the

study (errors in dose preparation and administration, errors in recording the times of doses and

blood samples, model mis-specification, and any changing model parameter values during the period of

the data analysis). In the same way, the clinical part of the noise gives you no information about

the assay error. They are quite independent of each other. For this reason it does not seem correct

to have the gamma affect the assay error polynomial in any way. Lambda appears better to us. Keep

them separate and simply add them.

However, none of this is really correct. The measurement error (the assay error polynomial) is

truly measurement noise. But the other clinical sources of noise mentioned above truly are not. They

are actually process noise, and should be described as noise in the differential equations that

describe the behavior of the drug in the patient. Measurement noise and process noise are quite

separate. No one except a few people have correctly described process noise. That requires

stochastic differential equations, and none of us do that yet. One group has, but I cannot recall

the reference.

You can have the most accurate radar in the world, and still you cannot predict the course of an

airplane as it gets blown about by the winds at various altitudes. That is because you can have

extremely accurate and precise radar measurements of the plane, but you know nothing about the winds

aloft and how they will affect the course of the plane.

This is why we currently like lambda. We are also working on describing the process noise

properly. When this is done, we will describe things not as 2 sources of measurement noise, but as

separate measurement and process noise. Then we will be doing better.

All the best,

Roger - On 30 Jul 2013 at 14:01:09, Hans Proost sent the message

Back to the Top

Dear Roger,

Thank you for your extensive reply. I fully agree with your statement that the assay error and the

noise/errors in the clinical part are quite independent, and that the assay error can be handled

best using the error function determined from the analytical procedure. Since the 'ultimate

solution' for the errors in the clinical part is not yet available, we must rely on more simplistic

approaches.

However, I still have a problem in using an error term for the clinical part that is independent of

the concentration. As I argued in my message, it is likely that, at least a major component of that

error, is proportional (or something close to that) to the concentration. This is not taken into

account in your lambda approach, if I understand you correctly. This may result in inappropriate

weighting of the data, and consequently, biased parameter estimates. So in my opinion, the lambda

approach is questionable, at least in data sets with a wide range of observed concentrations.

More important, we seem to agree that this is an important issue, and that more work should be done.

best regards,

Hans

Johannes H. Proost

Dept. of Pharmacokinetics, Toxicology and Targeting

University Centre for Pharmacy

Antonius Deusinglaan 1

9713 AV Groningen, The Netherlands

Email: j.h.proost.at.rug.nl

Want to post a follow-up message on this topic?
If this link does not work with your browser send a follow-up message to PharmPK@lists.ucdenver.edu with "Phoenix combination error model" as the subject |

Copyright 1995-2014 David W. A. Bourne (david@boomer.org)