/*                              SAS Program

/*                              PCIC_SAS Program
Paired Comparisons Information Criteria
	The program implements the information theoretic approach to paired comparisons among sample means presented in the following references:
Dayton, C. M. (2003) Information criteria for pairwise comparisons. Psychological Methods, 8, 61-71.
Dayton, C. M. (2001) SUBSET: Best subsets using information criteria. Journal of Statistical Software, Vol. 06, Issue 02, April.
Dayton, C.M. (1998) "Information Criteria for the Paired-Comparisons Problem." American Statistician, 52, 144-151.

Akaike AIC and Schwarz BIC statistics are computed for all ordered subsets of independent means for both homogeneous and heterogeneous models. For heterogeneous cases, variances are pooled in the same pattern as means. Data can be from a SAS dataset or can be imported from an Excel spreadsheet or other database program. The user must input or import the data and identify the grouping and dependent variables in a manner
comparable to the following examples: */

/*Exemplary setup using SAS dataset
/* data 'c:\sasdata\nels.sas7bdat';
/* set nels1;                             /* Input data from a SAS data set*/ /*
/* if race>5 then delete;                 /* Race is grouping variable*/ /*
/* if MATH>=99.98 then delete;            /* Delete missing observations*/ */

/*Exemplary setup using Excel dataset 
/* Proc import datafile='c:\sas_pcic\ufcollege1.xls' out =nels3 ;
/* data nels1;
/* set nels3;                            /* Input data from an external Excel  file */ /*
/* if race>5 then delete;                /* Race is grouping variable */ /*
/* if MATH>=99.98 then delete;           /* Delete missing observations */

/* The example below assumes that raw data are imported from an Excel spreadsheet named NELSMATH.xls in the indicated sub-directory */
/* When running the example, substitute the name of your own sub-directory */

%let prnum = 10 ;   /* Set number of smallest AIC or BIC values to print */
%let k=5;           /* Set is number of groups =k  */

Proc import datafile='C:\Documents and Settings\Administrator\My Documents\SAS_PCIC\NELSMath.xls' out =nels3 ;
data nels1;
set nels3;                            /* Input data from an external Excel file */
if race>5 then delete;                /* Race is grouping variable */
if MATH>=99.98 then delete;           /* Delete missing observations */
run;
proc sort data=nels1;
by race;
proc means noprint data=nels1;
output out=stats mean=mean stddev=sd N=n;  
var math;                                /* Math is dependent variable for analysis*/
by race;                                 /* Produce mean, standard deviation, sample sizes groups */
run;  
/*beginning of analysis*/
data data;
set stats;
varunb=sd**2;
varmle=varunb*(n-1)/n;
sum=mean*n;
ss=varmle*n;
proc print data=data;
var race _freq_ mean sd n varunb varmle sum ss;
sum n;
sum ss;
title "Summary Table - group means in original order";
run;
proc sort data=data out= matrix1;     /* Rank order groups. Partitionings are based upon the ranked means and not the means in original order */
by mean;
data matrix;
set matrix1;
proc print;
var race _freq_ mean sd n varunb varmle sum ss;
title "Summary Table - group means in rank order";
run;
title"AIC and BIC for Homogeneous Case";
proc iml;                              /* Create  matrix index that is used to obtain the  
                                                patterns of subgroups */
b=2**(&k-1);              /* Transform from decimal base to binary digits and store in matrix X */
x= j(b,&k,0);
do m=0 to 2**(&k-1)-1;
	y=m;
	x[m+1,1]=1;
	do l=2 to &k;
		x[m+1,l]=floor(y/2**(&k-l));
		y=y-x[m+1,l]*2**(&k-l);
	end;
end;
n=j(2**(&k-1),&k-1,0);                        
 do i=1 to 2**(&k-1);
  do j=1 to &k-1;
if x[i,j]^=x[i,j+1] then n[i,j]=1;           /*Calculate degrees of freedom for
                                               homogeneous variance case*/
	   c=2*(n[,+]+2);
		end;
	 end;
