RE: How to calculate confidential interval of a plasma conc. curve

From: "Piotrovskij, Vladimir [JanBe]" <VPIOTROV@janbe.jnj.com> Subject: RE: How to calculate confidential interval of a plasma conc. curve Date: Mon, 7 Dec 1998 15:08:32 +0100 Dear Takeshi, You can use NONMEM to perform Monte-Carlo simulations needed for constructing confidence intervals for predictions, but you have to create first a suitable data set that is cumbersome. If you have access to S-PLUS you can obtain the plots you want in a few seconds. Below you will find an S-PLUS script which does the job for a simple model you mention in your mail. I used the following arbitrary parameter values. Fixed effects were: THETA(1)=5, THETA(2)=.2, THETA(3)=100. SD for random effects were: SD(ETA(1))=SD(ETA(2)).2, SD(EPS(1))=.1. The time interval for simulations was 0 to 24 h. WT was a normal random variable with mean=70 and SD=4. If you source the script to S-PLUS you will get a plot with the following lines: a model prediction for a 70-kg subject; 5% and 95% confidence intervals for typical values (referred to variability in WT); 5% and 95% confidence intervals for individual concentrations (residual error excluded); 5% and 95% confidence intervals for individual concentrations (residual error included). The script can be easily updated for more complex models. Apparently, other packages can be applied as well. BTW, using partial derivatives to obtain confidence intervals is good for low variability. If inter- or intraindividual variability or SE for parameter estimates exceeds, say, 30% you may get highly biased CI since Wald approximation which underlies this approach usually fails in this situation. Hope this helps Vladimir #*************** S-PLUS script ************** n _ 200 dose _ 100 # Fixed effect parameters TH1 _ 5 TH2 _ .2 TH3 _ 100 # Random effects ETA1 _ rnorm(n,sd=.2) ETA2 _ rnorm(n,sd=.2) EPS1 _ rnorm(n,sd=.1) # Covariate WT _ rnorm(n, mean=70,sd=4) # Typical values TVCL _ TH1 + TH2*WT TVV _ TH3 # Individual values CL _ TVCL*(1+ETA1) V _ TVV*(1+ETA2) # Conc predictions x _ seq(0,24,.5) time _ t(matrix(rep(x,n),ncol=n)) PRED.70 _ dose/TVV*exp(-(TH1 + TH2*70)/TVV*x) PRED _ dose/TVV*exp(-TVCL/TVV*time) IPRED _ dose/V*exp(-CL/V*time) conc _ IPRED*(1+EPS1) # Plotting xlm _ range(x) ylm _ range(conc) plot(1,1,log="y",type="n",xlab="Time",ylab="Conc",xlim=xlm,ylim=ylm) lines(x,PRED.70,lwd=4) y _ apply(PRED,2,quantile, probs=c(.05,.95)) lines(x,y[1,],lwd=2) ; lines(x,y[2,],lwd=2) y _ apply(IPRED,2,quantile, probs=c(.05,.95)) lines(x,y[1,],lwd=.5) ; lines(x,y[2,],lwd=.5) y _ apply(conc,2,quantile, probs=c(.05,.95)) lines(x,y[1,],lwd=.5) ; lines(x,y[2,],lwd=.5)