R:Статистическая проверка принадлежности нормальному распределения
 
|
Материал «R:Статистическая проверка принадлежности нормальному распределения», созданный автором Артём Клевцов, публикуется на условиях лицензии Creative Commons «Attribution» («Атрибуция») 4.0 Всемирная. | |
|
Предположение о принадлежности случайной величины нормальному закону распределения лежит в основе многих статистических методов и критериев. В ряде случаев соблюдение данного требования является критичным для применения того или иного метода.
На практике мы встречаемся с двумя вариантами задач по проверке принадлежности распределения нормальному закону: для одномерного и многомерного распределения.
 |
Перед использованием функций из пакетов их необходимо предварительно установить и загрузить: КодR <syntaxhighlight lang="r">> install.packages(pkgs = "pkgname") > library(package = "pkgname")</syntaxhighlight> |
Содержание
Одномерное нормальное распределение
В качестве [math]H_0[/math] для всех нижеприведённых критериев является предположение, что «случайная величина [math]X[/math] распределена нормально».
Для демонстрации работы функций, реализующих различные критерии проверки принадлежности распределения нормальному закону сгенерируем вектор случайных чисел, имеющих стандартное нормальное распределение:
<syntaxhighlight lang="r">> x <- rnorm(n = 1000) </syntaxhighlight>
Статистические критерии
В R реализовано множество критериев проверки соответствия распределения нормальному закону.
Сравнительная таблица реализации критериев в пакетах
| Критерии | stats
|
nortest
|
moments
|
fBasics
|
tseries
|
lawstat
|
|---|---|---|---|---|---|---|
| Критерий Шапиро - Уилка | + | - | - | + | - | - |
| Критерий Колмогорова - Смирнова | + | - | - | + | - | - |
| Критерий Андерсона - Дарлинга | - | + | - | + | - | - |
| Критерий Крамера - фон Мизеса | - | + | - | + | - | - |
| Критерий Лиллиефорса | - | + | - | + | - | - |
| Критерий [math]\chi^2[/math] Пирсона | - | + | - | + | - | - |
| Критерий Шапиро - Франчия | - | + | - | + | - | - |
| Критерий Д'Агостино | - | - | + | + | - | - |
| Критерий Бонетта – Сайера | - | - | + | - | - | - |
| Критерий Жарка - Бера | - | - | + | + | + | + |
Пакет fBasics содержит также функцию normalTest(), которая является «обёрктой» для ряда функций из того же пакета. Необходимый критерий можно задать с помощью аргумента method. Доступны следующие критерии:
-
sw- критерий Шапиро - Уилка -
jb- критерий Жарка-Бера -
ks- критерий Колмогорова - Смирнова -
da- критерий Д'Агостино -
ad- критерий Андерсона - Дарлинга.
Пример вызова данной функции:
<syntaxhighlight lang="r">> normalTest(x, method = "sw") Title: Shapiro - Wilk Normality Test Test Results: STATISTIC: W: 0.9831 P VALUE: 0.2301 Description: Fri Feb 14 19:59:59 2014 by user: </syntaxhighlight>
Пакет TeachingDemos содержит функцию SnowsPenultimateNormalityTest(), реализующую неописанный в литературе критерий. Данная функция возвращает только уровень статистической значимости, свидетельствующий об отклонения распределения от нормального закона.
Таблица вызова функций в пакетах
| Критерии | stats
|
nortest
|
moments
|
fBasics
|
tseries
|
lawstat
|
|---|---|---|---|---|---|---|
| Критерий Шапиро - Уилка | shapiro.test |
- | - | shapiroTest |
- | - |
| Критерий Колмогорова - Смирнова | ks.test[1] |
- | - | ksnormTest |
- | - |
| Критерий Андерсона - Дарлинга | - | ad.test |
- | adTest |
- | - |
| Критерий Крамера - фон Мизеса | - | cvm.test |
- | cvmTest |
- | - |
| Критерий Лиллиефорса | - | lillie.test |
- | lillieTest |
- | - |
| Критерий [math]\chi^2[/math] Пирсона | - | pearson.test |
- | pchiTest |
- | - |
| Критерий Шапиро - Франчия | - | sf.test |
- | sfTest |
- | - |
| Критерий Д'Агостино | - | - | agostino.test |
dagoTest |
- | - |
| Критерий Бонетта – Сайера | - | - | bonett.test |
- | - | - |
| Критерий Жарка - Бера | - | - | jarque.test |
jarqueberaTest |
jarque.bera.test |
rjb.test
|
Маленькие хитрости
Применение функций к нескольким переменным
С помощью apply-функций можно последовательно применить функцию к вектору, списку или массиву. Прежде чем всего нам необходимо сформировать таблицу данных. С помощью функции replicate() сгенерируем 10 переменных, имеющих стандартное нормальное распределение, которые объединяются в класс data.frame.
<syntaxhighlight lang="r">> DF <- data.frame(replicate(n = 10, rnorm(n = 100))) </syntaxhighlight>
Структура сгенерированной таблицы выглядит следующим образом:
<syntaxhighlight lang="r">> str(DF) 'data.frame': 100 obs. of 10 variables: $ X1 : num 1.051 1.08 -0.477 -1.396 3.423 ... $ X2 : num -0.602 2.29 -0.758 -1.615 -0.364 ... $ X3 : num 0.0559 -1.0117 0.5242 0.4105 -0.3191 ... $ X4 : num -0.0965 0.2006 0.29 0.7702 -0.0182 ... $ X5 : num -0.7074 -1.6111 0.3478 0.2504 0.0609 ... $ X6 : num -1.432 0.535 -0.932 0.581 -1.606 ... $ X7 : num -1.42407 -0.31827 -2.04648 -0.19856 0.00301 ... $ X8 : num 0.511 0.192 0.467 -1.308 2.496 ... $ X9 : num -0.8508 0.4481 -0.2828 -0.5464 0.0605 ... $ X10: num 1.421 0.408 1.254 -0.956 -1.91 ... </syntaxhighlight>
Для решения поставленной задачи можно воспользоваться функцией sapply(). Но прежде, нам необходимо немного отформатировать формат вывода результатов нашей функции: нам нужно извлечь значения критерия и его уровень значимости, т.к. результат функции shapiro.test() содержит также информацию, которая не подлежит включению в итоговую таблицу, например, информация об используемом методе (критерии) и уточнение характера альтернативной гипотезы. Вывод результатов тест Шапиро - Уилка выглядит следующим образом:
<syntaxhighlight lang="r">> shapiro.test(x) Shapiro-Wilk normality test data: x W = 0.9903, p-value = 0.6882 </syntaxhighlight>
Структура результата применения функции shapiro.test() представлена ниже:
<syntaxhighlight lang="r">> str(shapiro.test(x)) List of 4 $ statistic: Named num 0.99 ..- attr(*, "names")= chr "W" $ p.value : num 0.688 $ method : chr "Shapiro-Wilk normality test" $ data.name: chr "x" - attr(*, "class")= chr "htest" </syntaxhighlight>
Как видим, помимо значений критерия и уровня значимости здесь содержится информация о применяемом методе. Мы можем отфильтровать вывод следующим образом:
<syntaxhighlight lang="r">> normTest <- function (x) { + res <- shapiro.test(x) + return(c(res$statistic, p.value = res$p.value)) + } </syntaxhighlight>
Результат теперь будет выглядеть следующим образом:
<syntaxhighlight lang="r">> normTest(x) W p.value 0.9903 0.6882 </syntaxhighlight>
Теперь можно использовать данную функцию при обработке столбцов нашей таблицы.
<syntaxhighlight lang="r">> t(sapply(DF, normTest)) W p.value X1 0.9831 0.2301 X2 0.9936 0.9213 X3 0.9800 0.1333 X4 0.9829 0.2219 X5 0.9874 0.4625 X6 0.9862 0.3874 X7 0.9839 0.2617 X8 0.9833 0.2360 X9 0.9915 0.7834 X10 0.9808 0.1531 </syntaxhighlight>
Того же результата можно добиться и с помощью функции lapply()[2]:
<syntaxhighlight lang="r">> do.call(rbind, lapply(DF, normTest)) W p.value X1 0.9831 0.2301 X2 0.9936 0.9213 X3 0.9800 0.1333 X4 0.9829 0.2219 X5 0.9874 0.4625 X6 0.9862 0.3874 X7 0.9839 0.2617 X8 0.9833 0.2360 X9 0.9915 0.7834 X10 0.9808 0.1531 </syntaxhighlight>
Применение функций к нескольким по группам
Добавим к нашей таблице группы испытуемых:
<syntaxhighlight lang="r">> DF$GRP <- factor(sample(LETTERS[1:3], size = 100, replace = TRUE))</syntaxhighlight>
Состав групп получился следующим:
<syntaxhighlight lang="r">> table(DF$GRP) A B C 38 25 37 </syntaxhighlight>
Рассчитаем значения критерия Шапиро - Уилка для первого столбца для каждоый группы испытуемых:
<syntaxhighlight lang="r">> do.call(rbind, tapply(DF$X1, DF$GRP, normTest)) W p.value A 0.9522 0.13281 B 0.9607 0.28697 C 0.9410 0.07256</syntaxhighlight>
Графические методы
Многие исследователи также используют графические методы для определения степени отклонения распределения от нормального закона. В R реализована возможность построения Q-Q и P-P графиков, гистограмм и кривых распределения плотности вероятностей.
Гистограмма
Гистограмма представляет собой графическое изображение зависимости частоты попадания элементов выборки от соответствующего интервала группировки. Построить гистограмму в R можно с помощью следующей команды:
<syntaxhighlight lang="r">> hist(x)</syntaxhighlight>
На гистограмме изображены абсолютные частоты. Также можно построить гистограмму, отражающую плотности вероятностей:
<syntaxhighlight lang="r">> hist(x, freq = FALSE)</syntaxhighlight>
График плотностей вероятности
Пакет stats
<syntaxhighlight lang="r">> plot(density(x))</syntaxhighlight>
Пакет car
<syntaxhighlight lang="r">> densityPlot(x)</syntaxhighlight>
Гистограммы с наложением графика плотностей вероятнотси
Пакет stats
<syntaxhighlight lang="r">> hist(x, freq = FALSE) > lines(density(x))</syntaxhighlight>
Теперь наложим на наш график кривую плотностей вероятности для нормального распределения:
<syntaxhighlight lang="r">> xfit <- seq(min(x), max(x), length = 100) > yfit <- dnorm(xfit, mean = mean(x), sd = sd(x)) > hist(x, freq = FALSE) > lines(density(x), col = "red") > lines(xfit, yfit, col = "blue") </syntaxhighlight>
Пакет gamlss
Более простой способ сравнение графиков плотностей вероятности представлен в функции histDist из пакета gamlss:
<syntaxhighlight lang="r">> histDist(x, family = "NO", density = TRUE) Family: c("NO", "Normal") Fitting method: "nlminb" Call: gamlssML(y = y, family = "NO", formula = x) Mu Coefficients: [1] -0.0462 Sigma Coefficients: [1] 0.023 Degrees of Freedom for the fit: 2 Residual Deg. of Freedom 998 Global Deviance: 2884 AIC: 2888 SBC: 2898 </syntaxhighlight>
С помощью аргумента family можно задать семейство распределений для подгонки и сравнения[3].
Q-Q график
Q-Q график (Q - квантиль) — это график, на котором квантили из двух распределений расположены относительно друг друга. Чем ближе точки на графике к диагональной прямой, тем ближе распределение исследуемой переменной к нормальному закону.
Построение квантильных графиков в R реализовано в нескольких пакетах.
Пакет stats
Построение Q–Q plot с помощью пакета stats выглядит следующим образом:
<syntaxhighlight lang="r">> qqnorm(x) > qqline(x)</syntaxhighlight>
Пакет QTLRel
Построение Q–Q plot с помощью пакета QTLRel выглядит следующим образом:
<syntaxhighlight lang="r">> qqPlot(x, x = "norm")</syntaxhighlight>
Пакет car
Альтернативный вариант реализован в функции qqPlot() из пакета car:
<syntaxhighlight lang="r">> qqPlot(x, distribution = "norm")</syntaxhighlight>
Пакет e1071
Построение Q-Q plot можно осуществить с помощью функции probplot из пакета e1071:
<syntaxhighlight lang="r">> probplot(x, qdist = qnorm)</syntaxhighlight>
Многомерное нормальное распределение
Перед началом обзора функций, реализующий критерии проверки многомерной нормальности, сгенерируем массив данных. Сделать это можно при помощью следующих функций
-
mvrnormиз пакетаMASS -
rmvnormиз пакетаmvtnorm -
rmnormиз пакетаmnormt
Вот пример кода, генерирующего массив данных, имеющих многомерное нормальное распределение:
<syntaxhighlight lang="r">> means <- c(0, 0, 0, 0) # средние для переменных > sigmas <- diag(length(means)) # ковариационная матрица > mx <- rmvnorm(100, mean = means, sigma = sigmas)</syntaxhighlight>
Примечания
- ↑ Для оценки нормальности вызов выглядит следующим образом:
ks.test(x, y = "pnorm"). - ↑ По результатам сравнения производительности, данный вариант оказался чуть быстрее предыдущего.
- ↑ Более подробную информацию о доступных семействах распределений можно получить с помощью команды
help("gamlss.family").
