/* This macro computes Krippendorff's alpha reliability estimate for judgments */;
/* made at any level of measurement, any number of judges, with or */;
/* without missing data.  The macro assumes the data file is set up */;
/* in a SAS data file with judges as the variables and the units being */;
/* judged in the rows.  The entries in the data matrix should be */;
/* the coding (quantified or numerically coded for nominal judgments) given */;
/* to the unit in that row by the judge in that column.  Once the macro is */;
/* activated (by running the command set below), the syntax is */;
/* */;
/* %KALPHA (data = dat, judges = judgelist, level = a, detail = b, boot = z); */;
/* */;
/* where 'dat' is the name of the SAS data file, */;
/* 'judgelist' is a list of variable names holding the names of the */;
/* judges, 'a' is the level of measurement (1 = nominal, 2 = ordinal, */;
/* 3 = interval, 4 = ratio), 'b' is set to 1 if you desire SPSS to print */;
/* the coincidence and delta matrices, and 'z' is the number of bootstrap */;
/* samples desired for inference;  z must be at least 1000 and is truncated to the */;
/* lowest 1000 entered (for exampole, 2300 is truncated to 2000) */;
/* The '/level' and '/detail' and '/boot' subcommands are */;
/* optional and default to 1,0, and 0, respectively, if omitted */;
/* */;
/* Missing data should be represented with a 'period' character */;
/* Units that are not coded by at least one judge are excluded from */;
/* the analysis.  The bootstrap estimates are saved into a file named 'bootstp' */;
/* */;
/* This version of the macro was uploaded on June 12, 2011  */;
/* */;
/* */;
/* Written by Andrew F. Hayes */;
/* School of Communication */;
/* The Ohio State University */;
/* hayes.338@osu.edu */;
/* http://www.afhayes.com */;

%macro kalpha (data = ,judges = ,level = 1,detail = 0,boot = 0);
proc iml;
use &data;
read all var{&judges} into dat;
nm={&judges};
btn = &boot;
expdis=1;
indx=0;
pmat=1;
if (&boot > 0) then;
   do;
   btn = floor(&boot/1000)*1000;
   end;
if (&boot > 0) then;
   do;
   if (btn = 0) then;
   do;
   print "Number of bootstraps must be at least 1000.";
   end;
   end;

rw = 1;
do i = 1 to nrow(dat);
good = 0;
do j = 1 to ncol(dat);
  if dat[i,j] ^= . then;
     do;
	   good = good + 1;
	 end;
end;
  if (good > 1) then;
     do;
	   dat[rw,]=dat[i,];
	   rw=rw+1;
	 end;
end;
dat = dat[1:(rw-1),];
nj = ncol(dat);
nobj = nrow(dat);
dat3 = dat;
allm = nobj*nj;
m = shape(dat,allm,1);
j = 0;
do i = 1 to nrow(m);
  if m[i,1] ^= . then;
    do;
	  j=j+1;
	  m[j,]=m[i,];
	end;
end;
m=m[1:j,1];
mss=nrow(m);
mss=allm-mss;
m=.//m;
do i = 2 to nrow(m);
   ix = m[i,1];
   do k = i to 2 by -1;
      if m[k-1,1] > ix then;
	     do;
		   m[k,1]=m[k-1,1];
		 end;
	  if m[k-1,1] <= ix then;
	     do;
         goto stpit;
         end;
   end;
   stpit:
   m[k,1]=ix;
end;
m=m[2:nrow(m),1];
m2=j(nrow(m),1,m[1,1]);
meq4=(m = m2);
yass = meq4[+]/nrow(m);

