R:Статистическая проверка принадлежности нормальному распределения — различия между версиями

Материал Psylab.info - энциклопедии психодиагностики
Перейти к: навигация, поиск
м
м (Сравнительная таблица реализации критериев в пакетах)
 
(не показано 30 промежуточных версий этого же участника)
Строка 9: Строка 9:
 
== Одномерное нормальное распределение ==
 
== Одномерное нормальное распределение ==
  
В качестве <math>H_0</math> для всех нижеприведённых критериев является предположение, что «случайная величина <math>X</math> распределена нормально».
+
Нулевой гипотезой (<math>H_0</math>) для всех нижеприведённых критериев является предположение, что «случайная величина <math>X</math> распределена нормально».
  
Для демонстрации работы функций, реализующий различные критерий проверки принадлежности распределения нормальному закону сгенерируем вектор случайных чисел, имеющих стандартное нормальное распределение:
+
Для демонстрации работы функций, реализующих различные критерии проверки принадлежности распределения нормальному закону сгенерируем вектор случайных чисел, имеющих стандартное нормальное распределение:
  
 
{{r-code|code=
 
{{r-code|code=
> x <- rnorm(n = 100)
+
<nowiki>> x <- rnorm(n = 1000)
 +
</nowiki>
 
}}
 
}}
  
 
=== Статистические критерии ===
 
=== Статистические критерии ===
  
==== Пакет <code>stats</code> ====
+
В R реализовано множество критериев проверки соответствия распределения нормальному закону.
  
В данном пакете реализованы две функции, которые позволяют осуществить проверку принадлежности распределения нормальному закону.
+
==== Сравнительная таблица реализации критериев в пакетах ====
  
* <code>shapiro.test</code> - критерий Шапиро - Уилка
+
{| class="wide wikitable sortable" style="text-align: center"
* <code>ks.test</code> - критерий Колмогорова - Смирнова<ref>Для оценки нормальности вызов выглядит следующим образом:<code>ks.test(x, y = "pnorm")</code></ref>
+
! Критерии !! {{r-package|stats|core=true}} !! {{r-package|nortest}} !! {{r-package|moments}} !! {{r-package|fBasics}} !! {{r-package|tseries}} !! {{r-package|lawstat}}
 +
|-
 +
| style="text-align: left" | Критерий Шапиро - Уилка || <code>shapiro.test</code> || - || - || <code>shapiroTest</code> || - || -
 +
|-
 +
| style="text-align: left" | Критерий Колмогорова - Смирнова || <code>ks.test</code><ref>Для оценки нормальности вызов выглядит следующим образом: <code>ks.test(x, y = "pnorm")</code>.</ref> || - || - || <code>ksnormTest</code> || - || -
 +
|-
 +
| style="text-align: left" | Критерий Андерсона - Дарлинга || - || <code>ad.test</code> || - || <code>adTest</code> || - || -
 +
|-
 +
| style="text-align: left" | Критерий Крамера - фон Мизеса || - || <code>cvm.test</code> || - || <code>cvmTest</code> || - || -
 +
|-
 +
| style="text-align: left" | Критерий Лиллиефорса || - || <code>lillie.test</code> || - || <code>lillieTest</code> || - || -
 +
|-
 +
| style="text-align: left" | Критерий <math>\chi^2</math> Пирсона || - || <code>pearson.test</code> || - || <code>pchiTest</code> || - || -
 +
|-
 +
| style="text-align: left" | Критерий Шапиро - Франчия || - || <code>sf.test</code> || - || <code>sfTest</code> || - || -
 +
|-
 +
| style="text-align: left" | Критерий Д'Агостино || - || - || <code>agostino.test</code> || <code>dagoTest</code> || - || -
 +
|-
 +
| style="text-align: left" | Критерий Бонетта – Сайера || - || - || <code>bonett.test</code> || - || - || -
 +
|-
 +
| style="text-align: left" | Критерий Жарка - Бера || - || - || <code>jarque.test</code> || <code>jarqueberaTest</code> || <code>jarque.bera.test</code> || <code>rjb.test</code>
 +
|}
  
Данные функции возвращают результат в виде S3-класса - <code>htest</code>.
+
Пакет {{r-package|fBasics}} содержит также функцию <code>normalTest()</code>, которая является «обёрктой» для ряда функций из того же пакета. Необходимый критерий можно задать с помощью аргумента <code>method</code>. Доступны следующие критерии:
 
+
==== Пакет <code>nortest</code> ====
+
 
+
В данный пакет входят следующие функции:
+
 
+
* <code>ad.test</code> - критерий Андерсона - Дарлинга
+
* <code>cvm.test</code> - критерий Крамера - фон Мизеса
+
* <code>lillie.test</code> - критерий Лиллиефорса
+
* <code>pearson.test</code> - критерий <math>\chi^2</math> Пирсона
+
* <code>sf.test</code> - критерий Шапиро - Франчия
+
 
