/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 6.4                                               */
/*   TITLE: Computations for Figures 6.7 & 6.8 - Standardized         */
/*          variances for two-factor designs                          */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Computations for Figures 6.7 & 6.8";
title2 "Standardized variances for two-factor designs";


/*
/  Create a grid of points over the unit square.
/---------------------------------------------------------------------*/
data Grid;
   do x1 = -1 to 1 by 0.05;
      do x2 = -1 to 1 by 0.05;
         output;
         end;
      end;
run;

/*
/  Create the designs---a 2*2 and a 3*3 factorial, respectively.
/---------------------------------------------------------------------*/
proc factex;
   factors x1-x2;
   output out=Design1;
proc factex;
   factors x1-x2 / nlev=3;
   output out=Design2;
run;

proc iml;

/*
/  Create the coefficients of the first-order model over the grid.
/---------------------------------------------------------------------*/
   use Grid; read all var ("x1":"x2") into grd;
   xx = j(nrow(grd),1) || grd;

/*
/  Create the variance matrix for the 2*2 design.
/---------------------------------------------------------------------*/
   use Design1; read all var ("x1":"x2") into des;
   n = nrow(des);
   F = j(n,1) || des;
   V = inv(F`*F);

/*
/  Compute standardized variance for each grid point.
/---------------------------------------------------------------------*/
   free d; do i = 1 to nrow(xx); d = d // n*xx[i,]*V*xx[i,]`; end;

/*
/  Put the standardized variances in a data set and print it.
/---------------------------------------------------------------------*/
   create StdVar1 var ("d");
   append from d;
data StdVar1; merge Grid StdVar1;
run;

title3 "2*2 design and first-order model (first ten obs)";
proc print data=StdVar1(obs=10);
run;
title3;

proc iml;

/*
/  Create the coefficients of the second-order model over the grid.
/---------------------------------------------------------------------*/
   use Grid;    read all var ("x1":"x2") into grd;
   xx = j(nrow(grd),1) || grd;
   do i = 1 to ncol(grd);               /* General code to add square */
      do j = i to ncol(grd);            /* and crossproduct terms     */
         xx = xx || grd[,i]#grd[,j];
         end;
      end;

/*
/  Create the variance matrix for the 3*3 design.
/---------------------------------------------------------------------*/
   use Design2; read all var ("x1":"x2") into des;
   n = nrow(des);
   F = j(n,1) || des;
   do i = 1 to ncol(des);
      do j = i to ncol(des);
         F = F || des[,i]#des[,j];
         end;
      end;
   V = inv(f`*f);

/*
/  Compute standardized variance for each grid point.
/---------------------------------------------------------------------*/
   free d; do i = 1 to nrow(xx); d = d // n*xx[i,]*V*xx[i,]`; end;

/*
/  Put the standardized variances in a data set and print it.
/---------------------------------------------------------------------*/
   create StdVar2 var ("d");
   append from d;
data StdVar2; merge Grid StdVar2;
run;
title3 "3*3 design and second-order model (first ten obs)";
proc print data=StdVar2(obs=10);
run;
title3;
title2;
title1;