/*-------------------------------------------------------------------*
 * KAPPA SAS      Compute Cohen's Kappa for agreement in a           *
 *                frequency table [Weighted/Unweighted]              *
 *-------------------------------------------------------------------*
 *  Author:  Michael Friendly            <friendly@yorku.ca>         *
 * Created:  5 Mar 1991 21:41:00            (c) 1991                 *
 * Revised:  13 Apr 1991 23:59:30                                    *
 * Version:  1.0                                                     *
 * Requires: SAS/IML                                                 *
 *-------------------------------------------------------------------*/
*-- Calculate Cohens Kappa for agreement;
proc iml;
start kappa(freq, w, kappa, var);
   /*--------------------------------------------------------------*
    |  Arguments:                                                  |
    |    freq  - a square matrix of frequencies, rater A x rater B |
    |    w     - weight matrix, same size as freq. w=identity for  |
    |            for unweighted kappa.                             |
    |    kappa - output kappa value                                |
    |    var   - variance of kappa                                 |
    *--------------------------------------------------------------*/
   row_sum = freq[,+];
   col_sum = freq[+,];
   n       = freq[+,+];
   obs_agr = sum( w # freq );
   exp_agr = col_sum * w * row_sum / n;
   kappa = ( obs_agr - exp_agr ) / ( n - exp_agr );
 
   pobs  = obs_agr / n;
   pexp  = exp_agr / n;
   Print 'Observed and Expected Agreement (under independence)',
       (obs_agr / n)[c='Observed']
       (exp_agr / n)[c='Expected'];
 
   *-- Estimated variance of kappa --;
   wbari = col_sum * w / n;
   wbarj = w * row_sum / n;
   do i = 1 to nrow(freq);
   do j = 1 to nrow(freq);
      wd = ( w[i,j] - wbari[i] - wbarj[j] )##2;
      end; end;
   var = (( col_sum * wd * row_sum / n ) - pexp##2 ) / n # (1 - pexp##2);
finish;
 
   print "Cohen's Kappa: Husbands and Wives Sexual Fun";
   freq =
     {  7     7     2      3,     /* Never fun  */
        2     8     3      7,     /* Sometimes  */
        1     5     4      9,     /* Very Often */
        2     8     9     14 };   /* Always fun */
 
   print "Observed Frequencies", f[format=5.];
   reset noname;
   r = nrow(freq);
   w = I(r);
   print 'Weights' , w[format=fract8.];
   run kappa(freq, w, kappa, var);
   print "Cohen's Kappa (Unweighted) for Agreement",
         kappa[colname='Kappa'] (sqrt(var))[c='Std. Error'];
 
   w = diag({1 1 1 1});
   do i=1 to r;
   do j=1 to r;
      w[i,j] = 1 / ( 1 + abs(i-j) );
      end; end;
   print 'Weights' , w[format=fract8.];
   run kappa(freq, w, kappa, var);
   print "Cohen's Kappa (Weighted) for Agreement",
         kappa[colname='Kappa'] (sqrt(var))[c='Std. Error'];
 
   do i=1 to r;
   do j=1 to r;
      w[i,j] = 1 -((i-j)##2 / (r-1)##2);
      end; end;
   print 'Fleiss-Cohen Weights' , w[format=fract8.];
   run kappa(freq, w, kappa, var);
   print "Cohen's Kappa (Weighted) for Agreement",
         kappa[colname='Kappa'] (sqrt(var))[c='Std. Error'];