+
Данные функции возвращают результат в виде S3-класса - <code>htest</code>.
+
 
+
==== Пакет <code>moments</code> ====
+
 
+
В данный пакет входят следующие функции:
+
 
+
* <code>agostino.test</code> - критерий Д'Агостино
+
* <code>bonett.test</code> - критерий Бонетта – Сайера
+
* <code>jarque.test</code> - критерий Жарка-Бера
+
 
+
Данные функции также возвращают результат в виде S3-класса - <code>htest</code>.
+
 
+
==== Пакет <code>fBasics</code> ====
+
 
+
В данном пакете не предлагается никакой оригинальной реализации критериев - код в основном заимствован из пакетов <code>stats</code>, <code>nortest</code>, <code>moments</code>. Данный пакет предлагает альтернативный вывод результатов в виде объекта S4-класса <code>fHTEST</code>, в том время как все предыдущие функции использовали S3-класс <code>htest</code>.
+
 
+
Функция <code>normalTest()</code> является «обёрктой» для ряда функций из того же пакета - <code>fBasics</code>. Задать необходимый критерий можно задать с помощью аргумента <code>method</code>. Доступны следующие критерии:
+
  
 
* <code>sw</code> - критерий Шапиро - Уилка
 
* <code>sw</code> - критерий Шапиро - Уилка
Строка 65: Строка 59:
  
 
{{r-code|code=
 
{{r-code|code=
> normalTest(x, method = "sw")
+
<nowiki>> normalTest(x, method = "sw")
  
 
Title:
 
Title:
Строка 78: Строка 72:
 
Description:
 
Description:
 
  Fri Feb 14 19:59:59 2014 by user:
 
  Fri Feb 14 19:59:59 2014 by user:
 +
</nowiki>
 
}}
 
}}
  
Помимо функции <code>normalTest()</code> данный пакет включает в себя следующие функции:
+
Пакет <code>lawstat</code> содержит также функцию <code>sj.test()</code>, которая является реализацией рабастного критерия нормальности, созданного на основа критерия Шапиро - Уилка.
  
* <code>shapiroTest</code> - критерий Шапиро - Уилка
+
Пакет <code>TeachingDemos</code> содержит функцию <code>SnowsPenultimateNormalityTest()</code>, реализующую неописанный в литературе критерий. Данная функция возвращает только уровень статистической значимости, свидетельствующий об отклонения распределения от нормального закона.
* <code>ksnormTest</code> - критерий Колмогорова - Смирнова<ref>Данная функция вызывает <code>ks.test(x, "pnorm")</code> для трёх альтернативных гипотез - двусторонней и двух односторонних.</ref>
+
* <code>jarqueberaTest</code> - критерий Жарка-Бера
+
* <code>dagoTest</code> - критерий Д'Агостино
+
* <code>adTest</code> - критерий Андерсона - Дарлинга
+
* <code>cvmTest</code> - критерий Крамера - фон Мизеса
+
* <code>lillieTest</code> - критерий Лиллиефорса
+
* <code>pchiTest</code> - критерий Пирсона
+
* <code>sfTest</code> - критерий Шапиро - Франчия
+
 
+
Данные функции также возвращают результат в виде S4-класса - <code>fHTEST</code>.
+
 
+
==== Пакет <code>TeachingDemos</code> ====
+
 
+
Данные пакет содержит только одну функцию, имеющую отношение к критериям проверки принадлежности распределения нормальному закону - <code>SnowsPenultimateNormalityTest()</code>. Данная функция возвращают результат в виде S3-класса - <code>htest</code>.
+
 
+
==== Пакет <code>tseries</code> ====
+
 
+
Данный пакет содержит только одну функцию, имеющую отношение к критериям проверки принадлежности распределения нормальному закону - <code>jarque.bera.test</code>, которая является реализацией критерия Жарка-Бера. Данная функция возвращают результат в виде S3-класса - <code>htest</code>.
+
 
+
==== Пакет <code>lawstat</code> ====
+
 
+
В данный пакет входят следующие функции:
+
 
+
* <code>rjb.test</code> - критерий Жарка-Бера
+
* <code>sj.test</code> - SJ-критерий
+
  
 
==== Маленькие хитрости ====
 
==== Маленькие хитрости ====
Строка 116: Строка 86:
  
 
{{r-code|code=
 
{{r-code|code=
> DF <- data.frame(replicate(n = 10, rnorm(n = 100)))
+
<nowiki>> DF <- data.frame(replicate(n = 10, rnorm(n = 100)))
 +
</nowiki>
 
}}
 
}}
  
