Статистически устойчивый анализ данных тест Манна-Уитни-Уилкоксона и Score функции

В настоящей статье развиваются идеи и распространяются методы, изложенные в прошлой публикации Статистически устойчивый анализ данных: тест Уилкоксона на случай двух выборок. Это простая, но широко используемая на практике модель, так как даже в более сложных ситуациях целевые показатели часто сопоставляются на двух уровнях.

Анализ модели о сдвиге параметров положения двух генеральных совокупностей начинается с описания свободной от распределения ранговой процедуры Манна-Уитни-Уилкоксона (Mann-Whitney-Wilcoxon, MWW), здесь строятся точечные и интервальные оценки для величины сдвига. Далее кратко описывается метод анализа, основанный на применении score функций и, с его помощью, также проверяется нулевая гипотеза о величине параметре сдвига. В заключение, модель для параметра положения формулируется в виде регрессионной задачи, решение которой также позволяет построить точечную и интервальную оценки для параметра сдвига.

Все изложенные в статье методы проиллюстрированы на сквозном примере, реализованном в виде алгоритмов на языке R.

1.Пустьи две непрерывные случайные величины: F(t) и f(t) обозначают функцию (cdf) и плотность (pdf) распределения случайной величины, а G(t) и g(t) обозначают соответственно функцию (cdf) и плотность (pdf) случайной величины. Будем говорить, чтои следуют модели для параметра положения (location model), если для некоторого параметра $\Delta$ , $-\infty<\Delta<\infty$ имеем

$G(t)=F(t-\Delta),~~~~g(t)=f(t-\Delta).$

Параметр $\Delta$ это сдвиг параметра положения случайных величини, например, это может быть разница между медианами или средними (в случае, если средние существуют). Отметим, что предложенная модель предполагает равенство параметров масштаба случайных величини.

2.Рассмотрим две независимые друг от друга выборки, извлеченные в соответствии с выбранной моделью. Пусть $X_1,\ldots,X_{n_1}$ случайная выборка из генеральной совокупности с законом распределения(с функциями cdf и pdf F(t) и f(t) соответственно), и $Y_1,\ldots,Y_{n_2}$ случайная выборка из генеральной совокупности с законом распределения(с функциями cdf и pdf $F(t-\Delta)$ и $f(t-\Delta)$ соответственно). Пусть также n=n_1+n_2 размер объединённой выборки $X_1,\ldots,X_{n_1},Y_1,\ldots,Y_{n_2}$ . Рассмотрим гипотезу

$H_0:\Delta=0,~~~H_a:\Delta\neq0.$

При этом односторонние альтернативы также допустимы.

Далее изложены методы проверки (точный и асимптотический) данной гипотезы, а также обсуждаются подходы для получения оценки параметра сдвига $\Delta$ : точечной оценки и доверительного интервала.

3.Упорядочим элементы выборки по возрастанию, то есть запишем её в виде вариационного ряда. Каждому элементу ряда поставим в соответствие его номер в ряду ранг. Если несколько элементов ряда совпадают по величине, то каждому из них присваивается ранг, равный среднему арифметическому их номеров. Последний элемент в ранжированной выборке изэлементов должен иметь ранг. Этот факт можно использовать при проверке правильности ранжирования.

> z <- c(12, 18, 11, 5, 11, 5, 11, 11)> rank(z)[1] 7.0 8.0 4.5 1.5 4.5 1.5 4.5 4.5

Пусть R(Y_i) означает ранг Y_i среди элементов объединенной выборки, то есть среди $X_1,\ldots,X_{n_1},Y_1,\ldots,Y_{n_2}$ . Тогда статистика теста Уилкоксона имеет вид

$T=\sum_{i=1}^{n_2}R(Y_i)$

На практике часто используется эквивалентнаяследующая статистика Манна-Уитни-Уилкоксона (Mann-Whitney-Wilcoxon, MWW). Рассмотрим набор из всех $n_1\cdot n_2$ разностей вида $\left\{Y_j-X_i\right\}$ и пусть T^+ число положительных таких разностей, то есть

