/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 19.6                                              */
/*   TITLE: SAS Task 19.1 - Augmenting a second-order design to third */
/*          order                                                     */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "SAS Task 19.1 - Augmenting a second-order design to third order";

title2 "Figure 19.1 - 13-trial D-optimum third-order design";

/*
/  Create 4x4 candidate grid over [-1,1]**2.
/---------------------------------------------------------------------*/
proc factex;
   factors x1-x2 / nlev=4;
   output out=Can3;
data Can3; set Can3;
   x1 = 2*(x1/3) - 1;
   x2 = 2*(x2/3) - 1;
   Source = "Can3  ";
run;

/*
/  Select 13 runs D-optimally.
/---------------------------------------------------------------------*/
proc optex data=Can3 seed=2 coding=orthcan noprint;
   model x1 x2 x1*x1 x1*x2 x2*x2 x1*x1*x1 x1*x1*x2 x1*x2*x2 x2*x2*x2;
   generate n=13 method=m_fedorov niter=1000 keep=10;
   output out=Order3;
run;


title2 "Augmenting a second-order design to a third-order design";

/*
/  The optimum 9-point, second-order design is the 3x3 grid on
/  [-1,1]**2.
/---------------------------------------------------------------------*/
proc factex;
   factors x1-x2 / nlev=3;
   output out=Order2;
data Order2; set Order2;
   Source = "Order2";
run;

/*
/  Use the AUGMENT= option in OPTEX to add points to Order2 for
/  fitting the third-order model.  Keep track of the D-efficiencies.
/  Although these are not the true D-efficiencies, as explained in
/  the text, they are proportional to the tru efficiencies, and thus
/  can be used to examine how the efficiency changes with design
/  size.
/---------------------------------------------------------------------*/
%macro Design(nMax);
   data AllEff;    if (0); run;
   data AllDesign; if (0); run;
   %do N=10 %to &nMax;
      proc optex data=Can3 coding=orthcan;
         model x1 x2 x1*x1 x1*x2 x2*x2 x1*x1*x1 x1*x1*x2 x1*x2*x2 x2*x2*x2;
         generate n=&N augment=Order2 method=m_fedorov niter=1000 keep=1;
         ods output Efficiencies=Eff;
         id Source;
         output out=Design(rename=(x1=x1_&N x2=x2_&N));
      data Eff; set Eff;
         N = &N;
         InfoPerRun = DCriterion / N;
      data AllEff;
	     set   AllEff    Eff;
      run;
      %if (&N = 10) %then %do;
         data AllDesign; set Design;
         run;
         %end;
      %else %do;
         data AllDesign; merge AllDesign Design; run;
         %end;
      %end;
%mend;
%Design(50);

/*
/  Print the efficiencies at different design sizes.
/---------------------------------------------------------------------*/
proc print data=AllEff noobs;
   var N DCriterion ACriterion GCriterion AvePredStdErr InfoPerRun;
run;