WilcxPowerContinuousNties.sas

/***********************************************************

* Program Name: WilcxPowerContinuousNTies.sas *

* Program dir : /udd/nhron/macrodevelopment/WilcxRank/ *

* Program Date: Dec., 2007 *

* *

* Note : This program is to estimate the power *

* for Continuous Random Variables with No *

* ties *

* Base on paper "Power and Sample Size Estimation for the *

* Wilcoxon Rank Sum Test" by B. Rosner and R. Glynn *

* Programmed by: Rong Chen and Marion McPhee *

***********************************************************/

%macro WilcxPowerContinuousNties(m =, /*sample size of one group */

n =, /*sample size of comparing group */

estimate=, /*y: use observed data to estimate muhat and theta */

dsn =, /*if estimate =y, then enter data set name for estimation */

var1 =, /* variable name for groups; users can specify any 2 distinct */

/* value for group number */

/* if there are more than 2 values, program will stop and post */

/* an error message */

var2 =, /* variable name for measurement */

alpha =, /*significance level */

theta =, /*user defined theta, if it is not estimated from data */

mu = /* effect size */);

/* if observe data set exist, check how many values that a group variable has */

%if &dsn^=%str() %then %do;

proc sql;

select count(distinct(&var1)) into: ngrp

from &dsn;

quit;

%if &ngrp>2 %then %goto out3; ;

%if &ngrp<2 %then %goto out4; ;

%end;

/* check if user enters theta, if not check to see if user enters mu */

%if &theta=%str() %then %do;

%if &mu^=%str() %then %do;

data _null_;

t=probnorm(&mu/sqrt(2));

call symput('theta', trim(left(t)));

run;

%end;

%else %if &mu=%str() %then %do;

%if %upcase(&estimate)^=Y %then %goto out1;

%end;

%if %upcase(&estimate)=Y %then %do;

%if &dsn=%str() %then %goto out2;

proc freq data=&dsn;

tables &var1/out=freqcnt;

run;

data _null_;

set freqcnt;

if _n_=1 then do;

call symput('grp1', trim(left(count)));

end;

else if _n_=2 then do;

call symput('grp2', trim(left(count)));

end;

run;

%put &grp1 &grp2;

proc rank data=&dsn out=rank;

var &var2;

ranks zrank;

run;

data probit;

set rank;

retain wstat 0;

if &var1=1 then wstat=wstat+zrank;

run;

data theta;

set probit end=eof;

if eof;

uformula=wstat-&grp1*(&grp1+1)/2;

ustat=&grp1*&grp2 - uformula;

theta=ustat/(&grp1*&grp2);

call symput('theta', trim(left(theta)));

run;

proc print data=theta;

var wstat uformula ustat theta;

run;

%put &theta;

%end;

%end;

data power;

z=probit(1-&alpha/2);

ptheta=probit(&theta);

vvar=(&theta*(1-&theta)+(&m+&n-2)*(probbnrm(ptheta,ptheta,0.5)-&theta**2))/(&m*&n);

a=abs(&theta-0.5)/sqrt(vvar)-z*sqrt(((&m+&n+1)/(12*&m*&n))/vvar);

b=-1*abs(&theta-0.5)/sqrt(vvar)-z*sqrt(((&m+&n+1)/(12*&m*&n))/vvar);

power1=probnorm(a);

power2=probnorm(b);

power=power1+power2;

m=&m;

n=&n;

theta=&theta;

output;

run;

title;

title "Power Estimation for Sample1 (&m) and Sample2 (&n) with theta=&theta";

proc print noobs;

var m n theta vvar a b z power1 power2 power;

run;

%goto out;

**** Error Message ***;

%out1: %err_msg(1);

%goto out;

%out2: %err_msg(2);

%goto out;

%out3: %err_msg(3);

%goto out;

%out4: %err_msg(4);

%goto out;

%out:

quit;

%mend WilcxPowerContinuousNties;

%macro err_msg(num);

%if &num=1 %then %do;

%put ****************ERROR IN SAS MACRO WilcxPowerContinuousNTies****************;

%put ****************Please define macro variables theta or mu *****************;

%end;

%if &num=2 %then %do;

%put ****************ERROR IN SAS MACRO WilcxPowerContinuousNTies****************;

%put ****************Please define data set name for estimating theta************;

%end;

%if &num=3 %then %do;

%put ****************ERROR IN SAS MACRO WilcxPowerContinuousNTies****************;

%put ****************Threre are more than 2 values in &var1 ****************;

%put ****************Please make sure there are 2 values in &var1*****************;

%end;

%if &num=4 %then %do;

%put ****************ERROR IN SAS MACRO WilcxPowerContinuousNTies****************;

%put ****************Threre is only 1 value in &var1 ****************;

%put ****************Please make sure there are 2 values in &var1*****************;

%end;

%mend err_msg;