if (yass <> 1) then;
   do;
   des=design(m);
   uniq=ncol(des);
   coinc=j(uniq,uniq,0);
   delta=coinc;
   map=j(uniq,1,0);
   do i = 1 to nrow(m);
     do j = 1 to uniq;
     if (des[i,j]=1) then;
	   do;
	     map[j,1]=m[i,1];
	   end;
    end;
  end;
  do i = 1 to nobj;
  do j = 1 to nj;
    if dat[i,j] ^= . then;
	   do;
		 do k = 1 to uniq;
		   if dat[i,j]=map[k,1] then;
			  do;
			    dat[i,j]=k;
				goto stpit2;
			  end;
		 end;
		 stpit2:
	   end;
  end;
  end;
  datms = (dat ^= .);
  mu = datms[,+];
  tmp=mu#(mu-1);
  nprs=tmp[+]*0.5;
  btalp=j(btn+1,1,0);
  do k = 1 to nobj;
    temp = j(uniq,uniq,0);
	do i = 1 to nj;
	  do j = 1 to nj;
	     if (dat[k,i] ^= .) then;
		    do;
			  if (dat[k,j] ^= .) then;
			    do;
				   if (i ^= j) then;
				     do;
					 temp[dat[k,i],dat[k,j]]=temp[dat[k,i],dat[k,j]]+(1/(mu[k,1]-1));
                     end;
				end;
			end;
		end;
	end;
	coinc = coinc + temp;
  end;

  q=shape(coinc, nrow(coinc)*ncol(coinc), 1);
  q=(q > 0);
  q=q[+];
  nc=coinc[,+];
  n=nc[+,];
  coinct=coinc;
  dmat = vecdiag(coinc);
  tmp= (dmat > 0);
  nzero = tmp[+];
  bootm = nprs;
  nx=(dmat/n)##bootm;
  nx=round(btn*nx[,+]);
  numone = 0;
   
  expect = coinc;
  do i = 1 to uniq;
    do j = 1 to uniq;
      if (i = j) then;
	     do;
		   expect[i,j]=nc[i,1]*(nc[j,1]-1)/(n-1);
		 end;
	  if (i ^= j) then;
	     do;
		   expect[i,j]=nc[i,1]*nc[j,1]/(n-1);
		 end;
	end;
  end;
  

  tst = q*25;
  tst = tst||(((nj-1)*n)*.5);
  bootm2 = tst[><];
  do z = 1 to (btn+1);
    if (z > 1) then;
	  do;
	   rand = ranuni(0);
       rand = ranuni(repeat(1,bootm2,1));
	   numsum = 0;
	   do i = 1 to bootm2;
	     do j = 2 to indx+1;
           if rand[i,1] <= pmat[j,1] then;
		      do;
			    if rand[i,1] >= pmat[(j-1),1] then;
				  numsum = numsum + pmat[j,2];
			  end;
		 end;
	   end;
	   alpha = 1 - (numsum*(1/(expdis*bootm2)));
	   if (alpha < -1) then;
	     do;
          alpha = -1;
	   end;
	   if (alpha = 1) then;
	     do;
		   if (nzero = 1) then;
		     do;
			   alpha = 0;
			end;
		end;
		if (alpha = 1) then;
		  do;
		    if (nzero > 1) then;
			  do;
			    numone = numone+1;
			end;
		end;
	  btalp[z,1]=alpha;
    end;
	
	if (z = 1) then;
	   do;
	     if (&level = 2) then;
		   do;
	         delta = j(uniq, uniq, 0);
             do i = 1 to uniq;
			   do j = i to uniq;
			     if (i ^=j) then;
				   do;
				     tmp=nc[i:j,1];
				     delta[i,j]=(tmp[+]-(nc[i,1]/2)-(nc[j,1]/2))##2;
					 delta[j,i]=delta[i,j];
				 end;
			   end;
			 end;
			 v={"Ordinal"};
			 if (z = 1) then;
			   do;
                 deltat=delta;
			 end;
          end;
       if (&level = 1) then;
	      do;
		    delta=1-I(uniq);
			v = {"Nominal"};
			deltat=delta;
	   end;
	end;
    if (&level = 3) then;
	  do;
	    do i = 1 to uniq;
		  do j = i to uniq;
		     if (i ^= j) then;
			    do;
				  delta[i,j]=(map[i,1]-map[j,1])##2;
				  delta[j,i]=delta[i,j];
			  end;
		   end;
	     end;
		 v = {"Interval"};
		 deltat = delta;
	end;
	if (&level = 4) then;
	  do;
	    do i = 1 to uniq;
		  do j = i to uniq;
		     if (i ^= j) then;
			    do;
				  delta[i,j]=((map[i,1]-map[j,1])/(map[i,1]+map[j,1]))##2;
				  delta[j,i]=delta[i,j];
			  end;
		   end;
	     end;
		 v = {"Ratio"};
		 deltat = delta;
	 end;
   tmp = delta#coinc;
   num = tmp[+];
   tmp = delta#expect;
   den = tmp[+];
   if (den > 0) then;
     do;
	   alp = 1-(num/den);
	   btalp[1,1]=alp;
	   expdis=den/n;
   end;
   pcoinc = 2*(coinc/n)-(diag(vecdiag(coinc))/n);
   temp=diag(vecdiag(coinc))/n;
   pmat=j((uniq+((uniq*(uniq-1))/2)),2,0);
   psum=0;
   ct=1;
   do i = 1 to uniq;
     do j = i to uniq;
	   psum=psum+pcoinc[j,i];
	   pmat[ct,1]=psum;
	   pmat[ct,2]=delta[j,i];
	   ct=ct+1;
	 end;
   end;
   indx=nrow(pmat);
   t3={0 0};
   pmat=t3//pmat;
  end;
  alpfirst = btalp[1,1];
  if (btn > 0) then;
    do;
      btalp[1,1]=-999;
  if (nx > 0) then;
    do;
	  if (nzero > 1) then;
	    do;
		  chk1 = 0;
		  chk2 = 0;
		    do i = 2 to nrwo(btalp);
			  if (nx >= numone) then;
			     do;
				   if (btalp[i,1] = 1) then;
				     do;
					   if (chk1 < numone) then;
					      do;
						  btalp[i,1]=0;
						  chk1=chk1+1;
						end;
					end;
				end;
				if (nx < numone) then;
				  do;
				    if (btalp[i,1] = 1) then;
					   do;
					      if (chk2 < nx) then;
						    do;
							btalp[i,1]=0;
							chk2=chk2+1;
							end;
					   end;
				  end;
				end;
		  end;   
	end;
  btpp = btalp[2:nrow(btalp),1];
  create bootstp from btpp [colname = 'kalpha'];
  append from btpp;
  do i = 2 to nrow(btalp);
  ix = btalp[i,1];
    do k = i to 2 by -1;
    if (btalp[k-1,1] > ix) then;
       do;
         btalp[k,1]=btalp[k-1,1];
       end;
    if (btalp[k-1,1] <= ix) then;
       do;
         goto stpit3;
       end;
    end;
    stpit3:
    btalp[k,1]=ix;
  end; 
  btalp=btalp[2:nrow(btalp),1];
  mn=btalp[+]/btn;
  low95 = int(0.025*btn);
  high95 = int(0.975*btn)+1;
  low95 = btalp[low95,1];
  high95 = btalp[high95,1];
  median = btalp[0.50*btn,1];
  q = {.9 0, .8 0, .7 0, 0.67 0, .6 0, .5 0};
  do i = 1 to 6;
    qcomp = (btalp < q[i,1]);
    qcomp = qcomp[+]/btn;
    q[i,2]=qcomp;
  end;
  end;
  print "Krippendorff's Alpha Reliability Estimate";
  if (btn = 0) then;
    do;
    res = alpfirst||nobj||nj||nprs;
	lab = {"Alpha" "Units" "Obsrvrs" "Pairs"};
  end;
  if (btn > 0) then;
    do;
	res = alpfirst||low95||high95||nobj||nj||nprs;
	lab = {"Alpha" "LL95CI" "UL95CI" "Units" "Observrs" "Pairs"};
  end;
  print res [rowname = v colname = lab format = 10.4];
  if (btn > 0) then;
    do;
      print "Probability (q) of failure to achieve an alpha of at least alphamin:";
      cn = {"Alphamin" "q"};
      print q [colname = cn format = 10.4];
      print "Number of Bootstrap Samples";
      print btn;
  end;
  print "Judges used in these computations";
  print nm;
  if (&detail = 1) then;
    do;
      print "==========================================================";
      print "Observed Coincidence Matrix";
      print coinct [format = 9.2];
	  print "Expected Coincidence Matrix";
	  print expect [format = 9.2];
	  print "Delta Matrix";
	  print deltat [format = 9.2];
  tmap = map`;
  print "Rows and columns correspond to the following unit values";
  print tmap [format = 9.2];
  end;
end;
if (yass = 1) then;
  do;
    print "ERROR: Input Reliability Data Matrix Exhibits No Variation";
end;
quit;
%mend kalpha;


data test;
  input a b c;
  cards;
  0 1 1
  1 0 1
  1 1 0 
  0 0 .
run;

%kalpha (data = test, judges = a b c, detail = 1, level = 1, boot = 2000);