Непараметрический метод сравнения выборок, По медиане и квартильному размаху Опции

[Диагностик]

Здравствуйте, уважаемые. Существуют две несвязанные выборки, по которым получены: n, Q1, Q2, Q3. Q2 - медиана. Других данных нет. Как по ним проверить гипотезу о значимости различия выборок? Спасибо.

[Диагностик]

Собственно практическая задача такая. Известны выборочные параметры распределения количества лейкоцитов для группы здоровых людей объёмом 54:
Q1=5.6; Me=6.8; Q3=9.3.
Тоже самое для группы больных людей объёмом 37:
Q1=7.0; Me=7.9; Q3=9.7.
Ничего нельзя сказать о значимом отличии групп по этим данным?

Сообщение отредактировал Диагностик - 11.04.2015 - 07:40

[nokh]

Собственно практическая задача такая. Известны выборочные параметры распределения количества лейкоцитов для группы здоровых людей объёмом 54:
Q1=5.6; Me=6.8; Q3=9.3.
Тоже самое для группы больных людей объёмом 37:
Q1=7.0; Me=7.9; Q3=9.7.
Ничего нельзя сказать о значимом отличии групп по этим данным?

Такие данные можно сравнить с помощью критерия Колмогорова, который находит различия в распределениях. Для этого нужно отложить значения квартилей на кривой накопленных частот (CDF). Максимальное значение разности D между этими кривыми (в нашем случае они достаточно грубые - ломаные) используется далее в расчёте статистики критерия: формулы см. http://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test в разделе Two-sample Kolmogorov-Smirnov test. Здесь возможны 2 ситуации:

а) Расположение кривых позволяет найти D (рис. а). Вычисления не выходят за пределы заштрихованной области. Максимальные значения высот этой этой фигуры (D) всегда или от Q1 вверх или от Q3 вниз, по крайней мере на тех фигурах, что я накидал за 3 мин было так. Видимо это можно доказать геометрически и алгоритмизировать поиск D. Я думал ситуация а) - это для случая перекрывающихся интервалов, но ваш пример показывает, что нет - и при перекрытии возможна ситуация б).

б) Расположение кривых позволяет найти не D, а его нижнюю границу. Т.к. минимумы и максимумы не известны, реальное D тоже не известно: зелёная стрелка упирается в пунктир. Такая ситуация всегда будет при неперекрывании межквартильных размахов. Для такой ситуации если бы были известны минимумы и максимумы, то вместо пунктира продолжилась бы ломаная CDF, что сделало бы определение D также однозначным.

Ваша ситуация попадает в категорию б), т.к. возможно, что реальное D будет больше - пунктирная зелёная стрелка с вопросом. Но за неимением большей информации - работаем с тем, что есть.
1) По значениям для группы здоровых (синие точки) находим уравнение прямой, проходящей через медиану Ме (6,8; 50) и Q3 (9,3; 75): у=10х-18.
2) Находим ординату точки пересечения прямой D c этой прямой в точке х=7. у=10х7-18=52. Находим D как 52-25=27, 27/100%=0,27.
3) Подставляем это значение в формулу и находим р интерполяцией по табличным значениям с википедии (интерполировал полиномом 4-ой степени). Для ваших данных р=0,081. Я трактую это как "тенденцию к различиям", хотя знаю что многие ругают такую формулировку.
Вот такой мой подход. В "материале и методах" можно писать, что использовали критерий Колмогорова. Нужно сказать, что реальное значение р, возможно, будет меньше, т.к. мы смогли найти только нижнюю границу D. Если реальное D будет больше нашего, значит и С(альфа) будет больше, а р - меньше. А может и не будет больше. Поэтому написал р <=0,081. В общем, критерий Колмогорова D_{(37; 54)}=0,27; р <=0,081.

И ещё ремарка по поводу "грубости" такого подхода. Вообще говоря, при вычислении CDF и далее статистики критерия выбор интервала группировки классов строго не регламентируется. Кто-то работает "на глазок", кто-то по формуле Стургеса, кто-то по EM-алгоритму. То, что мы берём такие аршинные интервалы как квартили, конечно не очень хорошо, но почему бы и нет, если это позволяет принимать решения. В данном случае решение о различии распределений.

PS. В подходе р2004r не понял откуда выбирать значения, если интервал min-max не определён. Был бы признателен за код с вашим вариантом р.

Сообщение отредактировал nokh - 12.04.2015 - 07:28

Эскизы прикрепленных изображений

 

[p2004r]

PS. В подходе р2004r не понял откуда выбирать значения, если интервал min-max не определён. Был бы признателен за код с вашим вариантом р.

1. Я сразу сказал, что берусь оценить доверительный интервал только для положения медианы. Подход точно такой же как и в вашем случае, часть выборки нам известно в каких границах лежит. Для этой части мы принимаем предположение (как и у вас) что всё линейно и делаем генерацию точек из униформного распределения с границами заданными квантиляли и медианой. Для получившейся ситуации считаем положение медианы, накапливая такие перевыборок получаем доверительный интервал .

