Множественные сравнения в R

После проведения дисперсионного анализа, события развиваются по трем направлениям:

  1. Анализ не показал различий между группами. Сравнение групп на этом закончено.
  2. Изучались две группы, различия между которыми были подтверждены статистически с помощью дисперсионного анализа. Сравнение групп на этом закончено.
  3. Групп больше двух, и дисперсионный анализ привел нас к выводу, что выборки взяты из разных генеральных совокупностей. Далее необходимо выяснить, какие именно группы выбиваются из общей колеи.

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

Задача 4.5

Автор предлагает нам решить задачу 3.2, но отвечать придется на новый вопрос: необходимо найти группы с одинаковой функцией легких.

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

Для вычисления t-статистики у нас уже есть готовая функция, которой не нужны исходные данные эксперимента, а достаточно уже имеющихся средних величин и стандартных отклонений. Эта функция описана в решении задачи 4.4. Нам остается только провести попарные сравнения и подгонку значений p.

Для начала следует подключить файл с функцией, проводящей тест Стьюдента:

# Папка с файлом, где определена функция tmeans(). Измените путь на ваш!!!
SOURCE_DIR <- "~/dev/R_code/glantz/r_2_sample/";
# Имя файла с функцией tmeans(). Измените имя файла, если это необходимо!!!
SOURCE_FILE <- "ex4_4.R";
# Поключить этот файл
source(paste(SOURCE_DIR,SOURCE_FILE,sep=""));


Потом задаем переменные, как показано ниже:

gname <- c("NSC", "NSD", "SL", "SM", "SB"); # имена групп
n     <- rep(200, times=5);                 # размеры выборок
means <- c(3.17, 2.72, 2.63, 2.29, 2.12);   # средние арифметические
sd    <- c(0.74, 0.71, 0.73, 0.70, 0.72);   # стандартные отклонения
p     <- numeric();          # вектор для значений p в тесте Стьюдента
p.adj <- numeric();          # вектор для значений p после поправки Бонферрони
# Задать имена элементам векторов
names(means) <- gname;
names(sd)    <- gname;
names(n)     <- gname;


Теперь попарно сравниваем все группы методом Стьюдента:

cb <- combn(gname, 2); # Комбинации из gname по 2 эелемента
l  <- length(cb[1,]);  # Количество пар групп
for (i in 1:l) {
p  <- c(p, tmeans(means[cb[,i]], sd[cb[,i]], n[cb[,i]], cb[,i])$p.value);
}
names(p) <- NULL; # Удалить имена элементов из вектора


В этой части кода я использовал функцию combn(). Она необходима для отбора из вектора gname всех возможных комбинаций по два элемента. Результат работы этой функции представляет собой матрицу, столбцы которой содержат нужные комбинации. После прохождения цикла переменная p содержит в себе результаты теста Стьюдента. Однако в ней остаются еще и имена элементов вектора, попавшие из переменных, содержащих исходные значения. С помощью names(p) = NULL мы избавляемся от них.

Далее вносим в результаты поправку Бонферрони и сохраняем итог в переменную res:

p.adj <- p.adjust(p, "bonferroni"); # Поправка Бонферрони
res <- cb[,which(p.adj > 0.05)]; # Результат в виде матрицы. Читать по столбцам!


Согласно нашим расчетам, группы NSD и SL не отличаются друг от друга, что говорит об одинаковом влиянии на функции легких как пассивного курения, так и выкуривания небольшого количества сигарет за день. Одной генеральной совокупности принадлежат, к тому же, и группы курильщиков с большим (SB) и средним (SM) количеством сигарет, выкуренных за сутки.

Файл с решением этой задачи приложен к статье.

Задача 4.6

Нам снова придется использовать исходные данные из задачи 3.2. Но если в прошлый раз мы проводили сравнение «многие ко многим», то теперь произведем сравнение «один ко многим». Необходимо определить статистическую значимость различий группы некурящих, работающих в помещении, где не курят (NSC), со всеми остальными группами (NSD, SL, SM, SB).

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

