/**********************************************************************/ /* Atkinson, Donev, and Tobias, "Optimum Experimental Designs" */ /* */ /* NAME: Program 17.4 */ /* TITLE: Tables 17.1 and 17.2 - D- and c-optimum designs for a */ /* compartmental model */ /* KEYS: */ /* DATA: */ /* */ /* Author: Randy Tobias */ /* History: */ /* Created...............................................27Jun2007 */ /**********************************************************************/ title1 "Tables 17.1 and 17.2 - D- and c-optimum designs for compartmental model"; /* / The experimental design of Fresen. /---------------------------------------------------------------------*/ data Theophylline; input Time Concentration; Time = round(60*Time)/60; datalines; 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 ; title2 "NLIN analysis of original design"; proc nlin data=Theophylline; parm t1=1 t2=1 t3=1; model Concentration = t3*(exp(-t1*Time) - exp(-t2*Time)); run; title2; title2 "NLMIXED analysis, with estimated functions of interest"; proc nlmixed data=Theophylline; parameters t1=1 t2=1 t3=1; model Concentration ~ normal(t3*(exp(-t1*Time) - exp(-t2*Time)),s2); Estimate "AUC" t3*(1/t1 - 1/t2); Estimate "MaxTime" (log(t2) - log(t1))/(t2 - t1); Estimate "MaxConc" t3*( exp(-t1*((log(t2) - log(t1))/(t2 - t1))) - exp(-t2*((log(t2) - log(t1))/(t2 - t1)))); ods output ParameterEstimates=Parm; run; title2; /* / Set assumed parameter values to the estimated ones. /---------------------------------------------------------------------*/ data _null_; set Parm; call symput(Parameter,trim(left(Estimate))); run; proc iml; /* / The Jacobian at given time points. /------------------------------------------------------------------*/ start MakeZ(Time); Z = &t3* Time#exp(-&t1*Time) || -&t3* Time#exp(-&t2*Time) || ( exp(-&t2*Time) - exp(-&t1*Time)); return(Z); finish; /* / First, define a general function for finding an optimum one- / factor design, depending on a function Crit(w,x) of the weights / and factor values. The approach is the same as that of Program / 9.1, which see for comments on the method. /------------------------------------------------------------------*/ 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; start CritW(w) global(x); return(Crit((w##2)/sum(w##2),x)); finish; start CritX(x) global(w); return(Crit(w,x`)); finish; start CritWX(wx); w = shape(wx,2,ncol(wx)/2)[1,]; w = (w##2)/sum(w##2); x = shape(wx,2,ncol(wx)/2)[2,]; return(Crit(w`,x`)); finish; start Optimize(nInt) global(x,w); x = 20*(1:nint)`/nint; winit = j(nint,1) / nint; call nlpqn(rc,wOp,"CritW",winit,1); wOp = ((wOp##2)/sum(wOp##2))`; xOp = x`; call Consolidate(xOp,wOp,1e-2,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,"CritWX",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; return (xOp` || wOp`); finish; /* / Critical functions for D- and c-optimality. /---------------------------------------------------------------------*/ start DCrit(w,x); Z = MakeZ(x); d = det(Z`*diag(w)*Z); return(llog(d)); finish; start cCrit(w,x) global(c); Z = MakeZ(x); Dw = diag(w); return(c`*inv(Z`*Dw*Z + 1e-5*i(ncol(Z)))*c); finish; /* / Set Crit <- DCrit and run Optimize to compute D-optimum design. /---------------------------------------------------------------------*/ start Crit(w,x); return(DCrit(w,x)); finish; xwDOp = Optimize(50); wOp = xwDOp[,2]`; xOp = xwDOp[,1]`; /* / Coefficients of linearized combinations for / * c1 - area under the curve, / * c2 - time to maximum concentration, and / * c3 - maximum concentration, / as defined in Section 17.5 of the text. /---------------------------------------------------------------------*/ c1 = -&t3/(&t1**2) // &t3/(&t2**2) // (1/&t1 - 1/&t2); 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; c3 = &t3*(-xmax*e1 + f*c2[1]) // &t3*( xmax*e2 + f*c2[2]) // e1-e2; /* / Set Crit <- -cCrit (in order to minimize variance) and run Optimize / for each linearized combination c1, c2, and c3 in turn, to compute / c-optimum designs. /---------------------------------------------------------------------*/ start Crit(w,x); return(-cCrit(w,x)); finish; options nonotes; c = c1; xwOp1 = Optimize(50); cDCrit1 = cCrit(xwOp1[,2],xwOp1[,1]); c = c2; xwOp2 = Optimize(50); cDCrit2 = cCrit(xwOp2[,2],xwOp2[,1]); c = c3; xwOp3 = Optimize(50); cDCrit3 = cCrit(xwOp3[,2],xwOp3[,1]); options notes; /* / Read in the original experimental design. /---------------------------------------------------------------------*/ use Theophylline; read all var {Time} into x18; w18 = j(nrow(x18),1)/nrow(x18); /* / Compute the four criteria (D and c for each of c1, c2, and c3) for / each of the five designs (D-optimum; c-optimum for each of c1, c2, / and c3; and the original design) ... /---------------------------------------------------------------------*/ 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); c = c1; d = d || ( (1/cCrit(xwDOp[,2],xwDOp[,1])) // (1/cCrit(xwOp1[,2],xwOp1[,1])) // (1/cCrit(xwOp2[,2],xwOp2[,1])) // (1/cCrit(xwOp3[,2],xwOp3[,1])) // (1/cCrit(w18 ,x18 ))); c = c2; d = d || ( (1/cCrit(xwDOp[,2],xwDOp[,1])) // (1/cCrit(xwOp1[,2],xwOp1[,1])) // (1/cCrit(xwOp2[,2],xwOp2[,1])) // (1/cCrit(xwOp3[,2],xwOp3[,1])) // (1/cCrit(w18 ,x18 ))); c = c3; d = d || ( (1/cCrit(xwDOp[,2],xwDOp[,1])) // (1/cCrit(xwOp1[,2],xwOp1[,1])) // (1/cCrit(xwOp2[,2],xwOp2[,1])) // (1/cCrit(xwOp3[,2],xwOp3[,1])) // (1/cCrit(w18 ,x18 ))); /* / ... and divide by each criterions maximum to compute efficiencies. /---------------------------------------------------------------------*/ do i = 1 to ncol(d); d[,i] = d[,i]/max(d[,i]); end; print "Table 17.1 - D- and c-optimum designs",, "Design" "Support Weight " "Crit ",, "D-opt " (1*xwDOp) [format=7.4] (DCrit(xwDOp[,2],xwDOp[,1])) [format=7.4],, "AUC " (1*xwOp1) [format=7.4] (1*cDCrit1) [format=7. ],, "tMax " (1*xwOp2) [format=7.4] (1*cDCrit2) [format=7.5],, "yMax " (1*xwOp3) [format=7.4] (1*cDCrit3) [format=7.4] ; print "Table 17.2 - Efficiencies of designs for compartmental models",, "Design" " D-Eff AUC-Eff tMax-Eff yMax-Eff ",, ({"D-opt " "AUC " "tMax " "yMax " "18 pnt "}`) (1*d) [format=percent8.2];