В предыдущей статье по среде статистических вычислений 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))
Значение критерия: -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 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(). Помимо этого, я рассказал о функциях распределения (плотность вероятности, кумулятивная функция распределения, квантили распределения и генератор случайных значений) и мы усовершенствовали навыки написания собственных функций.

Файлы для задач 4.1, 4.3 и 4.4 приложены к статье, кодировка - UTF8.