/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 15.1                                              */
/*   TITLE: Example 15.1 - Second-order block designs with variable   */
/*          block sizes                                               */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Example 15.1 - Second-order block designs with variable block sizes";

%macro BlockedDesign(n1,n2,n3,Replication);

options nonotes;

/*
/  Create candidate points (-1,0,1)**m and blocks with the given
/  configuration.
/---------------------------------------------------------------------*/
proc factex;
   factors x1 x2 / nlev=3;
   output out=Candidates;
data Blocks;
   Block=1; do i=1 to &n1; output; end;
   Block=2; do i=1 to &n2; output; end;
   Block=3; do i=1 to &n3; output; end;
run;


ods listing close;

/*
/  Find the D-optimal design with OPTEX and save the efficiencies.
/  If replication is disallowed, we use the INIT=CHAIN option to
/  force the candidates to be the initial design, and NOEXCHANGE to
/  force only interchange steps, making sure the final design is a
/  reordering of the 9 candidates.
/---------------------------------------------------------------------*/
%if ("&Replication" = "Yes") %then
   %let BlockOpt =;
%else
   %let BlockOpt = %quote(init=chain noexchange);
proc optex data=Candidates coding=none;
   model x1|x2@2 x1*x1 x2*x2;
   block design=Blocks &BlockOpt;
   class Block;
   model Block;
   output out=opDesign;
run;

ods listing;

/*
/  Compute D-criteria as given in the text and save to a data set.
/---------------------------------------------------------------------*/
proc iml;
   use opDesign;
   read all var ("Block");
   read all var {"x1" "x2"} into xop;
   Z = design(Block);
   F = xop || xop[,1]#xop[,1] || xop[,1]#xop[,2] || xop[,2]#xop[,2];
   X = Z || F;
   N = nrow(Z);
   detM  = det(X`*X/N);
   detMb = det(X`*X/N)/det(Z`*Z);
   create d var {"detM" "detMb"};
   data = detM || detMb;
   append;
data d; set d; Blocks = "&n1:&n2:&n3"; Replication = "&Replication";
run;
%mend;

/*
/  There are three possible ways to configure 9 runs in 3 blocks of
/  size at least 2, and we can either allow replication or force the
/  (-1,0,1)**2 candidates to be simply blocked. 
/---------------------------------------------------------------------*/
data Results; if (0); run;
%BlockedDesign(3,3,3,Yes); data Results; set Results d; run;
%BlockedDesign(2,2,5,Yes); data Results; set Results d; run;
%BlockedDesign(2,3,4,Yes); data Results; set Results d; run;

%BlockedDesign(3,3,3,No);  data Results; set Results d; run;
%BlockedDesign(2,2,5,No);  data Results; set Results d; run;
%BlockedDesign(2,3,4,No);  data Results; set Results d; run;

proc sort data=Results;
   by descending Replication descending detMb;
proc print data=Results noobs;
   by descending Replication;
   var Blocks detM detMb;
   format detM   8.7;
   format detMb 10.9;
run;

title1;