Сравнение двух выборок в R

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

В этой статье будет приведено решение первых четырех задач из 4й главы книги С. Гланца «Медико-биологическая статистика».

Задача 4.1.

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

Располагая данными из условия, мы можем вычислить t-критерий (иногда называемый t-статистикой), а затем сравнить с критическим значением. Итак, напишем функцию:

###############################################################
#
# Функция для вычисления t-критерия Стьюдента. Можно сравнивать выборки разного
# размера.
# Расчеты по С.Гланц Медико-биологическая статистика, 1999.
# Функция возвращает значение t-критерия.
#
# Аргументы:
# m - вектор, содержащий средние арифметические групп
# s - вектор, содержащий стандартные отклонения
# n - вектор, содержащий численность групп (размеры выборок)
#
###############################################################
t <- function(m, s, n){

  # Провести проверку длины векторов. Должно быть не более 2.
  if  (length(m)==2 && length(s)==2 && length(n)==2) {

    # Вычислить число степеней свободы
    df <- n[1]+n[2]-2;

    # Вычислить объединенную оценку дисперсии
    SS <- (((n[1]-1)*s[1]^2)+((n[2]-1)*s[2]^2))/df;

    # Рассчитать t-статистику
    T <- (m[1]-m[2])/sqrt(SS/n[1]+SS/n[2]);

    cat("Значение t-критерия: ", T, "\n");
    cat("Число степеней свободы: ", df, "\n");

  }
  else {stop("Сравнивать можно только 2 группы!!!\n")}
}


Как видите, здесь реализован способ расчетов, описанный у С. Гланца, – предварительно проверяется длина векторов. Поскольку метод Стьюдента рассчитан только на две группы, то и значений в векторах должно быть ровно два. Если окажется, что это не так, работа функции t() прекратится, и на экран будет выдано сообщение, записанное в функцию stop(), которая предназначена специально для таких ситуаций.

Для вывода результатов на экран использована функция cat(). А для того, чтобы следующие сообщения показывались с новой строки, я добавил «\n».

Осталось подключить скрипт с этой функцией и можно рассчитать t-критерий:

> source("ex4_1.R")
> t(c(23,45), c(2,7), c(9,16))
Значение критерия:&nbsp;&nbsp;-9.14324
Число степеней свободы:  23
> source("ttest.R")
> m<-c(76.8, 91.4)
> s<-c(13.8, 19.6)
> n<-c(9, 16)
> t(m,s,n)
Значение t-критерия:  -1.968728
Число степеней свободы:  23
> m<-c(2210, 2830)
> s<-c(1200, 1130)
> t(m,s,n)
Значение t-критерия:  -1.288502
Число степеней свободы:  23


Теперь немного теории. Одно из удобств использования R – это наличие большого набора статистических таблиц. Для каждого типа распределения доступно четыре функции: плотность вероятности, кумулятивная функция распределения, квантили распределения и генератор случайных значений из данного распределения.

Плотность вероятности непрерывного распределения показывает вероятность получить значение, близкое к x. Все функции, вычисляющие эту вероятность, начинаются с буквы d. Для распределения t-статистики это dt(). График этой функции для нашей задачи можно построить так:

> curve(dt(x,23),-3,3)


Кумулятивная функция – это вероятность получить величину x или меньшую в данном распределении. Кумулятивные функции для каждого распределения начинаются с буквы p. В нашем случае это pt(). График этой функции можно построить так:

> curve(pt(x,23),-3,3)


Квантильная функция выполняет действия, обратные кумулятивной функции. Она вычисляет величину x при заданной вероятности p. Другими словами, мы задаем ей p-й квантиль и получаем ответ, что с заданной вероятностью p может быть получена величина, равная или меньше чем x. Построим ее график:

> curve(qt(x,23),-3,3)


Если сравнить этот график с предыдущим, то станет ясно, что здесь просто поменялись местами оси x и y.

Для решения задачи нам понадобится квантильная функция распределения qt(p, df), где p – квантиль (или вероятность), df – число степеней свободы.

Как вы знаете, чтобы считать различия между группами, значимыми при двустороннем варианте критерия Стьюдента, величина критерия t должна быть больше или меньше (при t<0) определенного критического значения. При уровне значимости 0,05 критическое значение t вычисляется так:

> qt(0.025, 23)


Думаю, не стоит объяснять, почему был задан 0,025 квантиль, а не 0,05.

Таким образом, задача решена. При уровне значимости 0,05 полученные значения t больше критического, поэтому различия между двумя группами статистически незначимы.

Задача 4.2.

В этой задаче мы располагаем не средними величинами, а первичными данными. Поэтому решать задачу будем несколько по-другому. Вместо вычисления t-статистики применим статистическую функцию t.test().

Как вы знаете, t вычисляется как частное от деления разности выборочных средних на стандартную ошибку разности выборочных средних.

