Проверка теста на равенство двух дисперсий требует предварительной проверки двух выборок на нормальное распределение. Тест очень чувствителен к отклонениям от нормальности, поэтому в случае отклонений следует воспользоваться другими тестами, указанными в конце параграфа.
Предположим, что мы имеем две выборки x и y, извлеченные из двух нормально распределенных независимых генеральных совокупностей.
Пусть
alpha <- 0.05
x <- c(5.90, 5.90, 5.90, 8.10, 3.90, 4.10, 4.10, 4.10, 6.10, 11.60, 11.60, 5.90, 9.20, 5.90, 5.90, 5.90, 4.10, 5.90, 5.90, 2.30, 9.20, 3.90, 3.90, 6.20, 0.30, 2.90, 3.80, 4.30, 4.80, 4.10, 4.10, 7.90, 1.20, 7.90, 4.10, 4.90, 2.80, 7.80, 3.90, 7.90, 3.30, 6.20, 2.30, 3.90, 0.30, 1.20, 3.90, 7.90, 4.10, 2.70)
y <- c(2.30, 5.90, 7.90, 4.10, 5.00, 5.90, 4.30, 10.60, 10.60, 2.30, 3.90, 3.20, 1.20, 4.10, 4.80, 3.40, 5.30, 4.70, 9.20, 4.10, 9.20, 8.10, 3.90, 1.20, 5.20, 4.80, 3.90, 4.80, 4.90, 1.20, 0.87, 0.84, 7.90, 3.90, 5.90, 6.40, 5.10, 10.60, 6.30, 9.20, 9.20, 11.60, 3.90, 10.60, 10.60, 4.10, 6.40, 3.90, 4.10, 1.80)
Для проведения данного теста при оговоренных допущениях используется F-статистика Фишера.
Сам тест в R можно провести в различных пакетах, например, при помощи стандартного встроенного var.test(). По умолчанию данный тест проводится для двусторонней критической области. Поскольку для теста Фишера используется правосторонняя критическая область, это необходимо указать в качестве аргумента функции как alternative = “greater”:
var.test(x,y, alternative = "greater")
##
## F test to compare two variances
##
## data: x and y
## F = 0.73562, num df = 49, denom df = 49, p-value = 0.857
## alternative hypothesis: true ratio of variances is greater than 1
## 95 percent confidence interval:
## 0.4576774 Inf
## sample estimates:
## ratio of variances
## 0.73562
По умолчанию тест проверяет равенство дисперсий (ratio=1), но при желании можно проверить превышение в любое интересующее число раз одной дисперсии над другой.
Результаты теста включают: F-наблюдаемое, p-value, то есть все, что необходимо для проверки нулевой гипотезы. Так как p-value=0.857 больше, чем 0,05, то H0 не отвергается. Значит, можно считать на уровне значимости 0,05, что гипотеза о равенстве дисперсий генеральных совокупностей не отвергается. Тест, кроме того, выдает выборочное отношение дисперсий двух выборок (здесь 0.73562) и его интервальную оценку (она включает 1 и поэтому гипотеза не отвергается).
В нашем случае при использовании правосторонней критической области данного критерия нужна вероятность правого хвоста распределения, т.е. P[X ≥x], поэтому задаем это в функции (по умолчанию она выдает квантиль, левый хвост, т.е. вероятность P[X <x]).
qf(alpha,49,49,lower.tail = FALSE)
## [1] 1.607289
Как видим, F-наблюдаемое=0.73562 < F-критическое=1.607289, следовательно, на уровне значимости 0,05 гипотеза о равенстве дисперсий генеральных совокупностей не отвергается.
Проверка гипотезы о значении генеральной средней 𝑯𝟎: 𝝁 = 𝝁𝟎 при неизвестной генеральной дисперсии σ2 требует предварительной проверки выборки на нормальное распределение. Тест очень чувствителен к отклонениям от нормальности, поэтому в случае отклонений следует воспользоваться другими тестами, указанными далее.
Провести этот тест в R можно с помощью стандартной встроенной функции t.test() (аналогичная функция TTestA есть в пакете DescTools), которая поддерживает проверку гипотезы о равенстве среднего константе. Данная функция имеет аргумент mu=, который отвечает за параметр проверяемого среднего значения гипотезы 𝐻0: 𝜇 = 𝜇0.
t.test(x, y = NULL,
alternative = c("two.sided", "less", "greater"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
##
## One Sample t-test
##
## data: x
## t = 14.39, df = 49, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 4.370564 5.789436
## sample estimates:
## mean of x
## 5.08
Пример: Найдем среднее введенной нами ранее выборки х. Предположим, что выборка x извлечена из нормально распределенной генеральной совокупности (это надо обязательно проверять). Проверим состав нашей переменной с помощью удобной функции проверки целостности данных glimpse из пакета dplyr.
library(dplyr)
##
## Присоединяю пакет: 'dplyr'
## Следующие объекты скрыты от 'package:stats':
##
## filter, lag
## Следующие объекты скрыты от 'package:base':
##
## intersect, setdiff, setequal, union
x <- c(5.90, 5.90, 5.90, 8.10, 3.90, 4.10, 4.10, 4.10, 6.10, 11.60, 11.60, 5.90, 9.20, 5.90, 5.90, 5.90, 4.10, 5.90, 5.90, 2.30, 9.20, 3.90, 3.90, 6.20, 0.30, 2.90, 3.80, 4.30, 4.80, 4.10, 4.10, 7.90, 1.20, 7.90, 4.10, 4.90, 2.80, 7.80, 3.90, 7.90, 3.30, 6.20, 2.30, 3.90, 0.30, 1.20, 3.90, 7.90, 4.10, 2.70)
glimpse(x)
## num [1:50] 5.9 5.9 5.9 8.1 3.9 4.1 4.1 4.1 6.1 11.6 ...
mean(x)
## [1] 5.08
Получаем результат, показывающий нам, что переменная числовая, состоящая из 50 значений, ее среднее по выборке равно 5.08.
Допустим, по нормативу среднее значение Х должно составлять 6 единиц. Таким образом, нам по имеющейся в нашем распоряжении выборке необходимо проверить нулевую гипотезу: 𝐻0: 𝜇 = 6.
Проведем тест на проверку равенства средней нормативу.
t.test(x,mu = 6)
##
## One Sample t-test
##
## data: x
## t = -2.606, df = 49, p-value = 0.0121
## alternative hypothesis: true mean is not equal to 6
## 95 percent confidence interval:
## 4.370564 5.789436
## sample estimates:
## mean of x
## 5.08
Видим, что p-value=0.0121 меньше, чем α=0,05. Значит, гипотеза 𝐻0: 𝜇 = 6 отвергается на уровне значимости 0,05 (средняя соответствует нормативу).
Наблюдаемое значение статистики критерия: t = -2.606.
По умолчанию нулевая гипотеза проверяется против двусторонней критической области, одностороннюю альтернативу можно указать, задав alternative = “less” (ЛКО) или “greater” (ПКО). Доверительный интервал для генеральной средней по умолчанию строится в рамках теста с надежностью 0,95, в опциях эту вероятность можно изменить. Видим, что он не содержит значение норматива 6.
В нашем случае при использовании двусторонней критической области данного критерия нужна вероятность P[|T|≥t], поэтому задаем в функции вероятность α/2 (по умолчанию она выдает квантиль уровня α, левый хвост, т.е. вероятность P[T <t]= α).
qt(alpha/2, 49)
## [1] -2.009575
Итак, |tнабл.| = 2.606 > t-крит.= 2.01. Значит, гипотеза 𝐻0: 𝜇 = 6 отвергается на уровне значимости 0,05 (средняя вероятнее всего НЕ соответствует нормативу).
После проведения теста на равенство дисперсий (п. 2.1), можно провести тест на равенство средних в двух (независимых) выборках, извлеченных из нормальных совокупностей. В нем используется t-статистика Стьюдента.
Выполнить этот тест можно с помощью уже знакомой нам по п. 2.2.1 функции t.test(). На входе функция получает две выборки — x и y — и два параметра: var.equal и paired. Первый из них отвечает за равенство дисперсий: var.equal = TRUE означает, что был проведён предыдущий тест (п. 2.1) и генеральные дисперсии неизвестны, но равны, вследствие чего будет использоваться t-тест Стьюдента.
В случае же, если гипотеза о равенстве дисперсий отверглась, следует поставить var.equal = FALSE, и будет использоваться t-тест Уэлча (Welch’s t-test) — устойчивый к неравным дисперсиям и объёмам выборок тест, статистика которого основана на выборочных дисперсиях, а не объединённой выборочной дисперсии, а степени свободы аппроксимируются уравнениями Уэлча- Саттертуэйта. Параметр paired = FALSE означает, что выборки независимы. В качестве результата функция t.test() также возвращает t-наблюдаемое, степени свободы df и p-value и, кроме того, оценки средних выборок:
t.test(x,y,var.equal = TRUE, paired = FALSE)
##
## Two Sample t-test
##
## data: x and y
## t = -0.70851, df = 98, p-value = 0.4803
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.4603018 0.6919018
## sample estimates:
## mean of x mean of y
## 5.0800 5.4642
Как можно видеть, p-value=0.4803 больше, чем α=0,05. Значит, гипотеза 𝑯𝟎: 𝝁𝒙 = 𝝁𝒚 не отвергается на уровне значимости 0,05 (средние двух совокупностей с большой долей вероятности равны). На это же указывает и 95%- ный доверительный интервал их разности (𝝁𝒙 − 𝝁𝒚), который содержит 0.
Наблюдаемое значение статистики критерия: t = -0.70851.
Выведем t-критическое:
qt(0.975, 98)
## [1] 1.984467
Итак, |tнабл.| = 0.70851 < t-крит.= 1.984, значит, гипотеза 𝑯𝟎: 𝝁𝒙 = 𝝁𝒚 не отвергается на уровне значимости 0,05 (средние двух совокупностей с большой долей вероятности равны).
Визуализировать различия (сходства) между средними можно, например, с помощью уже упомянутых ящичковых диаграмм (boxplots).
Для построения этого графика организуем данные в виде таблицы:
# Создаем таблицу данных из двух столбцов
my_data <- data.frame(
group = rep(c("x", "y"), each = 50),
something = c(x, y)
)
print(my_data)
## group something
## 1 x 5.90
## 2 x 5.90
## 3 x 5.90
## 4 x 8.10
## 5 x 3.90
## 6 x 4.10
## 7 x 4.10
## 8 x 4.10
## 9 x 6.10
## 10 x 11.60
## 11 x 11.60
## 12 x 5.90
## 13 x 9.20
## 14 x 5.90
## 15 x 5.90
## 16 x 5.90
## 17 x 4.10
## 18 x 5.90
## 19 x 5.90
## 20 x 2.30
## 21 x 9.20
## 22 x 3.90
## 23 x 3.90
## 24 x 6.20
## 25 x 0.30
## 26 x 2.90
## 27 x 3.80
## 28 x 4.30
## 29 x 4.80
## 30 x 4.10
## 31 x 4.10
## 32 x 7.90
## 33 x 1.20
## 34 x 7.90
## 35 x 4.10
## 36 x 4.90
## 37 x 2.80
## 38 x 7.80
## 39 x 3.90
## 40 x 7.90
## 41 x 3.30
## 42 x 6.20
## 43 x 2.30
## 44 x 3.90
## 45 x 0.30
## 46 x 1.20
## 47 x 3.90
## 48 x 7.90
## 49 x 4.10
## 50 x 2.70
## 51 y 2.30
## 52 y 5.90
## 53 y 7.90
## 54 y 4.10
## 55 y 5.00
## 56 y 5.90
## 57 y 4.30
## 58 y 10.60
## 59 y 10.60
## 60 y 2.30
## 61 y 3.90
## 62 y 3.20
## 63 y 1.20
## 64 y 4.10
## 65 y 4.80
## 66 y 3.40
## 67 y 5.30
## 68 y 4.70
## 69 y 9.20
## 70 y 4.10
## 71 y 9.20
## 72 y 8.10
## 73 y 3.90
## 74 y 1.20
## 75 y 5.20
## 76 y 4.80
## 77 y 3.90
## 78 y 4.80
## 79 y 4.90
## 80 y 1.20
## 81 y 0.87
## 82 y 0.84
## 83 y 7.90
## 84 y 3.90
## 85 y 5.90
## 86 y 6.40
## 87 y 5.10
## 88 y 10.60
## 89 y 6.30
## 90 y 9.20
## 91 y 9.20
## 92 y 11.60
## 93 y 3.90
## 94 y 10.60
## 95 y 10.60
## 96 y 4.10
## 97 y 6.40
## 98 y 3.90
## 99 y 4.10
## 100 y 1.80
Теперь можем построить для двух групп (наших двух переменных x и y) ящичковые диаграммы. Для этого понадобится пакет ggpubr, подключим его (верхнее меню: Tools - Install Packages -> ggpubr). И подгрузим его для работы в библиотеку (напоминаем, все пакеты ставятся один раз, а подгружаются каждый раз при запуске файла), вызвав из него затем нужную нам функцию построения ящичковой диаграммы ggboxplot:
library("ggpubr")
## Загрузка требуемого пакета: ggplot2
ggboxplot(my_data, x = "group", y = "something",
color = "group", palette = c("#00AFBB", "#E7B800"),
ylab = "something", xlab = "Groups")
Рис. 4. Ящичковые диаграммы (boxplots) для переменных Х и У, построенные с помощью функции ggboxplot пакета ggpubr
Задача из ДКР2: Медицинское обследование показало, что из 300 случайно отобранных по картотеке детей города 60 болеют гриппом. В предположении, что генеральная совокупность подчиняется биномиальному закону распределения, проверить на уровне значимости 0,05 следующие гипотезы: а) Соответствуют ли полученные данные утверждению медицинского департамента, что заболеваемость среди детей в городе составляет 15%?
n=300 m=60 𝐻0: 𝑝 = 𝑝0 = 0.15 𝐻1: 𝑝 ≠ 𝑝0 (проверяем процент заболеваемости среди детей в городе)
prop.test(60, 300, p = 0.15)
##
## 1-sample proportions test with continuity correction
##
## data: 60 out of 300, null probability 0.15
## X-squared = 5.4967, df = 1, p-value = 0.01905
## alternative hypothesis: true p is not equal to 0.15
## 95 percent confidence interval:
## 0.1571501 0.2507122
## sample estimates:
## p
## 0.2
Надо отметить, что результаты без поправки Йетса не сильно отличаются:
prop.test(60, 300, p = 0.15, correct = FALSE)
##
## 1-sample proportions test without continuity correction
##
## data: 60 out of 300, null probability 0.15
## X-squared = 5.8824, df = 1, p-value = 0.01529
## alternative hypothesis: true p is not equal to 0.15
## 95 percent confidence interval:
## 0.1586569 0.2489289
## sample estimates:
## p
## 0.2
И конечно, в этом случае применима и первая функция с точным биномиальным тестом:
binom.test(60, 300, p = 0.15)
##
## Exact binomial test
##
## data: 60 and 300
## number of successes = 60, number of trials = 300, p-value = 0.01879
## alternative hypothesis: true probability of success is not equal to 0.15
## 95 percent confidence interval:
## 0.1562313 0.2498044
## sample estimates:
## probability of success
## 0.2
Все три теста дали близкие результаты - p-value = 0.01905 (с поправкой Йетса) - 0.01529 все меньше уровня значимости α=0,05 – следовательно, гипотеза 𝐻0: 𝑝 = 𝑝0 = 0.15 отвергается с вероятностью ошибки α, утверждение продюсера нельзя считать справедливым. (95% интервальная оценка вероятности р также не включает в себя заявленную долю 0.15).
Задача из ДКР2: Медицинское обследование показало, что из 300 случайно отобранных по картотеке детей города 60 болеют гриппом. В предположении, что генеральная совокупность подчиняется биномиальному закону распределения, проверить на уровне значимости 0,05 следующие гипотезы: б) Сравнить заболеваемость детей гриппом в этом и соседнем городе, где из 500 случайно отобранных детей 125 болеют гриппом.
n1=300 m1=60 n2=500 m2=125
𝑯𝟎: 𝒑𝟏 = 𝒑2
prop.test(x = c(60, 125), n = c(300, 500),correct = FALSE)
##
## 2-sample test for equality of proportions without continuity correction
##
## data: c(60, 125) out of c(300, 500)
## X-squared = 2.6368, df = 1, p-value = 0.1044
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## -0.10907051 0.00907051
## sample estimates:
## prop 1 prop 2
## 0.20 0.25
qchisq(0.95,2)
## [1] 5.991465
Хи-квадрат наблюдаемое = 2.6368 < хи-квадрат критического=5.991465, p-value = 0.1044 > α=0.05, следовательно, нулевая гипотеза об однородности вероятностей не отвергается с вероятностью ошибки 0.05, с большей долей вероятности заболеваемость детей в городах совпадает.
Можно также для проверки гипотезы использовать модификацию хи-квадрат критерия – стандартной встроенной функции chisq.test.
chisq.test(x, y = NULL, correct = TRUE,
p = rep(1/length(x), length(x)), rescale.p = FALSE,
simulate.p.value = FALSE, B = 2000)
##
## Chi-squared test for given probabilities
##
## data: x
## X-squared = 60.106, df = 49, p-value = 0.1328
Построим эмпирическую плотность наших данных. По ее виду, как и по виду гистограммы можно сделать вывод о предполагаемом теоретическом законе распределения.
library("ggpubr")
ggdensity(x,
main = "Эмпирическая плотность Коэффициент соотношения заёмных и собственных средств 100 предприятий",
xlab = "Коэффициент соотношения заёмных и собственных средств 100 предприятий")
Рис. 5. Эмпирическая плотность распределения объема основных фондов, построенная с помощью функции ggdensity пакета ggpubr
Видим (рис. 5), что плотность очень похожа на кривую Гаусса нормального распределения, поэтому проверять гипотезу о принадлежности изучаемой случайной величины нормальному распределению надо в первую очередь.
Построим график типа квантиль-квантиль (Q-Q plot) для проверки соответствия квантилей выборки квантилям нормального распределения. Если точки эмпирического распределения близки к прямой под углом 45 градусов на графике – это говорит в пользу нормального распределения генеральной совокупности.
ggqqplot(x)
## Warning: The following aesthetics were dropped during statistical transformation: sample
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: sample
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
Рис.6. Квантильный график (Q-Q plot) - эмпирическое распределение квантилей объема основных фондов и квантилей нормального распределения, построенный с помощью функции ggdensity пакета ggpubr
Видим, что Q-Q plot (рис. 6) также позволяет выдвинуть гипотезу о нормальном распределении изучаемой совокупности.
Один из первых и наиболее распространенных критериев согласия,применяемый как для дискретных, так и для непрерывных распределений при условии достаточно большого объема выборки (не менее 50 наблюдений, в нашем случае речь идет о 100 наблюдениях). Основан на сравнении эмпирических и теоретических частот предполагаемого распределения. Как свидетельствуют результаты научных исследований, показывает достаточно неплохие результаты по сравнению с другими критериями, при этом очень прост и понятен в реализации.
x <- c(5.90, 5.90, 5.90, 8.10, 3.90, 4.10, 4.10, 4.10, 6.10, 11.60, 11.60, 5.90, 9.20, 5.90, 5.90, 5.90, 4.10, 5.90, 5.90, 2.30, 9.20, 3.90, 3.90, 6.20, 0.30, 2.90, 3.80, 4.30, 4.80, 4.10, 4.10, 7.90, 1.20, 7.90, 4.10, 4.90, 2.80, 7.80, 3.90, 7.90, 3.30, 6.20, 2.30, 3.90, 0.30, 1.20, 3.90, 7.90, 4.10, 2.70)
y <- c(2.30, 5.90, 7.90, 4.10, 5.00, 5.90, 4.30, 10.60, 10.60, 2.30, 3.90, 3.20, 1.20, 4.10, 4.80, 3.40, 5.30, 4.70, 9.20, 4.10, 9.20, 8.10, 3.90, 1.20, 5.20, 4.80, 3.90, 4.80, 4.90, 1.20, 0.87, 0.84, 7.90, 3.90, 5.90, 6.40, 5.10, 10.60, 6.30, 9.20, 9.20, 11.60, 3.90, 10.60, 10.60, 4.10, 6.40, 3.90, 4.10, 1.80)
result = cor.test(x, y, method = "pearson")
print(result)
##
## Pearson's product-moment correlation
##
## data: x and y
## t = -0.90365, df = 48, p-value = 0.3707
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3935156 0.1545781
## sample estimates:
## cor
## -0.1293345
Получаем, что коэффициент корелляции пирсона равен -0.1293345, корелляционная связь отрицательная и слабая.
Так как p-value = 0.3707 > 0.05, значит, что H0 не отвергается на уровне значимости. То есть, сгенерированная величина подчиняется нормальному распределению.
library(DescTools)
JarqueBeraTest(x, robust = FALSE)
##
## Jarque Bera Test
##
## data: x
## X-squared = 2.6681, df = 2, p-value = 0.2634
library(DescTools)
JarqueBeraTest(y, robust = FALSE)
##
## Jarque Bera Test
##
## data: y
## X-squared = 2.6908, df = 2, p-value = 0.2604
В обоих случаях р-value < 0.9, что говорит о том, что гипотеза о нормальности данных отвергается на уровне значимости 0,05.
library(nortest)
set.seed(1)
x <- rnorm(100, 0, 1)
ad.test(x)
##
## Anderson-Darling normality test
##
## data: x
## A = 0.16021, p-value = 0.9471
Нулевая гипотеза для теста AD состоит в том, что данные действительно следуют нормальному распределению. Таким образом, если наше p-значение для теста ниже нашего уровня значимости (обычно выбирают 0,10, 0,05 и 0,01), то мы можем отклонить нулевую гипотезу и сделать вывод, что у нас есть достаточно доказательств, чтобы сказать, что наши данные не соответствуют нормальному. распределение.
В этом случае наше p-значение равно 0,9471. Поскольку это не ниже нашего уровня значимости (скажем, 0,05), у нас нет достаточных доказательств, чтобы отвергнуть нулевую гипотезу. Можно с уверенностью сказать, что наши данные следуют нормальному распределению.
heights=unique(x)
length(heights)
## [1] 100
m=mean(heights)
s=sd(heights)
ks.test(heights,pnorm,mean=m,sd=s)
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: heights
## D = 0.047014, p-value = 0.9799
## alternative hypothesis: two-sided
library(nortest)
lillie.test(heights)
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: heights
## D = 0.047014, p-value = 0.8479
В этом случае наше p-значение равно 0,8479. Поскольку это не ниже нашего уровня значимости (скажем, 0,05), у нас нет достаточных доказательств, чтобы отвергнуть нулевую гипотезу. Можно с уверенностью сказать, что наши данные следуют нормальному распределению.
library(moments)
set.seed(1234)
x = rnorm(100)
skewness(x)
## [1] 0.5961274
agostino.test(x)
##
## D'Agostino skewness test
##
## data: x
## skew = 0.59613, z = 2.43285, p-value = 0.01498
## alternative hypothesis: data have a skewness