data ryan;
infile 'd:\kurt\ryan\ryan.dat';
input v1-v3;
run;
proc iml;

start main;
print "--------------------------------------------------------------------";
print "Ryan and Joiner (1994) data";
print "---------------------------------------------------------";
use ryan;
read all var _num_ into xmatr;
x=xmatr;
n=nrow(x);
p=ncol(x);
nc=19;

print "Sample size n=" n;

*see the magnitude of the variables in x1 and x2;
a={1,1,0};
y1=x(|,1|); y2=x(|,2|); y3=(x*a);

**order the cases according to observed marginals;
do i=1 to n-1;
        do j=i+1 to n;
                if y1(|i|)>y1(|j|) then do;
                        tt=y1(|i|);
                        y1(|i|)=y1(|j|);
                        y1(|j|)=tt;
                end;
        end;
end;

do i=1 to n-1;
        do j=i+1 to n;
                if y2(|i|)>y2(|j|) then do;
                        tt=y2(|i|);
                        y2(|i|)=y2(|j|);
                        y2(|j|)=tt;
                end;
        end;
end;

do i=1 to n-1;
        do j=i+1 to n;
                if y3(|i|)>y3(|j|) then do;
                        tt=y3(|i|);
                        y3(|i|)=y3(|j|);
                        y3(|j|)=tt;
                end;
        end;
end;

*centralize data and abs;
*xm1=(y1(|14|)+y1(|15|))/2; 
xm1=j(1,n,1)*y1/n;
print "x1 mean=" xm1;
xabs1=abs(x(|,1|)-j(n,1,xm1)); print "xabs1=" xabs1;

*xm2=(y2(|14|)+y2(|15|))/2; 
xm2=j(1,n,1)*y2/n;
print "x2 mean=" xm2;
xabs2=abs(x(|,2|)-j(n,1,xm2)); print "xabs2=" xabs2;

*xm12=(y3(|14|)+y3(|15|))/2; 
xm12=j(1,n,1)*y3/n;
print "x1+x2 mean=" xm12;

xabs12=abs(x(|,1|)+x(|,2|)-j(n,1,xm12)); print "xabs12=" xabs12;

x13=xabs1||x(|,3|); 
x23=xabs2||x(|,3|); 
x123=xabs12||x(|,3|); 

*order the cases according to abs(x-mid(x));
do i=1 to n-1;
        do j=i+1 to n;
                if x13(|i,1|)>x13(|j,1|) then do;
                        tt=x13(|i,|);
                        x13(|i,|)=x13(|j,|);
                        x13(|j,|)=tt;
                end;
        end;
end;


do i=1 to n-1;
        do j=i+1 to n;
                if x23(|i,1|)>x23(|j,1|) then do;
                        tt=x23(|i,|);
                        x23(|i,|)=x23(|j,|);
                        x23(|j,|)=tt;
                end;
        end;
end;

do i=1 to n-1;
        do j=i+1 to n;
                if x123(|i,1|)>x123(|j,1|) then do;
                        tt=x123(|i,|);
                        x123(|i,|)=x123(|j,|);
                        x123(|j,|)=tt;
                end;
        end;
end;


xn13=(1:n)`||x13;
print "new order according to x1=";
print xn13;
w1=0; 
do i=1 to n;
if x13(|i,2|)=-99 then w1=w1+i;
end;

xn23=(1:n)`||x23;
print "new order according to x2=";
print xn23;
w2=0; 
do i=1 to n;
if x23(|i,2|)=-99 then w2=w1+i;
end;


xn123=(1:n)`||x123;
print "new order according to x1+x2=";
print xn123;
w3=0; 
do i=1 to n;
if x123(|i,2|)=-99 then w3=w3+i;
end;

print w1;
print w2;
print w3;

n=9; m=19;
ew=n*(n+m+1)/2;
varw=n*m*(n+m+1)/12;
z1=(w1-ew)/sqrt(varw);
z2=(w2-ew)/sqrt(varw);
z3=(w3-ew)/sqrt(varw);
print z1;
print z2;
print z3;

finish;
run main;