%macro testpow;
%window start
#5 @4 "POWER CALCULATIONS FOR GLM CROSSOVER"
#8 @7 "DO YOU WANT PATIENT IN THE MODEL? (Y/N)" 
       @50 _pt_ 1 attr=underline
#10 @7 "DO YOU WANT PERIOD IN THE MODEL? (Y/N)"
        @50 _p_ 1 attr=underline
#12 @7 "DO YOU WANT CARRYOVER IN THE MODEL? (Y/N)"
        @50 _c_ 1 attr=underline
#14 @7 "STATE THE NUMBER OF PERIODS" @40 p 2 attr=underline
;                                          /* number of periods  */

%window getdes
#5 @4 "ENTER YOUR DESIGN SEQUENCE FOR SEQUENCE WITH TREATMENTS"
#6 @4 "CODED AS 1,2 ETC AND NUMBER OF REPLICATES LAST IN
       EACH ROW."
#10 @4 _in1 40 attr=underline
#11 @4 _in2 40 attr=underline
#12 @4 _in3 40 attr=underline
#13 @4 _in4 40 attr=underline
#14 @4 _in5 40 attr=underline
#15 @4 _in6 40 attr=underline
#16 @4 _in7 40 attr=underline
#17 @4 _in8 40 attr=underline
;

%window hypothes
#5 @4 "STATE YOUR HYPOTHESIS CONCERNING THE CONTRASTS IN TERMS"
#6 @4 "OF STD. (RETURN) WILL MEAN NO POWER CALCULATION."
#7 @4 hypo 40 attr=underline               /* hypothesis vector  */
#8 @4 "GIVE THE REQUIRED SIGNIFICANCE LEVEL OR JUST (RETURN)."
#9 @4 alpha 5 attr=underline               /* significance level */
;

%window dim
#5 @4 "DO YOU WANT A GRAPH OF POWER AS A FCN OF N (Y/N)."
#6 @4 GRAPH 1 attr=underline
#8 @4 "IF YOU WANT TO FIND N STATE THE REQUIRED POWER"
#9 @4 "GIVE THAT VALUE. NO VALUE WILL MEAN NO CALCULATION."
#10 @4 VALUE 5 attr=underline              /* required power     */
;

%window contr
#5 @4 "WHICH TYPE OF CONTRASTS DO YOU WANT TO WORK WITH?"
#7 @4 "A) RELATIVE A REFERENCE, I.E. 2-1, 3-1 ETC."
#9 @4 "B) CONSECUTIVE DIFFERENCES, I.E. 2-1,  3-2 ETC."
#11 @4 "C) CONTRASTS BETWEEN ALL PAIRS."
#12 @4 "ANSWER A, B OR C."
#13 @4 CTR 1 attr=underline
;

%put  _____________________________________________________;
%put | Experimental Release 2 of Macro TESTPOW            |;
%put | Author: P.Broberg Astra Arcus                      |;
%put | November 30 1995                                   |;
%put  ---------------------------------------------------- ;
title 'TEST OF CROSSOVER DESIGN';
OPTIONS  NODATE pageno=1;
%display start;
%display getdes;
%display contr;
%display hypothes;
%display dim;

%let ns=8;                  /* Maximum number of sequences.      */
%let lcount= ;
%do %until (&lcount ne);
      %if %bquote(&&_in&ns) ne %then
            %let lcount=&ns.;
      %else
            %let ns=%eval(&ns-1);
%end;

FOOTNOTE ' ';

data inmat;        /* Input the sequences, the sequence numbers, */
                   /* and the number of patients.                */
   keep per1-per&p repete seq;
array in $40 in1-in&ns;
array per per1-per&p;
   %do i=1 %to &ns;
           in&i="&&_in&i";
   %end;

do i=1 to &ns;
   do j=1 to &p;
          per[j]=input(scan(in[i],j),8.);
   end;
   repete=input(scan(in[i],&p+1),8.);
   seq=i;
   output;
end;
run;

proc print data=inmat;
run;

proc means data=inmat noprint sum;
   var repete;
   output out=npat sum=patients;
run;

data _null_;
   set npat;
   call symput('sekv',_freq_);
   call symput('n',patients);
run;

data patient;                         /* Create patient numbers. */
   do i=1 to &n;
      do j=1 to &p;
         pat=i;
         output;
      end;
   end;
