/**********************************************************************/ /* Atkinson, Donev, and Tobias, "Optimum Experimental Designs" */ /* */ /* NAME: Program 11.5 */ /* TITLE: Table 11.7 - Optimal quadratic designs */ /* KEYS: */ /* DATA: */ /* NOTES: The program uses OPTEX and compares the results to those */ /* given in the first edition of the text. */ /* */ /* Author: Randy Tobias */ /* History: */ /* Created...............................................27Jun2007 */ /**********************************************************************/ title1 "Table 11.7 - Optimal quadratic designs"; %macro Quadratic(m,N); options nonotes; /* / Create a macro variable that holds the quadratic model. /---------------------------------------------------------------------*/ /* / Linear effects and cross-products /------------------------------------------------------------------*/ %let model = x1; %do j = 2 %to &m; %let model = &model|x&j; %end; %let model = &model@2; /* / Quadratic effects /------------------------------------------------------------------*/ %do j = 1 %to &m; %let model = &model x&j*x&j; %end; /* / Create candidate points (-1,0,1)**m and add "standard order" IDs. /---------------------------------------------------------------------*/ proc factex; factors x1-x&m / nlev=3; output out=Candidates; proc sort data=Candidates; by x&m-x1; data Candidates; set Candidates; id = _N_; run; /* / Pull IDs for the points listed in Atkinson and Donev, Table 11.6, / and transform them into points in (-1,0,1)**m. /---------------------------------------------------------------------*/ data adDesign; set AD_11_6(where=((m=&m) & (N=&N))); array x{&m}; %do i=1 %to &m; x{&i} = mod(int((id-1)/(3**(&i-1))),3)-1; %end; run; ods listing close; /* / Find the D-optimal design with OPTEX and save the efficiencies. /---------------------------------------------------------------------*/ proc optex data=Candidates coding=orthcan; model &model; generate n=&N method=m_fedorov niter=10000 keep=1; output out=opDesign; id id; run; /* / Evaluate the designs and save the efficiencies. /---------------------------------------------------------------------*/ proc optex data=Candidates coding=none; model &model; generate method=sequential initdesign=adDesign; ods output Efficiencies=adEff; data adEff; set adEff; Source = "AD Table 11.6"; run; proc optex data=Candidates coding=none; model &model; generate method=sequential initdesign=opDesign; ods output Efficiencies=opEff; data opEff; set opEff; Source = "OPTEX"; run; ods listing; proc factex; factors x1-x&m / nlev=5; output out=Grid %do i = 1 %to &m; x&i nvals=(-1 to 1 by .5) %end; ; run; proc iml; start MakeZ(x); Z = j(nrow(x),1); do i = 1 to ncol(x); Z = Z || x[,i]; end; do i = 1 to ncol(x)-1; do j = i+1 to ncol(x); Z = Z || x[,i]#x[,j]; end; end; do i = 1 to ncol(x); Z = Z || x[,i]#x[,i]; end; return(Z); finish; use Candidates; read all var ("x1":"x&m") into x; Z = MakeZ(x); start DCritW(w) global(Z); Dw = diag((w##2)/sum(w##2)); return(log(det(Z`*Dw*Z)**(1/ncol(Z)))); finish; winit = j(nrow(x),1) / nrow(x); call nlpqn(rc,w,"DCritW",winit,1); call symput('DOpt',trim(left(char(exp(DCritW(w)),20)))); use Grid; read all var ("x1":"x&m") into xg; Zg = MakeZ(xg ); use adDesign; read all var ("x1":"x&m") into xad; Zad = MakeZ(xad); use opDesign; read all var ("x1":"x&m") into xop; Zop = MakeZ(xop); opV = inv(Zop`*Zop/nrow(Zop)); adV = inv(Zad`*Zad/nrow(Zad)); free dx; do i = 1 to nrow(Zg); zgi = Zg[i,]; dx = dx // (zgi*adV*zgi` || zgi*opV*zgi`); end; dAve = dx[:,]`; create dAve var ("dAve"); append from dAve; /* / Combine the two efficiencies and report. /---------------------------------------------------------------------*/ data Eff; set adEff opEff; m = &m; N = &N; GEff = (GCriterion/100)**2; p = 1 + 2*&m + &m*(&m-1)/2; DEff = exp(DCriterion/p - log(&N))/&DOpt; dAveC = &N*AvePredStdErr**2; data Eff; merge Eff dAve; proc print data=Eff noobs; var Source DEff GEff dAveC dAve; title1 "Comparing Atkinson & Donev Table 11.6 with OPTEX for m = &m, N = &N"; format DEff GEff 6.4; format dAveC dAve 5.2; run; title1; data Table11_8; set Table11_8 Eff; run; options notes; %mend; data AD_11_6; input m N @@; do i = 1 to N; input id @@; output; end; drop i; datalines; 2 6 1 3 4 7 8 9 2 7 1 3 4 7 8 9 5 2 8 1 3 4 7 8 9 5 6 3 10 27 19 3 8 16 21 1 23 15 11 3 14 7 3 19 21 9 25 1 5 15 11 27 13 17 23 3 15 27 19 3 8 16 21 1 23 15 11 9 25 4 7 2 3 16 27 19 3 8 16 21 1 23 15 11 9 25 4 7 2 6 3 18 27 19 3 8 16 21 1 23 15 11 9 25 4 7 2 6 21 19 3 20 27 19 3 8 16 21 1 23 15 11 9 25 4 7 2 6 21 19 18 26 4 15 54 7 38 78 22 6 57 73 58 79 21 26 1 63 18 4 18 10 66 9 41 55 25 21 79 30 34 46 81 60 4 62 74 2 27 4 24 72 36 54 3 55 21 25 73 1 18 40 38 61 59 75 80 7 19 57 6 79 78 8 27 4 25 5 75 69 49 21 57 1 38 79 62 55 9 25 54 34 19 3 73 70 26 63 18 7 81 58 4 27 5 75 69 49 21 57 1 38 79 62 55 9 25 54 34 19 3 73 70 26 63 18 7 81 58 24 10 5 21 189 21 81 1 46 217 57 241 169 61 225 25 9 165 64 181 237 157 204 125 77 5 26 55 157 126 70 223 235 209 14 19 187 219 171 63 7 27 3 102 243 163 137 75 52 193 186 80 174 5 27 169 114 216 44 229 218 21 237 73 61 165 81 57 4 225 25 9 185 181 241 91 108 190 2 180 155 55 ; data Table11_8; if (0); run; %Quadratic(2, 6); %Quadratic(2, 7); %Quadratic(2, 8); /* %Quadratic(2, 9); In 11.8, not 11.6 */ %Quadratic(3,10); %Quadratic(3,14); %Quadratic(3,15); %Quadratic(3,16); %Quadratic(3,18); %Quadratic(3,20); %Quadratic(4,15); %Quadratic(4,18); %Quadratic(4,24); %Quadratic(4,25); %Quadratic(4,27); %Quadratic(5,21); %Quadratic(5,26); %Quadratic(5,27); proc print data=Table11_8 noobs; title2 "OPTEX Designs"; where (Source = "OPTEX"); var m N DEff GEff dAve; format DEff GEff 6.4; format dAveC dAve 5.2; run; title1; proc print data=Table11_8 noobs; title2 "AD Table 11.6"; where (Source = "AD Table 11.6"); var m N DEff GEff dAve; format DEff GEff 6.4; format dAveC dAve 5.2; run; title2; title1;