Дисперсионный анализ в R

Настало время перейти от описания выборок к их сравнению, а точнее к оценке статистической значимости различий между ними. Займемся дисперсионным анализом. Дисперсионный анализ применяют для изучения влияния качественных признаков на количественную переменную. В англоязычной литературе он обозначается как 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)
> # Вектор стандартных отклонений
> # Здесь будут значения трех выборок
> 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<-(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 позволяет одним действием умножить каждое значение вектора на скаляр. И опять задача решена.

Задача 3.5

В этой задаче надо сравнить две одинаковые по численности группы животных. Исходные данные представлены в виде средних и стандартных ошибок среднего.

Как сгенерировать выборку, мы уже знаем. Но каждый раз писать несколько лишних строк кода только лишь для того, чтобы получить ряд цифр, – это непродуктивный путь работы. Есть более удобный и правильный вариант. Можно написать свою функцию, сохранить ее в файле с расширением «.R» и использовать при необходимости.

Команды, сохраненные во внешнем файле, могут быть выполнены с помощью функции source(). Аргумент функции – путь к файлу. Если указать только имя файла, то R будет искать его в рабочем каталоге. По умолчанию, рабочий каталог – это домашняя папка пользователя в Linux или папка «Мои документы» в Windows. Чтобы узнать путь к текущему рабочему каталогу, надо использовать функцию getwd(), для смены каталога – функцию setwd(). Пользователям Linux можно запускать R из той папки, которую надо сделать рабочей.

$ cd /путь/к/рабочему/каталогу
$ R


Итак, создадим файл sampleGen.R в рабочем каталоге и запишем в него свою функцию.

###############################################################
#
# Функция sampleGen генерирует выборки с заданным размером группы,
# средним арифметическим, стандартным отклонением
# (или стандартной ошибкой среднего).
# Функция возвращает объект типа data.frame, состоящий из 2х столбцов.
#
# Аргументы:
# gval - вектор, содержащий размеры групп
# means - вектор со средними арифметическими
# devs - вектор, содержащий стандартные отклонения или ошибки среднего
# sd - логическая переменная, указывающая содержит ли devs стандартные
#      отклонения. По умолчанию - TRUE.
#
###############################################################

sampleGen<-function(gval, means, devs, sd=TRUE) {
  X<-numeric();
  if (sd==FALSE) {
    devs<-devs*sqrt(gval)
  }
  for (i in 1:length(gval)) {
    D<-rnorm(gval[i], 0, 1);
    D<-(D-mean(D))/sd(D);
    x<-c(means[i]+D*devs[i]);
    names(x)<-rep(i, gval[i]);
    X<-c(X, x);
  }
  stack(X);
}


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

name <- function(Arg1, Arg2, ...) expression


Expression – это выражение на языке R, которое использует аргументы Arg1, Arg2 и т. д. для вычисления значения, возвращаемого функцией. Как правило, команды expression группируются скобками {}.

В этой функции я использовал новый оператор. Это оператор условия if. Формат записи:

if (условие) {выражение}


Обратите внимание, что условие равенства обозначено как «==», а не «=». Дело в том, что в языке R знак равенства может использоваться как оператор присваивания, равносильный «<-».

Еще одно новшество – это функция stack(). Данная функция создает из вектора таблицу данных, содержащую два столбца – значения и факторы. Факторы задаются через индексы, как это сделано в функции sampleGen. Можно эту функцию применить и к таблице, содержащей более двух столбцов с данными. Тогда все столбцы объединятся в один, а их имена станут факторами:

> x<-1:10
> y<-11:20
> z<-21:30
> s<-data.frame(x,y,z)
> s<-stack(s)


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

> source("sampleGen.R")
> bpart<-sampleGen(gval=c(36, 36), means=c(51.4, 59.4), devs=c(3.2, 3.9), sd=FALSE)
> model<-aov(values ~ ind, data=bpart)
> summary(model)


И опять новые функции. Функция aov(), как и lm() в одной из прошлых задач, выполняет подгонку статистической модели к данным. Функция summary(), как и anova(), выводит уже известную вам таблицу (т. н. ANOVA table). Более подробно на статистическом моделировании я сейчас останавливаться не буду. Это тема будет освещена позже.

Задача 3.6.

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

Сначала попробуем построить график для получения наглядной информации о разбросе значений:

> nurse<-sampleGen(gval=rep(16, 6), means=c(49.9, 51.2, 57.3, 46.4, 43.9, 65.2), devs=c(14.3, 13.4, 14.9, 14.7, 16.5, 20.5))
> plot(values ~ ind, data=nurse, main="Выраженность синдрома опустошенности у медсестер", xlab="группы", ylab="баллы")


Наверное, и без объяснений все понятно. Сначала были созданы выборки, а затем был построен график. Из графика видно, что последняя группа достаточно далеко отстоит от остальных. Возможно, дисперсионный анализ покажет, что выборки и на самом деле взяты из разных совокупностей.

> model<-aov(values ~ ind, data=bpart)
> summary(model)


Нулевую гипотезу отвергаем. Задача решена.

Задача 3.7.

В этой задаче даны выборки разного размера, среднее и стандартная ошибка среднего.

Нарисуем график другого типа:

> gval<-c(30, 13, 20, 20)
> means=c(15, 15, 9, 7)
> devs=c(1, 2, 2, 1)
> miocard<-sampleGen(gval, means, devs, sd=FALSE)
> stripchart(miocard$values~miocard$ind,
     vert=T,
     main="Вес пораженного участка миокарда\n(в % от веса левого желудочка)",
     xlab="Группы",
     ylab="%",
     group.names=c("Контроль", "Дофамин низ", "Дофамин выс", "Нитропруссид"),
     ylim=c(0, 30))
> arrows(1:4,xbar+2*devs,1:4,xbar-2*devs,angle=90,code=3)


Функция stripchart() выводит маркеры, соответствующие значениям выборки. Затем на этом же графике мы рисуем указатели, показывающие интервал шириной в 4 стандартные ошибки среднего. Если все эти указатели пересекаются, то очень высока вероятность правильности нулевой гипотезы.

Теперь проверим предположение о том, что выборки взяты из разных совокупностей:

> model<-aov(values ~ ind, data=miocard)
> summary(model)


Все встало на свои места, задача решена.

Задача 3.8.

По условию задачи имеется пять групп разной численности. Каждая группа охарактеризована средним и стандартным отклонением. Используем для решения задачи новую функцию oneway.test(). Особенность этой функции - нетребовательность к равенству дисперсий. В остальном ничего нового:

> gval<-c(15, 37, 31, 13, 10)
> means<-c(257, 196, 221, 280, 310)
> devs<-c(159, 359, 340, 263, 95)
> sd<-TRUE
> PlatCell<-sampleGen(gval, means, devs, sd)
> oneway.test(values ~ ind, data=PlatCell, var.equal=T)


Аргумент var.equal указывает на то, что дисперсии равны. А по условию задачи это именно так.

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

Во второй половине статьи (задачи 3.5. – 3.8.)  я рассказал о построении графиков. Немного расширил знания о настройке их параметров. Обратите внимание, что аргументы, использованные в функции stripchart() могут применяться и для других графических функций (например, plot(), boxplot()). Помимо графических, были показаны новые функции для проведения дисперсионного анализа. Наиболее важная часть статьи – описание создания собственных функций и скриптов, что значительно упрощает работу при статобработке результатов своих исследований.

06.05.2019 / 106 / Загрузок: 0 / Ришат Габидуллин / | Теги: ANOVA, статистика, язык R
Всего комментариев: 0
avatar
SixSigmaOnline.ru © 2009-2019            Хостинг от uWeb