Краткий рассказ про пакет MVN
Минутка теории
Допустим, у нас есть некоторое совместное распределение n переменных – и нам необходимо проверить, является ли оно нормальным. Решить эту задачу просто нам мешает один маленький факт – из нормальности многомерного распределения следует нормальность распределения каждой переменной в отдельности, но в обратную сторону это работает только при случае независимости компонентов распределения, что на практике не выполняется почти никогда. Поэтому приходится что-то изобретать.
Схема проверки статистической гипотезы о нормальности многомерного распределения идентична соответствующей для одномерного случая, только в ней используются другие тесты.
Тест Мардиа (оригинальная работа: K. V. Mardia. Measures of multivariate skewness and kurtosis with applications. Biometrika, 57(3): 519–530, 1970) основан на вычислении эксцесса и асимметрии многомерного распределения по формулам
![n - количество наблюдений, р – переменных n - количество наблюдений, р – переменных](https://habrastorage.org/getpro/habr/upload_files/a4b/459/05f/a4b45905ff4387d2998fa62fd3ffac45.png)
При этом m – это расстояние Маланхобиса между i-м и j-м наблюдениями
![S - ковариационная матрицы S - ковариационная матрицы](https://habrastorage.org/getpro/habr/upload_files/a0f/3b7/937/a0f3b7937b544a1af4a89bda34815c9a.png)
В такой трактовке рассчитанная величина асимметрии, умноженная на n/6, распределена по закону Хи-квадрат с p(p+1)(p+2)/6 степенями свободы, а величина эксцесса распределена по закону нормального распределения со средним p(p+2) и отклонением 8p(p+2)/n
Тест Хенце-Циклера (базовая работа: N. Henze and B. Zirkler. A class of invariant consistent tests for multivariate normality. Communications in Statistics - Theory and Methods, 19(10):3595–3617, 1990.) основан на следующей формуле расчета статистического критерия:
![D – расстояние Маланхобиса, β – параметр D – расстояние Маланхобиса, β – параметр](https://habrastorage.org/getpro/habr/upload_files/ad2/a81/d40/ad2a81d407f50e4d8df93493a180b659.png)
![](https://habrastorage.org/getpro/habr/upload_files/ba3/7b0/d9d/ba37b0d9d70806f68d2a5be91b3e604e.png)
Значения критерия распределены по логнормальному закону с параметрами
![](https://habrastorage.org/getpro/habr/upload_files/0e4/d87/1f6/0e4d871f6d756859df726a45f7e7d345.png)
Тест Ройстона основан на идее теста Шапиро-Уилкса. Значение статистического критерия рассчитывается по формуле
![](https://habrastorage.org/getpro/habr/upload_files/515/de7/880/515de78808abaef801d2bd5fdc16ddce.png)
Его величина распределяется по закону Хи-квадрат с количеством степеней свободы, равным е. Цепочка расчетов следующая:
![Wj – значение статистики Шапиро-Уилка для j-ой переменной, r – коэффициент корреляции Wj – значение статистики Шапиро-Уилка для j-ой переменной, r – коэффициент корреляции](https://habrastorage.org/getpro/habr/upload_files/845/f04/663/845f04663cfa00b2e9d97c92df6eca88.png)
Тест Дорника-Хансена (оригинальная работа: Doornik, J. A., and H. Hansen. 2008. An omnibus test for univariate and multivariate normality. Oxford Bulletin of Economics and Statistics 70: 927–939.) основан на преобразовании многомерных наблюдений и вычислении эксцесса и асимметрии для одномерной переменной.
Преобразование проводиться по формуле
![Первая матрица этого произведения – центрированная матрица исходных данных Х
Вторая матрица – диагональная матрица, в которой элементы равны S-1/2 для отдельной переменной
Третья матрица – матрица собственных векторов корреляционной матрицы С
Четвертая матрица – диагональная матрица собственных значений матрицы С Первая матрица этого произведения – центрированная матрица исходных данных Х
Вторая матрица – диагональная матрица, в которой элементы равны S-1/2 для отдельной переменной
Третья матрица – матрица собственных векторов корреляционной матрицы С
Четвертая матрица – диагональная матрица собственных значений матрицы С](https://habrastorage.org/getpro/habr/upload_files/c33/e3f/711/c33e3f7111c697a75299a11c4b61a5de.png)
Далее рассчитывается эксцесс и асимметрия для каждой переменной в новой матрице.
Значения асимметрии (b1) и эксцесса (b2) распределены не по нормальному закону. Для их трансформации применяются следующие преобразования:
![](https://habrastorage.org/getpro/habr/upload_files/071/b0b/26b/071b0b26bc0d323495b928b334ddb52c.png)
Полученные значения z1 и z2 объединяются в вектора Z1 и Z2, а рассчитанная величина статистики распределена по закону Хи-квадрат с числом степеней свободы, равном 2k
![Формула расчета тестовой статистики Формула расчета тестовой статистики](https://habrastorage.org/getpro/habr/upload_files/a90/8dd/7ea/a908dd7ea22280d9a1de4aa28e5a10a3.png)
Тест Е-статистики (тест Шекели-Риццо, базовая работа: G.J. Szekely, M.L. Rizzo. A new test for multivariate normality / Journal of Multivariate Analysis 93 (2005) 58–80) подразумевает вычисление тестовой статистики с помощью разложения в ряд Тейлора:
![n – количество наблюдений, d – количество переменных n – количество наблюдений, d – количество переменных](https://habrastorage.org/getpro/habr/upload_files/bc8/52f/c57/bc852fc578d5395fb4d1baee330dbfd8.png)
y – центрированная матрица наблюдений, получаемая путем постолбцовых пребразований как
![](https://habrastorage.org/getpro/habr/upload_files/0b3/c35/e68/0b3c35e685c94416fac93b802700089a.png)
Методология
Для примера выберем базу данных “Crime” из пакета plm, и возьмем оттуда три переменных:
prbpris – вероятность тюремного заключения
avgsen – средний срок заключения, дней
pctymle – доля в населении мужчин в возрасте 15-24 лет
Из этих трех переменных соберем две базы данных – с двумя и тремя переменными:
library(MVN)
library(tidyverse)
library(plm)
data("Crime")
glimpse(Crime)
ggplot(Crime, aes(x=Crime$avgsen)) + geom_density()
# Crime$prbpris - точно, avgsen - 70/30, pctymle - 50/50
Data_1 <- Crime[,c(6,7)]
Data_2 <- Crime[,c(6,7,24)]
Расчеты и описание
Базовая функция расчетов – функция mvn со следующими параметрами:
data - База данных (в виде матрицы или датафрейма)
subset - Факторная группировочная переменная
mvnTest - Определяет статистический тест, которым проводится проверка
desc - Логическая переменная. Если она равна истине, выводятся описательные статистики
univariateTest - Определяет статистический тест, которым проводится проверка на нормальность отдельных переменных
univariatePlot - Определяет вид выводимого одномерного графика нормальности
multivariatePlot - Определяет вид графика ошибок
multivariateOutlierMethod - Выбирает метод определения выбросов
Проверим наши данные на нормальность с помощью классического теста Мардиа
![Результат один – NO в графах «Result» и «Normality» говорят нам о том, что нельзя принять гипотезу как о нормальном многомерном распределении как всей совокупности переменных, так и нормальности распределения каждой переменной в отдельности. Результат один – NO в графах «Result» и «Normality» говорят нам о том, что нельзя принять гипотезу как о нормальном многомерном распределении как всей совокупности переменных, так и нормальности распределения каждой переменной в отдельности.](https://habrastorage.org/getpro/habr/upload_files/fe8/67b/c26/fe867bc269e2481ee828376b365b8b07.png)
Можно глазами посмотреть на Q-Q график
mvn(data = Data_1, mvnTest = "mardia", multivariatePlot = "qq")
![](https://habrastorage.org/getpro/habr/upload_files/106/25d/d2c/10625dd2ceb7b4db17393ee1d84c730a.png)
А также на двумерный график распределения
mvn(data = Data_1, mvnTest = "mardia", multivariatePlot = "qq")
![](https://habrastorage.org/getpro/habr/upload_files/de8/44e/e24/de844ee244d94306c12c97269a1ac30a.png)
И на двумерный контурный график
mvn(data = Data_1, mvnTest = "energy", multivariatePlot = "contour")
![Цифры - соответствующие квантили Цифры - соответствующие квантили](https://habrastorage.org/getpro/habr/upload_files/16e/9ae/3ed/16e9ae3edcdd215b8b7bb759b4de44d3.png)
Можно также вывести, например, Q-Q график по каждой переменной в отдельности
mvn(data = Data_1, mvnTest = "mardia", univariatePlot = "qqplot")
Особенно интересна возможность, предоставляемая переменной subset. Если есть группировочная переменная, есть возможность проверить многомерные / одномерные нормальности в зависимости от разных ее значений:
![В нашем примере гипотеза о многомерной нормальности не подтверждается для любого из рассматриваемых регионов, но переменная «вероятность тюремного заключения» распределена по нормальному закону для западного и центрального района. В нашем примере гипотеза о многомерной нормальности не подтверждается для любого из рассматриваемых регионов, но переменная «вероятность тюремного заключения» распределена по нормальному закону для западного и центрального района.](https://habrastorage.org/getpro/habr/upload_files/dd3/9b7/ad1/dd39b7ad1e8a8b82fc6c504bb8018a62.png)
Это основы функционала пакета MVN. Все материалы доступны на https://github.com/acheremuhin/Multivariate_normal
victor_1212
Привет Артем, что у вас за терминал или принтер стоит (ностальгию почувствовал увидев последнюю картинку) ?
И для каких (функциональных) применений пакет предназначен ?
acheremuhin Автор
Добрый день, Виктор. Что касается ностальгии по последней картинке - это результат симбиоза PrintScreen и Point. Насчет же функционального применения - любая область, в которой вам нужно проверить наличие многомерного нормального распределения. Как правило, в предобработке данных, поскольку некоторые методы анализа в качестве условия требуют именно это условие.
victor_1212
> любая область, в которой ...
Спасибо Артем, насколько полезным этот пакет покажет себя именно на конкретных проблемах вот вопрос, достаточно знаю теорию вероятностей чтобы понимать вас с полуслова, но правильная интерпретация результатов это своего рода искусство, и это для меня самое интересное :)