Настало время перейти от описания выборок к их сравнению, а точнее к оценке статистической значимости различий между ними. Займемся дисперсионным анализом. Дисперсионный анализ применяют для изучения влияния качественных признаков на количественную переменную. В англоязычной литературе он обозначается как ANOVA (Analisis of Variances). Есть еще MANOVA (Multiple Analysis of Variance). Это — многомерный дисперсионный анализ, но сейчас мы им заниматься не будем.

Дисперсионный анализ был разработан еще в 20-х годах прошлого века Р. Фишером. Позднее на его основе были созданы и другие критерии значимости. Для проведения этого анализа, как вы знаете, требуется выполнение ряда условий, в том числе нормальность распределения, и равенство дисперсий. Однако существуют и непараметрические аналоги дисперсионного анализа, но в третьей главе изучаемой нами книги подразумевается, что все данные соответствуют требованиям однофакторного дисперсионного анализа.

А теперь перейдем к решению задач.

Задача 3.1.

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

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

Чтобы узнать какие пакеты у вас уже установлены наберите

> library()

Вы увидите большой список пакетов. В каждом пакете хранится набор функций, который сможете использовать после подключения этого пакета. Для подключения, просто наберите library(имя_пакета). Например:

> library(boot)

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

Но как узнать какой пакет содержит нужную нам функцию? Просто зайдите на сайт CRAN (The Comprehensive R Archive Network) и в разделе Search начните поиск интересующей вас функции. Кроме этого, искать можно по задачам, которые стоят перед вами и по названию пакета. На сайте есть соответствующие разделы. Учтите, что вся информация представлена на английском языке.

Итак, для решения задач из этой главы нам потребуется пакет HH. Чтобы установить его наберите:

> install.packages("HH")

Теперь вам будет предложено выбрать зеркало CRAN откуда R скачает и установит пакет.

После установки подключаем пакет.

> library("HH")

Теперь начнем работу. Чтобы решить задачу достаточно ввести одну команду:

> anova.mean(c("Gel", "Control"), c(21, 21), c(8.5, 13.9), c(4.7, 4.1), ylabel="hours")

В результате выполнения функции получаем таблицу. Способ представления данных, а точнее терминология, в ней отличается от той, что используется при расчете критерия Фишера F в книге Гланца. Несмотря на это, сразу видно, что последнее значение в первой строке и есть заветное значение p, из которого можно сделать вывод о статистической значимости различий. F value тоже не вызывает сомнений, а в остальном сейчас разберемся.

В первой строке указаны межгрупповые число степеней свободы (Df), сумма квадратичных отклонений (Sum of Sq) и среднее квадратичных отклонений (дисперсия, Mean Sq), а в нижней внутригрупповые. Если и сейчас не понятно, то советую вспомнить как вычисляется дисперсия. Дисперсия — это отношение суммы квадратичных отклонений к числу степеней свободы.

Вот и все. Задача решена.

Задача 3.2.

Как решить это задание вам уже понятно.

> anova.mean(c(NSC, NSD, SL, SM, SB), rep(200, times=5), c(3.17, 2.72, 2.63, 2.29, 2.12), c(0.74, 0.71, 0.73, 0.70, 0.72), ylabel="MOSSV")

Вектор (NSC, NSD, SL, SM, SB) — это названия групп. Новой для вас является функция rep(). Эта функция повторяет значение 200 пять раз подряд. Использовал я ее только для того, чтобы сократить количество написанных цифр. И снова задача решена.

А теперь еще один интересный трюк с функциями в языке R. Допустим вас не интересует вся таблица, выводимая функцией anova.mean() и вы хотите получить конкретные ее элементы. Делать это проще, записав результат работы функции в переменную.

> anova.table<-anova.mean(c(NSC, NSD, SL, SM, SB), rep(200, times=5), c(3.17, 2.72, 2.63, 2.29, 2.12), c(0.74, 0.71, 0.73, 0.70, 0.72), ylabel="MOSSV")

Для выбора подмножества элементов вектора или матрицы язык R имеет систему индексов. Индексный вектор добавляется в квадратных скобках сразу после имени объекта. Например, создадим вектор x.

> x<-1:10

А теперь выберем из него 3, 4 и 5 элементы.

> x[3:5]

Или 1, 6 и 9 элементы

> x[с(1, 6, 9)]

В случае с матрицами, в квадратных скобках можно указать несколько векторов (количество зависит от размерности матрицы). Наша переменная anova.table как раз и содержит две строки и пять столбцов. Поэтому для вывода первой строки надо написать следующее.

> anova.table[1,]

А для выбора первого столбца пишем так.

> anova.table[,1]

Если мы хотим вывести только значение p, выбираем пятый элемент из первой строки.

> anova.table[1,5]

Но это еще не все. Индекс может быть как числом, так и строкой. Это возможно в том случае, когда заданы имена строк или столбцов. Посмотрим, какие имена есть в нашей переменной.

> names(anova.table)

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

> anova.table[,"Mean Sq"]

