/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 15.6                                              */
/*   TITLE: Section 15.5.5 - Db-optimum design for clinical trial     */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*   NOTES: Took ~2.5 minutes to run on server-class PC in May, 2007. */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Section 15.5.5 - Db-optimum design for clinical trial";

%macro temp; options nonotes; ods listing close;

%let Start = %sysfunc(datetime());

   %let Lower =  0;
   %let Upper = 10;

   %let nDesign = 0;

   data AllBlockSizes; if (0);
   %do  n7=&Lower %to 30;                       %if (&n7<=&Upper) %then %do;
   %do  n6=&n7    %to 30- &n7;                  %if (&n6<=&Upper) %then %do;
   %do  n5=&n6    %to 30-(&n7+&n6);             %if (&n5<=&Upper) %then %do;
   %do  n4=&n5    %to 30-(&n7+&n6+&n5);         %if (&n4<=&Upper) %then %do;
   %do  n3=&n4    %to 30-(&n7+&n6+&n5+&n4);     %if (&n3<=&Upper) %then %do;
   %do  n2=&n3    %to 30-(&n7+&n6+&n5+&n4+&n3); %if (&n2<=&Upper) %then %do;
   %let n1=  %eval(30-(&n7+&n6+&n5+&n4+&n3+&n2));
   %if ((&n2<=&n1) & (&n1<=&Upper)) %then %do;
      %put &n1 &n2 &n3 &n4 &n5 &n6 &n7;
   /*
   /  Create candidates for quantitative factors.
   /---------------------------------------------------------*/
   data Candidates;
      do Dose = -1 to 1;
         do Admin = -1 to 1;
            output;
            end;
         end;

   /*
   /  Create blocks.
   /---------------------------------------------------------*/
   data Hospitals;
      Hospital = 1; do i = 1 to &n1; output; end;
      Hospital = 2; do i = 1 to &n2; output; end;
      Hospital = 3; do i = 1 to &n3; output; end;
      Hospital = 4; do i = 1 to &n4; output; end;
      Hospital = 5; do i = 1 to &n5; output; end;
      Hospital = 6; do i = 1 to &n6; output; end;
      Hospital = 7; do i = 1 to &n7; output; end;

   /*
   /  Select quantitative runs D_beta optimal for a response
   /  surface model with the given blocks.
   /---------------------------------------------------------*/
   proc optex data=Candidates coding=orthcan;
      class Admin;
      model Dose|Admin Dose*Dose Admin*Admin;
      block design=Hospitals niter=100 keep=1;
      class Hospital;
      model Hospital;
      ods output BlockDesignEfficiencies=BDEff;
   run;
      data BDeff; set BDEff;
         n1 = &n1; n2 = &n2; n3 = &n3; n4 = &n4;
         n5 = &n5; n6 = &n6; n7 = &n7;
      data AllBlockSizes; set AllBlockSizes BDEff;
      run;

      %let nDesign = %eval(&nDesign+1);
/*
      %if ((&nDesign/10)*10 = &nDesign) %then %do;
         ods listing;
         proc sort data=AllBlockSizes; 
            by descending DCriterion;
         data BDEff; set AllBlockSizes(obs=1) BDEff; drop DesignNumber;
         proc print data=BDEff;
         run;
         ods listing close;
         %end;
*/
      %end;
      %end; %end;
      %end; %end;
      %end; %end;
      %end; %end;
      %end; %end;
      %end; %end;
ods listing; options notes;

%let Stop = %sysfunc(datetime());
%put Time = %sysevalf(&Stop - &Start);

%mend;
%temp;
proc sort data=AllBlockSizes; 
   by descending DCriterion;
run;

title2 "Best 10 designs";
proc print data=AllBlockSizes(obs=10);
run;
title2;

title1;