Строка 122: Строка 93:
  
 
{{r-code|code=
 
{{r-code|code=
> str(DF)
+
<nowiki>> str(DF)
 
'data.frame': 100 obs. of  10 variables:
 
'data.frame': 100 obs. of  10 variables:
 
  $ X1 : num  1.051 1.08 -0.477 -1.396 3.423 ...
 
  $ X1 : num  1.051 1.08 -0.477 -1.396 3.423 ...
Строка 134: Строка 105:
 
  $ X9 : num  -0.8508 0.4481 -0.2828 -0.5464 0.0605 ...
 
  $ X9 : num  -0.8508 0.4481 -0.2828 -0.5464 0.0605 ...
 
  $ X10: num  1.421 0.408 1.254 -0.956 -1.91 ...
 
  $ X10: num  1.421 0.408 1.254 -0.956 -1.91 ...
 +
</nowiki>
 
}}
 
}}
  
Для решения поставленной задачи можно воспользоваться функцией <code>sapply()</code>. Но прежде, нам необходимо немного отформатировать формат вывода результатов нашей функции: нам нужно извлечь значения критерия и его уровень значимости, т.к. результат функции <code>shapiro.test()</code> содержит также информацию, которая не подлежит включению в таблицу.
+
Для решения поставленной задачи можно воспользоваться функцией <code>sapply()</code>. Но прежде, нам необходимо немного отформатировать формат вывода результатов нашей функции: нам нужно извлечь значения критерия и его уровень значимости, т.к. результат функции <code>shapiro.test()</code> содержит также информацию, которая не подлежит включению в итоговую таблицу, например, информация об используемом методе (критерии) и уточнение характера альтернативной гипотезы. Вывод результатов тест Шапиро - Уилка выглядит следующим образом:
  
 
{{r-code|code=
 
{{r-code|code=
> shapiro.test(x)
+
<nowiki>> shapiro.test(x)
  
 
Shapiro-Wilk normality test
 
Shapiro-Wilk normality test
Строка 145: Строка 117:
 
data:  x
 
data:  x
 
W = 0.9903, p-value = 0.6882
 
W = 0.9903, p-value = 0.6882
 +
</nowiki>
 
}}
 
}}
  
Строка 150: Строка 123:
  
 
{{r-code|code=
 
{{r-code|code=
> str(shapiro.test(x))
+
<nowiki>> str(shapiro.test(x))
 
List of 4
 
List of 4
 
  $ statistic: Named num 0.99
 
  $ statistic: Named num 0.99
Строка 158: Строка 131:
 
  $ data.name: chr "x"
 
  $ data.name: chr "x"
 
  - attr(*, "class")= chr "htest"
 
  - attr(*, "class")= chr "htest"
 +
</nowiki>
 
}}
 
}}
  
Строка 163: Строка 137:
  
 
{{r-code|code=
 
{{r-code|code=
> normTest <- function (x) {
+
<nowiki>> normTest <- function (x) {
 
+    res <- shapiro.test(x)
 
+    res <- shapiro.test(x)
 
+    return(c(res$statistic, p.value = res$p.value))
 
+    return(c(res$statistic, p.value = res$p.value))
 
+ }
 
+ }
 +
</nowiki>
 
}}
 
}}
  
Строка 172: Строка 147:
  
 
{{r-code|code=
 
{{r-code|code=
> normTest(x)
+
<nowiki>> normTest(x)
 
       W p.value  
 
       W p.value  
 
  0.9903  0.6882
 
  0.9903  0.6882
 +
</nowiki>
 
}}
 
}}
  
Строка 180: Строка 156:
  
 
{{r-code|code=
 
{{r-code|code=
> t(sapply(DF, normTest))
+
<nowiki>> t(sapply(DF, normTest))
 
         W p.value
 
         W p.value
 
X1  0.9831  0.2301
 
X1  0.9831  0.2301
Строка 192: Строка 168:
 
X9  0.9915  0.7834
 
X9  0.9915  0.7834
 
X10 0.9808  0.1531
 
X10 0.9808  0.1531
 +
</nowiki>
 
}}
 
}}
  
Строка 197: Строка 174:
  
 
{{r-code|code=
 
{{r-code|code=
> do.call(rbind, lapply(DF, normTest))
+
<nowiki>> do.call(rbind, lapply(DF, normTest))
 
         W p.value
 
         W p.value
 
X1  0.9831  0.2301
 
X1  0.9831  0.2301
Строка 209: Строка 186:
 
X9  0.9915  0.7834
 
X9  0.9915  0.7834
 
X10 0.9808  0.1531
 
X10 0.9808  0.1531
 +
</nowiki>
 
}}
 
}}
  
Строка 216: Строка 194:
  
 
{{r-code|code=
 
{{r-code|code=
> DF$GRP <- factor(sample(LETTERS[1:3], size = 100, replace = TRUE))
+
<nowiki>> DF$GRP <- factor(sample(LETTERS[1:3], size = 100, replace = TRUE))</nowiki>
 
}}
 
}}
  
