/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 9.3                                               */
/*   TITLE: Example 9.3 - Quadratic regression: Variance function     */
/*                        d(x,xi) for optimum exact designs           */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Computations for Figure 9.2";
title2 "Variance function d(x,xi) for optimum exact designs";


/*
/  First, try to use OPTEX to find the optimal exact D-optimum design.
/---------------------------------------------------------------------*/
%let N = 4;

data Candidates;
   do x = -1 to 1 by 0.01;
      output;
      end;
proc optex data=Candidates coding=orthcan;
   model x x*x;
   generate n=&N method=m_fedorov niter=1000 keep=10;
   output out=DOpDesign;
run;

/*
/  The D-best 4-run design will replicate any run of the D/G-optimum
/  3-run design {-1,0,1}.  OPTEX won't necessarily choose to replicate
/  the centerpoint ...
/---------------------------------------------------------------------*/
title3 "OPTEX-chosen D-optimum &N-point design";
proc freq data=DOpDesign;
   table x / nocum;
run;
title3;

/*
/  ... so we force it.
/---------------------------------------------------------------------*/
data DOpDesign;
   do x = -1 to 1;
      output;
      if (x = 0) then output;
      end;
run;
title3 "Another choice of D-optimum &N-point design";
proc freq data=DOpDesign;
   table x / nocum;
run;
title3;


/*
/  Give the points of the D-optimal design arbitrary response values
/  and append a grid of other points with missing response values.
/---------------------------------------------------------------------*/
data DOpDesign; set DOpDesign;
   y = rannor(1);
data Grid;
   do x = -1 to 1 by 0.01;
      output;
      end;
data Grid; set Grid DOpDesign;
run;

/*
/  Use GLM to compute the standard error for prediction over all
/  points and also save the root MSE.  Then compute the standardized
/  variance by dividing out the MSE from the prediction variance and 
/  scaling by the design size.
/---------------------------------------------------------------------*/
title3 "Computing standardized prediction variances indirectly";
proc glm data=Grid;
   model y = x x*x;
   output out=pGrid stdp=stdp;
   ods output FitStatistics=Fit;
data dD; set pGrid(where=(y = .));
   if (_N_ = 1) then set Fit;
   dD = &N*(stdp/RootMSE)**2;
   keep x dD;
run;
title3;






/*
/  Exact G-optimum designs are difficult to find, since they involve
/  concentric optimization: finding the points which minimize the
/  maximum standardized variance over the design region.  The
/  following approach comes quite close to the theoretical optimum,
/  but it is sensitive to the parameters of the optimization.
/---------------------------------------------------------------------*/
options nonotes;
proc iml;
   start d(x) global(V);             /* Define the standardized var   */
      Fx = j(nrow(x),1) || x || x#x; /* for given variance matrix as  */
      free d;                        /* a function of the prediction  */
      do i = 1 to nrow(x);           /* points.                       */
         d = d // Fx[i,]*V*Fx[i,]`;
         end;
      return(d);
   finish;

   %let xLo = -1;                    /* Define a function to maximize */
   %let xHi =  1;                    /* the prediction variance for a */
   %let Inc = 10;                    /* given design.                 */
   xGrid = &xLo + (&xHi - &xLo)*(0:&Inc)`/&Inc;
   start dmax(xDesign) global(V,xGrid);
      xd = shape(xDesign,nrow(xDesign)*ncol(xDesign),1);
      F = j(nrow(xd),1) || xd || xd#xd;
      V  = inv(F`*F);
      d = d(xGrid);
      call nlpqn(rc,xmax,"d",xGrid[d[<:>]],1,{-1,1});
      return(d(xmax));
   finish;

/*
/  Find the design which minimizes the maximum variance.
/---------------------------------------------------------------------*/
   xDesignInit = {-1 -0.3 0.3 1}`;
   con = shape(-1,1,4) // shape(1,1,4);
   call nlpqn(rc,xDesignOpt,"dmax",xDesignInit,0,con);
   print "Numerically-determined G-optimum",,
         xDesignOpt (dmax(xDesignOpt)) [colname={"Maximum d"}];

/*
/  Find the design which minimizes the maximum variance.
/---------------------------------------------------------------------*/
   a = sqrt(sqrt(5)-2);
   xGOpt = (-1 || -a || a || 1);
   print "Theoretical G-optimum",,
          xGOpt (dmax(xGOpt)) [colname={"Maximum d"}];

create GOpDesign var {"x"};
xGOpt = xGOpt`;
append from xGOpt;


/*
/  Give the points of the G-optimal design arbitrary response values
/  and append a grid of other points with missing response values.
/---------------------------------------------------------------------*/
data GOpDesign; set GOpDesign;
   y = rannor(1);
data Grid;
   do x = -1 to 1 by 0.01;
      output;
      end;
data Grid; set Grid GOpDesign;
run;

/*
/  Use GLM to compute the standard error for prediction over all
/  points and also save the root MSE.  Then compute the standardized
/  variance by dividing out the MSE from the prediction variance and 
/  scaling by the design size.
/---------------------------------------------------------------------*/
title3 "Computing standardized prediction variances indirectly";
proc glm data=Grid;
   model y = x x*x;
   output out=pGrid stdp=stdp;
   ods output FitStatistics=Fit;
data dG; set pGrid(where=(y = .));
   if (_N_ = 1) then set Fit;
   dG = &N*(stdp/RootMSE)**2;
   keep x dG;
run;
title3;

/*
/  Merge the results for the D- and G-optimum designs and print.
/---------------------------------------------------------------------*/
data d; merge dD dG;
proc print data=d;
run;

title2;
title1;