/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 15.5                                              */
/*   TITLE: SAS Task 5.3 - Second-order designs with random effects   */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "SAS Task 15.3 - Second-order designs with random effects";

data Design1; input x1 x2 @@; cards;
-1 -1   1 -1  -1  1    0  0  -1  1   1  1    0 -1  -1  0   1  1
-1 -1   1  0   0  1   -1 -1   0  0   1 -1    1 -1  -1  1   1  1
;
data Design2; input x1 x2 @@; cards;
-1 -1   0  0   1  1    0 -1  -1  1   1  1   -1  1   0  0   1 -1
-1  1   0  1   1 -1   -1  0   0 -1   1  0   -1 -1   0  1   1  1
;
data Design3; input x1 x2 @@; cards;
-1  0   0 -1   1  1   -1  1   0  0   1 -1   -1 -1   0  1   1  0
-1  0   0  1   1 -1   -1 -1   0  0   1  1    0 -1  -1  1   1  0
;
data Design4; input x1 x2 @@; cards;
-1 -1   1 -1  -1  1    0  0  -1  1   1  1   -1 -1   0 -1   1  0
 0 -1  -1  0   1  1   -1 -1   0  0   1 -1   -1  1   0  1   1  0
;

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;

%macro DEval(EtaR,EtaC);
options nonotes;

   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;

   data Eff; if (0); run;
   %do iDesign = 1 %to 4;
      ods listing close;
      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;
      quit;
      ods listing;

      data iEff; set iEff;
         Design = &iDesign;
      data Eff; set Eff iEff;
      run;
      %end;

proc sort data=Eff; by descending DCriterion;
data _null_; set Eff(obs=1);
   call symput('DMax',trim(left(DCriterion)));
data Eff; set Eff;
   DCriterion = 100*(DCriterion/&DMax);
run;

title2 "D-efficiencies for (EtaR,EtaC) = (&EtaR,&EtaC)";
title3 "Relative to best exact design";
proc print data=Eff noobs;
   var Design DCriterion;
run;
title3;
title2;
%mend;

%DEval( 8.00, 0.05);
%DEval(20.00, 0.15);
%DEval(13.00, 1.20);
%DEval( 0.02, 0.30);