homoaicmatrix=j(2**(&k-1),1,0);
homobicmatrix=j(2**(&k-1),1,0);                
                                                      /* Compute AIC and BIC for Homogeneous case */
do i=1 to b;
a=j(&k);                                              /* a is set of matrices for different 
                                                                 patterns of subgroups for homogeneous case*/
   do c=1 to &k;
     exit=&k-c;
       do r=1 to exit;
if x[i,c]^= x[i,c+r] then do;                   /* First set all elements of matrix to 0
                                                 	     for every row if x, if the previous
                                                                 number is different from the one after, 
                                                                 keep the element zero. If same, set 1 and 
                                                                  keep checking the one after it. */ 
           		   do j=c+r to &k;
		     a[c,j]=0;
			end;
			r=exit;
		   end;
			 end;
			 end;
			   do p=1 to &k;
			      do q=1 to p-1;
				     a[p,q]=a[q,p];                /* Transform a from upper 								triangular to symmetric matrix */
					 end;
					 end;
			 Use matrix;
			 read ALL var{n}into mm;
			 read all var{varmle} into nn;
			 read all var{sum} into dd;
			 read all var{mean} into mean1;
			 kk=a*mm;
			 ll=a*dd;
			 o=ll#(1/kk);
			 pp=(o-mean1)##2;
			 rr=sum(mm);
			 qq=(sum((mm#nn)))/rr;
			 bb=sum(mm#nn);
			 ss=qq*rr;
			 tt=(sum(mm#pp))+ss;
			uu=((-1/2)*rr)*((log(2*3.1415926))+log(qq))-(tt/(2*qq));
			vv=(-2)*uu;
			AIC=vv+2*(n[i,+]+2);                   /*Compute AIC for 										homogeneous case*/
			homoaicmatrix[i,1]=AIC;
		    BIC=vv+log(mm[+,])*(n[i,+]+2);         /*Compute BIC for 										homogeneous case*/                   
			homobicmatrix[i,1]=BIC;
			 end;
			 r1=rank(homoaicmatrix);               /*Rank AIC and BIC*/
			 r2=rank(homobicmatrix);
	         group=j(2**(&k-1),&k,1);
do i=1 to 2**(&k-1);
  do j=1 to &k-1;
       if x[i,j]^=x[i,j+1] then group[i,j+1]=group[i,j]+1;  /* Create the ordered subsets*/
	   else group[i,j+1]=group[i,j];
	  		end;
	 end;
	 aic_homog1=r1||homoaicmatrix||group;     /*Print output for homogeneous case*/
	 bic_homog1=r2||homobicmatrix||group;
	 if &k<=3 then do;
	 print "Rank of AIC, value of AIC and ordered subsets for homogeneous variance case:",aic_homog1;
	 print "Rank of BIC, value of BIC and ordered subsets for homogeneous variance case:",bic_homog1;
     end;
     else 
     do;
     aic_homog=j(&prnum,&k+2,0);     /* Calculate the smallest prnum AIC value and 
					associated ordered subsets for homogeneous case */ 
     do i=1 to 2**(&k-1);
	   do j=1 to &k+2;
	     r= aic_homog1[i,1];
	     if r<=&prnum then aic_homog[r,j]=aic_homog1[i,j];
		 end;
		 end;
     bic_homog=j(&prnum,&k+2,0);     /*calculate the smallest prnum BIC value and 
	 			  associated ordered subsets for homogeneous case */ 
     do i=1 to 2**(&k-1);
	   do j=1 to &k+2;
	     r= bic_homog1[i,1];
	     if r<=&prnum then bic_homog[r,j]=bic_homog1[i,j];
		 end;
		 end;
	 print "Rank of AIC, value of AIC and ordered subsets for homogeneous variance case:",aic_homog;
	 print "Rank of BIC, value of BIC and ordered subsets for homogeneous variance case:",bic_homog;
    end;
title"AIC and BIC for Heterogeneous Case";
/*AIC and BIC For heterogeneous case */
proc iml;
b=2**(&k-1);
x= j(b,&k,0);
do m=0 to 2**(&k-1)-1;
	y=m;
	x[m+1,1]=1;
	do l=2 to &k;
		x[m+1,l]=floor(y/2**(&k-l));
		y=y-x[m+1,l]*2**(&k-l);
	end;
end;
n=j(2**(&k-1),&k-1,0);
do i=1 to 2**(&k-1);
  do j=1 to &k-1;
       if x[i,j]^=x[i,j+1] then n[i,j]=1;
	   c=2*(2*(n[,+]+1));
		end;
	 end;
heteroaicmatrix=j(2**(&k-1),1,0);
heterobicmatrix=j(2**(&k-1),1,0);
do i=1 to b;
a=j(&k);
   do c=1 to &k;
     exit=&k-c;
       do r=1 to exit;
	      if x[i,c]^= x[i,c+r] then do;
           
		   do j=c+r to &k;
		     a[c,j]=0;
			end;
			r=exit;
		   end;
			 end;
			 end;
			   do p=1 to &k;
			      do q=1 to p-1;
				     a[p,q]=a[q,p];
					 end;
					 end;
			 /*print a;*/
			 Use matrix;
			 read ALL var{n}into mm;
			 read all var{varmle} into nn;
			 read all var{sum} into dd;
			 read all var{mean} into mean1;
			 read all var{ss} into ss;
			 kk=a*mm;
			 ll=a*dd;
			 o=ll#(1/kk);
			 pp=mm#((o-mean1)##2);
		             h1=a*pp;
			 h2=a*ss;
			 h3=log((h1+h2)#(1/kk));
			 h4=0.5*(sum(mm#h3));
			 rr=sum(mm);
			uu=((-1/2)*rr)*(log(2*3.1415926)+1)-h4;
			vv=(-2)*uu;
			AIC=vv+2*(2*(n[i,+]+1));
			/*print AIC;*/
			heteroaicmatrix[i,1]=AIC;
		           BIC=vv+log(mm[+,])*(2*(n[i,+]+1));
			/*print BIC;*/
			heterobicmatrix[i,1]=BIC;
			 end;
			 group=j(2**(&k-1),&k,1);      /*Calculate degrees of freedom for 								heterogeneous variance case*/
 do i=1 to 2**(&k-1);
  do j=1 to &k-1;
       if x[i,j]^=x[i,j+1] then group[i,j+1]=group[i,j]+1; 
	   else group[i,j+1]=group[i,j];
		end;
	 end;
			r3=rank(heteroaicmatrix);
			r4=rank(heterobicmatrix);
            aic_heterog1=r3||heteroaicmatrix||group;
			bic_heterog1=r4||heterobicmatrix||group;
			aic_heterog=j(&prnum,&k+2,0);
			if &k<=3 then do;
 print "Rank of AIC, value of AIC and ordered subsets for patterned heterogeneous variance case:",aic_heterog1;
 print "Rank of BIC, value of BIC and ordered subsets for patterned heterogeneous variance case:",bic_heterog1;
end;
else
do;
     					/* Calculate the smallest prnum AIC value and 
					associated ordered subsets for heterogeneous case */ 
     do i=1 to 2**(&k-1);
	   do j=1 to &k+2;
	     r= aic_heterog1[i,1];
	     if r<=&prnum then aic_heterog[r,j]=aic_heterog1[i,j];
		 end;
		 end;
bic_heterog=j(&prnum,&k+2,0);     /* Calculate the smallest prnum BIC value and 
				  associated ordered subsets for heterogeneous case */ 
     do i=1 to 2**(&k-1);
	   do j=1 to &k+2;
	     r= bic_heterog1[i,1];
	     if r<=&prnum then bic_heterog[r,j]=bic_heterog1[i,j];
		 end;
		 end;
 print "Rank of AIC, value of AIC and ordered subsets for patterned heterogeneous variance case:",aic_heterog;
 print "Rank of BIC, value of BIC and ordered subsets for patterned heterogeneous variance case:",bic_heterog;
  end;                           
quit;