Строка 222: Строка 200:
  
 
{{r-code|code=
 
{{r-code|code=
> table(DF$GRP)
+
<nowiki>> table(DF$GRP)
 
  A  B  C  
 
  A  B  C  
38 25 37  
+
38 25 37 </nowiki>
 
}}
 
}}
  
Строка 230: Строка 208:
  
 
{{r-code|code=
 
{{r-code|code=
> do.call(rbind, tapply(DF$X1, DF$GRP, normTest))
+
<nowiki>> do.call(rbind, tapply(DF$X1, DF$GRP, normTest))
 
       W p.value
 
       W p.value
 
A 0.9522 0.13281
 
A 0.9522 0.13281
 
B 0.9607 0.28697
 
B 0.9607 0.28697
C 0.9410 0.07256
+
C 0.9410 0.07256</nowiki>
 
}}
 
}}
  
 
=== Графические методы ===
 
=== Графические методы ===
  
Многие исследователи также используют графические методы для определения степени отклонения распределения от нормального закона. В R реализована возможность построения Q-Q и P-P графиков, гистограмм и кривых распределения плотности вероятностей.
+
Многие исследователи также используют графические методы для определения степени отклонения распределения от нормального закона. В R реализована возможность построения Q-Q графиков, гистограмм и кривых распределения плотностей вероятности.
  
==== Пакет <code>stats</code> ====
+
==== Гистограмма ====
 +
 
 +
Гистограмма представляет собой графическое изображение зависимости частоты попадания элементов выборки от соответствующего интервала группировки. Построить гистограмму в R можно с помощью следующей команды:
 +
 
 +
{{r-code|code=
 +
<nowiki>> hist(x)</nowiki>}}
 +
 
 +
[[Файл:Graphics-hist.svg|400px|центр]]
 +
 
 +
На гистограмме изображены абсолютные частоты. Также можно построить гистограмму, отражающую плотности вероятностей:
 +
 
 +
{{r-code|code=
 +
<nowiki>> hist(x, freq = FALSE)</nowiki>}}
 +
 
 +
[[Файл:Graphics-hist-probs.svg|400px|центр]]
 +
 
 +
==== График плотностей вероятности ====
 +
 
 +
===== Пакет <code>stats</code> =====
 +
 
 +
{{r-code|code=
 +
<nowiki>> plot(density(x))</nowiki>}}
 +
 
 +
[[Файл:Stats-density.svg|400px|центр]]
 +
 
 +
===== Пакет <code>car</code> =====
 +
 
 +
{{r-code|code=
 +
<nowiki>> densityPlot(x)</nowiki>}}
 +
 
 +
[[Файл:Car-densityPlot.svg|400px|центр]]
 +
 
 +
==== Гистограммы с наложением графика плотностей вероятнотси ====
 +
 
 +
===== Пакет <code>stats</code> =====
 +
 
 +
{{r-code|code=
 +
<nowiki>> hist(x, freq = FALSE)
 +
> lines(density(x))</nowiki>}}
 +
 
 +
[[Файл:Stats-hist-density.svg|400px|центр]]
 +
 
 +
Теперь наложим на наш график кривую плотностей вероятности для нормального распределения:
 +
 
 +
{{r-code|code=
 +
<nowiki>> xfit <- seq(min(x), max(x), length = 100) # Координаты по оси X
 +
> yfit <- dnorm(xfit, mean = mean(x), sd = sd(x)) # Вычисление координат по оси Y
 +
> hist(x, freq = FALSE)
 +
> lines(density(x), col = "red") # Накладываем кривую плотностей вероятности
 +
> lines(xfit, yfit, col = "blue") # Накладываем «нормальную» кривую
 +
</nowiki>}}
 +
 
 +
[[Файл:Stats-density-compare.svg|400px|центр]]
 +
 
 +
===== Пакет <code>gamlss</code> =====
 +
 
 +
Более простой способ сравнение графиков плотностей вероятности представлен в функции <code>histDist</code> из пакета <code>gamlss</code>:
 +
 
 +
{{r-code|code=
 +
<nowiki>> 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
 +
</nowiki>
 +
}}
 +
 
 +
[[Файл:Gamlss-histdist.svg|400px|центр]]
 +
 
 +
С помощью аргумента <code>family</code> можно задать семейство распределений для подгонки и сравнения<ref>Более подробную информацию о доступных семействах распределений можно получить с помощью команды <code>help("gamlss.family")</code>.</ref>.
 +
 
 +
==== Q-Q график ====
 +
 
 +
Q-Q график (Q - квантиль) — это график, на котором квантили из двух распределений расположены относительно друг друга. Чем ближе точки на графике к диагональной прямой, тем ближе распределение исследуемой переменной к нормальному закону.
 +
 
 +
Построение квантильных графиков в R реализовано в нескольких пакетах.
 +
 
 +