Стандартная ошибка разности выборочных средних (SEDM) может вычисляться двумя способами. В книге С. Гланца описан вариант, при котором необходимо, чтобы дисперсии групп были одинаковы, но можно использовать альтернативный вариант: метод Уэлша (Welch), который не требует равенства дисперсий.

Чтобы подогнать статистику Welch под распределение t-критерия Стьюдента, число степеней свободы вычисляется несколько другим способом и, в отличие от «классического» метода, не является целым. В среде R для функции t.test() по умолчанию используется метод Welch.

Итак, создадим переменные, содержащие данные контрольной и экспериментальной выборок:

P<-scan(what=numeric(), n=11)
N<-scan(what=numeric(), n=11)


Функция scan() принимает значения из файла или стандартного ввода (клавиатуры) и записывает их в переменные. Это значит, что после ее запуска вам надо просто вводить значения переменной, при этом каждое значение вводится с новой строки. Это очень удобная функция для тех случаев, когда у вас есть небольшое количество данных, не записанных в файл, доступный для импорта в R.

Теперь немного об аргументах scan(). Аргумент what задает тип значений, которые будут содержаться в переменной, а n – количество этих значений. Это не единственные доступные аргументы – об остальных можете узнать в справке по функции. Следует отметить, что число значений (n) можно не указывать – в этом случае для прекращения ввода данных необходимо дважды нажать Enter.

Проведем сравнение групп:

> t.test(P,N)

      Welch Two Sample t-test

data:  P and N
t = 3.1405, df = 18.567, p-value = 0.005498
alternative hypothesis: true difference in means methods/html/is.html">is not equal to 0
95 percent confidence interval:
  8.795494 44.113597
sample estimates:
mean of x mean of y
128.9091  102.4545


Вывод функции очень прост для понимания. После названия теста нам показано, какие группы были протестированы:

data:  P and N


Далее показаны: t-статистика, число степеней свободы и точное значение p. Как я уже говорил, число степеней свободы отличается от того, что указано в ответах к задачам в книге, поскольку там используется метод Стьюдента, а не Welch:

t = 3.1405, df = 18.567, p-value = 0.005498


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

alternative hypothesis: true difference in means is not equal to 0


Следующие цифры – доверительный интервал. Если вам необходим 99% доверительный интервал, можно использовать аргумент conf.level=0.99. О доверительных интервалах подробнее будет сказано в статье, которую я напишу к задачам седьмой главы.

95 percent confidence interval:
8.795494 44.113597


Показаны средние для каждой выборки:

sample estimates:
mean of x mean of y
128.9091 102.4545


Из полученных результатов можно сделать вывод, что нифедипин снижает артериальное давление.

Задача 4.3

Эта задача похожа на предыдущую, но решим мы ее несколько иначе. Будем использовать метод Стьюдента, а не Уэлша, и другим способом зададим исходные данные.

Во всех предыдущих задачах ввод данных осуществлялся непосредственно в строке R. Сейчас я покажу, как можно импортировать в переменную таблицы из *csv файла. Файлы этого типа – это обычные текстовые файлы, в которых разделение столбцов задается каким-либо символом (я использую «;»), а строки таблицы соответствую строкам в файле. Создавать такие таблицы можно в любом текстовом редакторе, но проще использовать табличный процессор вроде MS Excel или Calc из пакета OpenOffice.org. Эти программы позволяют сохранить страницу в виде *.csv файла.

Поэтому сейчас я создам таблицу с исходными данными в Calc, а затем сохраню как *.csv. В таблице должно быть два столбца: столбец, содержащий значения зависимых переменных, и столбец с независимыми переменными. Теперь запишем таблицу в переменную:

> data<-read.table("ex4_3.csv", sep=";", dec=",", header=TRUE)


Функция read.table() импортирует таблицу из файла в переменную data. Тип переменной – data.frame. Первый аргумент – это имя файла. Если файл лежит в рабочей директории, то полный путь указывать необязательно. Параметр sep задает разделитель полей, dec – указывает разделитель в десятичных дробях. Этот аргумент следует задавать каждый раз, когда десятичный разделитель в исходном файле не является точкой, как это принято в среде R. И последним указан аргумент header. Если задать ему значение TRUE, то первая строка таблицы будет расцениваться как заголовок и будет использована для задания имен столбцам в переменной.

Теперь перейдем к следующему этапу расчетов. Для того, чтобы провести сравнение групп по «классическому» методу Стьюдента, необходимо удостовериться в равенстве дисперсий:

> var.test(var ~ class, data)


Функция var.test() проверяет гипотезу равенства дисперсий, а точнее то, что их соотношение равно 1. В аргументах задано имя переменной data и т. н. формула модели – это специальный синтаксис для задания статистических моделей. Он используется во многих функциях, в том числе и в тех, которые не применяются для работы со статистическими моделями. Параметры var и class – это имена столбцов, содержащих значения зависимой и независимой переменных соответственно. Поскольку p-value = 0,4153, считаем, что дисперсии равны. Переходим к тесту Стьюдента:

