/**********************************************************************/ /* Atkinson, Donev, and Tobias, "Optimum Experimental Designs" */ /* */ /* NAME: Program 18.4 */ /* TITLE: Tables 18.8-18.11 - Bayesian D- and c-optimum designs */ /* for compartmental model */ /* KEYS: */ /* DATA: */ /* NOTES: Uses the QUAD subroutine in IML to perform double */ /* integration over the prior distribution for the */ /* parameters. Took ~15 minutes to run on server-class PC */ /* in May, 2007. */ /* */ /* Author: Randy Tobias */ /* History: */ /* Created...............................................27Jun2007 */ /**********************************************************************/ title1 "Example 18.3 - Compartmental model: Bayesian D- and c-optimum designs"; /* / Create data set with original data. /---------------------------------------------------------------------*/ data Theophylline; input Time Concentration; Time = round(60*Time)/60; cards; 0.166 10.1 0.333 14.8 0.5 19.9 0.666 22.1 1 20.8 1.5 20.3 2 19.7 2.5 18.9 3 17.3 4 16.1 5 15.0 6 14.2 8 13.2 10 12.3 12 10.8 24 6.5 30 4.6 48 1.7 ; %macro Tables; options nonotes; proc iml; /* / The approach is much the same as we've used in many other programs, / nonlinear optimization of first w for a grid of x values, and then / of sqrt(w) and x jointly. /---------------------------------------------------------------------*/ start MakeZ(Time,t1,t2,t3); emt1Time = exp(-t1*Time); emt2Time = exp(-t2*Time); Z = t3* Time#emt1Time || -t3* Time#emt2Time || ( emt2Time - emt1Time); return(Z); finish; start Consolidate(x,w,ex,ew); x = x[loc(w > ew)]`; w = w[loc(w > ew)]`; Group = shape(0,1,ncol(x)); do i = 1 to ncol(x); if (^Group[i]) then do; Group[i] = max(Group) + 1; do j = i+1 to ncol(x); if (abs(x[i] - x[j]) < ex) then Group[j] = Group[i]; end; end; end; xNew = shape(.,1,max(Group)); wNew = shape(.,1,max(Group)); do i = 1 to max(Group); xNew[i] = (x[loc(Group = i)])[:]; wNew[i] = (w[loc(Group = i)])[+]; end; r = rank(xNew); x = xNew; w = wNew; x[r] = xNew; w[r] = wNew; finish; start llog(x); if (x > 0) then return(log(x)); else return(-9999); finish; /* / Define critical value functions for D-optimality. The fundamental / DCrit(w,x) is a double-integral. /---------------------------------------------------------------------*/ start DCritFun1(t1) global(xFun,wFun,t2Fun); Z = MakeZ(xFun,t1,t2Fun,&Truet3); Mi = Z`*diag(wFun)*Z; return(llog(det(Mi))); finish; start DCritFun2(t2) global(t2Fun); t2Fun = t2; call quad(z,"DCritFun1",(&t1Lo || &t1Hi)) eps=1e-4 cycles=1; return(z/(&t1Hi-&t1Lo)); finish; start DCrit(w,x) global(xFun,wFun); xFun = x; wFun = w; call quad(z,"DCritFun2",(&t2Lo || &t2Hi)) eps=1e-4 cycles=1; return(z/(&t2Hi-&t2Lo)); finish; start DCritW(w) global(x); return(DCrit((w##2)/sum(w##2),x)); finish; start DCritX(x) global(w); return(DCrit(w,x)); finish; start DCritWX(wx); w = shape(wx,2,ncol(wx)/2)[1,]`; w = (w##2)/sum(w##2); x = shape(wx,2,ncol(wx)/2)[2,]`; return(DCrit(w,x)); finish; /* / Compute the D-optimum design. /---------------------------------------------------------------------*/ nint = 10; x = 20*2**(-(0:nint))`; winit = j(nrow(x),1) / nrow(x); call nlpqn(rc,wOp,"DCritW",winit,1); wOp = ((wOp##2)/sum(wOp##2)); xOp = x`; call Consolidate(xOp,wOp,0,1e-5); i = 0; diff = 1; do while ((i < 100) & (diff > 1e-6)); i = i + 1; last = (xOp` || wOp`); wxinit = sqrt(wOp) || xOp; con = shape(.,1,ncol(xOp)) || shape(0,ncol(xOp),1)`; call nlpqn(rc,wxOp,"DCritWX",wxinit,1,con); wOp = shape(wxOp,2,ncol(wxOp)/2)[1,]; wOp = (wOp##2)/sum(wOp##2); xOp = shape(wxOp,2,ncol(wxOp)/2)[2,]; call Consolidate(xOp,wOp,1e-2,2e-3); if (ncol(xOp) ^= nrow(last)) then diff = 1; else diff = max(abs(last - (xOp` || wOp`))); end; xwDOp = xOp` || wOp`; start cDCritFun1(t1) global(cType,xFun,wFun,t2Fun); t2 = t2Fun; t3 = &Truet3; if (cType = 1) then do; c = -t3/(t1**2) // t3/(t2**2) // (1/t1 - 1/t2); end; else if (cType = 2) then do; a = t2 - t1; b = log(t2/t1); c = (b-a/t1)/(a**2) // (a/t2-b)/(a**2) // 0; end; else do; a = t2 - t1; b = log(t2/t1); c2 = (b-a/t1)/(a**2) // (a/t2-b)/(a**2) // 0; xmax = b/a; e1 = exp(-t1*xmax); e2 = exp(-t2*xmax); f = t2*e2 - t1*e1; c = t3*(-xmax*e1 + f*c2[1]) // t3*( xmax*e2 + f*c2[2]) // e1-e2; end; Z = MakeZ(xFun,t1,t2,t3); Mi = Z`*diag(wFun)*Z + 1e-5*i(ncol(Z)); return(c`*inv(Mi)*c); finish; start cDCritFun2(t2) global(t2Fun); t2Fun = t2; call quad(z,"cDCritFun1",(&t1Lo || &t1Hi)) eps=1e-4 cycles=1; return(z/(&t1Hi-&t1Lo)); finish; start cDCrit(w,x) global(xFun,wFun); xFun = x; wFun = w; call quad(z,"cDCritFun2",(&t2Lo || &t2Hi)) eps=1e-4 cycles=1; return(z/(&t2Hi-&t2Lo)); finish; start cDCritW(w) global(x); return(cDCrit((w##2)/sum(w##2),x`)); finish; start cDCritX(x) global(w); return(cDCrit(w,x`)); finish; start cDCritWX(wx); w = shape(wx,2,ncol(wx)/2)[1,]; w = (w##2)/sum(w##2); x = shape(wx,2,ncol(wx)/2)[2,]; return(cDCrit(w`,x`)); finish; start cDOpt(cIn,nInt) global(cType,x,w); cType = cIn; x = 20*2**(-(0:nint)); winit = j(ncol(x),1) / ncol(x); call nlpqn(rc,wOp,"cDCritW",winit,0); wOp = ((wOp##2)/sum(wOp##2)); xOp = x`; call Consolidate(xOp,wOp,0,1e-5); i = 0; diff = 1; do while ((i < 100) & (diff > 1e-6)); i = i + 1; last = (xOp` || wOp`); wxinit = sqrt(wOp) || xOp; con = shape(.,1,ncol(xOp)) || shape(0,ncol(xOp),1)`; call nlpqn(rc,wxOp,"cDCritWX",wxinit,0,con); wOp = shape(wxOp,2,ncol(wxOp)/2)[1,]; wOp = (wOp##2)/sum(wOp##2); xOp = shape(wxOp,2,ncol(wxOp)/2)[2,]; call Consolidate(xOp,wOp,1e-2,2e-3); if (ncol(xOp) ^= nrow(last)) then diff = 1; else diff = max(abs(last - (xOp` || wOp`))); end; return (xOp` || wOp`); print (xOp` || wOp`); finish; options nonotes; xwOp1 = cDOpt(1,10); cDCrit1 = cDCrit(xwOp1[,2],xwOp1[,1]); xwOp2 = cDOpt(2,10); cDCrit2 = cDCrit(xwOp2[,2],xwOp2[,1]); xwOp3 = cDOpt(3,10); cDCrit3 = cDCrit(xwOp3[,2],xwOp3[,1]); use Theophylline; read all var {Time} into x18; w18 = j(nrow(x18),1)/nrow(x18); d = exp(DCrit(xwDOp[,2],xwDOp[,1])/3) // exp(DCrit(xwOp1[,2],xwOp1[,1])/3) // exp(DCrit(xwOp2[,2],xwOp2[,1])/3) // exp(DCrit(xwOp3[,2],xwOp3[,1])/3) // exp(DCrit(w18 ,x18 )/3); cType = 1; d = d || ( (1/cDCrit(xwDOp[,2],xwDOp[,1])) // (1/cDCrit(xwOp1[,2],xwOp1[,1])) // (1/cDCrit(xwOp2[,2],xwOp2[,1])) // (1/cDCrit(xwOp3[,2],xwOp3[,1])) // (1/cDCrit(w18 ,x18 ))); cType = 2; d = d || ( (1/cDCrit(xwDOp[,2],xwDOp[,1])) // (1/cDCrit(xwOp1[,2],xwOp1[,1])) // (1/cDCrit(xwOp2[,2],xwOp2[,1])) // (1/cDCrit(xwOp3[,2],xwOp3[,1])) // (1/cDCrit(w18 ,x18 ))); cType = 3; d = d || ( (1/cDCrit(xwDOp[,2],xwDOp[,1])) // (1/cDCrit(xwOp1[,2],xwOp1[,1])) // (1/cDCrit(xwOp2[,2],xwOp2[,1])) // (1/cDCrit(xwOp3[,2],xwOp3[,1])) // (1/cDCrit(w18 ,x18 ))); do i = 1 to ncol(d); d[,i] = d[,i]/max(d[,i]); end; print "Designs",, "Crit" " Time" " Weight" " Crit",, "DOpt" (1*xwDOp) [format=7.4] (DCrit(xwDOp[,2],xwDOp[,1])) [format=8.4],, "AUC " (1*xwOp1) [format=7.4] (1*cDCrit1) [format=8.1],, "tmax" (1*xwOp2) [format=7.4] (1*cDCrit2) [format=8.6],, "yMax" (1*xwOp3) [format=7.4] (1*cDCrit3) [format=8.4] ; name = {"DOpt" "AUC " "tMax" "yMax" "18-pnt"}; Efficiencies = d; print Efficiencies [rowname=name colname=name format=percent8.1]; %mend; %let Truet1 = 0.05884; %let Truet2 = 4.298; %let Truet3 = 21.80; %let t1Lo = (&Truet1-0.01); %let t1Hi = (&Truet1+0.01); %let t2Lo = (&Truet2-1); %let t2Hi = (&Truet2+1); title2 "Prior distribution I"; %Tables; title2; %let t1Lo = (&Truet1-0.04); %let t1Hi = (&Truet1+0.04); %let t2Lo = (&Truet2-4); %let t2Hi = (&Truet2+4); title2 "Prior distribution II"; %Tables; title2; title1;