===== Пакет <code>stats</code> =====
  
 
Построение Q–Q plot с помощью пакета <code>stats</code> выглядит следующим образом:
 
Построение Q–Q plot с помощью пакета <code>stats</code> выглядит следующим образом:
  
 
{{r-code|code=
 
{{r-code|code=
> qqnorm(x)
+
<nowiki>> qqnorm(x)
> qqline(x)
+
> qqline(x)</nowiki>
 
}}
 
}}
  
 
[[Файл:Stats-qqnorm.svg|400px|центр]]
 
[[Файл:Stats-qqnorm.svg|400px|центр]]
  
==== Пакет <code>QTLRel</code> ====
+
===== Пакет <code>QTLRel</code> =====
  
 
Построение Q–Q plot с помощью пакета <code>QTLRel</code> выглядит следующим образом:
 
Построение Q–Q plot с помощью пакета <code>QTLRel</code> выглядит следующим образом:
  
 
{{r-code|code=
 
{{r-code|code=
> qqPlot(x, x = "norm")
+
<nowiki>> qqPlot(x, x = "norm")</nowiki>
 
}}
 
}}
  
 
[[Файл:Qtlrel-qqplot.svg|400px|центр]]
 
[[Файл:Qtlrel-qqplot.svg|400px|центр]]
  
==== Пакет <code>car</code> ====
+
===== Пакет <code>car</code> =====
  
 
Альтернативный вариант реализован в функции <code>qqPlot()</code> из пакета <code>car</code>:
 
Альтернативный вариант реализован в функции <code>qqPlot()</code> из пакета <code>car</code>:
  
 
{{r-code|code=
 
{{r-code|code=
> qqPlot(x, distribution = "norm")
+
<nowiki>> qqPlot(x, distribution = "norm")</nowiki>
 
}}
 
}}
  
[[Файл:Car-qqplot.svg|400px|центр]]
+
[[Файл:Car-qqPlot.svg|400px|центр]]
  
==== Пакет <code>e1071</code> ====
+
===== Пакет <code>e1071</code> =====
  
Построение P-P plot можно осуществить с помощью функции <code>probplot</code> из пакета <code>e1071</code>:
+
Построение Q-Q plot можно осуществить с помощью функции <code>probplot</code> из пакета <code>e1071</code>:
  
 
{{r-code|code=
 
{{r-code|code=
> probplot(x, qdist = qnorm)
+
<nowiki>> probplot(x, qdist = qnorm)</nowiki>
 
}}
 
}}
  
 
[[Файл:E1071-probplot.svg|400px|центр]]
 
[[Файл:E1071-probplot.svg|400px|центр]]
  
==== Пакет <code>gamlss</code> ====
+
== Многомерное нормальное распределение ==
  
Ещё один интересный способ графического анализа представлен функцией <code>histDist</code> из пакета <code>gamlss</code>:
+
Перед началом обзора функций, реализующий критерии проверки многомерной нормальности, сгенерируем массив данных. Сделать это можно при помощью следующих функций
 +
* <code>mvrnorm</code> из пакета <code>MASS</code>
 +
* <code>rmvnorm</code> из пакета <code>mvtnorm</code>
 +
* <code>rmnorm</code> из пакета <code>mnormt</code>
 +
 
 +
Вот пример кода, генерирующего массив данных, имеющих многомерное нормальное распределение:
  
 
{{r-code|code=
 
{{r-code|code=
> histDist(x, family = "NO", density = TRUE)
+
> means <- c(0, 0, 0, 0) # средние для переменных
 +
> sigmas <- diag(length(means)) # ковариационная матрица
 +
> mx <- rmvnorm(100, mean = means, sigma = sigmas)
 +
}}
  
Family:  c("NO", "Normal")  
+
Пакет {{r-package|mvnormtest}} реализует модификацию критерия Шапиро - Уилка для многомерных данных - функция <code>mshapiro.test()</code><ref>В качестве аргумента необходимо передать транспонированную матрицу: <code>mshapiro.test(t(mx))</code>.</ref>.
Fitting method: "nlminb"
+
  
Call: gamlssML(y = y, family = "NO", formula = x)  
+
Пакет {{r-package|ICS}} предлагает реализацию критериев эксцесса и асимметрии для многомерных данных: <code>mvnorm.kur.test()</code>, <code>mvnorm.skew.test()</code>.
  
Mu Coefficients:
+
Пакет {{r-package|energy}} реализует E-статистики для сравнения распределений. Критерия для проверки гипотезы о соответствия распределения многомерной переменной многомерному нормальному распределению предлагается функция <code>mvnorm.etest()</code><ref>Для вычисления уровня значимости критерия используется метод бутстрепа (bootstrap). Число итераций для бутстрепа можно задать с помощью аргумента <code>R</code>.</ref>.
[1]  -0.2273
+
Sigma Coefficients:
+
[1]  0.09813
+
  
Degrees of Freedom for the fit: 2 Residual Deg. of Freedom  98
+
== Ссылки ==
Global Deviance:    303.414
+
            AIC:    307.414
+
            SBC:    312.624
+
}}
+
  
