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.pdf
Bayesian 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] <- 10*exp(-k*time[j]) } k ~dnorm(10,20) # this is mean, precision where precision = 1/var sigma2 <- 0.04 tau <- 1/(sigma2) } list( time=c(1, 2, 3), data=c(3.87, 1.66, 0.44) ) # data list( k=2) # initial estimate of k (I tried others with the same result) When the prior variance was set to 4 Khat was estimated at 0.95 When the prior variance was set to 1 Khat was estimated at 0.96 When the prior variance was set to 0.1 Khat was estimated at 1.2 When the prior variance was set to 0.05 Khat was estimated at 10 These are similar to the NONMEM and SAS results when variance was 4. In addition, there appeared to be a rapid change in Khat occurred with small changes in variance (results not shown). 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 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] <- 10*exp(-k*time[j]) } omega2 <- 0.095 #sharp change here io<- 1/(omega2) k ~dnorm(10,io) } list(time=c(1, 2, 3),data=c(3.87, 1.66, 0.44)) list(k=2) Yaning Wang, PhD Pharmacometrician OCPB, FDA _______________________________________________________ From: "Nick Holford" n.holford@auckland.ac.nz Subject: Re: [NMusers] posthoc step Date: Thu, December 9, 2004 10:25 pm Jerry, I enjoyed reading through your posting to nmusers with the analytical solution for the linear model. I had not appreciated how much the precision of the observations dominated the solution. This made me realize that the model I had been using to test this with NONMEM differs from yours because I was estimating both K *and* the residual error variance. When I fix the residual error to the prior estimate of 0.04 (which is the more usual thing to do for this kind of Bayesian estimation) then NONMEM gives an estimate of Khat of 0.944 using both the DATA and PRIOR methods. This is now in good agreement with the estimates from SAS/MAP (Jerry) and WinBUGS (Steve Duffull). Apologies to all for the confusion created by my original posting of results from NONMEM using a different model. I am attaching an updated PDF showing the change in Khat with the prior variance of K. As noted by Yaning using WinBUGS there is quite a sharp transition from Khat ~10 to Khat ~1 with prior variance of K going from 0.094 to 0.095. NONMEM differs in having a transition when going from 0.3 to 0.4. I've double checked these NONMEM values - they are variances not SDs. Based on the use of a comparable model it seems that NONMEM gets somewhat similar results to SAS and WinBUGS when estimation is performed using the DATA or the PRIOR methods. It is still unexplained why the MAXEVAL=0 method gives very different results (Khat 9.78 with prior variance of K = 4) Nick 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 _______________________________________________________