/**********************************************************************/ /* Atkinson, Donev, and Tobias, "Optimum Experimental Designs" */ /* */ /* NAME: Program 15.2 */ /* TITLE: Example 15.2 - Second-order designs with random effects */ /* KEYS: */ /* DATA: */ /* NOTES: Took ~8.5 minutes to run on server-class PC in May, 2007. */ /* */ /* Author: Randy Tobias */ /* History: */ /* Created...............................................27Jun2007 */ /**********************************************************************/ title1 "Example 15.2 - Second-order designs with random effects"; %macro Example15_2; %let Start = %sysfunc(datetime()); ods listing close; options nonotes; %let nRInc = 20; %let nCInc = 20; %let lEtaRLo = 2 ; %let lEtaRHi = -2 ; %let lEtaCLo = -1.5; %let lEtaCHi = 0.5; /* / Create treatment and block candidates. /---------------------------------------------------------------------*/ data Candidates; do x1 = -1 to 1; do x2 = -1 to 1; output; end; end; data RowCol; do Row = 1 to 2; do Col = 1 to 3; do i = 1 to 3; output; end; end; end; run; %let nDesign = 0; /* / Now, for each (EtaR,EtaC) pair, ... /---------------------------------------------------------------------*/ %do iEtaR = 1 %to &nRInc; %do iEtaC = 1 %to &nCInc; %let EtaR = %sysevalf(10**( &lEtaRLo + (&iEtaR - 1)/(&nRInc - 1)*(&lEtaRHi - &lEtaRLo)) ); %let EtaC = %sysevalf(10**( &lEtaCLo + (&iEtaC - 1)/(&nCInc - 1)*(&lEtaCHi - &lEtaCLo)) ); /* / ... create the plot variance matrix for this pair and store it in a / data set named Var, ... /---------------------------------------------------------------------*/ proc iml; use RowCol; read all var {Row Col}; ZR = design(Row); ZC = design(Col); s2r = &etaR; s2c = &etaC; s2e = 1; V = s2r * Zr*Zr` + s2c * Zc*Zc` + s2e*i(nrow(Zr)); create Var from V; append from V; /* / ... then run OPTEX/BLOCK COVAR= to find the optimal (x1,x2) design / for this plot variance matrix. Use ODS to grab the calculated / efficiencies, rounding them to get rid of numerical fuzz. /---------------------------------------------------------------------*/ proc optex data=Candidates coding=orthcan; model x1 x2 x1*x2 x1*x1 x2*x2; block covar=Var var=(col:) niter=1000 keep=1; ods output BlockDesignEfficiencies=Eff; output out=DesignB; data Eff; set Eff; Design = &nDesign+1; DCriterion = round(DCriterion,0.0001); ACriterion = round(ACriterion,0.0001); run; /* / Compute the efficiencies with the current (EtaR,EtaC) for all the / other best design designs found so far, too, ... /---------------------------------------------------------------------*/ %do iDesign = 1 %to &nDesign; proc optex data=Candidates coding=orthcan; model x1 x2 x1*x2 x1*x1 x2*x2; generate initdesign=Design&iDesign method=sequential; block covar=Var var=(col:) init=chain noexchange niter=0; ods output BlockDesignEfficiencies=iEff; data iEff; set iEff; Design = &iDesign; DCriterion = round(DCriterion,0.0001); ACriterion = round(ACriterion,0.0001); data Eff; set Eff iEff; run; %end; /* / ... and then find the best of all of them. /---------------------------------------------------------------------*/ proc sort data=Eff; by descending DCriterion descending ACriterion Design; data _null_; set Eff(obs=1); call symput('Best',trim(left(Design))); run; /* / If the best is the new one, save it. /---------------------------------------------------------------------*/ %if (&Best = &nDesign+1) %then %do; %let nDesign = &Best; data Design&nDesign; set DesignB; run; %end; %put 1: &iEtaR &iEtaC &Best; %end; %end; /* / Now, having found a bunch of designs each of which is optimal for / some (EtaR,EtaC), we run through the loop again but just evaluating / these designs, to make sure we've really got the right ones as / "best" for each pair. /---------------------------------------------------------------------*/ data All; if (0); run; %do iEtaR = 1 %to &nRInc; %do iEtaC = 1 %to &nCInc; %let EtaR = %sysevalf(10**( &lEtaRLo + (&iEtaR - 1)/(&nRInc - 1)*(&lEtaRHi - &lEtaRLo)) ); %let EtaC = %sysevalf(10**( &lEtaCLo + (&iEtaC - 1)/(&nCInc - 1)*(&lEtaCHi - &lEtaCLo)) ); proc iml; use RowCol; read all var {Row Col}; ZR = design(Row); ZC = design(Col); s2r = &etaR; s2c = &etaC; s2e = 1; V = s2r * Zr*Zr` + s2c * Zc*Zc` + s2e*i(nrow(Zr)); create Var from V; append from V; Z0 = j(nrow(Zr),1); P0 = Z0*ginv(Z0`*Z0)*Z0`; Pr = Zr*ginv(Zr`*Zr)*Zr`; Pr = Pr - P0; Pc = Zc*ginv(Zc`*Zc)*Zc`; Pc = Pc - P0; Pi = i(nrow(Zr)) - P0 - Pr - Pc; e0 = ((P0*V)#(1/P0))[:]; er = ((Pr*V)#(1/Pr))[:]; ec = ((Pc*V)#(1/Pc))[:]; ei = ((Pi*V)#(1/Pi))[:]; e = e0 || er || ec || ei; create E var {"e0" "er" "ec" "ei"}; append from e; data Eff; if (0); run; %do iDesign = 1 %to &nDesign; proc optex data=Candidates coding=orthcan; model x1 x2 x1*x2 x1*x1 x2*x2; generate initdesign=Design&iDesign method=sequential; block covar=Var var=(col:) init=chain noexchange niter=0; ods output BlockDesignEfficiencies=iEff; data iEff; set iEff; Design = &iDesign; DCriterion = round(DCriterion,0.0001); ACriterion = round(ACriterion,0.0001); data Eff; set Eff iEff; run; %end; proc sort data=Eff; by descending DCriterion descending ACriterion Design; data _null_; set Eff(obs=1); call symput('Best',trim(left(Design))); run; data This; set E; iEtaR = &iEtaR; EtaR = &EtaR; iEtaC = &iEtaC; EtaC = &EtaC; Best = &Best; data All; set All This; run; %put 2: &iEtaR &iEtaC &Best; %end; %end; ods listing; %do iDesign = 1 %to &nDesign; data d; merge RowCol Design&iDesign; run; title2 "Design &iDesign"; proc print data=d noobs; run; title2; %end; %let pagesize = %sysfunc(getoption(pagesize)); options ps=40; title2 "Regions of optimality"; data All; set All; iEtaR = -iEtaR; proc plot data=All; plot iEtaR*iEtaC = Best; run; title2; options ps=&pagesize; %let Stop = %sysfunc(datetime()); %put Time = %sysevalf(&Stop - &Start); options notes; %mend; %Example15_2; title1;