Установка пакетов и рабочей директории

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

# install.packages(c('foreign', 'pwr'))

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

library(foreign)
library(pwr)

Устанавливаем рабочую директорию: через меню или командой: Session -> Set working directory -> Choose directory

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

Открываем файл psyhol_15.csv. Он в формате sav (формат, используемый в SPSS), поэтому его можно открыть с помощью функции read.spss из пакета foreign

df <- read.spss("psyhol_15.sav", use.value.labels = T, to.data.frame = T, use.missings = T)
str(df)
## 'data.frame':    37 obs. of  9 variables:
##  $ n           : num  1 2 3 4 5 6 7 8 9 10 ...
##  $ group       : num  143 142 142 141 141 141 143 141 141 143 ...
##  $ sex         : Factor w/ 2 levels "girl","boy": 1 1 1 2 1 1 1 1 1 1 ...
##  $ height      : num  165 156 167 180 175 166 170 165 159 173 ...
##  $ weight      : num  60 50 56 76 61 58 60 53 49 63 ...
##  $ EGE_math    : num  70 68 56 80 82 68 NA NA NA 73 ...
##  $ mark_hist_ps: num  8 8 8 7 7 7 5 3 3 5 ...
##  $ habitation  : Factor w/ 2 levels "Moscow","other": 2 2 1 2 1 1 2 2 2 1 ...
##  $ angle       : num  30 48 60 50 49 50 40 50 90 50 ...
##  - attr(*, "codepage")= int 65001

Тест Уилкоксона

t-тест устойчив к небольшим отклонениям от нормальности распределения, однако если отклонение слишком большие или распределение несимметрично, следует использовать непараметрический аналог - ранговый критерий Уилкоксона.

Посмотрим на распределение оценок по истории психологии ваших старших коллеги.

table(df$mark_hist_ps)
## 
##  3  4  5  6  7  8 
##  2  6 10  6  5  8
barplot(table(df$mark_hist_ps))

Оно явно отличается от нормального и не является симметричным.

Тест Уилкоксона для одной выборки

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

wilcox.test(df$mark_hist_ps, mu=6)
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  df$mark_hist_ps
## V = 220, p-value = 0.5806
## alternative hypothesis: true location is not equal to 6

p-value = 0.5806, т.е. больше 0.05, поэтому нет оснований отклонить H0 об отсутствии различий между средней оценкой по истории психологии и оценкой 6. Попробуйте сравнить с другими оценками.

Тест чувствителен к наличию значений с одинаковыми рангами - невозможно вычислить точное значение p-value, но это не страшно при большом размере выборки.

Тест Уилкоксона для двух независимых выборок

С помощью теста Уилкоксона можно сравнивать и две независимые выборки. Сравним средние оценки по истории психологии между двумя группами (т.к. этот тест может сравнивать только 2 группы). Для этого сначала отберём студентов только двух групп.

two_groups <- subset(df, group==142 | group==143)
wilcox.test(two_groups$mark_hist_ps ~ two_groups$group)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  two_groups$mark_hist_ps by two_groups$group
## W = 113, p-value = 0.05409
## alternative hypothesis: true location shift is not equal to 0

p-value = 0.05409, т.е. больше 0.05, следовательно нет основания отклонять нулевую гипотезу об отсутствии различий между группами.

Тест Уилкоксона для двух зависимых выборок

С помощью теста Уилкоксона можно сравнивать две связанные выборки. Допустим есть оценки какого-то явления по 10-балльной шкале до какого-то события, и после

set.seed(1234)
pre <- round(runif(50, min = 1, max = 10), 0)
pre
##  [1]  2  7  6  7  9  7  1  3  7  6  7  6  4  9  4  9  4  3  3  3  4  4  2
## [24]  1  3  8  6  9  8  1  5  3  4  6  3  8  3  3 10  8  6  7  4  7  4  6
## [47]  7  5  3  8
post <- round(runif(50, min = 1, max = 10), 0)
post
##  [1] 2 4 7 6 2 6 5 8 3 9 9 1 4 1 3 7 4 6 1 6 2 9 1 8 2 6 4 2 4 7 9 5 2 6 3
## [36] 9 5 4 2 9 2 9 2 2 2 6 4 1 4 8

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

wilcox.test(pre, post, paired = TRUE)
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  pre and post
## V = 567.5, p-value = 0.2541
## alternative hypothesis: true location shift is not equal to 0

p-value = 0.2541, т.е. больше 0.05, следовательно нет оснований отклонять нулевую гипотезу об отсутствии различий между средними оценками. Следовательно оценка после события не изменилась.

Сравнение более двух групп: тест Краскела-Уоллиса и тест Фридмана

Часто возникает потребность сравнить более двух групп. Для этого задачи чаще всего используют дисперсионный анализ (ANOVA), но он требует, чтобы распределение зависимой переменной было похоже на нормальное распределение. При существенном отклонении распределения зависимой переменой от нормального следует использовать непараметрические тесты. Если группы независимы, лучше всего подойдет тест Краскела-Уоллиса (Kruskal-Wallis test). Если группы зависимы (например, это результаты повторных измерений или данные эксперимента с блочным дизайном), то нужно использовать тест Фридмана (Friedman test).

Сравним рост в трёх группах с помощью теста Краскела-Уоллиса.

kruskal.test(height ~ group, data=df)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  height by group
## Kruskal-Wallis chi-squared = 1.1262, df = 2, p-value = 0.5694

Про попарные сравнения см. Главу 7 в книге “R в действии” (Kabacoff).

Тест Фридмана вызывается так:

friedman.test(y ~ A | B, data)

где y – это числовая переменная отклика, A – группирующая переменная, B – блочная переменная, в которой содержится информация о парных наблюдениях.

Скачаем данные прямо из интернета:

demo3 <- read.csv("http://www.ats.ucla.edu/stat/data/demo3.csv")

В данных содержаться значения трёх измерений пульса у одних и тех же людей. Переменная id - идентификатор человека. Переменная time - идентификатор номера измерения.

friedman.test(pulse ~ time|id, demo3)
## 
##  Friedman rank sum test
## 
## data:  pulse and time and id
## Friedman chi-squared = 16, df = 2, p-value = 0.0003355

Посмотрим медианы переменной pulse в каждом из трёх измерений.

aggregate(demo3$pulse, list(demo3$time), median)
##   Group.1  x
## 1       1 45
## 2       2 34
## 3       3 19

Анализ мощности

Допусти вы хотели изчучать некоторый эффект. Вы спланировали эксперимент, в котором планируется сравнить средние в двух независимых выборках с помощью t-теста, и теперь хотите понять, сколько человек должно принять участие в неё, чтобы обеспечить достаточную мощность. Минимальной достаточно мощностью принято считать мощность в 80%. Другими словами, вероятность зафиксировать изучаемый эффект, если он действительно существует, должны быть не ниже 0.80. Чтобы рассчитать необходимый объём выборки, можно воспользоваться pwr, в котором есть функция pwr.t.test (есть и для других дизайнов).

pwr.t.test(d=0.5,
           power=.80,
           sig.level=.05,
           type="two.sample",
           alternative="two.sided")
## 
##      Two-sample t test power calculation 
## 
##               n = 63.76561
##               d = 0.5
##       sig.level = 0.05
##           power = 0.8
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

Анализ мощности показал, что для того, чтобы с 80% мощностью зафиксировать эффект (d = 0.5), необходимо, чтобы в каждой группе было по 64 человека.