FOCUS AREAS


Contents | SAS Program | PDF

   /*--------------------------------------------------------------------*/
   /* Copyright (c) 2014 SAS Institute Inc. Cary NC 27512  USA           */
   /*                                                                    */
   /*    NAME: RBDEVAL                                                   */
   /* PURPOSE: Evaluate treatment designs with nested blocks             */
   /* SUPPORT: Randy.Tobias@SAS.com          UPDATE: July 24, 2015       */
   /*                                                                    */
   /*--------------------------------------------------------------------*/
   /* INPUT:                                                             */
   /*    Design - name of data set with the design to be evaluated    .  */
   /*    vName  - name of treatment or variety variable                  */
   /*    rName  - name of replicate or first-level block variable        */
   /*    bName  - name of blocks-within-replicate or second-level        */
   /*             block variable                                         */
   /*    nCheck - non-nested check designs, 0 by default                 */
   /*                                                                    */
   /* ALGORITHM:                                                         */
   /* 1) Levels and joint levels are counted to determine whether the    */
   /*    block structure is balanced and orthogonal.                     */
   /* 2) OPTEX is used to compute the D- and A-efficiency of the design  */
   /*    with respect to first- and second-level blocking only,          */
   /*    respectively.  For balanced block structures, the block design  */
   /*    D-efficiency for each situation is also computed.               */
   /* 3) If a non-zero value for NCHECK= is specified, OPTEX is used to  */
   /*    find the best design for each of these two situations, with     */
   /*    ITER=&NCHECK, and the ratios of the design's efficiency         */
   /*    measures to those of these best designs are computed.           */
   /* 4) Finally, the resulting efficiencies are displayed---just the    */
   /*    the raw ones computed in Step (2) by default, or additionally   */
   /*    the efficiencies relative to OPTEX if NCHECK= is specified.     */
   /*                                                                    */
   /* Note that both the block design efficiency computed for balanced   */
   /* block structures in (2) and the relative efficiencies computed     */
   /* optionally in (3) get at the question of how good the design is    */
   /* with what you can truly achieve given the blocking structure. The  */
   /* block design efficiency is the ratio of the D criterion for the    */
   /* design to that of a balanced incomplete block design, if it        */
   /* exists.                                                            */
   /*--------------------------------------------------------------------*/
   
   %macro RBDEval(Design, vName, rName, bName, nCheck=0);
   %local _nopt; %let _nopt = %sysfunc(getoption(notes));
   options nonotes;
   
   ods select none;

   %if (^%sysfunc(exist(&Design))) %then %do;
      %put ERROR: Design data set &Design does not exist.;
      %goto exit;
      %end;
      
   %let NoVName = 0;
   %let NoRName = 0;
   %let NoBName = 0;
   data _null_;
      dsid = open("&Design");
      if (varnum(dsid,"&vName") = 0) then call symput("NoVName","1");
      if (varnum(dsid,"&rName") = 0) then call symput("NoRName","1");
      if (varnum(dsid,"&bName") = 0) then call symput("NoBName","1");
   run;
     
   %if (&NoVName) %then %do;
      %put ERROR: Design data set &Design does not contain variable &vName..;
      %goto exit;
      %end;
   %if (&NoRName) %then %do;
      %put ERROR: Design data set &Design does not contain variable &rName..;
      %goto exit;
      %end;
   %if (&NoBName) %then %do;
      %put ERROR: Design data set &Design does not contain variable &bName..;
      %goto exit;
      %end;
      

   /*
   /  Determine number of levels of variables and evaluate whether the
   /  block structure is balanced and orthogonal.
   /---------------------------------------------------------------------*/
      proc freq data=&Design noprint;
         table &vName        / out=_fV;
         table &rName        / out=_fR;
         table &rName*&bName / out=_fRB;
      run;
   
      proc summary data=_fV;
         var COUNT;
         output out=_vNum N=vNum;
      proc summary data=_fR;
         var COUNT;
         output out=_rNum N=rNum   std=rStd;
      proc summary data=_fRB;
         var COUNT;
         output out=_rbNum N=rbNum std=rbStd mean=rbAvg;
      run;
   
      data _null_; set &Design;
         call symput('N'   ,trim(left(_N_  )));
      data _null_; set _vNum;
         call symput('vNum',trim(left(vNum)));
      data _null_; set _rNum;
         call symput('rNum',trim(left(rNum)));
         call symput('rStd',trim(left(rStd)));
         call symput('rBal',trim(left(rStd < 1e-8)));
      data _null_; set _rbNum;
         call symput('rbNum',trim(left(rbNum)));
         call symput('rbStd',trim(left(rbStd)));
         call symput('rbBal',trim(left(rbStd < 1e-8)));
         call symput('rbAvg',trim(left(rbAvg)));
      run;
   
      %let BalBlock = %eval(&rBal & &rbBal);
   
      %if (&BalBlock) %then %do;
         %let pNum = &rbAvg;
         %let bNum = %eval(&rbNum / &rNum);
   
         %put Block structure: (&rNum * &bNum * &pNum);
         %end;
   
   /*
   /  Create the variety candidate set, the block design, and a candidate
   /  set for variety and replicate only,
   /---------------------------------------------------------------------*/
      data RBDTmtCan;
         do &vName = 1 to &vNum;
            output;
            end;
      data RBDBlkDesign; set &Design;
         keep &rName &bName;
      data RBDRepTmtCan;
         do &vName = 1 to &vNum;
         do &rName = 1 to &rNum;
            output;
            end; end;
      run;
   
      %let l1 = %length(&rName);
      %let l2 = %length(&bName);
      %let lEval = %eval(&l1 + &l2 + 2);
   
   
   /*
   /  Collect design's efficiency for 
   /     1) just treatments and replicates, ignoring blocks-within-
   /        replicates, and
   /     2) just treatments and blocks, ignoring replicates.
   /  If block structure is balanced, use the appropriate syntax for
   /  each, in order to get the BD efficiency; otherwise, just use the
   /  the block structure data set.
   /---------------------------------------------------------------------*/
      data Eff; if (0); run;
      %if (&BalBlock) %then %do;
         proc optex data=RBDTmtCan coding=orthcan;
            class &vName;
            model &vName;
            generate initdesign=&Design method=sequential;
            block structure=(&rNum)%eval(&bNum*&pNum) init=chain noexchange niter=0;
            ods output BlockDesignEfficiencies=iEff;
         quit;
   
         data iEff; set iEff;
            length Evaluation $ &lEval;
            Evaluation = "&rName";
         data Eff; set Eff iEff;
         run;
   
         proc optex data=RBDTmtCan coding=orthcan;
            class &vName;
            model &vName;
            generate initdesign=&Design method=sequential;
            block structure=(%eval(&rNum*&bNum))&pNum init=chain noexchange niter=0;
            ods output BlockDesignEfficiencies=iEff;
         quit;
   
         data iEff; set iEff;
            length Evaluation $ &lEval;
            Evaluation = "&bName(&rName)";
         data Eff; set Eff iEff;
         run;
   
         %let last = BCriterion;
         %let ASpace = 1;
         %end;
      %else %do;
         proc optex data=RBDTmtCan coding=orthcan;
            class &vName;
            model &vName;
            generate initdesign=&Design method=sequential;
            block design=RBDBlkDesign init=chain noexchange niter=0;
            class &rName;
            model &rName;
            ods output BlockDesignEfficiencies=iEff;
         quit;
   
         data iEff; set iEff;
            length Evaluation $ &lEval;
            Evaluation = "&rName";
         data Eff; set Eff iEff;
         run;
   
         proc optex data=RBDTmtCan coding=orthcan;
            class &vName;
            model &vName;
            generate initdesign=&Design method=sequential;
            block design=RBDBlkDesign init=chain noexchange niter=0;
            class &rName &bName;
            model &rName &bName &rName*&bName;
            ods output BlockDesignEfficiencies=iEff;
         quit;
   
         data iEff; set iEff;
            length Evaluation $ &lEval;
            Evaluation = "&bName(&rName)";
         data Eff; set Eff iEff;
         run;
   
         %let last = ACriterion;
         %let ASpace = 4;
         %end;
   
   /*
   /  If optimal checking is requested, go ahead and find the designs
   /  that are optimal with respect to the replicate and blocks-within-
   /  replicates models used previously for evaluation, and compute the
   /  given designs efficiencies relative to those.
   /---------------------------------------------------------------------*/
      %if (&nCheck) %then %do;
         data Check; if (0); run;
   
         %if (&BalBlock) %then %do;
            proc optex data=RBDTmtCan coding=orthcan;
               class &vName;
               model &vName;
               block structure=(&rNum)%eval(&bNum*&pNum) niter=&nCheck keep=1;
               ods output BlockDesignEfficiencies=iCheck;
            quit;
            data Check; set Check iCheck;
            run;
   
            proc optex data=RBDTmtCan coding=orthcan;
               class &vName;
               model &vName;
               generate initdesign=&Design method=sequential;
               block structure=(%eval(&rNum*&bNum))&pNum niter=&nCheck keep=1;
               ods output BlockDesignEfficiencies=iCheck;
            quit;
            data Check; set Check iCheck;
            run;
   
            data Eff; merge Eff Check(rename=(DCriterion=CheckDCriterion
                                              ACriterion=CheckACriterion
                                              BCriterion=CheckBCriterion));
               rDCriterion = 100*DCriterion/CheckDCriterion;
               rACriterion = 100*ACriterion/CheckACriterion;
               rBCriterion = 100*BCriterion/CheckBCriterion;
               format rDCriterion rACriterion rBCriterion 8.4;
            run;
            %end;
         %else %do;
            proc optex data=RBDTmtCan coding=orthcan;
               class &vName;
               model &vName;
               block design=RBDBlkDesign niter=&nCheck keep=1;
               class &rName;
               model &rName;
               ods output BlockDesignEfficiencies=iCheck;
            quit;
            data Check; set Check iCheck;
            run;
   
            proc optex data=RBDTmtCan coding=orthcan;
               class &vName;
               model &vName;
               block design=RBDBlkDesign niter=&nCheck keep=1;
               class &rName &bName;
               model &rName &bName &rName*&bName;
               ods output BlockDesignEfficiencies=iCheck;
            quit;
            data Check; set Check iCheck;
            run;
   
            data Eff; merge Eff Check(rename=(DCriterion=CheckDCriterion
                                              ACriterion=CheckACriterion));
               rDCriterion = 100*DCriterion/CheckDCriterion;
               rACriterion = 100*ACriterion/CheckACriterion;
               format rDCriterion rACriterion 8.4;
            run;
            %end;
   
         %end;
   
   ods select all;
   
   proc template;
      define table RBDEval;
   
         column Evaluation  DCriterion  ACriterion  BCriterion
                           rDCriterion rACriterion rBCriterion;
   
         define header h;
            text "&vname Efficiency for &rName and &bName-within-&rName";
            space=1;
         end;
   
         define Evaluation;
            id header="Evaluation" style=RowHeader;
         end;
   
         define DCriterion;
            header="D-Efficiency" format=8.4 space=1;
         end;
   
         define ACriterion;
            header="A-Efficiency" format=8.4 space=&ASpace;
         end;
   
         define BCriterion;
            header="D Relative to Upper Bound" format=8.4 space=4;
         end;
   
         %if (&nCheck) %then %do;
            define header hRaw;
               text "Raw";
               start=DCriterion end=&Last expand='-';
            end;
        %end;
   
         define rDCriterion;
            header="D-Efficiency" format=8.4 space=1;
         end;
   
         define rACriterion;
            header="A-Efficiency" format=8.4 space=&ASpace;
         end;
   
         define rBCriterion;
            header="D Relative to Upper Bound" format=8.4;
         end;
   
         define header hRel;
            text "Relative";
            start=rDCriterion end=r&Last expand='-';
         end;
   
      end;
   run;
   
   proc sgrender data=Eff template="RBDEval";
   run;

