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

Скачайте данные для этого занятия:

experiment.csv

Установим необходимые для занятия пакеты (это действие необходимо, только если пакеты ещё не установлены на ваш компьютер) (для установки выполните следующую строку без символа # вначале).

# install.packages(c('car', 'ggplot2', 'reshape2', 'dplyr', 'compute.es'))

Подгружаем загруженные пакеты для их использования в текущей сессии работы R.

library(car)
library(ggplot2)
library(reshape2)
library(dplyr)
library(compute.es)

Устанавливаем рабочую директорию: через меню или командой

#setwd('...')  
#(вместо '...' укажите путь к папке, в которой хранится файл)

Описание примера эксперимента

Загружаем файл с данными “experiment.csv”. В файле содержаться данные выдуманного эксперимента. Доктор Пилюлькин изобрёл таблетки, стимулирующие когнитивные способности. Чуя их большой потенциальный успех на рынке, доктор Пилюлькин пошёл в фармакологическую компанию, чтобы догоровиться о массовом производстве и продаже своих таблеток. Фармакологическая компания заинтересовалась таблетками, однако несмотря на сильное жаление на них заработать, компания была ответственной, поэтому решила сначала проверить эффективность таблеток. Проверку решили провести двойным слепым экспериментальным методом. 200 добровольцев, согласившихся испытать экспериментальные таблетки на себе, прошли тестирование на интеллект, а потом случайным образом были разделены на две группы. Одной группе дали экспериментальные таблетки, второй - таблетки-пустышки (витаминки). Ни участники эксперимента, ни экспериментатор не знали, какую именно таблетку принял каждый их них. После некоторого времени, в течение которого таблетки должны были подействовать, все 200 человек снова прошли тестирование на интеллект. В файле содержаться три переменные: group - указывает на то, в какой группе был испытуемый (в экспериментальной или в контрольной), pre - результаты тестирования интеллекта до экспериментального воздействия, post - результаты тестирования после экспериментального воздействия.

Загрузка данных

data <- read.table("experiment.csv")

Смотрим, всё ли верно загрузилось.

str(data) 
## 'data.frame':    200 obs. of  3 variables:
##  $ group: Factor w/ 2 levels "control","experiment": 1 1 1 1 1 1 1 1 1 1 ...
##  $ pre  : int  102 78 105 103 111 96 91 129 110 122 ...
##  $ post : int  83 100 95 90 126 96 108 57 121 95 ...

Поиск нехарактерных значений

Построим диаграммы boxplot для обеих групп по первому замеру.

ggplot(data, aes(group, pre)) +
  geom_boxplot()+
  theme_bw() + 
  theme(panel.grid.major.x = element_blank())


В контрольной группе выявлено одно значение, сильно отличающееся от остальных. В этой группе оказался человек, уровень IQ которого ниже 60 баллов. Это низкое значение, однако в принципе такое значение IQ возможно. Пока не будем его удалять, однако будем помнить о нём.

Проверка на нормальность

Проверим распределение наших переменных на нормальность с помощью критиерия Шапиро-Вилка и графически.


### Формальный критерий Проверим на нормальность распределение IQ баллов первичного замера в контрольной группе.

shapiro.test(subset(data, group=="control")$pre) # обратите внимание на фрагмент этой строки ...(subset(data, group=="control")$pre). В данном случае мы совершили отбор значений, относящихся к контрольной группе прямо внутри функции shapiro.test(), что сэкномило нам одну строку кода. 
## 
##  Shapiro-Wilk normality test
## 
## data:  subset(data, group == "control")$pre
## W = 0.98571, p-value = 0.357

Чтобы было понятно, мы только что сделали сразу две операции (одна внутри другой) вместо двух отдельных. А могли бы сделать как раньше: 1) сначала отобрать людей, относящихся к контрольной группы

control_group <- subset(data, group=="control")
  1. и только потом проверить статистическую гипотезу. Результат получается точно такой же.
shapiro.test(control_group$pre)
## 
##  Shapiro-Wilk normality test
## 
## data:  control_group$pre
## W = 0.98571, p-value = 0.357

Сделаем тоже самое для экспериментальной группы

shapiro.test(subset(data, group=="experiment")$pre)
## 
##  Shapiro-Wilk normality test
## 
## data:  subset(data, group == "experiment")$pre
## W = 0.97556, p-value = 0.05974

Во обоих случаях значения p-value превышает конвенциональное значение 0.05. Это означает, что если нулевая гипотеза (H0: Распределение переменной не отличается от нормального) верна, то в обоих случаях вероятность получить вычисленные значения критерия Shapiro-Wilk (W) выше 0.05 (или выше 5%). Следовательно, нет оснований отклонить нулевую гипотезу. Наши переменные следует рассматривать как нормально распределённые.

