%MACRO SurvQtest(DATA1=, TIME1=,
                 DATA2=, TIME2=,
                 P=, Q=0.25, ALPHA = 0.05,
                 OUT = _QTEST, PRINTOP=1);

%LET errorflg = 0;
%if &DATA1 =  %then %do;
    %put ERROR - Data set <DATA1> not defined;
    %LET errorflg = 1;
    %end;
%if &DATA2 =  %then %do;
    %put ERROR - Data set <DATA2> not defined;
    %LET errorflg = 1;
    %end;
%if &Time1 =  %then %do;
    %put ERROR - Variable <Time1> not defined;
    %LET errorflg = 1;
    %end;
%if &Time2 =  %then %do;
    %put ERROR - Data set <Time2> not defined;
    %LET errorflg = 1;
    %end;
%if &p =  %then %do;
    %put ERROR - Probability <p> not defined;
    %LET errorflg = 1;
    %end;
%if &errorflg = 1 %then %do;
    data _NULL_;
    error 'ERROR - detected in the input data to the macro <SurvQtest>.';
    %goto exit;
    %end;

proc lifetest data=&data1 outsurv=out1(keep=x _censor_ survival) noprint;
time x*c(0);
run;
proc lifetest data=&data2 outsurv=out2(keep=x _censor_ survival) noprint;
time x*c(0);
run;
data out1;
set out1;
if x>0 & _censor_=0 ;
run;
data out2;
set out2;
if x>0 and _censor_=0;
run;
proc means min data=out1 noprint; 
var survival;
output out=out3(keep=mymin3) min=mymin3;
run;
proc means min data=out2 noprint; 
var survival;
output out=out4(keep=mymin4) min=mymin4;
run;
data mytmp(keep=m);
merge out3 out4;
m = 1-max(mymin3,mymin4);
run;
data mytmp;
set mytmp;
if m < &p  then do;
    put 'ERROR: Probability <P> is not appropriate. Try another one.';
   _error_=1;
end;
run;

proc sort data=&data1 out=out1;  
by &Time1;    
run;
proc sort data=&data2 out=out2; 
by &Time2;    
run;

proc iml;
    
use out1;
read all into x;

use out2;
read all into y;

n1 = nrow(x);
n2 = nrow(y);


xtmp1 = j(n1,2,0);
do i=1 to n1;
    xtmp1[i,1] = (x[i,2]=1);       
    xtmp1[i,2] = sum((x[,1]>=x[i,1]));
end;
xtmp2 = j(n1,1,0);
do i=1 to n1;
    xtmp2[i] = 1-xtmp1[i,1]/xtmp1[i,2];
end;
xtmp3 = j(n1,1,1);
do i=1 to n1;
   do k=1 to i;
       xtmp3[i] = xtmp3[i]*xtmp2[k];    
   end;
end;

xtmp4 = j(n1,1,1)-xtmp3;                                       
xtmp5 = j(n1,1,xtmp4[1]);
do i=2 to n1;
    xtmp5[i] = xtmp4[i]-xtmp4[i-1];                         
end;
    
index1 = min(loc((xtmp4>=&p)>0));
xxi = x[index1,1];                                    
xhats = xtmp3[index1];                                      

xden = (xtmp1[,2]-xtmp1[,1])#xtmp1[,2];
xtmp6 = j(n1,1,0);
do i=1 to n1;
    if xden[i]^=0 then xtmp6[i] = xtmp1[i,1]/xden[i];