%exit:
   options &_nopt;
   %mend;

data Wines;
   do Wine = 1 to 50;
      output;
   end;
run;

data Setup;
   do Subject = 1 to 20;
      do Session = 1 to  3;
         if (Session < 3) then nPlot = 17;
         else                  nPlot = 16;
         do Plot = 1 to nPlot;
            output;
         end;
      end;
   end;
   drop nPlot;
run;

proc optex data=Wines coding=orthcan seed=16899;
   class Wine;
   model Wine;
   block design=Setup;
   class Subject Session;
   model Subject, Session(Subject) / prior=0,10;
   output out=Design;
run;

%RBDEval(Design,Wine,Subject,Session,nCheck=10);

data Candidates;
   do Treatment = 1 to 4;
      output;
   end;
run;

data Design;
   do Replicate = 1 to 3;
      do Block = 1 to 2;
         do Plot = 1 to 2;
            input Treatment @@;
            output;
         end;
      end;
   end;
datalines;
1 2   3 4
1 3   2 4
1 4   2 3
;

%RBDEval(Design,Treatment,Replicate,Block,nCheck=10);


data Candidates;
   do Treatment = 1 to 15;
      output;
   end;
run;

data BlockStructure;
   do Replicate = 1 to 7;
      do Block = 1 to 5;
         do Plot = 1 to 3;
            output;
         end;
      end;
   end;