keep pat;

data period;                          /* Create period numbers.  */
   do i=1 to &n;
      do j=1 to &p;
         per=j;
         output;
      end;
   end;
keep per;

data behandl;
   set inmat;                         /* Lay out the treatments. */
array peri {*} per1-per&p;
      do i=1 to repete;
         do j=1 to &p;
            treat=peri(j);
            output;
         end;
      end;
keep treat seq;

proc freq data=behandl noprint;
   table treat / out=ntreat;
run;

proc contents data=ntreat noprint out=conts(keep=nobs);
run;

data _null_;
   set conts;
   by nobs;

      /* b is the number of treatments.                          */
   if last.nobs then call symput('b',nobs);
run;                           

data carry;
   set period;
   set behandl;
   ca=lag(treat);
   if per=1 then co=&b+1;
   else co=ca;                 /* Add the carryover information. */
   output;
keep co;

data struc2;
   merge patient behandl period carry;
xi=rannor(0);
run;

%if %upcase(&_pt_) ne Y and %upcase(&_p_) ne Y and 
    %upcase(&_c_) ne Y %then %let model=treat;
%if %upcase(&_pt_) ne Y and %upcase(&_p_) ne Y and 
    %upcase(&_c_) eq Y %then %let model=treat co;
%if %upcase(&_pt_) ne Y and %upcase(&_p_) eq Y and 
    %upcase(&_c_) eq Y %then %let model=per treat co;
%if %upcase(&_pt_) eq Y and %upcase(&_p_) eq Y and 
    %upcase(&_c_) eq Y %then %do;
                                %let model=pat per treat co;
                             %end;
%if %upcase(&_pt_) eq Y and %upcase(&_p_) ne Y and 
    %upcase(&_c_) eq Y %then %do;
                                %let model=pat treat co;
                             %end;
%if %upcase(&_pt_) eq Y and %upcase(&_p_) ne Y and 
    %upcase(&_c_) ne Y %then %do;
                                %let model=pat treat;
                             %end;
%if %upcase(&_pt_) ne Y and %upcase(&_p_) eq Y and 
    %upcase(&_c_) ne Y %then %let model=per treat;
%if %upcase(&_pt_) eq Y and %upcase(&_p_) eq Y and 
    %upcase(&_c_) ne Y %then %do;
                                %let model=pat per treat;
                             %end;

   /* model will hold the model requested                        */
%put model Y=&model.%str(;);       
                                   
   /* The design matrix is output in the data.                   */
proc glmmod noprint outparm=parms outdesign=design;

      /* Set design.                                             */
   class &model. ;      
   model xi=&model.  ;
run;

proc contents data=design noprint out=out;
run;

proc means data=out noprint max;
   var varnum;
   output out=obsn max=nobs;
run;

data _null_;
   set obsn;
   call symput('dim',nobs);
run;

proc iml;
   reset noprint;
   use design;
   read all var {

   /* The macro columns creates a reference to the appropriate   */
   /* columns in the design matrix output by PROC GLMMOD.        */
%macro columns;
   %do i=1 %to %eval(&dim.-1);
       col&i.
   %end;
%mend;
%columns }

into x;     
xpx=t(x)*x;
r=1:%eval(&dim.-1);
z=sweep(xpx,r);
reset print;
print 'Contrast definition';

   /* The macro ell defines the matrix of contrast vectors,      */
   /* all against a reference, and has to fit the model          */
   /* statement used; it has to contain those factors included   */
   /* in the model.                                              */
%macro ell;
   %do i=1 %to %eval(&b.-1);   
   0
   %if %upcase(&_pt_) eq Y %then
      %do k=1 %to &n.;           
      0
      %end;
   %if %upcase(&_p_) eq Y %then
      %do k=1 %to &p.;
      0
      %end;
   -1
   %do j=1 %to %eval(&b.-1);
      %if &i.=&j. %then 1 ;
      %else 0 ;
   %end;
   %if %upcase(&_c_) eq Y %then
      %do k=1 %to %eval(&b.+1);
      0
      %end;
   %if &i. ne %eval(&b.-1) %then , ;
   %end;
%mend;

   /* The macro ell2 defines the matrix of contrast vectors,     */
   /* all consecutive contrasts.                                 */