$T^+=\#_{i,j}\{Y_j-X_i>0\}.$

В этом случае, верно равенство

$T^+=T-\frac{n_2(n_2+1)}{2}.$

Очевидно, что исследуемую двустороннюю гипотезу H_0 следует отвергнуть при малых и больших значениях статистики. При условии $H_0:\Delta=0$ обе выборки получены из одной и той же генеральной совокупности, поэтому любые наборы рангов равновероятны (например, вероятность того, что подмножество рангов n_2 относится к случайной величинеравна $1/C_n^{n_2}$ ). Следовательно, распределение статистикипри нулевой гипотезе H_0 не зависит (свободно) от распределения генеральной совокупности и точное значение p-value рассчитывается на основе распределения(точном распределении рангов).

4.Оценкой параметра сдвига $\Delta$ в тесте Манна-Уитни-Уилкоксона является медиана выборки, составленной из всех $N_d=n_1\cdot n_2$ парных разностей вида (оценка Ходжеса-Лемана (Hodges-Lehmann))

$\hat{\Delta}_W=\mbox{med}_{i,j}\{Y_j-X_i\}.$

Пусть $D_{(1)}<D_{(2)}<\cdots<D_{(N_d)}$ вариационный ряд для таких разностей. Если доверительная вероятность равна $1-\alpha$ и квантиль уровня $\alpha/2$ распределения T^+ , то есть $\alpha/2=P_{H_0}[T^+\leq c]$ , то интервал $\left(D_{(c+1)},D_{(N_d-c)}\right)$ является $(1-\alpha)100\%$ доверительным интервалом для $\Delta$ . Асимптотическим значение дляявляется величина

$c=\frac{n_1n_2-1}{2}-z_{\alpha/2}\sqrt{\frac{n_1n_2(n+1)}{12},}$

которая округляется до ближайшего целого и в общем случае довольно близка к фактическому значению.

5.Для примера применим тест Манна-Уитни-Уилкоксона к двум выборкам из t-распределения cстепенями свободы и значением параметра сдвига $\Delta=8$ .

> x <- round(rt(11, 5) * 10 + 42, 1)> y <- round(rt(9, 5) * 10 + 50, 1)> x [1] 76.6 41.0 59.3 34.9 29.1 45.0 42.6 31.1 32.4 52.5 47.9> y [1] 58.3 47.2 40.1 45.8 62.0 58.7 64.8 48.1 49.5> wilcox.test(y, x, exact = TRUE, conf.int = TRUE, conf.level = 0.95)Wilcoxon rank sum exact testdata:  y and xW = 72, p-value = 0.09518alternative hypothesis: true location shift is not equal to 095 percent confidence interval: -1.0 18.4sample estimates:difference in location                   10.4

Здесь: T^+=72 с p-value 0.09518 доверительный интервал (-1,18.4) включает истинное значение $\Delta=8$ , точечная оценка равна $\hat{\Delta}_W=10.4$ , доверительная вероятность $95\%$ . По умолчанию p-value основано на точном распределении статистики T^+ для малых выборок (n<50) без хвостов. В случае выбора аргумента exact = FALSE и correct = FALSE (без корректировки на непрерывность) используется асимптотический метод, изложенный далее. В этом случае значение p-value равно 0.08738 .

> wilcox.test(y, x, exact = FALSE, correct = FALSE)Wilcoxon rank sum testdata:  y and xW = 72, p-value = 0.08738alternative hypothesis: true location shift is not equal to 0

6.Рассмотрим функцию $a_\varphi(i)=\varphi(i/(n+1))$ , где функция $\varphi(u)$ (score функция) определена на интервале (0,1) и

$\int_0^1\varphi(u)du=0,~~~\int_0^1\varphi^2(u)du=1.$