Приступим к вычислениям. Начало как и в прошлой задаче:

# Папка с файлом, где определена функция sampleGen(). Измените путь на ваш!!! "
SOURCE_DIR <- "~/dev/R_code/glantz/anova/"
# Имя файла с функцией sampleGen(). Измените имя файла, если это необходимо!!! "
SOURCE_FILE <- "sampleGen.R"
# Поключить этот файл "
source(paste(SOURCE_DIR,SOURCE_FILE,sep=""));


Подключим библиотеку multcomp. Она нам понадобится в дальнейших вычислениях.

library("multcomp");


Задаем переменные:

gval   <- rep(200, times=5);                  # размеры выборок
gnames <- c("NSC", "NSD", "SL", "SM", "SB"); # имена групп
means  <- c(3.17, 2.72, 2.63, 2.29, 2.12);    # средние величины
devs   <- c(0.74, 0.71, 0.73, 0.70, 0.72);    # стандарные отклонения


Генерируем выборку:

smoke <- sampleGen(gval, gnames, means, devs, round = FALSE);


И завершаем вычисления:

amod      <- aov(var ~ class, data = smoke);
amod_glht <- glht(amod, linfct = mcp(class = "Dunnett"));
res       <- summary(amod_glht);


В этом участке кода используется функция glht() из пакета «multcomp». Данная функция умеет многое, но сейчас нас интересует ее способность выполнять множественные сравнения для параметрических моделей. В качестве аргументов я передал ей объект amod, который является результатом работы функции подгонки модели aov(), а так же указал, что множественное сравнение следует проводить только с контрольной группой. Контрольной будет считаться та группа, которая идет первой в таблице данных для модели. В последней строке мы получаем обобщенный результат и заносим его в переменную res.

Согласно проведенным расчетам, группа некурящих, работающих в помещении, где не курят, значительно отличается от всех остальных групп, в том числе и от пассивных курильщиков.

Файл с кодом R для этой задачи вы найдете в конце статьи.

Задача 4.7

Эта задача основана на условиях задачи 3.3. Решение ее можно точно таким же способом, как и предыдущую. Значения переменных берем из задания 3.3:

gval   <- rep(70, times=3);                  # размеры выборок 
gnames <- c("NS", "Run", "Mar"); # имена групп 
means  <- c(43.3, 58.0, 64.8);   # средние величины 
devs   <- c(14.2, 17.7, 14.3);   # стандарные отклонения


Но поскольку сравнивать необходимо все группы попарно, то еще одна строчка в коде будет отличаться. Вот она:

amod_glht <- glht(amod, linfct = mcp(class = "Tukey"));


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

Задача 4.12

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

Начало наших вычислений не нуждается в дополнительных комментариях:

# Папка с файлом, где определена функция sampleGen(). Измените путь на ваш!!! "
SOURCE_DIR <- "~/dev/R_code/glantz/anova/"
# Имя файла с функцией sampleGen(). Измените имя файла, если это необходимо!!!
SOURCE_FILE <- "sampleGen.R";
# Поключить этот файл
source(paste(SOURCE_DIR,SOURCE_FILE,sep=""));
gval   <- rep(16, 6);                                   # размеры выборок
gnames <- c("Sur1", "Ther1", "Sur2", "Ther2", "Sur3", "Ther3"); # имена групп
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);        # стандарные отклонения 
# генерировать выборку
nurse  <- sampleGen(gval, gnames, means, devs, round = FALSE);


Далее мы формируем объект-модель:

amod <- aov(var ~ class, data = nurse);


И заносим в переменную res результат работы функции:

res <- TukeyHSD(amod, "class", ordered = TRUE);


Вот и все. Думаю, структура вывода результатов этой функции будет понятна и без объяснений.

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