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
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
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/
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
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
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
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
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.
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
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
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/
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
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
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
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
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.
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
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
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
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.
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
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
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.
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
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
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
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
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
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
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.
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
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
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)