/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 16.5                                              */
/*   TITLE: Fig. 16.5. Example 16.4: second-order model for a         */
/*          three-component mixture with a qualitative factor at      */
/*          three levels                                              */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*   NOTES: Uses the ADX macros to create simplex lattice designs.    */
/*          Also addresses SAS Tasks 16.8 and 16.9.                   */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Second-order mixture designs with qualitative factors";


/*
/  Change this to order = 10 to explore SAS Task 16.9.  Note that if
/  you increase the order and thereby the number of candidates, you
/  probably also want to decrease the number of tries NITER for OPTEX,
/  to allow the procedure to run in a manageable amount of time,
/  especially for A-optimality.
/---------------------------------------------------------------------*/
%let order =    2;
%let niter = 1000;

/*
/  Create data sets holding the three designs in Figure 16.5.
/---------------------------------------------------------------------*/
data fig1605a; input A x1 x2 x3 @@; cards;
1 0   0   1      1 0   1   0      1 0.5 0.5 0 
2 1   0   0      2 0   0   1      2 0   0.5 0.5
3 1   0   0      3 0.5 0   0.5    3 0   0.5 0.5    3 0   1   0
;
data fig1605b; input A x1 x2 x3 @@; cards;
1 0   0   1      1 0   1   0
2 1   0   0      2 0   0   1      2 0   0.5 0.5    2 0.5 0.5 0
3 1   0   0      3 0.5 0   0.5    3 0   0.5 0.5    3 0   1   0
;
data fig1605c; input A x1 x2 x3 @@; cards;
1 0   0   1      1 0   1   0
2 1   0   0      2 0   0   1      2 0   0.5 0.5
3 1   0   0      3 0.5 0   0.5    3 0   0.5 0.5    3 0   1   0    3 0.5 0.5 0
;
%adxgen;
%adxmix;

/*
/  Create a candidate set of all vertices and centroids, within each
/  level of the qualitative factor.
/---------------------------------------------------------------------*/
%adxsld(Grid,x1 x2 x3,&order);
data Candidates; set Grid;
   do a = 1 to 3; output; end;
run;

/*
/  Compute OPTEX efficiencies for the D-optimum and A-optimum designs,
/  and also for the three designs depicted in Fig. 16.5.
/---------------------------------------------------------------------*/
ods listing close;
proc optex data=Candidates coding=orthcan;
   class a;
   model a x1|x2 x1*x1 x1*x2 x2*x2;
   id x3;
   generate n=10 method=m_fedorov niter=&niter keep=1 criterion=d;
   ods output Efficiencies=EffDOpt;
   output out=opDesign;
run;
   generate n=10 method=m_fedorov niter=&niter keep=1 criterion=a;
   ods output Efficiencies=EffAOpt;
run;
   generate initdesign=fig1605a method=sequential;
   ods output Efficiencies=Effa;
run;
   generate initdesign=fig1605b method=sequential;
   ods output Efficiencies=Effb;
run;
   generate initdesign=fig1605c method=sequential;
   ods output Efficiencies=Effc;
run;
ods listing;

/*
/  Collect efficiencies and round to get rid of numerical fuzz.
/---------------------------------------------------------------------*/
data EffDOpt; set EffDOpt; Source = "OPTEX D";
data EffAOpt; set EffAOpt; Source = "OPTEX A";
data Effa;    set Effa;    Source = "  (a)  ";
data Effb;    set Effb;    Source = "  (b)  ";
data Effc;    set Effc;    Source = "  (c)  ";
data Eff; set EffDOpt EffAOpt Effa Effb Effc;
   DCriterion = round(DCriterion,1e-8);
   ACriterion = round(ACriterion,1e-8);
   GCriterion = round(GCriterion,1e-8);
proc sort data=Eff;
   by descending DCriterion
      descending ACriterion
      descending GCriterion
      descending AvePredStdErr
      descending Source;
run;

/*
/  All designs are D- and G-equivalent, but designs (c) and (d) are
/  apparently A-optimum as well.
/---------------------------------------------------------------------*/
title2 "Efficiencies for designs in Figure 16.5";
proc print data=Eff noobs;
   var Source DCriterion ACriterion GCriterion AvePredStdErr;
run;
title2;

title2 "D-optimum design";
proc sort data=opDesign;
   by  a x1 x2 x3;
proc print data=opDesign;
   var a x1 x2 x3;
run;
title2;

/*
/  For order=2, consists of all candidate points with arbitrary
/  points replicated.  For order = 2*k where k > 1, doesn't have this
/  structure.
/---------------------------------------------------------------------*/
title2 "Replication factors for D-optimum design";
proc freq data=opDesign noprint;
   table x1*x2*x3 / list out=fDesign(rename=(COUNT=nRep) drop=PERCENT);
proc print data=fDesign noobs;
run;

title2;
title1;