[[Файл:gamlss-histdist.svg|400px|центр]]
+
* Juergen Gross and bug fixes by Uwe Ligges (2012). nortest: Tests for Normality. R package version 1.0-2.
 
+
*: http://CRAN.R-project.org/package=nortest
С помощью аргумента <code>family</code> можно задать семейство распределений для подгонки и сравнения<ref>Более подробную информацию о доступных семействах распределений можно получить с помощью команды <code>help("gamlss.family")</code>.</ref>.
+
* Lukasz Komsta and Frederick Novomestky (2012). moments: Moments, cumulants, skewness, kurtosis and related tests. R package version 0.13.
 
+
*: http://CRAN.R-project.org/package=moments
== Многомерное нормальное распределение ==
+
* Diethelm Wuertz, Rmetrics core team members, uses code builtin from the following R contributed packages: gmm from Pierre Chauss, gld from Robert King, gss from Chong Gu, nortest from Juergen Gross, HyperbolicDist from David Scott, sandwich from Thomas Lumley, Achim Zeileis, fortran/C code from Kersti Aas and akima from Albrecht Gebhardt (2013). fBasics: Rmetrics - Markets and Basic Statistics. R package version 3010.86.
 +
*: http://CRAN.R-project.org/package=fBasics
 +
* Adrian Trapletti and Kurt Hornik (2013). tseries: Time Series Analysis and Computational Finance. R package version 0.10-32.
 +
*: http://CRAN.R-project.org/package=tseries
 +
* Joseph L. Gastwirth; Yulia R. Gel <ygl@math.uwaterloo.ca>; W. L. Wallace Hui <wlwhui@uwaterloo.ca>; Vyacheslav Lyubchich <vlyubchich@uwaterloo.ca>; Weiwen Miao <miao@macalester.edu>; Kimihiro Noguchi <kinoguchi@ucdavis.edu> (2013). lawstat: An R package for  biostatistics, public policy, and law. R package version 2.4.1.
 +
*: http://CRAN.R-project.org/package=lawstat
 +
* John Fox and Sanford Weisberg (2013). car: Companion to Applied Regression. R package version 2.0-19/r346.
 +
*: http://R-Forge.R-project.org/projects/car/
 +
* Mikis Stasinopoulos, Bob Rigby with contributions from Calliope Akantziliotou and Vlasios Voudouris (2014). gamlss: Generalised Additive Models for Location Scale and Shape. R package version 4.2-7.
 +
*: http://CRAN.R-project.org/package=gamlss
 +
* Riyan Cheng (2013). QTLRel: Tools for mapping of quantitative traits of genetically related individuals and calculating identity coefficients from a pedigree. R package version 0.2-14.
 +
*: http://CRAN.R-project.org/package=QTLRel
 +
* David Meyer, Evgenia Dimitriadou, Kurt Hornik, Andreas Weingessel and Friedrich Leisch (2014). e1071: Misc Functions of the Department of Statistics (e1071), TU Wien. R package version 1.6-2.
 +
*: http://CRAN.R-project.org/package=e1071
  
 
== Примечания ==
 
== Примечания ==

Текущая версия на 08:55, 7 декабря 2014

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

На практике мы встречаемся с двумя вариантами задач по проверке принадлежности распределения нормальному закону: для одномерного и многомерного распределения.


Одномерное нормальное распределение

Нулевой гипотезой ([math]H_0[/math]) для всех нижеприведённых критериев является предположение, что «случайная величина [math]X[/math] распределена нормально».

Для демонстрации работы функций, реализующих различные критерии проверки принадлежности распределения нормальному закону сгенерируем вектор случайных чисел, имеющих стандартное нормальное распределение:

КодR

<syntaxhighlight lang="r">> x <- rnorm(n = 1000) </syntaxhighlight>

Статистические критерии

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

Сравнительная таблица реализации критериев в пакетах

Критерии 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

Пакет fBasics содержит также функцию normalTest(), которая является «обёрктой» для ряда функций из того же пакета. Необходимый критерий можно задать с помощью аргумента method. Доступны следующие критерии:

  • sw - критерий Шапиро - Уилка
  • jb - критерий Жарка-Бера
  • ks - критерий Колмогорова - Смирнова
  • da - критерий Д'Агостино
  • ad - критерий Андерсона - Дарлинга.

Пример вызова данной функции:

КодR

<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>

Пакет lawstat содержит также функцию sj.test(), которая является реализацией рабастного критерия нормальности, созданного на основа критерия Шапиро - Уилка.

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