Задание: Проверьте на нормальность распределение IQ баллов второго замера в контрольной и экспериментальной группах.


Попробуем использовать графические методы

qqPlot(subset(data, group=="control")$pre, xlab="Квантили нормального распределения", ylab="Наблюдаемые квантили")

qqPlot(subset(data, group=="experiment")$pre, xlab="Квантили нормального распределения", ylab="Наблюдаемые квантили")

Все точки лежат на линии или в пределах доверительного интервала, поэтому графический метода также говорит о том, что переменные следует рассматривать как нормально распределённые.


Задание: Проверьте графически нормальность распределения баллов IQ при втором замере в обеих группах


Проверка эквивалентности экспериментальной и контрольной групп

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

Внимание! Одним из базовых предположений, лежищих в основе t-теста Стьюдента в случае сравнения двух выборок, является предположение о равенстве дисперсий анализируемой переменной в этих группах. Поэтому необходимо сначала проверить выполнение этого предположения. Это можно сделать с помощью теста Ливиня. В R это можно сделать с помощью функции leveneTest(pre ~ group, data=data). Если тест показывает, что дисперсии равны, то можно использовать обычный t-тест Стьюдента, если тест показывает, что дисперсии не равны, то необходимо использовать t-тест в модификации Уэлча (Welch Two Sample t-test). Однако симулиционные исследования показывают, что в такой двухшаговой процедуре нет большого смысла, и вместо неё можно сразу использовать t-тест в модификации Уэлча (Welch Two Sample t-test), т.к. он не подвержен смещению оценок в случае, если дисперсии равны.

При выполнении двухвыборочного t-теста функция t.test() по умолчанию принимает, что дисперсии сравниваемых совокупностей не равны, и, как следствие, выполняет t-тест в модификации Уэлча. Если необходимо снять настройку по умолчанию, включающую t-тест Уэлча, необходимо воспользоваться аргументом var.equal = TRUE:

