data kappa;
N=27;
R=3;
K=4;
do sub=1 to N;
 do cat=1 to R;
        input count @@;
        output;
      end; end; cards;
1 3 .
4 . .
1 . 3
. 3 1
. . 4
. 2 2
2 . 2
3 . 1
3 1 .
. 1 3
1 . 3
4 . .
. 1 3
. 4 .
3 . 1
1 1 2
. 1 3
1 3 .
. . 4
1 2 1
4 . .
1 3 .
. 3 1
. . 4
2 2 .
. 4 .
3 . 1
;
run;
data kappa;
set kappa;
count2=count**2;
*count3=count2*count;
countag=count2-count;
run;
title "calculate phat for each rater";
proc means data=kappa noprint;
var countag;
id N R K;
by sub;
output out=sumag sum=sumag;
run;
data sumag;
set sumag;
phat=sumag/(K*(K-1));
run;
title "calculate pbar";
proc means data=sumag noprint;
var phat;
id N;
output out=pbar mean=pbar;
run;


title "calculate pe";
proc freq data=kappa order=data  noprint;
      weight count;
	  tables cat/out=col_total;
	  *test kappa ;
run;
data col_total;
set col_total;
perc2=(percent/100)**2;
perc3=(percent/100)**3;
run;
proc means data=col_total noprint;
var perc2;
output out=pe sum=pe;
run;
title "calculate estimate of KG";
data result1;
merge pbar pe;
drop _type_ _freq_;
KGhat=(pbar-pe)/(1-pe);
run;
title "calculate variance of KG and test statistics";
proc means data=col_total noprint;
var perc3;
output out=sump3 sum=sump3;
run;
data result0;
merge result1 sump3;
drop _type_ _freq_;
run;
data result;
merge sumag result0;
by N;
drop _type_ _freq_;
varKG=(pe-(2*K-3)*pe**2+2*(K-2)*sump3)/((1-pe)**2)*2/(N*K*(K-1));
z=KGhat/sqrt(varKG);
pval=probnorm(-1*abs(z))+1-probnorm(abs(z));
run;
data finalresult;
set result;
keep sub N K R phat pbar pe KGhat varKG z pval;
run;
title "display final Kappa test result";
/*proc print data=finalresult;
run;

title "display p-value";*/
data finalresult;
set finalresult;
if (pval<0.001) then result="p value < 0.001";
else result=pval;
run;
proc print data=finalresult;
run;