Re: posthoc step

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/