Subject:[NMusers] Obtaining gradients for the computation of prediction intervals
Date:11:38 AM 9/17/02 -0400

Hi everyone; 

      Aarons et al. (Br J Clin Pharmacol (1989) 28:305)) describes a method of calculating a prediction
	  interval for assessing how well a model predicts a validation data set. He first computes a
	  standardized prediction error for each concentration (SPE) which is 
      the SD(Observed) is obtained using the following expression (inserted as a graphic)

Can anyone tell us how to obtain the partial derivatives (gradients) which we need to compute this standard deviation? Thanks, Mike Fossler, Clinical PK, Modeling & Simulation Jagadeesh Aluri, PK Operations GlaxoSmithKline ___________________________________ From:Leonid Gibiansky Subject:Re: [NMusers] Obtaining gradients for the computation of prediction intervals Date:Tue, 17 Sep 2002 13:10:16 -0400 Michael, Unfortunately, I do not have the paper that you cited, and cannot find it on the web. But I am not sure that you can estimate SD(observed) from anything except observations. If you have identical dosing and sampling history, SD(observed) at each time point can be directly computed from the set of observations. The formula that you mentioned is based on a particular model and looks more like a SD(predicted). Concentration at any time point can be expressed as a function of thetas and etas + error term (using previously developed population model ). Then the formula can be used to compute SD (partial derivatives should be computed at etas=0). Leonid ___________________________________ From:"Gastonguay, Marc" Subject:RE: [NMusers] Obtaining gradients for the computation of prediction intervals Date:Tue, 17 Sep 2002 23:09:14 -0400 Mike, I recall that the SD you are looking for has been discussed previously in the context of prediction intervals, although I think the nmusers archive does not have the complete thread (see: ). The gradients you are looking for are stored in the NONMEM commons in the G and HH matrices for etas and epsilons, respectively. You can access these with verbatim code or presumably by using an INFN subroutine. A while ago, Tom Ludden and I worked-out some code for extracting the SD of the observation and using it in a 95% prediction interval (see below). This one works when OMEGA and SIGMA are diagonal (no covariances). A word of caution... These days, I typically create prediction intervals via simulation, as described in the archives, so this code has not been used since NONMEM IV. Also recognize that this is my implementation of the code, and any errors are mine, not Tom's. Best regards, Marc Gastonguay ; THIS IS AN NMTRAN CONTROL STREAM FOR GENERATING 95% PREDICTION INTERVALS ; OF THE TYPICAL POPULATION PREDICTION. ; USE FINAL MODEL WITH FINAL PARAMETER ESTIMATES AS INITIAL ESTIMATES AND ; MAXEVAL=0. ; USE DATA FROM A "TYPICAL" INDIVIDUAL - YOU ONLY NEED 1 INDIVIDUAL ; THE PHENOBARBITAL MODEL IS USED AS AN EXAMPLE HERE. $PROB PHENOBARBITAL POPULATION PK MODEL $INPUT ID TIME AMT WT APGR DV EVID MDV $DATA pheno $SUBROUTINES ADVAN1 $ABBREVIATED COMRES=20 ; RESERVES SPACE IN THE GLOBAL COMMON NMPRD4. YOU NEED TO RESERVE AS ; MANY SPACES (OR MORE) AS YOU HAVE ERROR TERMS. $PK ; BEGIN VERBATIM CODE " FIRST " COMMON/ROCM6/THETAF(20), OMEGAF(10,10), SIGMAF(10,10) " DOUBLE PRECISION THETAF, OMEGAF, SIGMAF "C YOU MUST DEFINE AS MANY COM() VARIABLES AS YOU HAVE ERROR TERMS. "C THE INTER- AND INTRA-INDIVIDUAL VARIANCE TERMS ARE THE DIAGONALS OF "C THE MATRICES OMEGAF AND SIGMAF, RESPECTIVELY. " COM(1)=OMEGAF(1,1) " COM(2)=OMEGAF(2,2) "C . "C . "C COM(n)=OMEGAF(n,n) "C . " COM(3)=SIGMAF(1,1) " COM(4)=SIGMAF(2,2) ; BACK TO ABBREVIATED CODE TVCL=THETA(1) CL=TVCL*DEXP(ETA(1)) TVV=THETA(2) V=TVV*DEXP(ETA(2)) K=CL/V S1=V $ERROR Y=F*(1+ERR(1))+ERR(2) ; THE FOLLOWING SECTION CALCULATES 95% PREDICTION INTERVALS OF THE PREDICTION. " LAST "C "C THE DATA FOR THE 95% PREDICTION INTERVAL WILL BE STORED IN FILE=NM.PI " OPEN(33, FILE = 'NM.PI') " IF (NEWIND.EQ.0)THEN " WRITE(33,3900) "ID","TIME","PRED","LOWER","UPPER" "3900 FORMAT(A7,4(A9)) " END IF "C "C THIS NEXT SECTION RETRIEVES VARIANCES WHICH ARE STORED IN COM(1-4). "C PARTIAL DERIVATIVES ARE ALSO RETRIEVED FROM G(I,J) AND HH(I,J). "C VETA1 IS THE VARIANCE IN PRED DUE TO ETA1. FOR EXAMPLE: "C VETA1=(PARTIAL DERIV OF PRED W.R.T. ETA1)**2*OMEGA(1), ETC... "C "C ADD VARIABLES BELOW AS NEEDED FOR THE APPROPRIATE # OF ETA & EPS "C " IF (ICALL.EQ.2)THEN " VETA1=G(1,1)**2*COM(1) " VETA2=G(2,1)**2*COM(2) "C . "C . "C VETAn=G(n,1)**2*COM(n) "C . " VEPS1=HH(1,1)**2*COM(3) " VEPS2=HH(2,1)**2*COM(4) "C "C THE TOTAL VARIANCE IS THE SUM OF EACH COMPONENT, SD=SQRT(TOTAL VAR). "C YOU MUST ADD VETA OR VEPS BELOW AS NEEDED. "C " SDPRED=SQRT(VETA1+VETA2+VEPS1+VEPS2) " WRITE(33,4000) ID,TIME,F,F-(1.96*SDPRED),F+(1.96*SDPRED) "C "4000 FORMAT(1X,F6.0,1X,F8.2,1X,3(F8.4,1X)) " " END IF $THETA 0.00548 1.4 $OMEGA 0.227 0.146 $SIGMA 0.0001 8 $EST MAXEVAL=0 ___________________________________ From:Alice I Nichols [] Subject:RE: [NMusers] Obtaining gradients for the computation of prediction intervals Date:Friday, September 20, 2002 11:24 AM Mike, As a followup regarding where to find these derivatives. Here's a way I have used to reliably identify partials related to this equation. If you look at the fortran code generated by NMTRAN for your model (fsubs.for) you can find the exact variable names for both partials w.r.t. etas and partials w.r.t. epsilons. There are comment lines identifying these equations. Since there are usually several equation lines throughout the code for each partial, I suggest you follow each equation thru all the code to correctly identify these variables. After identifing these parameters you can revise your control stream to add these parameters to your output file. As Marc Gastonguay mentioned however most people will do simulations and calculate SD based on this. Alice [Alice I Nichols, PhD Hawthorne Research and Consulting, INC 132 Hawthorn Rd King of Prussia, PA 19406 PH: 610-878-9112 FX: 610-878-9113] ___________________________________ From:"Alice I Nichols" Subject: [NMusers] Obtaining gradients for the computation of prediction intervals Date:Fri, 20 Sep 2002 11:34:13 -0700 Mike as a followup to my previous note the fortran subroutine to look at is for the error subroutine written by nmtran for your control stream. Alice. [Alice I Nichols, PhD Hawthorne Research and Consulting, INC 132 Hawthorn Rd King of Prussia, PA 19406 PH: 610-878-9112 FX: 610-878-9113] ___________________________________