Каппа Коэна (Cohen's Kappa)

*В 1997 году Дэвид Николс (David Nichols) из SPSS подготовил синтаксис для вычисления каппы Коэна (Cohen's Kappa),
*соответствующей стандартной ошибки, z-статистики и её уровня значимости.
*Этот синтаксис основан на предложенном им коде и использует его для вычисления
*перечисленных выше статистик. Исходный синтаксис был расширен мной после изучения материалов, представленных
*на ежегодном съезде Ассоциации Southwest Educational Research Association Джейсоном Кингом (Jason E. King) в Колледже Baylor
*College of Medicine. Работа называется "Software Solutions for Obtaining a Kappa-Type Statistic for Use with
*Multiple Raters" ("Программная реализация вычисления каппа-подобных статистик для множественных
*экспертных оценок"). Данный синтаксис группирует выдачу по 4 разделам. Первый раздел включает меру согласованности
*экспертных оценок (Pa), число объектов, число экспертов, число категорий оценок. Второй раздел включает статистику каппа (k), 
*стандартную ошибку, статистику z и её значимость p на основе исходного синтаксиса, а также добавляет верхнюю и нижнюю границы
*доверительного интервала на уровне доверия 95%. Третий раздел представляет исправленную стандартную ошибку по формуле Флейса
*(Fleiss) с пересчётом статистик z, p и доверительного интервала. Последний раздел включает следующую информацию по отдельным
*категориям оценок: каппа (k), стандартная ошибка, z, p и доверительный интервал.
*.
*Исходные данные интерпретируются таким образом, что в строках содержатся объекты/наблюдения, а в столбцах - оценки каждого
*эксперта. Значение в каждой ячейке - это категория (оценка), приписанная отдельным экспертом отдельному наблюдению/объекту.
*После кода определения макроса есть команда, запускающая макрос. При вызове макроса  переменные с оценками экспертов должны
*быть расположены последовательно в файле данных; в команде вызова они перечисляются так же, как и в примере в конце файла.
*Макрос протестирован с 20 экспертами, 20 категориями и 2000 наблюдениями. Если у вас больший масштаб, возможно, потребуется
*скорректировать параметр mxloops. Синтаксис подставляет -99999.0 в качестве пропущенного значения при выдаче для
*неиспользованных категорий оценок.

* Автор: Brian G. Dates, 23.02.2006.

data list list /rater1 rater2 rater3 .
begin data
1 1 1
2 1 2
2 2 2
2 1 1
2 1 2
2 1 2
2 2 1
2 1 2
2 1 2
2 1 1
2 1 3
4 4 3
4 6 4
5 5 5
4 4 4
6 6 6
4 4 3
2 1 6
2 2 3
1 3 1
end data .
preserve.
set printback=off mprint=off.
save outfile='k__tmp1.sav'.

define cohensk (vars=!charend('/')).
set mxloops=100000.
count ms__=!vars (missing).
select if ms__=0.

*В этом разделе формируется матрица x на основе исходных данных, матрица y с числом строк, равным числу объектов
*и числом столбцов, равным числу категорий, затем для в матрице y фиксируется число ответов на категорию
*для каждого из объектов. 

matrix.
get x /var=!vars.
compute c=mmax(x).
compute y=make(nrow(x),c,0).
loop i=1 to nrow(x).
loop j=1 to ncol(x).
loop k=1 to c.
do if x(i,j)=k.
compute y(i,k)=y(i,k)+1.
end if.
end loop.
end loop.
end loop.

*В этом разделе вычисляется сводная информация, каппа и связанные с ней статистики.

compute pe=msum((csum(y)/msum(y))&**2).
compute k=ncol(x).
compute n=nrow(y) .
compute r=ncol(y) .
compute pa=mssq(y)/(nrow(y)*k*(k-1))-(1/(k-1)).
compute ck=(pa-pe)/(1-pe).
compute num=2*(pe-(2*k-3)*(pe**2)+2*(k-2)*msum((csum(y)/msum(y))&**3)).
compute den=nrow(y)*k*(k-1)*((1-pe)**2).
compute ase=sqrt(num/den).
compute z=ck/ase.
compute sig=1-chicdf(z**2,1) .
compute ckul=ck+1.96*ase .
compute ckll=ck-1.96*ase .

*В этом разделе вычисляется альтернативная (исправленная) оценка стандартной ошибки (по формуле Флейса)
*и связанные с ней статистики.

compute nm=sqrt(n*k*(k-1)) .
compute vectora=csum(y)/msum(csum(y)) .
compute vectorb=1-csum(y)/msum(csum(y)) .
compute vectorc=1-2*(csum(y)/msum(csum(y))) .
compute vectord=vectora&*vectorb.
compute vectore=vectora&*vectorb&*vectorc .
compute e=msum(vectord) .
compute f =msum(vectore) .
compute fse=(sqrt(2)/(e*nm))*(sqrt(e**2-f)) .
compute fsez=ck/fse .
compute fsesig=1-chicdf(fsez**2,1) .
compute fseu=ck+1.96*fse .
compute fsel=ck-1.96*fse .


* Этот раздел рассчитывает каппу для отдельных категорий. Каждая статистика (например, k или стандартная ошибка)
* рассчитывается как вектор. Затем вектора собираются в матрицу из 6 статистик. В ходе вычислений значение
* -99999.0 интерпретируется как пропуск.

compute matz=k-y .
compute mata=y&*matz .
compute vectorf=csum(mata)+(.0001) .
compute vectorg=vectord*(n*k*(k-1))+(.0001) .
compute vectorh=1-(vectorf&/vectorg) .
compute vectori=(1+(2*(k-1)*(csum(y)/msum(csum(y)))))&**2 .
compute vectorj=(2*(k-1))*vectord .
compute vectork=(n*k*(k-1)**2)*vectord+(.0001) .
compute vectorse=sqrt((vectori+vectorj)&/vectork) .
compute vectorz=vectorh&/(vectorse+.0001).
compute vectorp=1-cdfnorm(vectorz) .
compute vectorll=vectorh-(1.96*vectorse) .
compute vectorul=vectorh+(1.96*vectorse) .
loop i=1 to ncol(vectorh) .
do if (vectorh(i)=0.00) .
compute vectorh(i)=-99999) .
end if .
end loop .
loop i=1 to ncol(vectorh) .
do if (vectorh(i)=-99999) .
compute vectorse(i)=-99999 .
end if .
end loop .
loop i=1 to ncol(vectorh) .
do if (vectorh(i)=-99999) .
compute vectorz(i)=-99999 .
end if .
end loop .
loop i=1 to ncol(vectorh) .
do if (vectorh(i)=-99999) .
compute vectorp(i)=-99999 .
end if .
end loop .
loop i=1 to ncol(vectorh) .
do if (vectorh(i)=-99999) .
compute vectorul(i)=-99999.
end if .
end loop .
loop i=1 to ncol(vectorh) .
do if (vectorh(i)=-99999) .
compute vectorll(i)=-99999 .
end if .
end loop .
compute ikstat={vectorh;vectorse;vectorz;vectorp;vectorll;vectorul} .
save ikstat /outfile='ikstat1.sav' .

*Следующий раздел сохраняет вычисленные данные и форматирует их к выдаче.

save {k,n,r,pa,ck,ase,z,sig,ckll,ckul,fse,fsez,fsesig,fseu,fsel} /outfile='k__tmp2.sav'
     /variables=k,n,r,pa,ck,ase,z,sig,ckll,ckul,fse,fsez,fsesig,fseu,fsel .
end matrix .
get file='k__tmp2.sav'.
formats k (f8.0) /n (f8.0) /r (f8.0) /pa (f8.4) /ck (f8.4) /ase (f8.4) /z (f8.4) /sig (f8.4) /ckul (f8.4) /ckll (f8.4)
 /fse(f8.4) /fsez (f8.4) /fsesig (f8.4) /fseu (f8.4) /fsel (f8.4) .
variable labels k 'Число экспертов' /n 'Число объектов' /r 'Число категорий' /pa '% согласованности мнений экспертов'
 /ck 'Каппа' /ase 'Станд. ошибка' /z 'z'/sig 'p' /ckul 'ВГ 95% ДИ' /ckll 'НГ 95% ДИ'
 /fse 'Станд. ошибка по Флейсу' /fsez 'z' /fsesig 'p' /fseu 'ВГ 95% ДИ' /fsel 'НГ 95% ДИ' .
report format=list automatic align(center)
  /variables=k,n,r,pa
  /title "Сводная информация" .
report format=list automatic align(center)
  /variables=ck ase z sig ckll ckul
  /title "Каппа Коэна".
report format=list automatic align(center)
  /variables=fse fsez fsesig fsel fseu
  /title "Каппа Коэна -- стандартная ошибка исправлена по Флейсу" .
get file='ikstat1.sav' .
flip .
delete variable case_lbl .
compute n1=$casenum .
formats n1 (f8.0) /var001 (f8.4) /var002 (f8.4)  /var003 (f8.4)  /var004 (f8.4)  /var005 (f8.4)  /var006 (f8.4)  .
variable labels n1 'Код категории' /var001 'k' /var002 'Станд. ошибка' /var003 'z' /var004 'p'
  /var005 'НГ 95% ДИ' /var006 'ВГ 95% ДИ' .
save outfile='ikstat.sav' .
report format=list automatic align(center)
  /variables=n1 var001 var002 var003 var004 var005 var006
  /Title "Статистики по отдельным категориям" .
!enddefine.
restore.
COHENSK VARS = rater1 to rater3.