Маленькие хитрости

Применение функций к нескольким переменным

С помощью apply-функций можно последовательно применить функцию к вектору, списку или массиву. Прежде чем всего нам необходимо сформировать таблицу данных. С помощью функции replicate() сгенерируем 10 переменных, имеющих стандартное нормальное распределение, которые объединяются в класс data.frame.

КодR

<syntaxhighlight lang="r">> DF <- data.frame(replicate(n = 10, rnorm(n = 100))) </syntaxhighlight>

Структура сгенерированной таблицы выглядит следующим образом:

КодR

<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() содержит также информацию, которая не подлежит включению в итоговую таблицу, например, информация об используемом методе (критерии) и уточнение характера альтернативной гипотезы. Вывод результатов тест Шапиро - Уилка выглядит следующим образом:

КодR

<syntaxhighlight lang="r">> shapiro.test(x) Shapiro-Wilk normality test data: x W = 0.9903, p-value = 0.6882 </syntaxhighlight>

Структура результата применения функции shapiro.test() представлена ниже:

КодR

<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>

Как видим, помимо значений критерия и уровня значимости здесь содержится информация о применяемом методе. Мы можем отфильтровать вывод следующим образом:

КодR

<syntaxhighlight lang="r">> normTest <- function (x) { + res <- shapiro.test(x) + return(c(res$statistic, p.value = res$p.value)) + } </syntaxhighlight>

Результат теперь будет выглядеть следующим образом:

КодR

<syntaxhighlight lang="r">> normTest(x) W p.value 0.9903 0.6882 </syntaxhighlight>

Теперь можно использовать данную функцию при обработке столбцов нашей таблицы.

КодR

<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]:

КодR

<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>

Применение функций к нескольким по группам

Добавим к нашей таблице группы испытуемых:

КодR

<syntaxhighlight lang="r">> DF$GRP <- factor(sample(LETTERS[1:3], size = 100, replace = TRUE))</syntaxhighlight>

Состав групп получился следующим:

КодR

<syntaxhighlight lang="r">> table(DF$GRP) A B C 38 25 37 </syntaxhighlight>

Рассчитаем значения критерия Шапиро - Уилка для первого столбца для каждоый группы испытуемых:

КодR

<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 графиков, гистограмм и кривых распределения плотностей вероятности.

Гистограмма

Гистограмма представляет собой графическое изображение зависимости частоты попадания элементов выборки от соответствующего интервала группировки. Построить гистограмму в R можно с помощью следующей команды:

КодR

<syntaxhighlight lang="r">> hist(x)</syntaxhighlight>

Graphics-hist.svg

На гистограмме изображены абсолютные частоты. Также можно построить гистограмму, отражающую плотности вероятностей:

КодR

<syntaxhighlight lang="r">> hist(x, freq = FALSE)</syntaxhighlight>

Graphics-hist-probs.svg

График плотностей вероятности

Пакет stats
КодR

<syntaxhighlight lang="r">> plot(density(x))</syntaxhighlight>

Stats-density.svg
Пакет car
КодR

<syntaxhighlight lang="r">> densityPlot(x)</syntaxhighlight>

Car-densityPlot.svg

Гистограммы с наложением графика плотностей вероятнотси

Пакет stats
КодR

<syntaxhighlight lang="r">> hist(x, freq = FALSE) > lines(density(x))</syntaxhighlight>

Stats-hist-density.svg

Теперь наложим на наш график кривую плотностей вероятности для нормального распределения:

КодR

<syntaxhighlight lang="r">> xfit <- seq(min(x), max(x), length = 100) # Координаты по оси X > yfit <- dnorm(xfit, mean = mean(x), sd = sd(x)) # Вычисление координат по оси Y > hist(x, freq = FALSE) > lines(density(x), col = "red") # Накладываем кривую плотностей вероятности > lines(xfit, yfit, col = "blue") # Накладываем «нормальную» кривую </syntaxhighlight>

Stats-density-compare.svg
Пакет gamlss

Более простой способ сравнение графиков плотностей вероятности представлен в функции histDist из пакета gamlss:

КодR

<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>

Gamlss-histdist.svg

С помощью аргумента family можно задать семейство распределений для подгонки и сравнения[3].

Q-Q график

Q-Q график (Q - квантиль) — это график, на котором квантили из двух распределений расположены относительно друг друга. Чем ближе точки на графике к диагональной прямой, тем ближе распределение исследуемой переменной к нормальному закону.

Построение квантильных графиков в R реализовано в нескольких пакетах.

Пакет stats

Построение Q–Q plot с помощью пакета stats выглядит следующим образом:

КодR

<syntaxhighlight lang="r">> qqnorm(x) > qqline(x)</syntaxhighlight>

Stats-qqnorm.svg
Пакет QTLRel

Построение Q–Q plot с помощью пакета QTLRel выглядит следующим образом:

КодR

<syntaxhighlight lang="r">> qqPlot(x, x = "norm")</syntaxhighlight>

Qtlrel-qqplot.svg
Пакет car

Альтернативный вариант реализован в функции qqPlot() из пакета car:

КодR

<syntaxhighlight lang="r">> qqPlot(x, distribution = "norm")</syntaxhighlight>

Car-qqPlot.svg
Пакет e1071

Построение Q-Q plot можно осуществить с помощью функции probplot из пакета e1071:

КодR

<syntaxhighlight lang="r">> probplot(x, qdist = qnorm)</syntaxhighlight>

E1071-probplot.svg

Многомерное нормальное распределение

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

  • mvrnorm из пакета MASS
  • rmvnorm из пакета mvtnorm
  • rmnorm из пакета mnormt

Вот пример кода, генерирующего массив данных, имеющих многомерное нормальное распределение:

КодR

<syntaxhighlight lang="r">> means <- c(0, 0, 0, 0) # средние для переменных > sigmas <- diag(length(means)) # ковариационная матрица > mx <- rmvnorm(100, mean = means, sigma = sigmas)</syntaxhighlight>

Пакет mvnormtest реализует модификацию критерия Шапиро - Уилка для многомерных данных - функция mshapiro.test()[4].

Пакет ICS предлагает реализацию критериев эксцесса и асимметрии для многомерных данных: mvnorm.kur.test(), mvnorm.skew.test().

Пакет energy реализует E-статистики для сравнения распределений. Критерия для проверки гипотезы о соответствия распределения многомерной переменной многомерному нормальному распределению предлагается функция mvnorm.etest()[5].

Ссылки

  • Juergen Gross and bug fixes by Uwe Ligges (2012). nortest: Tests for Normality. R package version 1.0-2.
    http://CRAN.R-project.org/package=nortest
  • Lukasz Komsta and Frederick Novomestky (2012). moments: Moments, cumulants, skewness, kurtosis and related tests. R package version 0.13.
    http://CRAN.R-project.org/package=moments
  • Diethelm Wuertz, Rmetrics core team members, uses code builtin from the following R contributed packages: gmm from Pierre Chauss, gld from Robert King, gss from Chong Gu, nortest from Juergen Gross, HyperbolicDist from David Scott, sandwich from Thomas Lumley, Achim Zeileis, fortran/C code from Kersti Aas and akima from Albrecht Gebhardt (2013). fBasics: Rmetrics - Markets and Basic Statistics. R package version 3010.86.
    http://CRAN.R-project.org/package=fBasics
  • Adrian Trapletti and Kurt Hornik (2013). tseries: Time Series Analysis and Computational Finance. R package version 0.10-32.
    http://CRAN.R-project.org/package=tseries
  • Joseph L. Gastwirth; Yulia R. Gel <ygl@math.uwaterloo.ca>; W. L. Wallace Hui <wlwhui@uwaterloo.ca>; Vyacheslav Lyubchich <vlyubchich@uwaterloo.ca>; Weiwen Miao <miao@macalester.edu>; Kimihiro Noguchi <kinoguchi@ucdavis.edu> (2013). lawstat: An R package for biostatistics, public policy, and law. R package version 2.4.1.
    http://CRAN.R-project.org/package=lawstat
  • John Fox and Sanford Weisberg (2013). car: Companion to Applied Regression. R package version 2.0-19/r346.
    http://R-Forge.R-project.org/projects/car/
  • Mikis Stasinopoulos, Bob Rigby with contributions from Calliope Akantziliotou and Vlasios Voudouris (2014). gamlss: Generalised Additive Models for Location Scale and Shape. R package version 4.2-7.
    http://CRAN.R-project.org/package=gamlss
  • Riyan Cheng (2013). QTLRel: Tools for mapping of quantitative traits of genetically related individuals and calculating identity coefficients from a pedigree. R package version 0.2-14.
    http://CRAN.R-project.org/package=QTLRel
  • David Meyer, Evgenia Dimitriadou, Kurt Hornik, Andreas Weingessel and Friedrich Leisch (2014). e1071: Misc Functions of the Department of Statistics (e1071), TU Wien. R package version 1.6-2.
    http://CRAN.R-project.org/package=e1071

Примечания

  1. Для оценки нормальности вызов выглядит следующим образом: ks.test(x, y = "pnorm").
  2. По результатам сравнения производительности, данный вариант оказался чуть быстрее предыдущего.
  3. Более подробную информацию о доступных семействах распределений можно получить с помощью команды help("gamlss.family").
  4. В качестве аргумента необходимо передать транспонированную матрицу: mshapiro.test(t(mx)).
  5. Для вычисления уровня значимости критерия используется метод бутстрепа (bootstrap). Число итераций для бутстрепа можно задать с помощью аргумента R.