end;
phi = (xhats**2)*sum((x[,1]<=xxi)#xtmp6);     


ytmp1 = j(n2,2,0);
do i=1 to n2;
    ytmp1[i,1] = (y[i,2]=1);
    ytmp1[i,2] = sum((y[,1]>=y[i,1]));
end;

ytmp2 = j(n2,1,0);
do i=1 to n2;
    ytmp2[i] = 1-ytmp1[i,1]/ytmp1[i,2];
end;
    
ytmp3 = j(n2,1,1);
do i=1 to n2;
   do k=1 to i;
       ytmp3[i] = ytmp3[i]*ytmp2[k];
   end;
end;

ytmp4 = j(n2,1,1)-ytmp3;
ytmp5 = j(n2,1,ytmp4[1]);
do i=2 to n2;
    ytmp5[i] = ytmp4[i]-ytmp4[i-1];
end;
    
index2 = min(loc((ytmp4>=&p)>0));
yxi = y[index2,1];
yhats = ytmp3[index2];

yden = (ytmp1[,2]-ytmp1[,1])#ytmp1[,2];
ytmp6 = j(n2,1,0);
do i=1 to n2;
    if yden[i]^=0 then ytmp6[i] = ytmp1[i,1]/yden[i];
end;
gamma = (yhats**2)*sum((y[,1]<=yxi)#ytmp6);                    

 
index3 = min(loc((xtmp4>=&q)>0));
xq_xi = x[index3,1];   
h1 = (n1**(-1/5))*4*xq_xi;
xcom = (((xxi - x[,1])/h1)+1 )#(-1 <= (xxi - x[,1])/h1 & (xxi - x[,1])/h1 <= 0)
      + (1-((xxi - x[,1])/h1))#( 0 < (xxi - x[,1])/h1 & (xxi - x[,1])/h1 <= 1 );

a1 = min(xxi+h1,max(x[,1]));
b1 = max(xxi,0);
c1 = min(xxi,max(x[,1]));
d1 = max(xxi-h1,0);
wf = (1/h1)*( (1+xxi/h1)*(a1-b1) - (a1**2-b1**2)/(2*h1) + (1-xxi/h1)*(c1-d1) + (c1**2-d1**2)/(2*h1) );
f = (1/wf)*(1/h1)*sum(xcom#xtmp5);


index4 = min(loc((ytmp4>=&q)>0));
yq_xi = y[index4,1];      
h2 = (n2**(-1/5))*4*yq_xi;           
ycom = (((yxi - y[,1])/h2)+1 )#(-1 <= (yxi - y[,1])/h2 & (yxi - y[,1])/h2 <= 0)
      + (1-((yxi - y[,1])/h2))#( 0 < (yxi - y[,1])/h2 & (yxi - y[,1])/h2 <= 1 );

a2 = min(yxi+h2,max(y[,1]));
b2 = max(yxi,0);
c2 = min(yxi,max(y[,1]));
d2 = max(yxi-h2,0);
wg = (1/h2)*( (1+yxi/h2)*(a2-b2) - (a2**2-b2**2)/(2*h2) + (1-yxi/h2)*(c2-d2) + (c2**2-d2**2)/(2*h2) );
g = (1/wg)*(1/h2)*sum(ycom#ytmp5);

beta = (xxi-yxi);
psi = phi/(f*f) + gamma/(g*g);
stat = (beta**2)/psi;

create mystat var{n1 n2 xxi yxi psi stat} ;
append;
close mystat;

quit;


data &out(keep=NUM_F NUM_G XI_F XI_G
               DIFF STD_ERR LOWER UPPER
               CHISQ PVALUE);
set mystat;
NUM_F = n1;
NUM_G = n2;    
XI_F = xxi;
XI_G = yxi;
DIFF = xxi - yxi;
z = probit(1-&alpha/2);  
STD_ERR = sqrt(psi);
LOWER = DIFF - z*sqrt(psi);
UPPER = DIFF + z*sqrt(psi);
CHISQ = stat;
PVALUE = 1-probchi(stat,1);
format xi_f xi_g diff std_err lower upper chisq pvalue 8.4;
run;


%if &printop = 0 %then %goto exit;

%if &printop = 1 %then %do;
proc print data=&out;
title "Testing for equality of &p quantile of two distributions";
run;
%end;

%exit:

%MEND SurvQtest;