Например, $\varphi_{Ns}(u)=\Phi^{-1}(u)$ , где $\Phi^{-1}(u)$ функция, обратная к функции cdf стандартной нормальной случайной величины N(0,1) . В этом случае преобразование $a_{Ns}$ посредством функции $\varphi_{Ns}$ (Normal score function) ставит в соответствие рангам элементов выборки значения, которые ожидаемы тем, как если бы данные были получены из нормального распределения. Заметим, что в данном определении термин normal score употребляется в смысле rankit, а не в смысле standard score или z-score. Наряду с функцией normal score, возможны другие виды score функций, например $\varphi_W(u)=\sqrt{12}[u-(1/2)]$ задаёт score функцию Уилкоксона.

Для выбранной модели для параметра положения определим следующую функцию переменной $\Delta$ :

$S_\varphi(\Delta)=\sum_{j=1}^{n_2}a\varphi[R(Y_j-\Delta)],$

где $a_\varphi$ некоторая score функция, $R(Y_j-\Delta)$ ранги $Y_j-\Delta$ среди $X_1,\ldots,X_{n_1}$ и $Y_1-\Delta,\ldots,Y_{n_2}-\Delta$ . В этом случае, статистику $S_\varphi=S_\varphi(0)$ можно использовать для проверки интересующей нас гипотезы:

$H_0:\Delta=0,~~~H_a:\Delta>0.$

При условии нулевой гипотезы H_0 , случайные величиныираспределены одинаково, отсюда следует, что распределение статистики $S_\varphi$ независимо (свободно) от общего распределения генеральной совокупности.

Хотя точное распределение $S_\varphi$ и может быть получено численно, обычно используют асимптотическое распределение. В условиях нулевой гипотезы H_0 распределение $S_\varphi(0)$ асимптотически нормально, с математическим ожиданием и дисперсией соответственно:

$E_{H_0}[S\varphi(0)]=0,~~~\sigma^2_\varphi=Var_{H_0}[S_\varphi(0)]=\frac{n_1n_2}{n(n-1)}\sum_{i=1}^na_\varphi^2(i).$

Это позволяет использовать стандартизованную статистику $z_\varphi=S_\varphi(0)/\sigma_\varphi$ чтобы отвергнуть гипотезу H_0 на уровне значимости $\alpha$ , если $z_\varphi\geq z_\alpha$ , где $z_\alpha$ квантиль уровня $(1-\alpha)$ стандартного нормального распределения. Двусторонняя и левосторонняя гипотезы проверяются аналогично.

7.В следующей R сессии для данных из рассматриваемого в статье примера представлен расчёт асимптотической статистики $z_\varphi$ и соответствующее ей значение p-value с использованием score функции Уилкоксона (пакет Rfit).

> x <- c(76.6, 41.0, 59.3, 34.9, 29.1, 45.0, 42.6, 31.1, 32.4, 52.5, 47.9)> y <- c(58.3, 47.2, 40.1, 45.8, 62.0, 58.7, 64.8, 48.1, 49.5)> # Объединим выборки x и y>   z = c(x, y)>   n1 = length(x)>   n2 = length(y)>   n = n1 + n2> # Выберем в качестве score функции функцию Уилкоксона>   scores = Rfit::wscores > # Найдём соответствующие score значения для рангов выборки z >   rs = rank(z)/(n + 1)>   asg = Rfit::getScores(scores, rs)> # Найдём значение статистики Sphi при нулевой гипотезе  >   Sphi = sum(asg[(n1 + 1):n])> # Найдём дисперсию Sphi>   asc = Rfit::getScores(scores, 1:n/(n + 1))>   varphi = ((n1 * n2)/(n * (n - 1))) * sum(asc^2)> # Найдём значения zphi и p-value >   zphi = Sphi/sqrt(varphi)>   alternative = "two.sided">   pvalue <-+     switch(+     alternative,+     two.sided = 2 * (1 - pnorm(abs(zphi))),+     less = pnorm(zphi),+     greater = 1 - pnorm(zphi)+   )> # Выведем результаты>   res <- list(Sphi = Sphi, statistic = zphi, p.value = pvalue)>   with(res, cat("statistic = ", statistic, ", p-value = ", p.value, "\n"))statistic =  1.709409 , p-value =  0.08737528

