/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 9.2                                               */
/*   TITLE: Example 9.2 - Quadratic regression: Variance function     */
/*                        d(x,xi) for optimum continuous design       */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*   NOTES: This program performs the calculation two different ways  */
/*          and compares the results: (1) directly with matrix code   */
/*          and (2) indirectly with SAS design and analysis PROCs.    */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Computations for Figure 9.1";
title2 "Variance function d(x,xi) for optimum continuous design";

/*
/  The direct approach.
/---------------------------------------------------------------------*/
proc iml;
   xOp = {-1 0 1}`;                        /* Define optimal design   */
   wOp = { 1 1 1}`/3;
   FOp = j(nrow(xOp),1) || xOp || xOp#xOp; /* and corresponding       */
   V   = inv(FOp`*diag(wOp)*FOp);          /* variance matrix.        */

   x = -1 + (1 - -1)*(0:200)`/200;         /* Prediction grid         */
   F = j(nrow(x),1) || x || x#x;
   free d;
   do i = 1 to nrow(x);                    /* Compute d = f`*V*f for  */
      d = d // F[i,]*V*F[i,]`;             /* each grid point.        */
      end;

   data = x || d;                          /* Save in a data set.     */
   create dMethod1 var {"x" "d"};
   append from data;


/*
/  The indirect approach: First, use OPTEX to find the optimal design.
/---------------------------------------------------------------------*/
%let N = 6;
title3 "Optimal &N-point design";

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=opDesign;
run;

proc freq data=OpDesign;
   table x / nocum;
run;

title3;


/*
/  Give the points of the optimal design arbitrary response values
/  and append a grid of other points with missing response values.
/---------------------------------------------------------------------*/
data opDesign; set opDesign;
   y = rannor(1);
data Grid;
   do x = -1 to 1 by 0.01;
      output;
      end;
data Grid; set Grid opDesign;
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 dMethod2; set pGrid(where=(y = .));
   if (_N_ = 1) then set Fit;
   d = &N*(stdp/RootMSE)**2;
   keep x d;
run;
title3;

title3 "Comparing different methods";
proc compare data=dMethod1 compare=dMethod2 method=relative brief;
   var d;
   id  x;
run;
title3;
title2;
title1;