t.test(data$pre ~ data$group, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
## 
##  Welch Two Sample t-test
## 
## data:  data$pre by data$group
## t = 0.12133, df = 195.32, p-value = 0.9036
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.813599  4.313599
## sample estimates:
##    mean in group control mean in group experiment 
##                   100.32                   100.07

Значения t = 0.12133, а p-value = 0.9036, т.е. если нулевая гипотеза (H0: средние значение в двух выборках равны) верна, то вероятность получить вычесленное значение критерия Уэлча (в данном случае t = 0.12133) равна 0.9036, что выше конвенционального значения 0.05. Следовательно, нет оснований отклонить нулевую гипотезу. Средние значение IQ в экспериментальной и контрольной группах следует рассматривать равными друг другу.

Рассчитаем также значение размера эффекта (size-effect). Для сравнения групп с помощью t-теста в качестве меры размера эффекта обычно используют значение d-Коэна (Cohen’s d). В R есть несколько пакетов, с помощью которых можно его расчитать. Воспользуемся пакетом compute.es. d-Коэна (Cohen’s d), как и другие варианты оценки размера эффекта можно рассчитать как на основе значения теста, который используется для сравнения средних (в нашем случае это t-тест), так и на основе информации о средних значениях, значениях стандартных отклонений и объёме групп. В первом случае нужно воспользоваться функцией tes(t, n.1, n.2), во втором - функцией mes(m.1, m.2, sd.1, sd.2, n.1, n.2).

tes(0.12133, 100, 100)
## Mean Differences ES: 
##  
##  d [ 95 %CI] = 0.02 [ -0.26 , 0.3 ] 
##   var(d) = 0.02 
##   p-value(d) = 0.9 
##   U3(d) = 50.68 % 
##   CLES(d) = 50.48 % 
##   Cliff's Delta = 0.01 
##  
##  g [ 95 %CI] = 0.02 [ -0.26 , 0.29 ] 
##   var(g) = 0.02 
##   p-value(g) = 0.9 
##   U3(g) = 50.68 % 
##   CLES(g) = 50.48 % 
##  
##  Correlation ES: 
##  
##  r [ 95 %CI] = 0.01 [ -0.13 , 0.15 ] 
##   var(r) = 0.01 
##   p-value(r) = 0.9 
##  
##  z [ 95 %CI] = 0.01 [ -0.13 , 0.15 ] 
##   var(z) = 0.01 
##   p-value(z) = 0.9 
##  
##  Odds Ratio ES: 
##  
##  OR [ 95 %CI] = 1.03 [ 0.62 , 1.71 ] 
##   p-value(OR) = 0.9 
##  
##  Log OR [ 95 %CI] = 0.03 [ -0.47 , 0.54 ] 
##   var(lOR) = 0.07 
##   p-value(Log OR) = 0.9 
##  
##  Other: 
##  
##  NNT = 206.68 
##  Total N = 200

В результаты выполнения функции были рассчитаны разные варианты оценки размера эффекта. Смотрим на результаты для d: значение d-Коэна = 0.02, границы его 95%-ного доверительного интервала = [-0.26, 0.3]. Для интерпретации сам Коэн предлагал следующие границы:

d = 0.2 - small effect,

d = 0.5 - medium effect,

d = 0.8 - large effect.

Рассчитаем d-Коэна с помощью функции mes(m.1, m.2, sd.1, sd.2, n.1, n.2).

mes(100.32, 100.07, 15.4, 13.69, 100, 100)
## Mean Differences ES: 
##  
##  d [ 95 %CI] = 0.02 [ -0.26 , 0.3 ] 
##   var(d) = 0.02 
##   p-value(d) = 0.9 
##   U3(d) = 50.68 % 
##   CLES(d) = 50.48 % 
##   Cliff's Delta = 0.01 
##  
##  g [ 95 %CI] = 0.02 [ -0.26 , 0.29 ] 
##   var(g) = 0.02 
##   p-value(g) = 0.9 
##   U3(g) = 50.68 % 
##   CLES(g) = 50.48 % 
##  
##  Correlation ES: 
##  
##  r [ 95 %CI] = 0.01 [ -0.13 , 0.15 ] 
##   var(r) = 0 
##   p-value(r) = 0.9 
##  
##  z [ 95 %CI] = 0.01 [ -0.13 , 0.15 ] 
##   var(z) = 0.01 
##   p-value(z) = 0.9 
##  
##  Odds Ratio ES: 
##  
##  OR [ 95 %CI] = 1.03 [ 0.62 , 1.71 ] 
##   p-value(OR) = 0.9 
##  
##  Log OR [ 95 %CI] = 0.03 [ -0.47 , 0.54 ] 
##   var(lOR) = 0.07 
##   p-value(Log OR) = 0.9 
##  
##  Other: 
##  
##  NNT = 206.68 
##  Total N = 200

Результат тот же самый.

Размер эффекта также можно рассчитать с помщью онлайн калькулятора, например, вот этого.

Если бы вы писали курсовую или статью, то вам необходимо было бы написать следующее:

Среднее значение баллов IQ при предварительном тестировании в контрольной группе (M = 100.32, SD = 15.4, n = 100) не отличается от среднего значения баллов IQ в экспериментальной группе (M = 100.07, SD = 13.69, n = 100). Различия между группами (M = 0.25, 95% CI = [-3.81;4.31]) анализировались с помощью t-теста Уэлча, t(195.32) = 0.12, p = 0.904, Cohen’s d = 0.02.

Таким образом, проверка показала, что до экспериментального воздействия выборки были эквивалентными, уровень IQ в них в среднем был одинаковым. Теперь может перейти к ответу на главный вопрос всего эксперимента: Увеличивают ли чудо-таблетки доктора Пилюлькина когнитивные способности, в частности, IQ.


Проверка эффекта таблеток

Для того, чтобы ответить на главный вопрос эксперимента: Увеличивают ли чудо-таблетки доктора Пилюлькина когнитивные способности, в частности, IQ, необходимо сравнить средние значения IQ, измеренного до экспериментального воздействия со средним значением IQ, измеренным после экспериментального воздействия в каждой из двух групп, т.е. до и после приёма чудо-таблеток экспериментальной группой и таблеток-пустышек контрольной. Если среднее значение IQ после экспериментального воздействия увеличится в экспериментальной группе, но не изменится в контрольной группе, то это будет будет свидетельствовать о том, что чудо-таблетки действительно эффективны и стимулируют когнитивные способности, в частности IQ. Проверим это.

В данном случай мы будем сравнивать два измерения, сделанные на одних и тех же людях, поэтому сравниваемые выборки будут связанными, или зависимыми, т.к. каждое значение из одной выборки можно однозначно сопоставить с одним значением из другой выборки. Поэтому необходимо использовать t-тест Стьюдента для двух связанных выборок: та же функция t.test, но с аргументом paired = TRUE

Сначала сравним IQ до и после экспериментального воздействия в экспериментальной группе

t.test(subset(data, group=="experiment")$pre, subset(data, group=="experiment")$post, paired = TRUE, conf.level = 0.95)
## 
##  Paired t-test
## 
## data:  subset(data, group == "experiment")$pre and subset(data, group == "experiment")$post
## t = -7.8378, df = 99, p-value = 5.334e-12
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -18.05806 -10.76194
## sample estimates:
## mean of the differences 
##                  -14.41

Значения t = -7.8378, а p-value = 5.334e-12, т.е. если нулевая гипотеза (H0: средние значение в двух выборках равны) верна, то вероятность получить вычесленное значение парного t-критерия Стьюдента (в данном случае t = -7.8378) равна 5.334e-12, что существенно ниже конвенционального значения 0.05. Следовательно, нужно отклонить нулевую гипотезу. Средние значение IQ до и после приёма чудо-таблеток в экспериментальной группе не равны друг другу.

Размер эффекта.

tes(-7.8378, 100, 100)
## Mean Differences ES: 
##  
##  d [ 95 %CI] = -1.11 [ -1.41 , -0.81 ] 
##   var(d) = 0.02 
##   p-value(d) = 0 
##   U3(d) = 13.38 % 
##   CLES(d) = 21.66 % 
##   Cliff's Delta = -0.57 
##  
##  g [ 95 %CI] = -1.1 [ -1.4 , -0.81 ] 
##   var(g) = 0.02 
##   p-value(g) = 0 
##   U3(g) = 13.47 % 
##   CLES(g) = 21.75 % 
##  
##  Correlation ES: 
##  
##  r [ 95 %CI] = 0.49 [ 0.37 , 0.59 ] 
##   var(r) = 0 
##   p-value(r) = 0 
##  
##  z [ 95 %CI] = 0.53 [ 0.39 , 0.67 ] 
##   var(z) = 0.01 
##   p-value(z) = 0 
##  
##  Odds Ratio ES: 
##  
##  OR [ 95 %CI] = 0.13 [ 0.08 , 0.23 ] 
##   p-value(OR) = 0 
##  
##  Log OR [ 95 %CI] = -2.01 [ -2.55 , -1.47 ] 
##   var(lOR) = 0.08 
##   p-value(Log OR) = 0 
##  
##  Other: 
##  
##  NNT = -5.73 
##  Total N = 200

Если бы вы писали курсовую или статью, то вам необходимо было бы написать следующее:

Среднее значение баллов IQ после экспериментального воздействия в экспериментальной группе (M = 114.48, SD = 15.01, n = 100) выше среднего значения баллов IQ до экспериментального воздействия (M = 100.07, SD = 13.69, n = 100). Различия между замерами (M = 14.41, 95% CI = [-18.06;-10.76]) анализировались с помощью парного t-теста Стьюдента, t(99) = -7.84, p < 0.001, Cohen’s d = -1.11.


Задание: Сравните средние уровни IQ до и после экспериментального воздействия в контрольной группе. Отличаются ли они?


Визуализация результатов эксперимента

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

Для этого необходимо немного преобразовать наши данные, переведя их из “широкого” (wide-format) в “длинный” (long-format) формат. Подробнее о различия можно прочитать, например, здесь: http://seananderson.ca/2013/10/19/reshape.html. Это можно сделать с помощью функции melt из пакета reshape2.

data_for_plot <- melt(data)
## Using group as id variables
data_for_plot <- rename(data_for_plot, measure=variable, IQ=value)

Теперь строим рисунок.

ggplot(data_for_plot, aes(x=measure, y=IQ, color=group)) + # задаём данные и переменные, на основе которых будет построем рисунок 
    stat_summary(fun.data = mean_cl_boot, geom = 'point', size = 5, position=position_dodge(0.3)) + # рисуем точки, обозначающие средние значения
    stat_summary(fun.data = mean_cl_boot, geom = 'errorbar', position=position_dodge(0.3)) + # рисуем доверительные интервалы
    theme_bw() # выбирем ч/б тему для рисунка

Готово. По рисунку видно, что средние значения IQ в экспериментальной и контрольной группах до экспериментального воздействия равны, а после воздействия среднее значение IQ в экспериментальной группе выросло, а в контрольной осталось на том же уровне. Это означает, что чудо-таблетки доктора Пилюлькина действительно стимулируют когнитивные способности, в частности IQ.

Но помните, что это выдуманный эксперимент и данные. На самом деле подобных таблеток не существует, поэтому чтобы успешно сдать экзамен по курсe матметодов, необходимо разбираться в статистике и много практиковаться в R.