%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;