2.

Код

# Восстановим сколько скорее всего попало при расчете в конкретный квантиль
> table(cut(1:54, quantile(1:54),5670252.lowest=T))

   [1,14.2] (14.2,27.5] (27.5,40.8]   (40.8,54]
         14          13          13          14

> table(cut(1:37, quantile(1:37),5670252.lowest=T))

[1,10] (10,19] (19,28] (28,37]
     10       9       9       9

# Перевыборка медианы происходит вот таким образом

> median(c(runif(9, min=7.0, max=7.9), runif(9, min=7.9, max=9.7)))
[1] 7.855721
> median(c(runif(9, min=7.0, max=7.9), runif(9, min=7.9, max=9.7)))
[1] 7.901809
> median(c(runif(9, min=7.0, max=7.9), runif(9, min=7.9, max=9.7)))
[1] 7.775273
> median(c(runif(9, min=7.0, max=7.9), runif(9, min=7.9, max=9.7)))
[1] 8.010139
> median(c(runif(9, min=7.0, max=7.9), runif(9, min=7.9, max=9.7)))
[1] 7.856393

Доверительный интервал для второй выборки

> quantile(replicate(10000, median(c(runif(9, min=7.0, max=7.9), runif(9, min=7.9, max=9.7)))), probs=c(0.025, 0.5, 0.975))
    2.5%      50%    97.5%
7.796376 7.927046 8.167319
> quantile(replicate(100000, median(c(runif(9, min=7.0, max=7.9), runif(9, min=7.9, max=9.7)))), probs=c(0.025, 0.5, 0.975))
    2.5%      50%    97.5%
7.795279 7.928610 8.165428

Для первой выборки

> quantile(replicate(10000, median(c(runif(13, min=5.6, max=6.8), runif(13, min=6.8, max=9.3)))), probs=c(0.025, 0.5, 0.975))
    2.5%      50%    97.5%
6.697030 6.828900 7.062315
> quantile(replicate(100000, median(c(runif(13, min=5.6, max=6.8), runif(13, min=6.8, max=9.3)))), probs=c(0.025, 0.5, 0.975))
    2.5%      50%    97.5%
6.698919 6.828830 7.073703

Медианы генсовокупностей из которых были извлечены выьорки не имеют шанса встретиться если извлечение выборок шло случайно.

Оценка получается более [s]широкой[\s] узкой чем просто бутстреп исходной выборки

Код

# модельная генсовокупность
> x<-rnorm(54, mean=6.8, sd=1.3)
> quantile(x)
      0%      25%      50%      75%     100%
2.716833 6.051314 6.823896 7.658319 9.626783
> quantile(replicate(100000, median(c(runif(13, min=6.051, max=6.823), runif(13, min=6.823, max=7.65)))), probs=c(0.025, 0.5, 0.975))
    2.5%      50%    97.5%
6.748091 6.824059 6.905382

# генсовокупность порождает вот такие варианты выборки в пределе
> quantile(replicate(10000, median(rnorm(54, mean=6.8, sd=1.3)) ), probs=c(0.025, 0.5, 0.975))
    2.5%      50%    97.5%
6.368884 6.799687 7.231649
> quantile(replicate(100000, median(rnorm(54, mean=6.8, sd=1.3)) ), probs=c(0.025, 0.5, 0.975))
    2.5%      50%    97.5%
6.371881 6.800581 7.228905
> quantile(replicate(100000, median(sample(x, )) ), probs=c(0.025, 0.5, 0.975))
x=        size=     replace=  prob=    

# настоящий бутстреп медианы выборки
> quantile(replicate(100000, median(sample(x, replace=T)) ), probs=c(0.025, 0.5, 0.975))
    2.5%      50%    97.5%
6.438816 6.823896 7.119538

И поправка плывет от размера sd (и скорее всего вида распределения), хотя и не зависит похоже от размера выборки. Похоже раз такая зависимость есть, то проще фитить в эти процентили какое то семейство распределений и считать уже по нему различия.

Код

> res.butstr <- sapply(20:100, function(n) {x<- rnorm(n, mean=6.8, sd=3.3); q<- quantile(x); qn<- table(q); quantile(replicate(10000, median(c(runif(qn[2], min=q[2], max=q[3]), runif(qn[3], min=q[3], max=q[4]))) ), probs=c(0.025, 0.5, 0.975)) / quantile(replicate(10000, median(sample(x, )) ), probs=c(0.025, 0.5, 0.975)) })
> plot(20:100, t(res.butstr)[,1], ylim=range(as.vector(res.butstr)))
> points(20:100, t(res.butstr)[,2], col="green")
> points(20:100, t(res.butstr)[,3], col="red")

Сообщение отредактировал p2004r - 12.04.2015 - 22:43