В первой части статьи Дисперсионный анализ в R я показал решение четырех задач из третьей главы книги Медико-биологическая статистика. В этой части будет продолжение. И как обычно, в статье будет немного новых функций языка 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, ...) expressionExpression — это выражение на языке 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)

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

> 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=nurse)
> 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 указывает на то, что дисперсии равны. А по условию задачи это именно так.

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