WilcxPowerWties_GroupData.sas

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

* Program Name: WilcxPowerWties_GroupData.sas *

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

* Program Date: Dec., 2007 *

* *

* Note : This program is to estimate the power *

* for grouped continuous random variables *

* with 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 WilcxPowerWties_GroupData(dsn=, /* input data */

prefix1=, /* prefix of variable from group1 (reference group):x need to supply all the time */

prefix2=, /* prefix of variable from comparison group:y */

px =, /* probability of reference group */

py =, /* probability of comparison group */

grp =, /* # of groups */

m=, /* sample size of one group */

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

alpha=, /* significance level */

mu= /* effect size */);

/* if mu is given then use mu to compute theta and theta */

%if &mu ne %then %do;

data _null_;

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

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

run;

%end;

/* estimate theta_star from the observed data */

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

/* read in data for estimating hx*/

data temp;

set &dsn;

total =round(sum(of &prefix1.1-&prefix1&grp),1);

hx0=-10.0;

hx%eval(&grp+1)=10.0;

run;

data temp;

set temp;

/* get Hx */

array x{*} &prefix1.1 - &prefix1&grp;

array px{*} &px.1 - &px&grp;

array cdf{*} cdf1- cdf&grp;

array hx{*} hx0 - hx%eval(&grp+1);

do j=1 to %eval(&grp);

cdf{j}=0;

do l=1 to j;

cdf{j}=cdf{j}+px{l}/100;

end;

jj=j+1;

if cdf{j} >1 then cdf{j}=0.99;

hx{jj}=probit(abs(cdf{j}));

end;

/* use the user defiend mu */

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

suma=0;

sumb=0;

do j=1 to &grp;

jj=j+1;

a=hx(jj-1)-sqrt(2)*probit(&theta);

b=hx(jj)-sqrt(2)*probit(&theta);

suma=sum(suma,probnorm(a)*probnorm(hx(jj)));

sumb=sum(sumb,probnorm(b)*probnorm(hx(jj-1)));

end;

theta_star=(1.0-suma+sumb)/2;

theta=θ

/* in this situation mu_est=mu */

mu_est=μ

%let mu_est=μ

%end;

/* calculation of theta_star */

/* use observed data to estimate theta_star */

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

/* get u function */

array y{*} &prefix2.1 - &prefix2&grp;

array py{*} &py.1 - &py&grp;

array uformat{*} u1-u&grp;

sumx=0;

sumy=0;

do j=1 to &grp;

uformat{j}=0;

do k=1 to &grp;

if x{j} > y{k} then uformat{j}=uformat{j}+0;

else if x{j}= y{k} then uformat{j}=uformat{j}+0.5*(px{j}/100)*(py{k}/100);

else if x{j}< y{k} then uformat{j}=uformat{j}+(px{j}/100)*(py{k}/100);

end;

end;

theta_star=sum(of u1-u&grp);

/* make it an macro for estimating theta */

theta0=theta_star;

diff=0.1;

count=1;

do while(diff > 0.00001);

suma=0;

sumb=0;

do j=1 to &grp;

jj=j+1;

a=hx(jj-1)-sqrt(2)*probit(theta0);

b=hx(jj)-sqrt(2)*probit(theta0);

suma=suma+probnorm(a)*probnorm(hx(jj));

sumb=sumb+probnorm(b)*probnorm(hx(jj-1));

end;

ftheta0=1.0-suma+sumb-2*theta_star;

/*calculation of dftheta0 */

z=sqrt(2)*probit(theta0);

yt=sqrt(2)*exp(0.5*(probit(theta0)**2));

sumc=0;

sumd=0;

do j=1 to &grp;

jj=j+1;

c=exp((-0.5)*((hx(jj-1)-z)**2));

d=exp((-0.5)*((hx(jj)-z)**2));

sumc=sumc+c*probnorm(hx(jj));

sumd=sumd+d*probnorm(hx(jj-1));

end;

dftheta0=yt*(sumc-sumd);

theta=theta0-(ftheta0/dftheta0);

diff=abs(theta-theta0);

/* iteration */

if diff>0.00001 then do;

theta0=theta;

end;

count+1;

if count>20 then do;

%err_msg(2);

end;

end;

/* get muest from theta */

mu_est=sqrt(2)*probit(theta);

call symput('muest',trim(left(mu_est)));

%end;

/*compute var(Uik_star), cov(uik1_star,uik2_star), cov(ui1k_star, ui2k_star) */

euik_star=0;

vuik_star1=0;

euik1uik2=0;

covuik1uik2=0;

eui1kui2k=0;

covui1kui2k=0;

part1=0;

part2=0;

part3=0;

do j=1 to &grp;

jj=j+1;

euik_star1=probnorm(hx(jj))-probnorm(hx(jj-1));

euik_star2=1-probnorm(hx(jj)-mu_est)+0.25*(probnorm(hx(jj)-mu_est)-probnorm(hx(jj-1)-mu_est));

euik_star=euik_star+euik_star1*euik_star2;

vuik_star1=vuik_star1+(euik_star1)*(1-0.5*(probnorm(hx(jj)-mu_est)+probnorm(hx(jj-1)-mu_est)));

euik1uik2=euik1uik2+euik_star1*(1-0.5*(probnorm(hx(jj)-mu_est)+probnorm(hx(jj-1)-mu_est)))**2;

eui1kui2k=eui1kui2k+(probnorm(hx(jj)-mu_est)-probnorm(hx(jj-1)-mu_est))*((probnorm(hx(jj))+probnorm(hx(jj-1)))/2)**2;

part1=part1+(euik_star1)*(1-probnorm(hx(jj)));

part2=part2+(probnorm(hx(jj))-probnorm(hx(jj-1)))**2;

part3=part3+(euik_star1)*((probnorm(hx(jj))+ probnorm(hx(jj-1)))**2-1)/4;

end;

/* var(Uik) */

vuik_star=euik_star-vuik_star1**2;

/* cov(Uik1,Uik2) */

covuik1uik2=euik1uik2-vuik_star1**2;

/* cov(Ui1k, Ui2k) */

covui1kui2k=eui1kui2k-vuik_star1**2;

/* var1thta_star */

var1_thetastar=(vuik_star+(&n-1)*covuik1uik2+(&m-1)*covui1kui2k)/(&m*&n);

/* var0theta_star */

var0_thetastar=(part1+(part2-1)/4+(&m+&n-2)*part3)/(&m*&n);

power1a=abs(theta_star-0.5)/sqrt(var1_thetastar)-1.96*sqrt(var0_thetastar/var1_thetastar);

power1b=-1.0*abs(theta_star-0.5)/sqrt(var1_thetastar)-1.96*sqrt(var0_thetastar/var1_thetastar);

power=probnorm(power1a)+probnorm(power1b);;

output;

run;

title;

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

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

proc print noobs;

var theta_star theta var1_thetastar var0_thetastar power;

run;

%end;

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

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

proc print noobs;

var theta_star ftheta0 dftheta0 theta count var1_thetastar var0_thetastar power;

run;

%end;

%goto out;

**** Error Message ***;

%out1: %err_msg(1);

%goto out;

%out:

%mend WilcxPowerWties_GroupData;

%macro err_msg(num);

%if &num=1 %then %do;

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

%put ****************Please define macro variables observed data set name *****************;

%end;

%else %if &num=2 %then %do;

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

%put ****************program goes through 20 iterations data does not converge ****************;

%end;

%mend err_msg;