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