1 function [kap,se,H,z,p0,SA,R]=kappa(d,c,arg3,w)
2 % KAPPA estimates Cohen's kappa coefficient
3 % and related statistics
5 % [...] = kappa(d1,d2);
6 % NaN's are handled as missing values and are ignored
7 % [...] = kappa(d1,d2,'notIgnoreNAN');
8 % NaN's are handled as just another Label.
9 % [kap,sd,H,z,ACC,sACC,MI] = kappa(...);
15 % kap Cohen's kappa coefficient point
16 % se standard error of the kappa estimate
17 % H Concordance matrix, i.e. confusion matrix
19 % ACC overall agreement (accuracy)
20 % sACC specific accuracy
21 % MI Mutual information or transfer information (in [bits])
22 % X is a struct containing all the fields above
23 % For two classes, a number of additional summary statistics including
24 % TPR, FPR, FDR, PPV, NPF, F1, dprime, Matthews Correlation coefficient (MCC) or
25 % Phi coefficient (PHI=MCC), Specificity and Sensitivity
26 % are provided. Note, the positive category must the larger label (in d and c), otherwise
27 % the confusion matrix becomes transposed and the summary statistics are messed up.
31 % [1] Cohen, J. (1960). A coefficient of agreement for nominal scales. Educational and Psychological Measurement, 20, 37-46.
32 % [2] J Bortz, GA Lienert (1998) Kurzgefasste Statistik f|r die klassische Forschung, Springer Berlin - Heidelberg.
33 % Kapitel 6: Uebereinstimmungsmasze fuer subjektive Merkmalsurteile. p. 265-270.
34 % [3] http://www.cmis.csiro.au/Fiona.Evans/personal/msc/html/chapter3.html
35 % [4] Kraemer, H. C. (1982). Kappa coefficient. In S. Kotz and N. L. Johnson (Eds.),
36 % Encyclopedia of Statistical Sciences. New York: John Wiley & Sons.
37 % [5] http://ourworld.compuserve.com/homepages/jsuebersax/kappa.htm
38 % [6] http://en.wikipedia.org/wiki/Receiver_operating_characteristic
40 % $Id: kappa.m 9608 2012-02-10 09:56:25Z schloegl $
41 % Copyright (c) 1997-2006,2008,2009,2011 by Alois Schloegl <alois.schloegl@gmail.com>
42 % This function is part of the NaN-toolbox
43 % http://pub.ist.ac.at/~schloegl/matlab/NaN/
45 % BioSig is free software: you can redistribute it and/or modify
46 % it under the terms of the GNU General Public License as published by
47 % the Free Software Foundation, either version 3 of the License, or
48 % (at your option) any later version.
50 % BioSig is distributed in the hope that it will be useful,
51 % but WITHOUT ANY WARRANTY; without even the implied warranty of
52 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
53 % GNU General Public License for more details.
55 % You should have received a copy of the GNU General Public License
56 % along with BioSig. If not, see <http://www.gnu.org/licenses/>.
63 if strcmpi(arg3,'notIgnoreNAN')
80 tmp(isnan(tmp)) = maxtmp+1;
81 [X.Label,i,j] = unique(tmp);
82 c = j(1+numel(d):end);
85 maxCLASS = kk - any(tmp>maxtmp);
89 % fprintf(2,'Warning KAPPA: some elements are NaN. These are handled as missing values and are ignored.\n');
90 % fprintf(2,'If NaN should be handled as just another label, use kappa(..,''notIgnoreNaN'').\n');
91 ix = find((c<=maxCLASS) & (d<=maxCLASS));
93 if ~isempty(w), w = w(ix); end;
96 X.Label(X.Label>maxtmp) = [];
98 X.Label(X.Label>maxtmp) = NaN;
102 H = full( sparse (d, c, 1, kk, kk) );
104 H = full( sparse (d, c, w, kk, kk) );
108 X.Label = 1:min(size(d));
109 H = d(X.Label,X.Label);
117 p0 = sum(diag(H))/N; %accuracy of observed agreement, overall agreement
118 %OA = sum(diag(H))/N);
123 SA = 2*diag(H)'./(p_i+pi_); % specific agreement
125 pe = (p_i*pi_')/(N*N); % estimate of change agreement
127 px = sum(p_i.*pi_.*(p_i+pi_))/(N*N*N);
130 kap = (p0-pe)/(1-pe);
131 sd = sqrt((pe+pe*pe-px)/(N*(1-pe*pe)));
134 se = sqrt((p0+pe*pe-px)/N)/(1-pe);
141 if ((1 < nargout) && (nargout<7))
147 pwi = sum(H,2)/N; % p(x_i)
148 pwj = sum(H,1)/N; % p(y_j)
149 pji = H./repmat(sum(H,2),1,size(H,2)); % p(y_j | x_i)
150 R = - sumskipnan(pwj.*log2(pwj)) + sumskipnan(pwi'*(pji.*log2(pji)));
152 if (nargout>1), return; end;
162 X.datatype = 'confusion';
165 % see http://en.wikipedia.org/wiki/Receiver_operating_characteristic
166 % Note that the confusion matrix used here is has positive values in
167 % the 2nd row and column, moreover the true values are indicated by
168 % rows (transposed). Thus, in summary H(1,1) and H(2,2) are exchanged
169 % as compared to the wikipedia article.
174 X.FNR = H(2,1) / sum(H(2,:));
175 X.FPR = H(1,2) / sum(H(1,:));
176 X.TPR = H(2,2) / sum(H(2,:));
177 X.PPV = H(2,2) / sum(H(:,2));
178 X.NPV = H(1,1) / sum(H(:,1));
179 X.FDR = H(1,2) / sum(H(:,2));
180 X.MCC = det(H) / sqrt(prod([sum(H), sum(H')]));
182 X.F1 = 2 * X.TP / (sum(H(2,:)) + sum(H(:,2)));
183 X.Sensitivity = X.TPR; %% hit rate, recall
184 X.Specificity = 1 - X.FPR;
186 X.dprime = norminv(X.TPR) - norminv(X.FDR);