Таким образом, результаты точного T^+=72 с p-value 0.0952 и асимптотического z_W=1.71 с p-value 0.0874 анализа крайне близки: нулевая гипотеза не отвергается на уровне значимости $5\%$ и существует значимая разница в сдвиге параметра положения при уровне значимости $10\%$ .

8.Cформулируем двухвыборочную задачу для параметра положения в виде регрессионной задачи. Пусть $\bar{Z}=(X_1,\ldots,X_{n_1},Y_1,\ldots,Y_{n_2})^T$ , $\bar{c}$ это $n\times1$ вектор с-ым элементомдля $1\leq i\leq n_1$ идля $n_1+1\leq i\leq n=n_1+n_2$ . Тогда модель для параметра положения можно записать в виде

$Z_i=\alpha+c_i\Delta+e_i,$

где $e_1,\ldots,e_n$ независимые, одинаково распределенные случайные величины с плотностью распределения f(t) . Таким образом, подобрав модель регрессии можно оценить параметр сдвига $\Delta$ . При использовании метода наименьших квадратов МНК-оценкой для $\Delta$ будет разность $\bar{Y}-\bar{X}$ . При использовании метода, основанного на рангах с использованием score функции Вилкоксона, оценкой для $\Delta$ будет уже рассмотренная выше оценка Ходжеса-Лемана медиана парных разностей.

В следующей R сессии представлены результаты ранговой регрессии для имеющихся данных.

> z = c(x, y)> ci <- c(rep(0, n1), rep(1, n2))> fit <- Rfit::rfit(z ~ ci, scores = Rfit::wscores)> coef(summary(fit))            Estimate Std. Error  t.value      p.value(Intercept)     41.8   4.400895 9.498068 1.960951e-08ci              10.4   5.720594 1.817993 8.574801e-02

Таблица с результатами, наряду с самой оценкой 10.4 содержит её стандартную ошибку 5.72 . Используя это, можно, например, построить приблизительный $95\%$ доверительный интервал для сдвига $\Delta$ , основанный на квантиле уровня 1-0.05/2 t-распределения с n-2 степенями свободы:

> conf.level <- 0.95> estse <- coef(summary(fit))[2, 1:2]> alpha <- 1 - conf.level> alternative = "two.sided"> tcvs <- switch(+   alternative,+   two.sided = qt(1 - alpha / 2, n - 2) * c(-1, 1),+   less = c(-Inf, qt(1 - alpha, n - 2)),+   greater = c(qt(alpha, n - 2), Inf)+ )> conf.int <- estse[1] + tcvs * estse[2]> cat(100 * conf.level, " percent confidence interval:\n", conf.int)95  percent confidence interval: -1.618522 22.41852

Построенные доверительный интервал (-1.62,22.42) несколько шире своего дискретного аналога (-1,18.4) , найденного ранее.

Мы рассмотрели два подхода (процедура Манна-Уитни-Уилкоксона и метод score функций) для проверки гипотезы о наличии сдвига у параметров положения двух генеральных совокупностей. Для величины сдвига построены точечная оценка и два типа доверительных интервала.

Работая в компании IT-аутсорса в качестве руководителя 3 линии поддержки, задумался, как автоматизировать подключение сотрудников по RDP, через VPN к серверам десятков клиентов.

Таблички с адресами, паролями и прочими настройками серверов, конечно, хорошо, но поиск клиента и вбивание адресов с аккаунтами занимает довольно существенное время.
Держать все подключения к VPN в Windows не самая лучшая идея, да и при переустан

Всем привет! Меня зовут Константин Измайлов, я руководитель направления Data Science в Delivery Club. Мы работаем над многочисленными интересными и сложными задачами: от формирования классических аналитических отчетов до построения рекомендательных моделей в ленте приложения.

Сегодня я расскажу пр

Представьте человека, который изучает алгоритмы. Чтобы понять как они работают, приходится изучать их код и представлять, как компьютер будет его выполнять. Это странно почему мы должны учиться думать как компьютер, вместо того, чтобы заставить его помогать на

Cодержание