%macro ell2;
   %do i=1 %to %eval(&b.-1);   
   0                           
   %if %upcase(&_pt_) eq Y %then
      %do k=1 %to &n.;
      0
      %end;
   %if %upcase(&_p_) eq Y %then
      %do k=1 %to &p.;
      0
      %end;
   %do j=1 %to %eval(&b.-1);
       %if &i.=&j. %then -1 1;
       %else 0 ;
   %end;
   %if %upcase(&_c_) eq Y %then
      %do k=1 %to %eval(&b.+1);
      0
      %end;
   %if &i. ne %eval(&b.-1) %then , ;
   %end;
%mend;

   /* The macro ell3 defines the matrix of contrast vectors, all */
   /* possible pairs.                                            */
%macro ell3;
   %do outer=1 %to  %eval(&b.-1);
      %do i=%eval(&outer.+1) %to &b.;
      0                               
      %if %upcase(&_pt_) eq Y %then
         %do k=1 %to &n.;
         0
         %end;
      %if %upcase(&_p_) eq Y %then
         %do k=1 %to &p.;
         0
         %end;
      %if %eval(&outer. gt 1) %then %do loop=2 %to %eval(&outer.);
      0 
      %end;
      -1
      %do j=%eval(&outer.+1) %to &b.;
         %if &i.=&j. %then 1 ;
         %else  0 ;
      %end;
      %if %upcase(&_c_) eq Y %then
         %do k=1 %to %eval(&b.+1);
         0
         %end;
      %if &i. ne &b. or &outer. ne %eval(&b.-1) %then , ;
      %end;
   %end;
%mend;

print 'Contrast vectors';
l={%if %upcase(&ctr.)=A %then %ell ;
   %else %if %upcase(&ctr.)=B %then %ell2 ;
   %else  %ell3; };

   /* l is the matrix defining the contrasts                     */
print 'Estimator covariance matrix';
cov=l*z*t(l);

   /* The standard error scales are found along the diagonal.    */
print 'Standard error scale';
scale=sqrt(diag(l*z*t(l)));   
                             
print 'The Average Standard Error';
average=trace(scale)/nrow(l) ;
print 'Used Estimates';
estimate=l*z*t(x);
%if %quote(&hypo) ne %then %do;    /* Calculation of matrix rank.*/
      rank1=round(trace(ginv(l)*l));
      rank2=round(trace(ginv(z)*z));
      edf=&n.*&p.-rank2;
      print 'Calculation of non-centrality';
      print 'given the following alternative';
      psi={&hypo.};
      delta2=psi*sweep(l*z*t(l))*t(psi);
      if edf<>0 then do;
      power=1-probf(finv(1-&alpha.,rank1,edf),
            rank1,edf,delta2);
   end;
%if %upcase(&_pt_.) eq Y %then step=&sekv.*(&p.-1);;
                                   /* Adjust increase in error   */
%if %upcase(&_pt_.) eq N %then step=&sekv.*&p.;;
                                   /* degrees of freedom as to   */
                                   /* whether a patient factor   */
                                   /* is present.                */
%if &value. ne %then %do;
   nn=&sekv.;m=1;d2=delta2;power=0;errordf=edf;
                                   /* Sample size determination. */
   lbl1: if power < &value. then do;
   m=m+1;
   print 'Iteration number' m;
   nn=nn+&sekv.;errordf=errordf+step;d2=delta2+d2;
   power=1-probf(finv(1-&alpha.,rank1,errordf),
         rank1,errordf,d2);
   end;
if power < &value. then goto lbl1;
print 'The required number of patients is'nn;
print power;
reset noprint;
%end;

%if %upcase(&graph.) eq Y %then %do;
reset noprint;
   n=0;
   power=0;


      /* At least 50 observations are created for the power graph. */
   create graf var {n power };         
                                     
      do i=2 to max(50,nn/%eval(&sekv.)+2);   
         n=&sekv.*i;d2=delta2*i;errdf=edf+step*(i-1);
         power=1-probf(finv(1-&alpha.,rank1,errdf),
            rank1,errdf,d2);
         append;
      end;
   close graf;
reset print;
%end;
run;
quit;

%if %upcase(&graph.) eq Y %then %do;
   symbol interpol=spline value=none;
   axis label=(h=3);

proc gplot data=graf;
   plot power*n;
run;

quit;
%end;
%end;
%mend testpow;