run;

proc optex data=Candidates coding=orthcan seed=492069001;
   class Treatment;
   model Treatment;
   blocks design=BlockStructure;
   class Replicate Block;
   model Replicate Block(Replicate);
   output out=Design;
   ods select BlockDesignEfficiencies;
run;


%RBDEval(Design,Treatment,Replicate,Block,nCheck=10);

proc freq data=Design;
   table Treatment*Replicate / norow nocol nopct nocum;
run;

proc optex data=Candidates coding=orthcan seed=607104001;
   class Treatment;
   model Treatment;
   blocks design=BlockStructure niter=10000 keep=10;
   class Replicate Block;
   model Replicate, Block(Replicate) / prior=0,100;
   output out=Design;
   ods select BlockDesignEfficiencies;
run;

%RBDEval(Design,Treatment,Replicate,Block,nCheck=10);

data Candidates;
   do Treatment = 1 to 4;
      output;
   end;
run;

data Design;
   do Replicate = 1 to 3;
      do Block = 1 to 2;
         do Plot = 1 to 2;
            input Treatment @@;
            output;
         end;
      end;
   end;
datalines;
1 2   3 4
1 3   2 4
1 4   2 3
;

%RBDEval(Design,Treatment,Replicate,Block);

%RBDEval(Design,Treatment,Replicate,Block,nCheck=10);