Еще раз обратите внимание, что перед запятой в квадратных скобках ничего не указано. Это значит, что надо взять все строки в столбце с именем Mean Sq.

 Задача 3.3.

Проводить дисперсионный анализ, когда есть уже готовые среднее и дисперсия — несколько странно. Хорошо, что С. Гланц уже подсчитал за нас эти параметры, но в реальной жизни такого не происходит. Мы всегда имеем дело с выборочными значениями. Конечно, можно вычислить среднее арифметическое и дисперсию, а затем передать их значения в функцию anova.mean(), но зачем делать лишние вычисления самому? Пусть этим займется компьютер. Тем более что язык R располагает для этого всем необходимым.

Но ведь у нас нет исходных данных для задач этой главы книги! Ничего страшного. Мы просто сгенерируем выборки с известными средними и дисперсиями. Таким образом, мы как бы проведем свой виртуальный эксперимент. По сути, это немного усложненный вариант решения задачи из статьи «R. Описательная статистика», когда мы программно бросали виртуальные игральные кости.

Описанный ниже способ формирования выборки с заданными параметрами, был взят мною из книги Crawley M.J. The R Book.

 В языке R есть функция rnorm(), которая генерирует заданное число случайных величин из нормального распределения с установленными параметрами. Например, чтобы получить десять случайных значений из нормального распределения со средним равным нулю и стандартным отклонением равным единице, надо записать следующее.

> rnorm(10, 0, 1)

Обратите внимание, что выборка получена из распределения с указанными параметрами, т. е. 0 и 1 — это параметры генеральной совокупности. Но у нас есть только выборочные средние и дисперсии, а значит эта функция нам не совсем подходит. Чтобы разобраться в ситуации, начнем решать задачу.

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

> devs<-rnorm(70, 0, 1)
> devs<-(devs-mean(devs))/sd(devs)

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

HDLP<-43.3+devs*14.2

Но это только одна выборка, а нам надо три. Чтобы не повторять три раза одно и то же воспользуемся средствами языка R.

> # Вектор средних
> Ms<-c(43.3, 58.0, 64.8)
> # Вектор стандартных отклонений
> SDs<-c(14.2, 17.7, 14.3)
> # Здесь будут значения трех выборок
> HDLP<-numeric()
> devs<-rnorm(70, 0, 1)
> devs<-(devs-mean(devs))/sd(devs)
> for (i in 1:3) {
+ # Создаем вектор значений
+ HDLP<-c(HDLP, Ms[i]+devs*SDs[i])
+ }
> # Таблица с данными
> THDLP<-data.frame(HDLP=HDLP, Groups=rep(c("NS", "Run", "Mar"), rep(70, 3)))

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

Далее идет оператор цикла for. Он задает переменной i значения от одного до трех и выполняет вычисления, заданные в фигурных скобках. Сами фигурные скобки в языке R обозначают группу команд. Более подробную информацию смотрите в справке (наберите ?''for'')

Последняя строка формирует таблицу данных, которая будет использоваться для дисперсионного анализа. Функция data.frame() создает объект из двух столбцов. В первом - значения выборок, во втором — обозначения каждой из групп («фактор» в терминологии языка R). Перед знаком равенства стоят названия столбцов, после — их содержимое.

Осталось провести сам однофакторный дисперсионный анализ.

> anova(lm(HDLP~Groups, THDLP))

Функция anova() в качестве параметра принимает объект, являющийся линейной регрессионной моделью. Функция lm() как раз и отвечает за подгонку линейной модели. Первый параметр функции — это формула модели: слева от тильды имя столбца с данными, справа — фактор или группирующая переменная. Второй параметр — таблица данных для анализа.

Задача 3.4.

У этой задачи есть одна особенность — вместо стандартных отклонений даны ошибки средних. Но поскольку стандартная ошибка средней это частное от деления стандарного отклонения и квадратного корня из численности выборки, то то проблем у нас не возникнет.

> # Вектрор выборочных средних
> m<-c(85.1, 83.5, 80.9, 72.6, 60, 73.5, 63.8)
> # Вектор станданых ошибок выборочных средних
> me<-c(0.3, 1, 0.6, 0.7, 1.3, 0.7, 2.6)
> # Вычислить стандартные отклонения
> sd<-me*sqrt(36)
> # Провести дисперсионный анализ
> anova.mean(c("0","15","30","50","75","75clean","150"), rep(36, 7), m, sd, ylabel="antibacterial")

Немного пояснений. Функция sqrt(), как вы уже догадались, вычисляет квадратный корень. Затем это значение умножается на стандартную ошибку среднего. Обратите внимание, что язык R позволяет одним действием умножить каждое значение вектора на скаляр. И опять задача решена.

Подведем итоги первой части статьи. Здесь я рассказал о системе пакетов среды R, о новом типе данных (data.frame) и генерации случайных и не случайных последовательностей. А так же об индексных векторах и циклической структуре for. Научил проводить однофакторный дисперсионный анализ двумя разными способами: основываясь на выборочных значениях и суммарных данных. Последние четыре задачи будут решены во второй части статьи.