Введение
Основные понятия и термины
Характеристика ИО как научной дисциплины
Этапы операционного исследования
- Постановка задачи
- Построени

Можно выделить ряд алгоритмов, которые являются базовыми и лежат в основе практически каждой строчки программ, написанных на языках высокого уровня. Хорошо иметь под руками классический многотомный труд Дональда Кнута "The Art of Computer Programming", там детально разобраны многие ба

Давайте представим, что вы параноик, и параноик вдвойне, когда дело касается многопоточности. Предположим, что вы делаете backend некого функционала приложения, а приложение переодически дергает на вашем серверы какие-то методы. Все вроде хорошо, но есть одно но. Что если ваш функционал напрямую зависит от каких-либо других данных, того же банального профиля например? Встает вопрос, как гарантировать то, что сценарий отработает име

Трансформеры за последние несколько лет штурмом захватили мир NLP, а сегодня они с успехом применяются в выходящих за рамки NLP приложениях. Они обладают такими возможностями благодаря модулю внимания, который схватывает отношения между всеми словами последовательностей. Но са

Пример расчётного доказательства в Lean

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

Зарабатывать продажей лекарств, которые заведомо не работают, не только аморально, но и не особо легко. Люди всё-таки обычно не хотят покупать препараты, неэффективность которых была доказана. А вот если вы сумели выдавить заветное p < 0.05 в пользу того, что акупунктура та

В этой статье мы рассмотрим, что такое классификатор, поговорим о мультиклассовой классификации с помощью нейронных сетей. Затем, ознакомившись с контекстом перейдем к основному топику поста к Log-Sum-Exp Trick. Напишем формулы и разберемся, как этот трюк помогает избежать переполнения чисел с плавающей то

В дополнение к открытым спутниковым данным, некоторые из которых перечислены в статье Общедоступные данные дистанционного зондирования Земли: как получить и использовать, существует и множество производных продуктов например, рельеф. Притом можно найти открытый рельеф разного пространственного разрешения, равно как и множество коммерческих, и появляется задача выбрать лучший продукт из доступных.

Призма Вельда-Бланделла

На рубеже четвертого и третьего тысячелетия до нашей эры на Земле возникли две первые цивилизации. В долине Нила после объединения верхнего и нижнего Египта образовалось

Зачем нужен Rust в R?

Первый вопрос, который должен возникнуть у читателя -- а зачем вообще использовать Rust вместе с R? Ответ довольно прост: Rust -- новый системный язык программирования, спроектированный специально для написания безопасного и легко распараллеливаемого кода. Rust довольно сложен в освоении (в сравнении с другими языками), но при этом предоставляет отли

Основной из проблем кластерного анализа практических данных является наличие выбросов. Большинство существующих методов кластеризации не учитывают их существование, из-за этого явно аномальные наблюдения включаются в состав каких-то кластеров, что может серьезно смещать их центры и влиять на качество классификации. Разумеется, можно сначала исходные данные проверить на выбросы, их отсеять и т.д., но тогда задача превратиться в двух

В этой статье, переводом которой мы решили поделиться специально к старту курса о Data Science

В ходе обсуждений возникла маленькая задачка построить динамику структуры кредитного портфеля (динамика кредитной карты, например). В качестве важной специфики необходимо применять метод FIFO для погашения займов. Т.е. при погашении первыми должны гаситься самые ранние займы. Это накладывает определенные требования на расчет статуса каждого отдельного займа и определения его даты погашения.

Ниже приведе

Проблематика

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

Какие есть варианты?

Сравнивая различные инструменты и таскменеджеры, выбор пал на одного из популярных игроков рынка - Trello. Однако, сразу после выбора программы, возникло новое препятствие. В бесплатной версии этого приложения

Анализ модели о сдвиге параметров положения двух генеральны

В настоящей статье развиваются идеи и распространяются методы, изложенные в прошлой публикации Статистически устойчивый анализа данных: Тест Уилкоксона на случай двух выборок. Это простая, но широко используемая на практике модель, так как даже в более сложных ситуациях целевые показатели часто сопоставляются на двух уровнях.