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

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

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

20.03.2021 21:27

ampopov 2 1200 Источник

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

Анализ модели о сдвиге параметров положения двух генеральных совокупностей начинается с описания свободной от распределения ранговой процедуры Манна-Уитни-Уилкоксона (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 test

data:  y and x
W = 72, p-value = 0.09518
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
 -1.0 18.4
sample 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 test

data:  y and x
W = 72, p-value = 0.08738
alternative 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-08
ci              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 функций) для проверки гипотезы о наличии сдвига у параметров положения двух генеральных совокупностей. Для величины сдвига построены точечная оценка и два типа доверительных интервала.

Комментарии (2)

excoder
21.03.2021 01:09
#22830322
Очень интересная и важная тема, спасибо за серию статей. Хотелось бы увидеть более живого изложения, с примерами – когда это применяют, в каких контекстах. В каком смысле вообще оно статистически устойчиво или неустойчиво. В бытность студентом на вес золота были преподаватели, которые цепляли сложными темами безвозвратно. Яркий пример – лекции Босса. Дисклэймер: я сам по образованию математик и применяю статистику на практике.
1. ampopov Автор
  21.03.2021 02:07
  #22830390
  Здравствуйте,
  
  Цель, которую я перед собой поставил при подготовке цикла статей о непараметрических методах статистики – это доступно, так чтобы это было понятно вдумчивому читателю с высшим техническим или экономическим образованием, разъяснить математические предпосылки, лежащие в основе тех или иных робастных статистических процедур.
  
  На мой взгляд, такого рода информации не хватает, хотя именно её наличие обуславливает обдуманное использование статистических методов на практике. В качестве конкретных примеров я привожу алгоритмы программ, изучив которые можно лучше понять, как получены те или иные статистические результаты и какова область их реального применения.
  
  Там, где это возможно я стараюсь приводить примеры использования рассматриваемых моделей в повседневной жизни. Например, в статье «Статистически устойчивый анализа данных: Тест Уилкоксона» однопараметрический тест Уилкоксона применяется для сравнительного ответа на вопрос о величине спроса в двух магазинах на основе данных о продажах.

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

Комментарии (2)

excoder

ampopov Автор