* Доверительный интервал для разности двух медиан по методу Ходжеса-Лимана (Hodges-Lehmann): Hodges JR and Lehmann EL (1963) 'Estimates of location based on rank tests' Annals of Mathematical Statistics 34, 598-611. * Улучшен алгоритм сортировки: 16 мая 2003 *. * Хотя пример приведён для выборок одинакового объёма, алгоритм будет работать и для выборок неравного объёма тоже. * Программа составлена: Marta Garcia-Granero. * Пример данных *. DATA LIST FREE /group(F4.0) wghtgain(F8.0). BEGIN DATA 1 175 1 149 1 132 1 187 1 218 1 123 1 151 1 248 1 200 1 206 1 219 1 179 1 234 1 206 2 142 2 214 2 311 2 249 2 337 2 176 2 262 2 211 2 302 2 216 2 195 2 236 2 253 2 199 END DATA. VAR LABEL group 'Экспериментальная группа' /wghtgain 'Урожай (фунтов)'. VALUE LABELS group 1 'Control' 2 'Vit. A'. ************************************************ * Программа расчёта доверительного интервала * ************************************************ * Если n1*n2>10000, измените параметр MXLOOPS соответствующим образом. * Данные предварительно сортируются по группирующей и зависимой переменной. * Случаи с пропущенными значениями временно удаляются перед использованием MATRIX *. SET MXLOOPS=10000. SORT CASES BY group (A) wghtgain (A) . TEMPORARY. SELECT IF ((NOT MISSING(group)) AND (NOT MISSING(wghtgain))). MATRIX. PRINT /TITLE 'ДОВЕРИТЕЛЬНЫЙ ИНТЕРВАЛ ХОДЖЕСА-ЛИМАНА ДЛЯ РАЗНОСТИ ДВУХ МЕДИАН'. PRINT /TITLE 'Допущение: распределения имеют одинаковую форму'. * Ввод данных: замените имена переменных на актуальные в вашем случае *. GET group /VAR=group. GET depvar /VAR=wghtgain. * Получим группирующие значения, подсчитаем объёмы выборок и разобъём данные на 2 группы *. COMPUTE totaln=NROW(group). COMPUTE ngrp1=group(1). COMPUTE ngrp2=group(totaln). PRINT {ngrp1;ngrp2} /TITLE='Группирующие значения' /RLABELS='Группа A=' 'Группа B=' /FORMAT='f5.0'. COMPUTE n=CSUM(DESIGN(group)). PRINT {T(n)} /FORMAT='f5.0' /TITLE='Объёмы выборок' /RLABELS=' N(a)=' ' N(b)='. DO IF RMIN(n) GT 1. /* Не считаем доверительные интервалы, если хотя бы в одной группе число наблюдений = 1. COMPUTE group1=MAKE(1,n(1),0). COMPUTE group2=MAKE(1,n(2),0). LOOP i=1 TO n(1). + COMPUTE group1(i)=depvar(i). END LOOP. LOOP i=1+n(1) TO totaln. + COMPUTE group2(i-n(1))=depvar(i). END LOOP. RELEASE ngrp1,ngrp2,group,depvar,totaln. * Пятичисельные итоги для обеих выборок *. COMPUTE min1=RMIN(group1). COMPUTE max1=RMAX(group1). COMPUTE pair=(TRUNC(n(1)/2) EQ n(1)/2). /* Проверяем, число n(1) нечётное (результат 0) или чётное (результат 1) (1) *. DO IF pair EQ 1. /* Квартили Тьюки для чётного числа объектов в выборке *. - COMPUTE lmedian=group1(n(1)/2). - COMPUTE umedian=group1(1+n(1)/2). - COMPUTE median1=(lmedian+umedian)/2. - RELEASE lmedian,umedian. - DO IF TRUNC(n(1)/4) NE (n(1)/4)). - COMPUTE q11=group1(1+(n(1)/2)/2). - COMPUTE q31=group1(1+(3*n(1)/2)/2). - ELSE. - COMPUTE lq11=group1(n(1)/4). - COMPUTE uq11=group1(1+n(1)/4). - COMPUTE lq31=group1(3*n(1)/4). - COMPUTE uq31=group1(1+3*n(1)/4). - COMPUTE q11=(lq11+uq11)/2. - COMPUTE q31=(lq31+uq31)/2. - RELEASE lq11,uq11,lq31,uq31. - END IF. ELSE. /* Квартили Тьюки для нечётного числа объектов в выборке *. - COMPUTE median1=group1((n(1)+1)/2). - DO IF TRUNC((1+n(1))/4) EQ (1+n(1))/4). - COMPUTE lq11=group1((n(1)+1)/4). - COMPUTE uq11=group1(1+(n(1)+1)/4). - COMPUTE lq31=group1(3*(n(1)+1)/4-1). - COMPUTE uq31=group1(3*(n(1)+1)/4). - COMPUTE q11=(lq11+uq11)/2. - COMPUTE q31=(lq31+uq31)/2. - RELEASE lq11,uq11,lq31,uq31. - ELSE. - COMPUTE q11=group1(((1+(n(1)+1)/2))/2). - COMPUTE q31=group1(((3*(n(1)+1)/2)-1)/2). - END IF. END IF. COMPUTE min2=RMIN(group2). COMPUTE max2=RMAX(group2). COMPUTE pair=(TRUNC(n(2)/2) EQ n(2)/2). /* Проверяем, число n(2) нечётное (результат 0) или чётное (результат 1) *. DO IF pair EQ 1. /* Квартили Тьюки для чётного числа объектов в выборке *. - COMPUTE lmedian=group2(n(2)/2). - COMPUTE umedian=group2(1+n(2)/2). - COMPUTE median2=(lmedian+umedian)/2. - RELEASE lmedian,umedian. - DO IF TRUNC(n(2)/4) NE (n(2)/4)). - COMPUTE q12=group2(1+(n(2)/2)/2). - COMPUTE q32=group2(1+(3*n(2)/2)/2). - ELSE. - COMPUTE lq12=group2(n(2)/4). - COMPUTE uq12=group2(1+n(2)/4). - COMPUTE lq32=group2(3*n(2)/4). - COMPUTE uq32=group2(1+3*n(2)/4). - COMPUTE q12=(lq12+uq12)/2. - COMPUTE q32=(lq32+uq32)/2. - RELEASE lq12,uq12,lq32,uq32. - END IF. ELSE. /* Квартили Тьюки для нечётного числа объектов в выборке *. - COMPUTE median2=group2((n(2)+1)/2). - DO IF TRUNC((1+n(2))/4) EQ (1+n(2))/4). - COMPUTE lq12=group2((n(2)+1)/4). - COMPUTE uq12=group2(1+(n(2)+1)/4). - COMPUTE lq32=group2(3*(n(2)+1)/4-1). - COMPUTE uq32=group2(3*(n(2)+1)/4). - COMPUTE q12=(lq12+uq12)/2. - COMPUTE q32=(lq32+uq32)/2. - RELEASE lq12,uq12,lq32,uq32. - ELSE. - COMPUTE q12=group2(((1+(n(2)+1)/2))/2). - COMPUTE q32=group2(((3*(n(2)+1)/2)-1)/2). - END IF. END IF. PRINT {median1,q11,q31,min1,max1,(q31-q11),(max1-min1);median2,q12,q32,min2,max2,(q32-q12),(max2-min2)} /FORMAT='f8.1' /TITLE='Пятичисельный итог и показатели разброса' /CLABELS='Медиана' 'Q1' 'Q3' 'Мин' 'Макс' 'ВтрКвРазм.' 'Размах' /RLABELS='Группа A' 'Группа B'. /*ВтрКвРазм - внутриквартильный размах*. RELEASE median1,q11,q31,min1,max1,median2,q12,q32,min2,max2. * Построим вектор попарных разностей *. COMPUTE n1n2=n(1)&*n(2). COMPUTE diff=MAKE(1,n1n2,0). LOOP i=1 TO n(1). - LOOP j=1 TO n(2). - COMPUTE diff(n(2)*(i-1)+j)=group1(i)-group2(j). - END LOOP. END LOOP. * Сортируем разности (алгоритм R Ristow и J Peck) *. COMPUTE sdiff = diff. COMPUTE sdiff(grade(diff)) = diff. RELEASE diff,group1,group2. * Считаем медиану разностей *. COMPUTE pair=(TRUNC(n1n2/2) EQ n1n2/2). /* Проверяем, произведение n1n2 нечётное (результат 0) или чётное (результат 1) *. DO IF pair EQ 0. /* Формула медианы для нечётных выборок. - COMPUTE median=sdiff((n1n2+1)/2). ELSE. /* Формула медианы для чётных выборок *. - COMPUTE median=(sdiff(n1n2/2)+sdiff(1+(n1n2/2)))/2. END IF. RELEASE pair. PRINT median /TITLE='Разность между медианами в генеральной совокупности * (A - B) ' /RLABELS='Тччн. оц.: ' /FORMAT='f8.1'. /*точечная оценка*. PRINT /TITLE='(*) Внимание! Это НЕ то же самое, что разность между выборочными медианами'. * Точные или асимптотические 95% и 99% доверительные интервалы *. DO IF ((n(1) LE 20) AND (n(2) LE 20)). /* Точные критические значения Мана-Уитни *. - COMPUTE d95={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0; 0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,2,2,2,2; 0,0,0,0,0,1,1,2,2,3,3,4,4,5,5,6,6,7,7,8; 0,0,0,0,1,2,3,4,4,5,6,7,8,9,10,11,11,12,13,13; 0,0,0,1,2,3,5,6,7,8,9,11,12,13,14,15,17,18,19,20; 0,0,1,2,3,5,6,8,10,11,13,14,16,17,19,21,22,24,25,27; 0,0,1,3,5,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34; 0,0,2,4,6,8,10,13,15,17,19,22,24,26,29,31,34,36,38,41; 0,0,2,4,7,10,12,15,17,20,23,26,28,31,34,37,39,42,45,48; 0,0,3,5,8,11,14,17,20,23,26,29,33,36,39,42,45,48,52,55; 0,0,3,6,9,13,16,19,23,26,30,33,37,40,44,47,51,55,58,62; 0,1,4,7,11,14,18,22,26,29,33,37,41,45,49,53,57,61,65,69; 0,1,4,8,12,16,20,24,28,33,37,41,45,50,54,59,63,67,72,76; 0,1,5,9,13,17,22,26,31,36,40,45,50,55,59,64,69,74,78,83; 0,1,5,10,14,19,24,29,34,39,44,49,54,59,64,70,75,80,85,90; 0,1,6,11,15,21,26,31,37,42,47,53,59,64,70,75,81,86,92,98; 0,2,6,11,17,22,28,34,39,45,51,57,63,69,75,81,87,93,99,105; 0,2,7,12,18,24,30,36,42,48,55,61,67,74,80,86,93,99,106,112; 0,2,7,13,19,25,32,38,45,52,58,65,72,78,85,92,99,106,113,119; 0,2,8,13,20,27,34,41,48,55,62,69,76,83,90,98,105,112,119,127}. - COMPUTE d99={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0; 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0; 0,0,0,0,0,0,0,0,0,0,0,1,1,1,2,2,2,2,3,3; 0,0,0,0,0,0,0,1,1,2,2,3,3,4,5,5,6,6,7,8; 0,0,0,0,0,1,1,2,3,4,5,6,7,7,8,9,10,11,12,13; 0,0,0,0,1,2,3,4,5,6,7,9,10,11,12,13,15,16,17,18; 0,0,0,0,1,3,4,6,7,9,10,12,13,15,16,18,19,21,22,24; 0,0,0,1,2,4,6,7,9,11,13,15,17,18,20,22,24,26,28,30; 0,0,0,1,3,5,7,9,11,13,16,18,20,22,24,27,29,31,33,36; 0,0,0,2,4,6,9,11,13,16,18,21,24,26,29,31,34,37,39,42; 0,0,0,2,5,7,10,13,16,18,21,24,27,30,33,36,39,42,45,48; 0,0,1,3,6,9,12,15,18,21,24,27,31,34,37,41,44,47,51,54; 0,0,1,3,7,10,13,17,20,24,27,31,34,38,42,45,49,53,57,60; 0,0,1,4,7,11,15,18,22,26,30,34,38,42,46,50,54,58,63,67; 0,0,2,5,8,12,16,20,24,29,33,37,42,46,51,55,60,64,69,73; 0,0,2,5,9,13,18,22,27,31,36,41,45,50,55,60,65,70,74,79; 0,0,2,6,10,15,19,24,29,34,39,44,49,54,60,65,70,75,81,86; 0,0,2,6,11,16,21,26,31,37,42,47,53,58,64,70,75,81,87,92; 0,0,3,7,12,17,22,28,33,39,45,51,57,63,69,74,81,87,93,99; 0,0,3,8,13,18,24,30,36,42,48,54,60,67,73,79,86,92,99,105}. - COMPUTE u95=d95(n(1),n(2)). - COMPUTE u99=d99(n(1),n(2)). - PRINT /TITLE='Точные доверительные интервалы для малых выборок (N(a) и N(b) =< 20)'. - DO IF u95 EQ 0. - PRINT /TITLE'Внимание: объёмы выборок слишком малы для надёжного расчёта 95% ДИ.'. - COMPUTE u95=1. /* Замена 0 на 1 для предотвращения ошибки вычисления. - END IF. - DO IF u99 EQ 0. - PRINT /TITLE='Внимание: объёмы выборок слишком малы для надёжного расчёта 99% ДИ.'. - COMPUTE u99=1. /* Замена 0 на 1 для предотвращения ошибки вычисления. - END IF. - RELEASE d95,d99. ELSE. /* Асимптотические критические значения Мана-Уитни *. - COMPUTE u95=TRUNC(n1n2/2-0.5-1.960*sqrt(n1n2*(n(1)+n(2)+1)/12)). - COMPUTE u99=TRUNC(n1n2/2-0.5-2.576*sqrt(n1n2*(n(1)+n(2)+1)/12)). - PRINT /TITLE'Asymptotic confidence intervals calculated (N(a) and/or N(b) > 20)'. END IF. COMPUTE low95=sdiff(u95). COMPUTE high95=sdiff(n1n2-u95). COMPUTE low99=sdiff(u99). COMPUTE high99=sdiff(n1n2-u99). PRINT {low95,high95;low99,high99} /FORMAT='f8.1' /TITLE='ДИ для разности двух медиан' /RLABELS='Ур. 95%' 'Ур. 99%' /CLABELS='НГ' 'ВГ'. ELSE. PRINT /TITLE='Объём по крайней мере одной выборки = 1. Доверительный интервал вычислить нельзя.'. END IF. END MATRIX.