/**********************************************************************/ /* Atkinson, Donev, and Tobias, "Optimum Experimental Designs" */ /* */ /* NAME: Program 20.3 */ /* TITLE: Example 20.3 & SAS Task 20.3 - Non-linear model for crop */ /* yield and plant density */ /* KEYS: */ /* DATA: */ /* NOTES: Using both large exact designs with OPTEX and nonlinear */ /* optimization with IML. Incorporates a two-factor version */ /* of the now-familiar template for finding an optimum */ /* continuous design by nonlinear optimization. */ /* */ /* Author: Randy Tobias */ /* History: */ /* Created...............................................27Jun2007 */ /**********************************************************************/ title1 "Example 20.3 - Non-linear model for crop yield and plant density"; /* / Table 20.2 - Preliminary data: Mean grain yield (g/m**2) for given / intra-row spacing (x1) and inter-row spacing (x2). /---------------------------------------------------------------------*/ data Beans; do x2 = 0.18 to 0.76 by 0.18; do x1 = 0.03 to 0.12 by 0.03; input Yield @@; logY = log(Yield); output; end; end; cards; 260.0 344.7 279.9 309.2 305.3 358.3 312.2 267.8 283.9 342.0 269.0 253.9 221.8 287.9 230.9 196.9 ; /* / Use preliminary to compute prior estimates of t1, y4, and t7. All / other parameter priors are set to zero. /---------------------------------------------------------------------*/ proc nlin data=Beans; parms t1 = 1 t4 = 1 t7 = 1; model logY = -(1/t7)*log(t1 + t4/(x1*x2)) - log(x1*x2); ods output ParameterEstimates=Parm; data _null_; set Parm; call symput(Parameter,trim(left(Estimate))); run; %put t1 = &t1; %let t2 = 0; %let t3 = 0; %put t4 = &t4; %let t5 = 0; %let t6 = 0; %put t7 = &t7; %macro N0(alpha,N); /* Inverting equation 20.5: */ %if (%sysevalf(&alpha > 0)) %then %do; /* prior N0 corresponding */ %sysevalf(((1-&alpha)/&alpha)*&N) /* to given alpha */ %end; %else %do; %sysevalf(1000*&N) %end; %mend; /* / Perform orthogonal scaling, to make the assumption of a common / prior variance tenable, as discussed in Section 20.4.2 of the text. /---------------------------------------------------------------------*/ data Can; do x1 = 0.15 to 0.8 by 0.01; do x2 = 0.03 to 0.2 by 0.01; f1 = -(&t7*(&t1+&t4/(x1*x2)))**(-1); /* Compute Jacobian */ f2 = f1/ x1 ; f3 = f1/ x2 ; f4 = f1/(x1*x2); f5 = f1/(x1*x1); f6 = f1/(x2*x2); f7 = ((&t7)**(-2))*log(&t1+&t4/(x1*x2)); output; end; end; proc reg data=Can noprint; /* Compute residuals */ model f2 f3 f5 f6 = f1 f4 f7; output out=rCan r = rf2 rf3 rf5 rf6; proc iml; /* Scale residuals */ use rCan; read all var { f1 f4 f7} into Fr; read all var {rf2 rf3 rf5 rf6} into Rs; Ws = diag(1/(Rs[<>,] - Rs[><,])); Fs = Rs*Ws; create Fs var {rf2 rf3 rf5 rf6}; append from Fs; data rJac; merge Can Fs; run; /* / Use OPTEX to generate large exact designs for given alpha. /---------------------------------------------------------------------*/ %macro Exact(N,alpha,ds,wgt); proc optex data=rJac coding=none noprint; model f1 f4 f7,rf2 rf3 rf5 rf6 / noint prior = 0,%N0(&alpha,&N); generate n=&N method=m_fedorov niter=100 keep=10; id x1 x2; output out=Design; proc freq data=Design noprint; table x1*x2 / list out=&ds(rename=(Percent=&wgt)); data &ds; set &ds; Alpha = α run; %mend; %Exact(1000,0 ,aFrq,Frq); %Exact(1000,1 ,bFrq,Frq); %Exact(1000,0.25,cFrq,Frq); %Exact( 7,1 ,Frq7,Frq); %Exact( 6,0.25,Frq6,Frq); data Grid; /* Jacobian over grid */ do ix1 = 1 to 14; x1 = 0.15 + (ix1-1)*(0.8-0.15)/(14-1); do ix2 = 1 to 14; x2 = 0.03 + (ix2-1)*(0.2-0.03)/(14-1); output; end; end; data gJac; set Grid; f1 = -(&t7*(&t1+&t4/(x1*x2)))**(-1); f2 = f1/ x1 ; f3 = f1/ x2 ; f4 = f1/(x1*x2); f5 = f1/(x1*x1); f6 = f1/(x2*x2); f7 = ((&t7)**(-2))*log(&t1+&t4/(x1*x2)); run; /* / Optimal weights for given alpha. /---------------------------------------------------------------------*/ options nonotes; proc iml; start MakeZ(x) global(B,Ws,Wr); x1 = x[,1]; x2 = x[,2]; ff1 = -(&t7*(&t1+&t4#(1/(x1#x2))))##(-1); ff2 = ff1/x1; ff3 = ff1/x2; ff4 = ff1/(x1#x2); ff5 = ff1/(x1#x1); ff6 = ff1/(x2#x2); ff7 = ((&t7)**(-2))*log(&t1+&t4#(1/(x1#x2))); Fr = ff1 || ff4 || ff7; Fs = ff2 || ff3 || ff5 || ff6; F = Fr*Wr || (Fs - Fr*B)*Ws; return(F); finish; start llog(x); if (x > 0) then return(log(x)); else return(-99999); finish; start DCritWF(w,F) global(alpha); Dw = diag((w##2)/sum(w##2)); if (alpha) then do; return(llog(det( (1-alpha)*diag((1:7)>3) + alpha *F`*Dw*F ))/ncol(F)); end; else do; return(llog(det(F[,1:3]`*Dw*F[,1:3]))/3); end; finish; start DCrit(w,x); return(DCritWF(w,MakeZ(x))); finish; start DCritW(w) global(F); return(DCritWF(w,F)); finish; start DCritX(xx) global(w); n = ncol(xx)*nrow(xx)/2; x1 = shape(xx[ 1:n ],n,1); x2 = shape(xx[n+(1:n)],n,1); F = MakeZ(x1 || x2); return(DCritWF(w,F)); finish; start DCritWX(wx); n = ncol(wx)*nrow(wx)/3; w = shape(wx[ 1:n ],n,1); x1 = shape(wx[ n+(1:n)],n,1); x2 = shape(wx[2*n+(1:n)],n,1); F = MakeZ(x1 || x2); return(DCritWF(w,F)); finish; /* / A multicolumn version of the function that is by now familiar. /------------------------------------------------------------------*/ if (inputn(symget('sysver'),'3.1') >= 9) then do; start sort(x,idx); sortx = j(nrow(x), ncol(x), .); do i = nrow(idx)*ncol(idx) to 1 by -1; sortx[rank(x[,idx[i]]),] = x[]; end; x = sortx; finish; end; start Consolidate(x,w,ex,ew); x = x[loc(w > ew),]; w = w[loc(w > ew)]; Group = shape(0,1,nrow(x)); do i = 1 to nrow(x); if (^Group[i]) then do; Group[i] = max(Group) + 1; do j = i+1 to nrow(x); if (max(abs(x[i,] - x[j,])) < ex) then Group[j] = Group[i]; end; end; end; xNew = shape(.,max(Group),ncol(x)); wNew = shape(.,max(Group),1); do i = 1 to max(Group); xNew[i,] = (x[loc(Group = i),])[:,]; wNew[i] = (w[loc(Group = i)])[+]; end; xw = xNew || wNew; call sort(xw,1:ncol(x)); w = xw[,ncol(x)+1]; x = xw[,1:ncol(x)]; finish; use gJac; read all var {x1 x2} into XCan; read all var {f1 f4 f7} into Fr; read all var {f2 f3 f5 f6} into Fs; B = inv(Fr`*Fr)*Fr`*Fs; Rs = Fs - Fr*B; Wr = diag(1/(Fr[<>,] - Fr[><,])); Ws = diag(1/(Rs[<>,] - Rs[><,])); start AlphaOpt; x = XCan; F = MakeZ(x); winit = j(nrow(X),1); call nlpqn(rc,w,"DCritW",winit,1); wOp = w; wOp = ((wOp##2)/sum(wOp##2))`; xOp = x; call Consolidate(xOp,wOp,0,1e-4); i = 0; Improvement = 1; do while ((i < 100) & (Improvement > 1e-4)); i = i + 1; last = (wOp // xOp[,1] // xOp[,2])`; n = nrow(xOp); con = (shape( . ,1,n) || shape(0.15,1,n) || shape(0.03,1,n)) // (shape( . ,1,n) || shape(0.8 ,1,n) || shape(0.2 ,1,n)); wxinit = (sqrt(wOp) // xOp[,1] // xOp[,2])`; call nlpqn(rc,wx,"DCritWX",wxinit,1,con); wOp = shape(wx[ 1:n ],n,1); xOp = shape(wx[ n+(1:n)],n,1) || shape(wx[2*n+(1:n)],n,1); wOp = (wOp##2)/sum(wOp##2); call Consolidate(xOp,wOp,1e-4,1e-3); next = (wOp // xOp[,1] // xOp[,2])`; if (ncol(next) ^= ncol(last)) then Improvement = 999; else Improvement = max(abs(last - next)); end; finish; alpha = 0; call AlphaOpt; data = wOp || xOp; create aWgt var {Wgt x1 x2}; append from data; alpha = 1; call AlphaOpt; data = wOp || xOp; create bWgt var {Wgt x1 x2}; append from data; f7Round = round(7*wOp); x7Round = xOp [loc(f7Round > 0),]; f7Round = f7Round[loc(f7Round > 0),]; use Frq7; read all var {Frq x1 x2}; f7Exact = round(7*Frq/100); x7Exact = x1 || x2; print "SAS Task 20.3 - Comparing exact optimum and rounded continuous designs", "Alpha = 1, N=7",, f7Round [format=1.] x7Round [format=5.3] ( exp(DCrit(sqrt(f7Round/7),x7Round)) / exp(DCrit(sqrt(wOp) ,xOp ))) [format=percent9.2], f7Exact [format=1.] x7Exact [format=5.3] ( exp(DCrit(sqrt(f7Exact/7),x7Exact)) / exp(DCrit(sqrt(wOp) ,xOp ))) [format=percent9.2]; alpha = 0.25; call AlphaOpt; data = wOp || xOp; create cWgt var {Wgt x1 x2}; append from data; f6Round = round(6*wOp); x6Round = xOp [loc(f6Round > 0),]; f6Round = f6Round[loc(f6Round > 0),]; use Frq6; read all var {Frq x1 x2}; f6Exact = round(6*Frq/100); x6Exact = x1 || x2; print "SAS Task 20.3 - Comparing exact optimum and rounded continuous designs", "Alpha = 0.25, N=6",, f6Round [format=1.] x6Round [format=5.3] ( exp(DCrit(sqrt(f6Round/6),x6Round)) / exp(DCrit(sqrt(wOp) ,xOp ))) [format=percent9.2], f6Exact [format=1.] x6Exact [format=5.3] ( exp(DCrit(sqrt(f6Exact/6),x6Exact)) / exp(DCrit(sqrt(wOp) ,xOp ))) [format=percent9.2]; data Fig2003a; set aFrq bFrq cFrq; Frq = Frq / 100; x1x2 = x1*x2; proc sort data=aWgt; by x1 x2; data aWgt; set aWgt; Alpha = 0; proc sort data=bWgt; by x1 x2; data bWgt; set bWgt; Alpha = 1; proc sort data=cWgt; by x1 x2; data cWgt; set cWgt; Alpha = 0.25; data Fig2003b; set aWgt bWgt cWgt; x1x2 = x1*x2; run; title2 "Optimal frequencies for 1000 runs"; proc print data=Fig2003a noobs; where (Alpha = 0); by notsorted Alpha; var x1 x2 x1x2 Frq; format x1 x2 Frq 5.3; run; title2; title2 "Optimal weights"; proc print data=Fig2003b noobs; where (Alpha = 0); by notsorted alpha; var x1 x2 x1x2 Wgt; format x1 x2 Wgt 5.3; run; title2; title2 "Optimal frequencies for 1000 runs"; proc print data=Fig2003a noobs; where (Alpha ^= 0); by notsorted Alpha; var x1 x2 Frq; format x1 x2 Frq 5.3; run; title2; title2 "Optimal weights"; proc print data=Fig2003b noobs; where (Alpha ^= 0); by notsorted alpha; var x1 x2 Wgt; format x1 x2 Wgt 5.3; run; title2; title1;