1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
*В 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.