From: "Pravin Jadhav" pravinj@gmail.com
Subject: [NMusers] posthoc step
Date: Mon, December 6, 2004 6:35 pm
Hello all,
This might be a naive question but I am puzzled by the requirement to
include POSTHOC command if individual parameters need to be obtained.
If I understand correctly, information on individual estimates is
available at the convergence in some form (OFV equation includes the
summation of difference between individual and population estimate
over all subjects). So what is the necessity of posthoc step? Could
some explain how does this empirical Bayesian step exactly work? or
direct me to some references where this principle is explained.
Also, do other equivalent softwares, specifically NLME implemented in
S-Plus, work on similar framework (meaning, population parameter
estimation followed by Bayesian step)? I apologize for question
concerning other software but I thought it might be relevant.
Any comments are appreciated.
Pravin
--
Pravin Jadhav
Graduate Student
Department of Pharmaceutics
MCV/Virginia Commonwealth University
DPE1/CDER/OCPB/Food and Drug Administration
Phone: (301) 594-5652
Fax: (301) 480-3212
_______________________________________________________
From: "Nitin Kaila" kail0004@umn.edu
Subject: Re: [NMusers] posthoc step
Date: Mon, December 6, 2004 9:05 pm
Hi Pravin: I am not sure if I understand your question correctly.
In the simplest analytical approximation i.e. the first order (FO) method
in NONMEM, all the etas are set to zero. Now what the POSTHOC option does
is that it use the estimated population parameters as priors for an
individual's data and computes the empirical bayesian post hoc
estimates. This happens at the end of a normal FO run with the POSTHOC option.
I hope this helps.
Regards,
Nitin
Nitin Kaila
Doctoral Candidate
308 Harvard St. SE
7-192 WDH
Minneapolis, MN 55455
Department of Experimental and Clinical Pharmacology
University of Minnesota
PH: 612-624-9683 (office)
_______________________________________________________
From: "Pravin Jadhav" pravinj@gmail.com
Subject: Re: [NMusers] posthoc step
Date: Tue, December 7, 2004 2:38 pm
Thank you all for replying to my query. I appreciate it.
I was not clear on these aspect. I was trying to recalculate the OFV
value using the NONMEM output, I got confused and was not sure about
the process.
Diane, I tried your suggestion. .....until one runs POSTHOC the IPRED
is the same as PRED for FO because eta isn't evaluated. (try it and
see)...I guess it is much more clear now.
However, I am still interested in specifics of the POSTHOC step. As
explained many times, it uses population estimates as priors to
calculate individual estimates. My question is- how? and where do I
get the details of so called empirical Bayesian step?. Is it using
MCMC algorithm like what WINBUGS uses?
Thanks.
Pravin
_______________________________________________________
From: "Nick Holford" n.holford@auckland.ac.nz
Subject: Re: [NMusers] posthoc step
Date: Tue, December 7, 2004 3:02 pm
Pravin,
As I understand the mysteries of NONMEM, it computes Maximum A Posteriori Bayesian
(MAPB) estimates from a prior parameter distribution using the population estimates
of THETA (the mean) and OMEGA (the uncertainty). This is done after each iteration
with FOCE and with the final estimates with FO and the POSTHOC option.
MAPB estimates can be obtained with an objective function like this:
SS = sum i=1 to Nobs [(Yhati - Yobsi)**2/SIGMA**2] + sum k=1 to NPAR [(THETAhatk -
THETAki)/OMEGAhatk**2)]
where Yhat is the prediction of the observation, Yobs, using individual estimates,
THETAki. THETAhatk is the population estimate of the kth THETA parameter. OMEGAhatk
is the estimate of OMEGA for the kth THETA parameter. I am sure this is a
oversimplification of what NONMEM actually does but this should give you an idea.
Sheiner LB, Beal SL. Bayesian individualization of pharmacokinetics: simple
implementation and comparison with non-Bayesian methods. J Pharm Sci
1982;71(12):1344-8
Nick
--
Nick Holford, Dept Pharmacology & Clinical Pharmacology
University of Auckland, 85 Park Rd, Private Bag 92019, Auckland, New Zealand
email:n.holford@auckland.ac.nz tel:+64(9)373-7599x86730 fax:373-7556
http://www.health.auckland.ac.nz/pharmacology/staff/nholford/
_______________________________________________________
From: "Bachman, William (MYD)" bachmanw@iconus.com
Subject: RE: [NMusers] posthoc step
Date: Tue, December 7, 2004 3:25 pm
On a much more simplistic level, one can consider the posthoc step in the
following manner:
1. the objective function value is a function of individual predictions,
observations, fixed and random effect parameters and distributions of the
parameters (which we minimize by "optimizing" the parameters).
3. we have the prior parameters and distributions of these parameters (and
the objective function value) from the first order analysis.
4. Bayes theorem allows us to "back calculate", the individual parameters
from the known individual observations, population parameters and
distributions (and objective function value). In other words, its like a
system of known and unknowns where we now rearrange the objective function
equation and solve for the individual parameters and predictions.
No apologies to the statisticians!
:)
Bill
_______________________________________________________
From: "Wang, Yaning" WangYA@cder.fda.gov
Subject: RE: [NMusers] posthoc step
Date: Tue, December 7, 2004 3:50 pm
Emperical Bayes Estimate is very typical in mixed effect modeling. Maybe you
should start with linear mixed effect modeling to get familiar with the
concept first before dealing with nonlinear mixed effect modeling where
various estimation methods may confuse you even further. Reading some basic
statistic books about mix effect modeling (Marie Davidian and David
Giltinan, or Geert Molenberghs and Geert Verbeke) and
Bayesian Analysis (Andrew Gelman, John B. Carlin, Hal S. Stern, and Donald
B. Rubin) will help.
Yaning Wang, PhD
Pharmacometrician
OCPB, FDA
_______________________________________________________
From: "Kowalski, Ken" Ken.Kowalski@pfizer.com
Subject: RE: [NMusers] posthoc step
Date: Tue, December 7, 2004 5:17 pm
Pravin and Nick,
I think Nick's objective function below is essentially correct when using an
additive error model and a diagonal Omega with the exception that (THETAhatk
- THETAki) in the second term of his expression for SS needs to be squared
and the i index he uses for the second term is for subjects while the i in
the first term indexes observations. If we assume additive errors for the
ETAs for ease of exposition such that THETAki = THETAhatk + ETAki and use j
as an index for observations then we can re-write the MAPB objective
function as:
SSi = sum j=1 to Nobsi [(Yobsij - Yhatij)^2/SIGMA^2] + sum k=1 to Npar
[ETAki^2/OMEGAhatk^2]
and we minimize SSi for subject i with respect to the ETAki (i.e., find the
values of ETAki, k=1,...,Npar that minimize SSi). I use an expression like
this to help illustrate the shrinkage estimation properties of MAPB
estimates. That is, the second term of the objective function will dominate
when Nobsi is small (sparse data) and if you take it to the extreme when
Nobsi=0 then the second term and hence SSi is minimized when ETAki=0 for all
k=1,...,Npar. That is, in absence of any data for a new individual the best
estimate is the typical individual prediction we obtain from the population
parameter estimates (THETAhatk) where ETAki=0 for all k.
Ken
_______________________________________________________
From: "Gastonguay, Marc" marcg@metrumrg.com
Subject: Re: [NMusers] posthoc step
Date: Tue, December 7, 2004 7:17 pm
Ken,
I'd like to suggest an additional point to consider regarding your
explanation about shrinkage estimation properties of MAP Bayes
estimates. As you point out, the number of data points does affect the
balance between the influences of the individual data and population
priors, but this is not the only consideration. The relative magnitude
of OMEGA vs SIGMA is also an important factor in determining if the MAP
Bayes estimate is dominated by the individual's data or the population
priors. If inter-individual variability (OMEGA) is very large relative
to measurement noise (SIGMA), a situation that is not uncommon in
population PKPD, the second term in the MAP Bayes OF is essentially
minimized regardless of the magnitude of ETAki. Under such
circumstances, it is possible to have a data-dominated individual MAP
Bayes estimate, even when sampling is sparse (but greater than zero).
We've all observed the situation where POSTHOC IPREDs from a base model
fit the individual data very well, even with sparse sampling. The
opposite is also theoretically possible (SIMGA>>OMEGA), and could lead
to prior-dominated MAP Bayes estimates even when sampling is extensive,
but this is rarely the case in population PKPD.
Marc
Marc R. Gastonguay
Metrum Research Group LLC
www.metrumrg.com
_______________________________________________________
From: jerry.nedelman@pharma.novartis.com
Subject: RE: [NMusers] posthoc step
Date: Tue, December 7, 2004 9:16 pm
Friends:
I did a little experiment to check the MAP calculations of NONMEM. Either
there's something wrong with my code, or something fishy is going on. I
hope somebody can help solve the mystery.
First I tried a simple nonlinear model:
y = 10*exp(-kt) + err, k~N(10,var=4), err~N(0,var=0.04)
I simulated 3 data points at t=1,2,3 with k=1. (This may not look like a very
good value of k considering the prior, but my experiment was just to check
arithmetic, and maybe the choice was sufficiently stressful to the methodology
to have found something.)
I wrote a SAS program to fit the model to the data in 2 ways:
1) Using ordinary least squares.
2) Using weighted least squares with a pseudo-observation for the parameter, which
is MAP as described in LB Sheiner and SL Beal, Bayesian individualization of
pharmacokinetics: Simple implementation and comparison with non-Bayesian
methods, J Pharm Sci, 71:1344-1348 (1982).
SAS/OLS gives khat = 0.9404, and SAS/MAP gives khat=0.9437.
Here is the SAS code:
data one;
******************************************;
* dat: an indicator for whether it's "real data" or an extra obs for the parameter.
* t: the t of y = exp(-kt)
* obs: observation (real data or theta)
* "Real data" generated by S-Plus: 10*exp(-c(1,2,3)) + 0.2*rnorm(3,0,1)
* wobs: weighted observation, obs/sqrt(var), used for weighted-least-squares.
* var = 0.04 for the real data, var = 4 for the prior.
******************************************;
input dat t obs wobs;
cards;
1 1 3.87 19.35
1 2 1.66 8.3
1 3 0.44 2.2
0 999 10 5
;
run;
**********************************;
title 'Ordinary least squares';
**********************************;
proc nlin data=one;
where dat=1;
model obs = 10*exp(-k*t);
parm k=0.1;
run;
**********************************;
title 'Weighted least squares, MAP';
**********************************;
proc nlin data=one;
model wobs = dat*( 10*exp(-k*t)/0.2 ) +
(1-dat)*( k/2 );
parm k=0.1;
output out=out1 pred=pred;
run;
Then I tried using NONMEM to find the MAP estimate by the posthoc step. The NONMEM/MAP
estimate of k was khat=9.78, pulling toward the prior much more than SAS/MAP.
Here are the control and data files:
$PROB BAYES TEST
$DATA OLSID.CSV IGNORE=#
$INPUT DAT ID T OBS=DV WT WOBS
$PRED
K = THETA(1) + ETA(1)
F = 10*EXP(-K*T)
IPRED = F
Y = F + ERR(1)
$ESTIMATE MAXEVALS=0 POSTHOC
$THETA 10
$SIGMA 0.04
$OMEGA 4
$TABLE K IPRED
Trying to understand why there would be such a discrepancy between the two approaches,
which I had expected to be identical, I got to wondering if it has to do with the
nonlinearity of the model. So I tried a linear model. The data looked like it might
be fairly well fit by y = 4 - kt with k=1. So I redid the above SAS and NONMEM steps
with this model, same data as before. The SAS/MAP and NONMEM/MAP estimates were
identical, khat = 1.1128.
Here are the relevant files:
data one;
******************************************;
* dat: an indicator for whether it's "real data" or an extra obs for the parameter.
* t: the t of y = exp(-kt)
* obs: observation (real data or theta)
* "Real data" generated by S-Plus: 10*exp(-c(1,2,3)) + 0.2*rnorm(3,0,1)
* wobs: weighted observation, obs/sqrt(var), used for weighted-least-squares.
* var = 0.04 for the real data, var = 4 for the prior.
******************************************;
input dat t obs wobs;
cards;
1 1 3.87 19.35
1 2 1.66 8.3
1 3 0.44 2.2
0 999 10 5
;
run;
**********************************;
title 'Weighted least squares, Bayesian, Linear Model';
**********************************;
proc nlin data=one;
model wobs = dat*( (4-k*t)/0.2 ) +
(1-dat)*( k/2 );
parm k=0.1;
output out=out1 pred=pred;
run;
$PROB BAYES TEST
$DATA OLSID.CSV IGNORE=#
$INPUT DAT ID T OBS=DV WT WOBS
$PRED
K = THETA(1) + ETA(1)
F = 4 - K*T
IPRED = F
Y = F + ERR(1)
$ESTIMATE MAXEVALS=0 POSTHOC
$THETA 10
$SIGMA 0.04
$OMEGA 4
$TABLE K IPRED FILE=LINEAR.FIT
Any ideas what might be wrong with the nonlinear case?
Thanks,
Jerry
Jerry R. Nedelman, Ph.D.
Director, Biostatistics and Statistical Programming
Novartis Pharmaceuticals
One Health Plaza
East Hanover, NJ 07936
jerry.nedelman@pharma.novartis.com
862-778-6730 phone
973-781-6498 fax
_______________________________________________________
From: "Pravin Jadhav" pravinj@gmail.com
Subject: Re: [NMusers] posthoc step
Date: Wed, December 8, 2004 12:49 am
Marc, Ken, Yaning, Bill, Nick and group:
Thank you very much. I appreciate your time. I learned something.
However, as Yaning suggested, it is time to go back to basics. The
only difference is now I kind of know (or may be not) where to start.
Apology if it was a really stupid question.
Pravin
_______________________________________________________
From: "Leonid Gibiansky" leonidg@metrumrg.com
Subject: RE: [NMusers] posthoc step
Date: Wed, December 8, 2004 8:27 am
Jerry,
Most likely, this is a problem with too extreme data and implementation
of POSTHOC in NONMEM.
I checked how solution depends on OMEGA. With OMEGA= 8.39 or higher, NONMEM
gives k=0.942, pretty close to what you would expect. But for OMEGA=8.38,
solution jumps to K=4.19; for OMEGA= 8.37, K= 9.04 and then goes up to 9.78
when you move OMEGA down to 4 (e.g, OMEGA=8, K= 9.27)
I also checked the probability of K=1 with mean=10 and various variances:
> pnorm(1, mean = 10, sd = sqrt(4))
[1] 3.397673e-006
> pnorm(1, mean = 10, sd = sqrt(8.4))
[1] 0.0009504466
Probability of the realization that you tried is about 3 / million. NONMEM
starts to give correct results at about p = 0.001. It this can be
generalized, then NONMEM POSTHOC fails at outliers with p < 0.001
Best regards,
Leonid
_______________________________________________________
From: "Kowalski, Ken" Ken.Kowalski@pfizer.com
Subject: RE: [NMusers] posthoc step
Date: Wed, December 8, 2004 9:05 am
Marc, NMusers,
I agree with your comments. I was not suggesting that the number of
observations is the only factor that determines the amount of shrinkage.
When OMEGA>>SIGMA, MAP Bayes will have less shrinkage even with sparse data
but there will still be some shrinkage. The degree of shrinkage will also
be influenced by the design (times at which samples are taken). For example
with sparse PK during a steady state dosing interval following oral
administration, there is more shrinkage for the ETAs on ka and V/F than
there is for CL/F since we often don't have rich information during the
absorption phase and V/F is better estimated during non-steady-state
conditions. So, with sparse data, depending on the placement of the time
points, there can still be a fair amount of shrinkage in one or more of the
ETAki.
I don't think you were implying this but just to be clear, the IPREDs can
fit the observations very well under sparse conditions even when there is
significant shrinkage in one or more of the individual parameter estimates.
With sparse data (e.g., fewer observations than number of parameters) if we
fit the individual model using WLS we will have over-parameterization as the
estimates will not be unique in that we can have an infinite set of
solutions that will result in essentially the same minimum value of the WLS
objective function. MAP Bayes estimates on the other hand, although they
shrink the estimates to mean (perhaps some parameters more than others) will
be unique.
A simple exercise one can do to assess the degree of shrinkage is to
calculate the sample variance for ETAki and compare that to the estimate of
OMEGAk from the population model fit. The sample variance should be smaller
and the greater the discrepancy the more shrinkage there is.
Ken
_______________________________________________________
From: "Nick Holford" n.holford@auckland.ac.nz
Subject: Re: [NMusers] posthoc step
Date: Wed, December 8, 2004 2:34 pm
Leonid,
I have also done something similar and found similar results.
Without relying on any explicit calculation it seems that if the prior for K is 10
with a SD of 2 and data of 3 observations simulated with K=1 that the posterior
estimate of K might be much closer to 10 than to 1. A value of 1 drawn from
N(10,SD=2) is quite unlikely (NORMDIST(1,10,2,TRUE) is 3E-6). The POSTHOC estimate
of 9.78 which is obtained using the NM-TRAN code and data Jerry supplied supports
this. The SAS MAP estimate of 0.9437 that Jerry reported seems unreasonable given
such a strong prior of SD=2 for K=10.
I have used NONMEM with both the pseudo-observation (DATA) method which Jerry
mentioned and the undocumented PRIOR method in NONMEM V (METHOD=ZERO and
METHOD=COND). The estimate of Khat for all methods remains stubbornly at 10 with SD
of 2. When the SD for the prior on K was increased to 9 then Khat changes abruptly
from 10 to 0.94. I was rather surprised not to find a gradual change in Khat as the
SD for the prior on K was increased. As Leonid shows below the transition from Khat
of ~10 to Khat ~ 1 happens over a very narrow range (between SD=8 and SD=9).
In addition to this sharp transition I also found 'spikes' of Khat dropping to ~1 at
certain values for the SD of K. These spikes were not present when I used NONMEM VI
with the DATA method of Bayesian estimation and the transition SD was at a lower
value (7.35). (see attached PDFs).
Bayesian Estimation with NONMEM V.pdfBayesian Estimation with NONMEM VI.pdf
Nick
--
Nick Holford, Dept Pharmacology & Clinical Pharmacology
University of Auckland, 85 Park Rd, Private Bag 92019, Auckland, New Zealand
email:n.holford@auckland.ac.nz tel:+64(9)373-7599x86730 fax:373-7556
http://www.health.auckland.ac.nz/pharmacology/staff/nholford/
_______________________________________________________
From: "Steve Duffull"
Subject: RE: [NMusers] posthoc step
Date: Wed, December 8, 2004 3:19 pm
Hi all
Nick wrote:
> When the SD for the prior on K was increased to 9 then Khat changes abruptly from
10 to 0.94.
I am not going to attempt to comment meaningfully on the findings to date - but want
to ask another question about the findings.
If the prior for K is N(10,x) and the data were computed based on K_i=1 then how
could NONMEM (or SAS...) have computed a Khat of < 1? I would have guessed that the
Bayesian solution would have wanted the value of Khat: 1 < Khat < 10. That Khat
was 0.94 would indicate that the prior would have to be noninformative, and
therefore a value of Khat < 1 is simply due to some estimation bias in the system. I
would have guessed that a N(10,x), where x lies between 4 and 10 is not
noninformative. So - why is this happening?
It is of course possible that the solution is obvious and I have missed part of the
previous thread which explains why this is so - in which case just ignore this
email.
Steve
========================================Stephen Duffull
School of Pharmacy
University of Queensland
Brisbane 4072
Australia
Tel +61 7 3365 8808
Fax +61 7 3365 1688
University Provider Number: 00025B
Email: sduffull@pharmacy.uq.edu.au
www: http://www.uq.edu.au/pharmacy/sduffull/duffull.htm
PFIM: http://www.uq.edu.au/pharmacy/sduffull/pfim.htm
MCMC PK example: http://www.uq.edu.au/pharmacy/sduffull/MCMC_eg.htm
========================================
_______________________________________________________
From: "Steve Duffull" sduffull@pharmacy.uq.edu.au
Subject: RE: [NMusers] posthoc step
Date: Wed, December 8, 2004 8:10 pm
Hi
To continue the story a little on Jerry's example.
I ran the model in WinBUGS - see code below.
model {
for (j in 1:3) {
data[j] ~ dnorm(model[j], tau)
model[j] http://www.uq.edu.au/pharmacy/sduffull/duffull.htm
PFIM: http://www.uq.edu.au/pharmacy/sduffull/pfim.htm
MCMC PK example: http://www.uq.edu.au/pharmacy/sduffull/MCMC_eg.htm
=========================================
_______________________________________________________
From: "Nick Holford" n.holford@auckland.ac.nz
Subject: Re: [NMusers] posthoc step
Date: Wed, December 8, 2004 10:32 pm
Steve,
I don't understand why you say:
> When the prior variance was set to 4 Khat was estimated at 0.95
> These are similar to the NONMEM and SAS results when variance was 4.
I wrote earlier:
> The estimate of Khat for all methods remains stubbornly at 10 with SD of 2.
i.e. the prior variance on K was 4.
and Leonid wrote:
> [K] goes up to 9.78 when you move OMEGA down to 4
I'ts not clear if Leonid means OMEGA is the variance of the prior or the SD of the
prior. But whichever convention one assumes it is evident that we both get estimates
of Khat very close to 10 and not 0.95.
It looks like WinBUGS and SAS seem to think Khat is dominated by the data (simulated
with K=1) whereas NONMEM seems to prefer the prior of 10 (SD=2). It appears to me
that NONMEM has a better sense of the prior on K because as Leonid pointed out there
is only 3 on one million chance that K=1.
Nick
--
Nick Holford, Dept Pharmacology & Clinical Pharmacology
University of Auckland, 85 Park Rd, Private Bag 92019, Auckland, New Zealand
email:n.holford@auckland.ac.nz tel:+64(9)373-7599x86730 fax:373-7556
http://www.health.auckland.ac.nz/pharmacology/staff/nholford/
_______________________________________________________
From: jerry.nedelman@pharma.novartis.com
Subject: RE: [NMusers] posthoc step
Date: Wed, December 8, 2004 11:08 pm
Friends: Thanks to all for the interesting insights. I guess the bottom
line is that something about how NONMEM handles the nonlinearity breaks
down when the data are unlikely relative to the prior. For real examples,
this might mean that posthoc estimates for outlying subjects are shrunk
substantially more than they should be, as Leonid pointed out. There seem
to be some issues with the PRIOR method, too, based on Nick's results.
Steve shows that WinBUGS manages things OK.
For those who thought the MAP estimate should have been close to the
prior mean 10 because the prior mass was concentrated away from the
sparse data, consider the linear case, where SAS and NONMEM agreed.
The same prior and data are used there. And there, too the MAP estimate
(1.1128) is close to the OLS estimate (1.1064) and far from the prior
mean. One can actually find the MAP estimate for the linear case
analytically. Let
omega = prior variance (= 4 in the example)
theta = prior mean (= 10 in the example)
k_ols = OLS estimate of k (= 1.1064 in the example)
v_ols = sampling variance of k_ols, i.e, the square of its standard
error (= 0.04/(1^2 + 2^2 + 3^2) = 1/350 in the example)
k_map = MAP estimate of k (= 1.1128 in the example)
Then
k_map = k_ols - v_ols*(k_ols - theta)/(v_ols + omega).
The amount of shrinkage is determined by
v_ols/(v_ols + omega) = (1/350) / ( (1/350) + 4 ) = 1/1401.
Thus, even though the data were generated by a very unlikely value of k relative
to the prior, the precision of the least squares estimate is so great, relative
to the prior, that it dominates in the blending of the prior and the data to
yield the posterior. The same thing should happen in the nonlinear case, and
does happen with SAS and WinBUGS, but something breaks down with NONMEM's
way of handling it.
Jerry
_______________________________________________________
From: "Wang, Yaning"
Subject: RE: [NMusers] posthoc step
Date: Thu, December 9, 2004 6:26 pm
Dear all:
It has been a very interesting topic. This discussion may lead to
something quite significant.
Jerry:
I think it is too early to say NONMEM is worse than SAS. When you
compared SAS/MAP estimate (khat=0.944) with NONMEM MAP estimate
(khat=9.78, independent of estimation methods), you were only
checking whether NONMEM MAP estimate was following the method of
weighted least squares with a pseudo-observation for the parameter.
This was actually tested in the linear case in your example very well.
NONMEM does seem to use this method to estimate khat. But it failed for
the nonliear case. Does this mean SAS is better in MAP estimation in
nonlinear mixed modeling? Maybe not (see the following SAS PROC NLMIXED
output) . If I applied your SAS code for OLS and WLS/MAP fitting in
NONMEM, I could get exact the same results (khat_OLS=0.94, khat_WLS=0.944).
It seems to me that the unreasonable NONMEM MAP estimate (9.78) may be due
to the linearization or integral approximation used in NONMEM. I was really
reluctent to think this way because those approximation methods are used to
estimate the 'real' parameters and those emperical Bayesian estimates of
random effects should not be so complicated. But when I wrote down the
linearized model for your nonlinear example and applied the same SAS (or
NONMEM) WLS/MAP code to this linear model, I got khat=9.82. Given the
following WinBUGS results (Nick, WinBugs results don't match NONMEM
results at all),
omega2 0.04 0.09 0.094 0.095 0.1 0.5 4
khat 9.992 9.984 9.984 1.165 1.146 0.9708 0.9472
there is clearly something wrong.
Then I tried SAS PROC NLMIXED to see whether SAS can do a better job. Surprisingly,
or as expected, if FIRO (first-order) method is used, khat=9.82, which is
identical to the result above based on the linearized model. If GAUSS, HARDY or
ISAMP method is used, khat=9.78, which is identical to original NONMEM MAP
result. I also used your linear case (4-kt) to make sure the SAS code is working
as I expected. Under linear model, SAS PROC NLMIXED is also doing the same thing
as weighted least squares with a pseudo-observation for the parameter. It seems
that NONMEM is implementing FO in a different way from SAS, at least on MAP
estimates, and achieves similar MAP estimates as those more computation intensive
methods in SAS.
But none of these nonlinear mixed effect modeling tools (SAS or NONMEM, someone
can try S+ and report the outcome here) is handling nonlinear models appropriately
for MAP estimates under this "sufficiently stressful" situation. Given this
observation, the impact of this surprising outcome may deserve more study.
Original model: y = d*exp(-kt) ( I use d here to replace 10 to avoid confusion
because the prior mean of k is also 10 accidentally)
Linearized model (first-order Taylor expansion around 10): y=d*exp(-10t){1-(k-10)*t}
SAS code:
proc nlin data=one;
model wobs = dat*( 10*exp(-10*t)*(1-(k-10)*t)/0.2 ) +(1-dat)*( k/2 );
parm k=0.1;
output out=out1 pred=pred;
run;
NONMEM code:
$PROB BAYES TEST
$DATA ../data/nlWLS.CSV IGNORE=#
$INPUT DAT ID T OBS WOBS=DV
$PRED
K = THETA(1)
F=DAT*10*EXP(-10*T)*(1-(K-10)*T)/0.2+(1-DAT)*(K/2)
IPRED = F
Y = F + ERR(1)
$ESTIMATE MAXEVALS=9999
$THETA 0.1
$OMEGA 0.04
$TABLE K IPRED FILE=WLS.FIT
SAS NLMIXED PROCEDURE
/*when use other methods, increase QPOINTS to 250*/
data oneb; set one; if dat=1;run;
proc nlmixed data=ONEB cov corr method=FIRO;
parms TVK=10 s2k=4 s2=0.04;
bounds 10 <=TVK <= 10, 0.04<=S2<=0.04;
K = TVK+ETAK;
F= 10*EXP(-K*T);
model OBS ~ normal(F,s2);
random etaK ~ normal([0],[4]) subject=DAT OUT=MAP;
PREDICT K OUT=KMAP;
run;
WinBugs code:
model {
for (j in 1:3) {
data[j] ~ dnorm(model[j], 25)
model[j] Bayesian Estimation with NONMEM V RUV FIXED.pdf
-- data method ---
$PROB BAYES TEST
$DATA data.csv IGNORE=#
$INPUT ID TIME OBS=DV DVID
$ESTIMATE MAXEVALS=9990
METHOD=COND SLOW
$THETA 10 ; K
$OMEGA 0 FIX ; PPVK
$SIGMA 0.04 FIX ; RUV
$SIGMA 4 FIX ; Kprior_uncertainty
$PRED
K = THETA(1) + ETA(1)
C = 10*EXP(-K*TIME)
IF (DVID.EQ.1) THEN ; prior
Y=THETA(1)+ERR(2)
ENDIF
IF (DVID.EQ.2) THEN ; obs
Y = C + ERR(1)
ENDIF
-- prior method ---
$PROB BAYES TEST
$DATA prior.csv IGNORE=#
$INPUT ID TIME OBS=DV
$ESTIMATE MAXEVALS=9990
METHOD=COND SLOW
$THETA 10 ; K
$OMEGA 0 FIX ; PPVK
$SIGMA 0.04 FIX; RUV
;prior THETA
$THETA 10 FIX ; Kprior
$OMEGA 4 FIX ; Kprior_uncertainty
$SUBR PRIOR=prior.for
$PRED
K = THETA(1) + ETA(1)
C = 10*EXP(-K*TIME)
Y = C + ERR(1)
-- maxeval=0 method ---
$PROB BAYES TEST
$DATA prior.csv IGNORE=#
$INPUT ID TIME OBS=DV
$ESTIMATE MAXEVALS=0
METHOD=COND SLOW
$THETA 10 ; K
$OMEGA 4 ; Kprior_uncertainty
$SIGMA 0.04 ; RUV
$PRED
K = THETA(1) + ETA(1)
C = 10*EXP(-K*TIME)
Y = C + ERR(1)
$TABLE ID TIME K
--- data.csv ---
#ID Time Obs DVID
1, 0, 10, 1
1, 1, 3.87, 2
1, 2, 1.66, 2
1, 3, 0.44, 2
--- prior.csv ---
#ID Time Obs
1, 1, 3.87
1, 2, 1.66
1, 3, 0.44
--- prior.for ---
SUBROUTINE PRIOR (ICALL,CNT,NTHP,NETP,NEPP)
DOUBLE PRECISION CNT
IF (ICALL.LE.1) THEN
NTHP=1
NETP=1
NEPP=1
ENDIF
CALL NWPRI(CNT)
RETURN
END
--
Nick Holford, Dept Pharmacology & Clinical Pharmacology
University of Auckland, 85 Park Rd, Private Bag 92019, Auckland, New Zealand
email:n.holford@auckland.ac.nz tel:+64(9)373-7599x86730 fax:373-7556
http://www.health.auckland.ac.nz/pharmacology/staff/nholford/
_______________________________________________________
From: "Ludden, Thomas (MYD)" luddent@iconus.com
Subject: RE: [NMusers] posthoc step
Date: Fri, December 10, 2004 11:48 am
Nick, Jerry, Leonid, Steve, Yaning et al.
The objective function value for this problem is rather ugly (see tabulation
below). It is generally flat around a K value of 10 but there is a minor
maximum at about K=7.2-7.25. Some gradient search algorithms will have
trouble with this. The Solver in Excel (obviously not the Gold Standard)
gives a K estimate of 9.78 with initial estimates of 7.25 up to 10. With
initial estimates of .1 to 7.2, the K estimate is 0.944. The initial
estimates for NONMEM's POSTHOC ETA search are not readily accessed so I have
not tested NONMEM itself but I have been in contact with Stuart Beal
regarding this question. Jerry, would you please try the SAS procedure with
an initial estimate of 10 with and without estimation of the residual
variance? It would be interesting to see if SAS's search procedure can get
past the minor maximum.
initial initial Final Final
K est OFV K est OFV
10 448.0646546 9.781 448.06
9 448.1637275 9.781 448.06
8 448.5035677 9.781 448.06
7.5 448.6452926 9.781 448.06
7.25 448.6697797 9.781 448.06
7.2 448.6687872 0.944 21.603
7.1 448.6595588
7 448.6393969 0.944 21.603
6.5 448.3096192
6 447.3663721 0.944 21.603
5.5 445.3349748
5 441.4403283
4.5 434.4249225
4 422.270891
3.5 401.8019287
3 368.1922873
2.5 314.6254852
2 233.1745355
1.5 121.66397
1 23.59852696 0.944 21.603
0.9 23.0073974
0.8 39.55006733
0.7 83.27305134
0.6 170.024752
0.5 325.154084
0.4 589.7760326
0.3 1031.4946
0.1 2973.924151 0.944 21.603
_______________________________________________________
From: jerry.nedelman@pharma.novartis.com
Subject: RE: [NMusers] posthoc step
Date: Sun, December 12, 2004 11:50 pm
Tom: SAS doesn't get past it either. That was with the Residual SD fixed.
When the Residual SD was estimated (starting always from 0.2), the estimate
of k was 10 for initial guesses from 10 to 3. But the estimated Residual SD
was weird then, essentially plus or minus infinity. When the initial guess
of k was 2 or 1, SAS failed to converge, although it stopped with k near 0.94.
Jerry
Results:
Initial k Final k SSE Res. SD
10 9.78 448.1 Fixed
7.25 9.78 448.1 Fixed
7.2 0.94 21.6 Fixed
Initial k Final k SSE Res. SD Initial Res. SD Final Res. SD
10 10 8.49E-14 Estimated 0.2 14526239
7 10 2.86E-19 Estimated 0.2 -7.92E+09
3 10 3.61E-20 Estimated 0.2 2.23E+10
2 0.9847 21.0218 Estimated 0.2 0.3671 Failure to converge
1 0.9578 21.7932 Estimated 0.2 0.1957 Failure to converge
Code:
data one;
******************************************;
* dat: an indicator for whether it's "real data" or an extra obs for the parameter.
* t: the t of y = exp(-kt)
* obs: observation (real data or theta)
* "Real data" generated by S-Plus: 10*exp(-c(1,2,3)) + 0.2*rnorm(3,0,1)
******************************************;
input dat t obs;
cards;
1 1 3.87
1 2 1.66
1 3 0.44
0 999 10
;
run;
**********************************;
title 'Weighted least squares, MAP, initial=10, fix residual sd';
**********************************;
proc nlin data=one;
y = dat*obs/0.2 + (1-dat)*obs/2;
model y = dat*( 10*exp(-k*t)/0.2 ) +
(1-dat)*( k/2 );
parm k=10;
run;
**********************************;
title 'Weighted least squares, MAP, initial=7.25, fix residual sd';
**********************************;
proc nlin data=one;
y = dat*obs/0.2 + (1-dat)*obs/2;
model y = dat*( 10*exp(-k*t)/0.2 ) +
(1-dat)*( k/2 );
parm k=7.25;
run;
**********************************;
title 'Weighted least squares, MAP, initial=7.2, fix residual sd';
**********************************;
proc nlin data=one;
y = dat*obs/0.2 + (1-dat)*obs/2;
model y = dat*( 10*exp(-k*t)/0.2 ) +
(1-dat)*( k/2 );
parm k=7.2;
run;
**********************************;
title 'Weighted least squares, MAP, initial=10, estimate residual sd';
**********************************;
proc nlin data=one;
y = dat*obs/res_sd + (1-dat)*obs/2;
model y = dat*( 10*exp(-k*t)/res_sd ) +
(1-dat)*( k/2 );
parm k=10 res_sd=0.2;
run;
**********************************;
title 'Weighted least squares, MAP, initial=7, estimate residual sd';
**********************************;
proc nlin data=one;
y = dat*obs/res_sd + (1-dat)*obs/2;
model y = dat*( 10*exp(-k*t)/res_sd ) +
(1-dat)*( k/2 );
parm k=7.25 res_sd=0.2;
run;
**********************************;
title 'Weighted least squares, MAP, initial=3, estimate residual variance';
**********************************;
proc nlin data=one;
y = dat*obs/res_sd + (1-dat)*obs/2;
model y = dat*( 10*exp(-k*t)/res_sd ) +
(1-dat)*( k/2 );
parm k=3 res_sd=0.2;
run;
**********************************;
title 'Weighted least squares, MAP, initial=2, estimate residual sd';
**********************************;
proc nlin data=one;
y = dat*obs/res_sd + (1-dat)*obs/2;
model y = dat*( 10*exp(-k*t)/res_sd ) +
(1-dat)*( k/2 );
parm k=2 res_sd=0.2;
run;
**********************************;
title 'Weighted least squares, MAP, initial=1, estimate residual sd';
**********************************;
proc nlin data=one;
y = dat*obs/res_sd + (1-dat)*obs/2;
model y = dat*( 10*exp(-k*t)/res_sd ) +
(1-dat)*( k/2 );
parm k=1 res_sd=0.2;
run;
_______________________________________________________
From: "Ludden, Thomas (MYD)" luddent@iconus.com
Subject: RE: [NMusers] posthoc step
Date: Mon, December 13, 2004 4:53 pm
Jerry,
Thanks very much for the SAS results. The multiple minima seems to be at
least part of the problem.
With an initial estimate of K=10 (ETA=0) and fixed residual error variance
the predicted values are all well below the prior err SD of 0.2. With K=10
the half-life is 0.0693. The first data point is at 1 time unit, > 10
halflives. The predicted values at T=1,2,3 are 0.000454, 2.06x10 ^ -8, 9x10
^ -13, respectively. It is my understanding that the NONMEM ETA search
begins with ETA=0 (THETA=10) and is in a region of the objective function
(OFV) where changes in the observation part of the OFV are very small.
[Subtracting a very small number, the prediction, from a much larger number,
the observation, results in very little decrease in this part of the OFV as
the K value is changed.] Moving the K estimate away from 10 toward 1 should
decrease the observation part of the OFV and increase the parameter part of
the OFV. In the region immediately around K=10 the objective function has a
minor minimum at K=9.78 but then as K is further decreased the changes are
slightly dominated by the prior and the total OFV increases (slightly) until
K is about 7.2 at which point the predictions are becoming large enough to
cause the OFV to begin to decrease again with the real minimum, 21.6 at
K=0.944. Unfortunately, I have not found a way to vary the initial estimate
for ETA when using MAXEVAL=0, POSTHOC with NONMEM so I cannot determine how
NONMEM handles the two minima under these conditions.
In this extreme case, the individual residuals (IRES) should have large
absolute values and would signal that this individual's data are poorly
described by the individual parameter estimates. This problem may occur only
rarely, since, in my experience, IPRED's and DV's almost always agree very
well.
The multiple minima in the OFV may be a function of the extreme conditions
of the example. As pointed out by Stuart Beal, "... posthoc is meant to be
used with realistic estimates of the population parameters; ones that are
commensurate with data." The T values of 1,2,3 are very high relative to
realistic sampling times given the prior of K~N(10,var=4). The majority of
subjects from the population described by the prior would have very low
concentrations, essentially "nondata" values given the residual error
variance of 0.04. T values of 0.1, 0.2, 0.3 would be more reasonable design
points for sampling from the prior. If I use Y values based on K=1 and
these more realistic T values, the empirical Bayes OFV appears to have a
single minimum at about K=1.011 over the range of K values from .1 to 10
(See tabulation below). Now NONMEM (MAXEVAL=0, POSTHOC) and Excel solver
both yield K=1.01 with an initial value of K=10. However, it is difficult
to know how to generalize this.
Thanks for bring this interesting problem to the attention of nmusers.
Users should examine IPRED VS DV or similar diagnostic plots to look for
substantial deviations from the line of unity that might indicate a problem
of this kind. I simulated a sample of 50 individuals from the prior and
then included the extreme individual in the data set. Plots of IPRED vs DV
revealed marked deviation from the line of unity for this individual.
Tom
initial est OFV
10 3082.1932
9.8 3032.9791
9.6 2982.4598
9.4 2930.6026
9.2 2877.3750
9 2822.7451
8 2527.5070
7.5 2365.3835
7.25 2280.5581
7.2 2263.2889
7.1 2228.4456
7 2193.1954
6.5 2010.8571
6.45 1992.0692
6 1818.5523
5.5 1616.8539
5 1406.8869
4.5 1190.5497
4 970.8134
3.5 752.1250
3 540.9518
2.5 346.5135
2 181.7666
1.5 64.7274
1.011 20.2260
1 20.2500
0.9 22.8657
0.8 30.0327
0.7 42.0970
0.6 59.4266
0.5 82.4128
0.4 111.4720957
0.3 147.0470862
0.1 239.6568697
_______________________________________________________
From: "Nick Holford" n.holford@auckland.ac.nz
Subject: Re: [NMusers] posthoc step
Date: Tue, December 14, 2004 3:56 pm
Hi,
Would it be fair to summarize the performance of the various methods (with fixed
residual SD) as follows?
1. NONMEM using either the DATA or the PRIOR methods and WinBugs (Steve, is this
correct?) get past the local minimum
2. NONMEM MAXEVAL=0 and SAS get stuck in the local minimum
Tom, thanks for discovering the local minimum and describing it. For the information
of those who (like me) can 'see' things better in graphical form this Excel
worksheet
http://www.health.auckland.ac.nz/pharmacology/staff/nholford/pkpd/Bayesian/BLS_solver.xls
illustrates how the OBJ value changes with differential values of K ('Ktry'). The
observation derived component (OLS) and the Bayesian derived component (BLS) of the
OBJ are graphed to show how they combine to produce the local and presumed global
minimum. If you have the Excel Solver installed you can test its performance e.g.
with initial estimates of Khat of 10 and 1.
Nick
--
Nick Holford, Dept Pharmacology & Clinical Pharmacology
University of Auckland, 85 Park Rd, Private Bag 92019, Auckland, New Zealand
email:n.holford@auckland.ac.nz tel:+64(9)373-7599x86730 fax:373-7556
http://www.health.auckland.ac.nz/pharmacology/staff/nholford/
_______________________________________________________
From: "Steve Duffull" sduffull@pharmacy.uq.edu.au
Subject: RE: [NMusers] posthoc step
Date: Wed, December 15, 2004 1:31 am
Nick
> 1. NONMEM using either the DATA or the PRIOR methods and
> WinBugs (Steve, is this correct?) get past the local minimum
I think the concept of local minima within the framework of MCMC with
variably informative priors a finite sampling space and a finite number
of samples is quite a complex issue. But I can say that I tried the
problem as discussed, starting from a variety of different initial
estimates for the chain and got the same results each time. There may
well be initial estimates that I did not try that for the limited number
of samples that I ran the chain (5000) may produce different answers.
So - I would not want to say that the MCMC process (for the number of
iterations that I used) in BUGS avoids the problem - so much as I did
not see it.
Steve
========================================Stephen Duffull
School of Pharmacy
University of Queensland
Brisbane 4072
Australia
Tel +61 7 3365 8808
Fax +61 7 3365 1688
University Provider Number: 00025B
Email: sduffull@pharmacy.uq.edu.au
www: http://www.uq.edu.au/pharmacy/sduffull/duffull.htm
PFIM: http://www.uq.edu.au/pharmacy/sduffull/pfim.htm
MCMC PK example: http://www.uq.edu.au/pharmacy/sduffull/MCMC_eg.htm
=======================================
_______________________________________________________
From: "Nick Holford" n.holford@auckland.ac.nz
Subject: Re: [NMusers] posthoc step
Date: Wed, December 15, 2004 1:47 am
Steve,
Let me ask the question more specifically -- if you used the same initial estimate
for K that was used for NONMEM, SAS, Excel Solver i.e. 10, and fixed the residual
error SD to 0.2 (no uncertainty!) with a prior on K of 10 and SD for the prior of 2
then what final estimate did you obtain? Did you get a number closer to 10 or closer
to 1?
Nick
--
Nick Holford, Dept Pharmacology & Clinical Pharmacology
University of Auckland, 85 Park Rd, Private Bag 92019, Auckland, New Zealand
email:n.holford@auckland.ac.nz tel:+64(9)373-7599x86730 fax:373-7556
http://www.health.auckland.ac.nz/pharmacology/staff/nholford/
_______________________________________________________
From: "Steve Duffull" sduffull@pharmacy.uq.edu.au
Subject: RE: [NMusers] posthoc step
Date: Wed, December 15, 2004 5:43 am
Nick
> Let me ask the question more specifically -- if you used the same initial estimate
for K that was used for NONMEM, SAS, Excel Solver i.e. 10, and fixed the residual
error SD to 0.2 (no uncertainty!) with a prior on K of 10 and SD for the prior of
2 then what final estimate did you obtain? Did you get a number closer to 10 or
closer to 1?
For the problem you describe above - I started 5 chains with initial estimates of K0
= 1, 5, 10, 20 & 50 and generated 1000 samples from the posterior distribution for
each chain. The mean of the posterior distribution of each chain (Khat) was 0.95.
Based on this scenario - WinBUGS did not experience local minima.
Steve
========================================Stephen Duffull
School of Pharmacy
University of Queensland
Brisbane 4072
Australia
Tel +61 7 3365 8808
Fax +61 7 3365 1688
University Provider Number: 00025B
Email: sduffull@pharmacy.uq.edu.au
www: http://www.uq.edu.au/pharmacy/sduffull/duffull.htm
PFIM: http://www.uq.edu.au/pharmacy/sduffull/pfim.htm
MCMC PK example: http://www.uq.edu.au/pharmacy/sduffull/MCMC_eg.htm
_______________________________________________________
From: "Ludden, Thomas (MYD)" luddent@iconus.com
Subject: RE: [NMusers] posthoc step
Date: Wed, December 15, 2004 5:32 pm
Hi Nick,
Thanks for providing the link to your Excel spreadsheet. It provides an
excellent visual display of the objective function.
I was finally able to find the place in NONMEM V where I could use a
debugger to change the initial estimate for the ETA search when MAXEVAL=0
POSTHOC was specified. Initial estimates of ETA from 0 to -2.750 (initial
est. of K = 10 to 7.25) yielded K=9.78. Initial estimates of ETA from
-2.768 to -9.9 (initial est. of K = 7.232 to 0.1) yielded K=0.944. Initial
estimates of ETA from -2.755 to -2.767 (initial est. of K = 7.245 to 7.233)
yielded values of K intermediate between 9.78 and 0.944. Over this very
small range of initial values, the step size for the ETA search is small and
the algorithm exceeds the maximum number of iterations that is allowed.
I was curious about how the "Data" Method missed the local minimum so I used
the debugger to trace the search for the "Data" method using the control
stream you provided, except that I delete the SLOW option. Starting with an
initial estimate of 10 the series of K values and the associated OFV's are
tabulated below. Note that the THETA search (different from the code used
for the ETA search) "jumps" well over both minima but then returns to the
region of the global minimum and rather quickly terminates at this minimum.
The ETA search assumes it is starting in the correct region whereas the
THETA search makes no such assumption. In this case, that is an obvious
advantage. The THETA search never detects the shallow local minimum near a
K value of 10.
K OFV
10.010 439.794
-20.000 3.2605X10^55
-5.000 2.6717X10^16
2.5 306.355
2.51 307.674
-5.000 2.6717X10^16
-1.25 4895847
0.625 134.856
0.623 125.167
0.802 30.721
0.979 14.1438
0.989 14.6477
0.962 13.5448
0.944 13.333
0.954 13.405
0.934 13.398
0.944 13.33 termination
Tom
_______________________________________________________
From: Vicente.Casabo@uv.es
Subject: RE: [NMusers] posthoc step
Date: Thu, December 16, 2004 7:57 am
Dear all
very interesting discussion. I have enjoyed it very much.
in Jerry's SAS code, if you are trying to estimate the Residual SD
as he did, values of plus minus infinity are not weird at all, are
obvious. should not you code a penalty function in that case?
I do not use SAS, but in the excel file from Nick you just have to
add the 2*LN of RUVsd to OLS (when you are estimating aleatory
effects you required extended least squares).
(If you recode the OLS in that way and add a restriction to RUVsd of
being >1*10-8 to avoid negative or zero values for RUV)
Using initial estimates of K from 10 to 1.8 solver converges to K=10
and SD=2.44
initial estimates K=0.1 and sd 0.2 lead to K 0.808 and sd0.4
initial estimates k=0.5 to 1.7 and sd=0.2 lead to K=0.94 and sd=0.12
initial estimates from K=1.8 and higher lead as I mentioned above
to K=9.9988 and SD=2.
Vicente G. Casabo
Associate Professor
University of Valencia
Spain
_______________________________________________________
From: "Nick Holford" n.holford@auckland.ac.nz
Subject: Re: [NMusers] posthoc step
Date: Thu, December 16, 2004 3:08 pm
Tom,
Thanks for the info on NONMEM's search methods. Base on this one example it looks
like one should not be too concerned about initial estimates for the usual
THETA/OMEGA/SIGMA search. It's a pity the ETA search is not so robust. Can you
confirm that the ETA search method used during MAXEVAL=0 is the same as that used by
METHOD=COND at each iteration?
Nick
--
Nick Holford, Dept Pharmacology & Clinical Pharmacology
University of Auckland, 85 Park Rd, Private Bag 92019, Auckland, New Zealand
email:n.holford@auckland.ac.nz tel:+64(9)373-7599x86730 fax:373-7556
http://www.health.auckland.ac.nz/pharmacology/staff/nholford/
_______________________________________________________
From: "Ludden, Thomas (MYD)" luddent@iconus.com
Subject: RE: [NMusers] posthoc step
Date: Thu, December 16, 2004 4:14 pm
Nick,
I believe there is just one ETA search algorithm in each version of NONMEM,
but let me confirm by tracing a simple example. I will let you know what I
find.
The ETA search in recent versions of NONMEM VI beta (since about June, 2003)
appears to be more stable than the one in NONMEM V. NONMEM VI beta does not
appear to produce the "spikes" in OFV due to inappropriate ETA estimates
that are sometimes seen with V. However, Jerry's example is a problem for
both versions.
Tom
_______________________________________________________
From: "Ludden, Thomas (MYD)" luddent@iconus.com
Subject: RE: [NMusers] posthoc step
Date: Mon, December 20, 2004 11:27 am
Nick et al.
NONMEM uses the same algorithm for all ETA searches. However, the details of
the algorithm differ between version V and version VI beta.
When MAXEVAL is not zero, the ETA search may benefit from different THETA,
OMEGA and SIGMA estimates as the overall search continues. As already
noted, rather small differences in the priors may alter the OFV shape. For
instance, using the same data (T=1,2,3,Y) as the example but somewhat
different priors, K~N(9.49,var=4.67), err~N(0,var=0.0371) instead of
K~N(10,var=4.), err~N(0,var=0.04), the Bayes OFV yields a shoulder over the
range of K=8 to K=9 but not a minimum. With this prior NONMEM V returns a K
value of 0.943 for the individual. As noted before, if one used a design
appropriate for the prior such that the data were (0.1*T,Y), NONMEM V
estimates K=1.01. In this case, the predicted Y values, based on the prior,
are above the err-SD. In fact, supplementing the previous data set with one
additional observation at T=0.1 provides a K estimate of about 1.
Best wishes to all for the holidays.
Tom
_______________________________________________________