data d1;

input a b c d;

*data from Siegel, S., & Castellan, N. J. (1988). Nonparametric statistics

for the behavioral sciences (p. 287). New York: McGraw-Hill.;

cards;

5 5 5 5

1 1 3 3

5 5 5 5

1 1 3 3

4 5 5 5

1 2 3 3

1 1 1 3

3 3 4 4

1 1 1 3

5 5 5 5

1 1 1 1

3 3 3 4

1 3 3 4

4 4 5 5

5 5 5 5

5 3 3 3

3 3 3 2

3 5 5 5

3 3 3 4

1 1 1 1

1 1 3 3

1 3 3 3

1 3 1 3

1 1 3 3

2 3 3 5

;

options nodate pagesize=60 linesize=80;

***********************************************************

* Each judge is a column

* Each row is an object that was judged

* The numbers in the columns are the category judgments

* assigned by that judge.

* All category numbers must be positive.

* Delete missing data prior to using IML, or else

* recode missing to be a unique positive number

* (e.g., 999 for all missing data).

***********************************************************;

proc iml;

* read data from dataset;

use d1;

read all into x;

*print x;

* find categories;

lb=min(x);

ub=max(x);

r=ub-lb+1;

hold1=j(r,2,0);

do l = lb to ub;

hold1[l,1]=l;

check=j(nrow(x),ncol(x),l);

hold1[l,2]=any(x=check);

end;

k = ncol(X); * k is the number of judges;

If k = 2 then do;

***********************************;

* find the 2-judge agreement matrix;

***********************************;

d=hold1[+,2]+1;

A=j(d,d,0);

m=1;

do l=lb to ub;

if hold1[l,2]=1 then do;

m=m+1;

A[m,1]=hold1[l,1];

A[1,m]=hold1[l,1];

end;

* fill the entries in the agreement matrix;

do l=1 to nrow(x);

do xc=1 to nrow(A);

if X[l,1] = A[xc,1] then p1=xc;

if X[l,2] = A[1,xc] then p2=xc;

end;

A[p1,p2]=A[p1,p2]+1;

end;

print "Agreement Matrix", A;

end;

****************************************************************;

*find contingency table (Siegel & Castellan, 1988, pp. 285-287);

****************************************************************;

numcat=hold1[+,2]; *numcat is the number of categories;

catlab=j(1,numcat,0); *catlab is the numerical labels of the categories;

b=1;

do a = 1 to nrow(hold1);

if hold1[a,2]=1 then do;

catlab[1,b]=hold1[a,1];

b=b+1;

end;

print "the judges used these catagories", catlab;

numobj=nrow(X); *numobj is the number of objects that were judged;

con=j(numobj,numcat,0);

do a = 1 to numobj;

do b = 1 to k;

do c = 1 to numcat;

if X[a,b]=catlab[1,c] then con[a,c]=con[a,c]+1;

end;

print "Contingency table",con;

* find kappa;

* ;

*find column totals;

totals=con[+,];

totn=totals[,+];

props=totals/totn;

PE=props*props`;

print "proportion of expected chance agreement", PE;

s=j(nrow(con),1,0);

do l=1 to nrow(con);

do m=1 to ncol(con);

if con[l,m] > 0 then s[l,1]=s[l,1]+con[l,m]*(con[l,m]-1);

end;

s=s*1/(k*(k-1));

PA = s[+,]/nrow(s);

print "proportion agreement", PA;

kappa=(PA-PE)/(1-PE);

*print kappa;

* find variance of kappa;

p3=props#props#props;

p3=p3[,+];

*print p3;

part1=2/(numobj*k*(k-1));

*print part1;

part2=(PE-(2*k-3)*(PE*PE)+2*(k-2)*p3)/((1-PE)*(1-PE));

*print part2;

var = part1*part2;

*print var;

ci95u=kappa+1.96*var;

ci95l=kappa-1.96*var;

print "95% lower bound, kappa, 95% upper bound", ci95l kappa ci95u;

quit;

run;