RE: posthoc step
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