RE: Obtaining gradients for the computation of prediction intervals
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:
http://www.cognigencorp.com/nonmem/nm/99dec071998.html ). 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