> t.test(var ~ class, data, var.equal=T)


Этот вариант использования функции t.test() отличается только заданием исходных данных в виде формулы и указанием на равенство дисперсий. Значение аргумента var.equal – T. Это не ошибка. Язык R позволяет сокращать TRUE и FALSE до одной буквы. В остальном эта функция вам уже знакома.

Так мы решили еще одну задачу. Ответ – диаметр коронарных сосудов не изменяется под действием нифедипина.

Задача 4.4

Здесь нам предложено решить задачи 3.1 и 3.5 используя критерий Стьюдента. Однако в этих задачах не даны исходные данные, а для функции t.test() требуются именно они, а не средние со стандартными отклонениями. Но выход из этой ситуации есть – напишем свою функцию.

Опыт создания простой функции для вычисления t-статистики у нас уже есть. Сейчас потребуется только усовершенствовать созданную ранее функцию t().

В новой функции, назовем ее tmeans(), будут вычисляться t-критерий, число степеней свободы и значение p. Все это будет выводиться на экран в формате, аналогичном выводу функции t.test(). И так же, как в t.test(), можно будет получить каждый элемент результатов расчета индивидуально – например, так: tmeans(...)$p.value. На данном этапе будет реализован только «классический» метод Стьюдента. Этого вполне достаточно для решения задачи.

###############################################################
#
# Функция для сравнения выборок методом Стьюлента. Можно сравнивать выборки
# разного размера.
# Расчеты по С.Гланц Медико-биологическая статистика, 1999.
# Функция возвращает значение t-критерия, число степеней свободы и p-значение.
# Результат - объект класса htest.
#
# Аргументы:
# m  - вектор, содержащий средние арифметические групп;
# s  - вектор, содержащий стандартные отклонения или ошибки среднего;
# n  - вектор, содержащий численность групп (размеры выборок);
# dn - вектор с именами групп.
# sd - логическая переменная, указывающая содержит ли s стандартные
#      отклонения. По умолчанию - TRUE.
#
###############################################################
tmeans <- function(m, s, n, dn, sd=TRUE) {

  # Провести проверку длины векторов. Должно быть не более 2.
  if  (length(m)==2 && length(s)==2 && length(n)==2) {

    if (sd==FALSE) { s<-s*sqrt(n) }
    # Вычислить число степеней свободы
    df <- n[1]+n[2]-2;
    # Вычислить объединенную оценку дисперсии
    SS <- (((n[1]-1)*s[1]^2)+((n[2]-1)*s[2]^2))/df;
    # Рассчитать t-статистику
    T <- (m[1]-m[2])/sqrt(SS/n[1]+SS/n[2]);
    # Рассчитать значение p
    pval <- 2 * pt(-abs(T), df);

    # Название метода
    method <- paste("Two Sample t-test");
    # Названия групп
    dname <- paste(dn[1], "and", dn[2]);

    names(T)  <- "t";
    names(df) <- "df";
    # Сформировать список с возвращаемыми значениями
    rval <- list(statistic = T,
                 parameter = df,
                 p.value   = pval,
                 method    = method,
                 data.name = dname);
    # Задать класс списку результатов
    class(rval) <- "htest";
    # Вывести результаты
    return(rval);
  }
  else {stop("Сравнивать можно только 2 группы!!!\n")}
}


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

Но если мы просто выведем такой список на экран, то читать его будет неудобно. Поэтому статистические функции, подобные t.test(), возвращают результат в виде объекта класса htest (hypotesis testing). Чтобы получить этот класс, необходимо выполнить два условия: задать элементам списка строго определенные имена и изменить значение атрибута class на «htest», после чего можно выводить результаты на экран.

Используя нашу новую функцию, решим задачи 3.1 и 3.5:

> source("ex4_4.R")
> m<-c(8.5, 13.9)
> s<-c(4.7, 4.1)
> n<-rep(21, 2)
> dn<-c("PE", "Control")
> tmeans(m, s, n, dn)

   Two Sample t-test

data:  PE and Control
t = -3.9676, df = 40, p-value = 0.0002931

> m<-c(51.4, 59.4)
> s<-c(3.2, 3.9)
> n<-rep(36, 2)
> dn<-c("THC", "Control")
> tmeans(m, s, n, dn, sd=FALSE)

   Two Sample t-test

data:  THC and Control
t = -1.5858, df = 70, p-value = 0.1173


Подведем итоги. В этой статье мы научились сравнивать две группы методом Стьюдента или Уэлша с помощью простой функции t.test(). Помимо этого, я рассказал о функциях распределения (плотность вероятности, кумулятивная функция распределения, квантили распределения и генератор случайных значений), и мы усовершенствовали навыки написания собственных функций.

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