• Выборочное распределение среднего значения
• Как устроено тестирование гипотез
• t-статистика
• Применение t-теста
• Статистическая значимость

Центральная предельная теорема (ЦПТ) говорит, что если мы возьмем достаточно большую выборку из генеральной совокупности, то среднее значение будет нормально распределено с параметрами \(\mu_{\bar x}\) и \(\sigma _{\bar{x}}\):

\[\bar X \sim N (\mu_{\bar x}, \sigma_{\bar x})\]

При чем \(\sigma_{\bar x} = \sigma/\sqrt{n}\).

Но самое важное, что при больших объемах выборки (\(N > 30\), или даже \(N > 100\)) это так, даже если \(x\) в генеральной совокупности не подчиняется нормальному распределению.

Давайте проверим на опыте, так ли это.

Представим, что данные об алмазах из датасета diamonds (пакет gplot2) — это генеральная совокупность.

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

Diamond Age by Steve Jurvetson on Flickr

1.Должна быть известна \(\sigma\) в генеральной совокупности.

2.Должны выполняться условия, при которых справедлива ЦПТ:

• Наблюдения в выборке должны быть независимы друг от друга.

• Большой объем выборки или нормальное распределение \(x\)

Если \(\sigma\) в генеральной совокупности не известна, ее можно оценить по выборочному стандартному отклонению \(s\).

\[\sigma / \sqrt{n} \approx s/\sqrt{n}\]

После стандартизации:

\[\frac{\bar X - \mu}{SE_{\bar x}} = \frac{\bar X - \mu}{s / \sqrt{n}} \sim t_{df = n - 1}\]

стандартизованное распределение выборочных средних подчиняется \(t\)-распределению с числом степеней свободы \(df = n - 1\)

Условия применимости

Выполняются условия, при которых справедлива ЦПТ:

• Наблюдения в выборке независимы друг от друга.
• Большой объем выборки и нет “выбросов” или нормальное распределение \(x\)

library(ggplot2)
data("diamonds")

# цена бриллиантов хорошего качества огранки
good <- diamonds$price[diamonds$cut == "Good"] 

.mean <- mean(good)                  # выборочное среднее
.n <- length(good)                   # объем выборки
SE <- sd(good)/ sqrt(.n)             # стандартная ошибка
t_crit <- qt(p = 0.975, df = .n - 1) # критич. зн. t для данного n и p = 0.95
err <- t_crit * SE                   # предел погрешности
# Границы доверительного интервала
.mean - err

## [1] 3825.819

.mean + err

## [1] 4031.909

Можем записать среднюю цену бриллиантов хорошей огранки и ее доверительный интервал: \(3928.9 \pm 103\)

theme_set(theme_bw())
ggplot(data = diamonds, aes(x = cut, y = price)) +
  stat_summary(geom = 'pointrange', fun.data = mean_cl_normal)

Посчитайте среднюю цену и доверительный интервал для бриллиантов с такими свойствами:

• качество огранки (cut) идеальное (Ideal)
• прозрачность (clarity) наивысшая (IF)

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

Вычислим границы доверительного интервала

ideal <- diamonds$price[diamonds$cut == 'Ideal' & diamonds$clarity == 'IF'] 

.mean <- mean(ideal)                  # выборочное среднее
.n <- length(ideal)                   # объем выборки
SE <- sd(ideal)/ sqrt(.n)             # стандартная ошибка
t_crit <- qt(p = 0.975, df = .n - 1) # критич. зн. t для данного n и p = 0.95
err <- t_crit * SE                   # предел погрешности
# Границы доверительного интервала
.mean - err

## [1] 2085.768

.mean + err

## [1] 2460.059

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

Способ 1 (некрасивый, требует доработки):

ggplot(data = diamonds, aes(x = cut, y = price, colour = clarity)) +
  stat_summary(geom = 'pointrange', fun.data = mean_cl_normal)

Способ 2 (удовлетворительный):

ggplot(data = diamonds, aes(x = cut, y = price, colour = clarity)) +
  stat_summary(geom = 'pointrange', fun.data = mean_cl_normal) +
  facet_wrap(~ clarity)

Осторожно! Группы могут быть неоднородны. Здесь не учтен размер и т.п.

Различия между выборками не всегда видны невооружённым глазом.

tres caracoles by Alberto Villen on Freeimages.com

Это первый шаг

tres caracoles by Alberto Villen on Freeimages.com

• Нулевая гипотеза \(H_0\) чаще всего формулируется как отсутствие различий между сравниваемыми объектами. Например: Улитки из обеих популяций одинакового размера

• Альтернативная гипотеза \(H_A\) формулируется как присутствие различий, она обратна нулевой гипотезе, т.е. включает все остальные случаи. Например: Улитки из обеих популяций разного размера.

Вне зависимости от нас, реальность может находиться в одном из двух состояний:

• \(H_0\) верна, улитки одинаковы
• \(H_0\) неверна, улитки различаются

После статистического теста мы принимаем решение о том, принять или отвергнуть \(H_0\). Но это решение не обязательно окажется верным. Возможно четыре исхода:

В мире где улитки одинаковы (\(H_0\) верна) мы можем:
- принять \(H_0\) (верное решение),
- отвергнуть \(H_0\) (ошибка).

Аналогично, в мире где улитки различаются (\(H_A\) верна), мы можем:
- принять \(H_0\) (ошибка),
- либо отвергнуть \(H_0\) (верное решение).

1. Формулируем нулевую и альтернативную гипотезы.

Гипотезы выражаются математически в виде тестовых статистик. На этом этапе мы делаем определенные допущения.

2. Проверяем условия применимости тестовой статистики.

3. По реальным данным вычисляем эмпирическое значение тестовой статистики.

Дальше мы должны ответить на вопрос:

Насколько вероятно получить такое или более экстремальное эмпирическое значение, если верна нулевая гипотеза \(H_0\)?

4. Строим теоретическое распределение тестовой статистики для случая, когда верна \(H_0\), и оцениваем по нему уровень значимости.

5. Решаем сохранить или отвергнуть \(H_0\).

Увы, мы не сможем узнать, какая гипотеза верна, но поймем, насколько с ней согласуются исходные данные.

Разберемся с одновыборочным \(t\)-тестом на вымышленном примере.

Представьте, что в одной статье сказано, что средняя плодовитость черепах определенного вида — 7 яиц в кладке.

В вашей выборке из \(33\) черепах — \(\bar x = 8.5\), \(s = 3.2\).

Отличается ли реальная плодовитость в обследованной вами популяции черепах от того, что указано в статье?

Gopher Tortoise by Judy Gallagher on Flickr

• \(H_0: \mu = \mu_0\) — Реальная средняя плодовитость черепах такая, как в статье.
• \(H_A: \mu \ne \mu_0\) — Средняя плодовитость отличается от того, что написано в статье.

\(\mu_0\) — это какое-то конкретное значение. В нашей задаче это — 7 яиц в кладке.

\[t = \cfrac{\bar x - \mu}{ s / \sqrt{n} }\]

Если выполняется ЦПТ, то одновыборочная \(t\)-статистика подчиняется \(t\)-распределению
с числом степеней свободы \(df = n - 1\).

Условия применимости:

• Наблюдения в выборке должны быть независимы друг от друга.
• Объем выборки достаточно велик или величины нормально распределены.

\[t = \cfrac{\bar x - \mu}{ s / \sqrt{n} }\] Средняя плодовитость в выборке из 33 черепах \(\bar x = 8.5\), стандартное отклонение \(s = 3.2\). В статье указана плодовитость \(7\).

\[t = \cfrac{8.5 - 7}{ 3.2 / \sqrt{33} } = 2.7\]

При \(H_0\) значение \(t\) будет близко к нулю.

Насколько необычны значения t меньше или больше \(\pm 2.7\)?

Уровень значимости — это вероятность получить значение \(t\) меньше или больше данного, если бы \(H_0\) была справедлива.

Можно вычислить значение \(p\)

2 * ( 1 - pt(2.7, df = 33 - 1) )

## [1] 0.01098686

Если бы плодовитость не отличалось от указанной в статье, получить \(t\) меньше или больше \(2.7\) можно было бы с вероятностью \(p = 0.011\).

Критический уровень значимости \(\alpha\) — это порог для принятия решений.

Обычно используют \(\alpha = 0.05\).

Если \(p \le \alpha\) — отвергаем \(H_0\) и принимаем \(H_A\).

Если \(p > \alpha\) — сохраняем \(H_0\), не можем ее отвергнуть.

Мы получили \(p < \alpha\), поэтому отвергаем \(H_0\) и принимаем \(H_A\).

Фактическая плодовитость статистически значимо отличается от указанной в статье.

• Правда ли, что \(p\) — вероятность того, что верна сама \(H_0\)?
  
• Нет! Значение \(p\) всегда считается при условии, что \(H_0\) верна. Но \(H_0\) может быть неверна.

• Правда ли, что \(p\) — это вероятность получить такое значение статистики при справедливой \(H_0\)?
• Нет! Вероятность вычисляется как площадь под участком кривой. Конкретное значение статистики — это точка и под ней нет площади.

• Правда ли, что если \(p > 0.05\), то различий между группами на самом деле нет?
  
• Нет! Это значит, что наблюдаемый эффект согласуется с нулевой гипотезой. Различия могут быть.

Как эти ошибки выглядят на графике? Как они взаимосвязаны?

Распределение статистики, когда справедлива \(H_0\), нам уже знакомо — его мы используем в тестах.

Но может быть справедлива \(H_A\) и ее тоже можно описать своим распределением.

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

\(\alpha\) (критический уровень значимости) — это вероятность ошибки I рода.

Если \(H_0\) справедлива, то при \(\alpha = 0.05\) мы отвергаем ее с 5% вероятностью.

Чтобы снизить вероятность таких ошибок, можно уменьшить \(\alpha\).

\(\beta\) — вероятность ошибки II рода.

Считается, что допустимо \(\beta \le 0.2\), но часто про нее забывают.

Если мы уменьшаем \(\alpha\) (график справа), то возрастает \(\beta\).

\(Power = 1 - \beta\) — мощность теста.

Хорошо, когда мощность не меньше \(0.8\).

Мощность теста зависит от величины наблюдаемого эффекта (от величины различий).

\(H_0: \mu_1 - \mu_2 = 0\) — средние значения не различаются в двух группах

\(H_A: \mu_1 - \mu_2 \ne 0\) — средние значения различаются

Т.е. нас интересует разность выборочных средних.

Ее ожидаемое значение при \(H_0\) будет 0.

t-тест в общем виде выглядит так

\[t=\frac{\text{Наблюдаемая величина - Ожидаемое значение}}{\text{Стандартная ошибка}}\]

Двухвыборочный t-тест используется для проверки значимости различий между средними

\[t=\frac{(\bar{x}_1 - \bar{x}_2) - (\mu_1 - \mu_2)}{SE_{\bar{x}_1 - \bar{x}_2}} \; = \; \frac{\bar{x}_1 - \bar{x}_2}{SE_{\bar{x}_1 - \bar{x}_2}}\]

\(SE_{\bar{x}_1 - \bar{x}_2}\) — стандартная ошибка разности двух средних, может рассчитываться по-разному

• t-тест Стьюдента — если считать, что дисперсии в группах равны
• t-тест Уэлча — если считать, что дисперсии могут быть разными

Если группы независимы и дисперсии в них равны, то по центральной предельной теореме

\[SE_{\bar{x}_1 - \bar{x}_2} = \sqrt{ \frac{\sigma^2}{n_{1}} + \frac{\sigma^2}{n_{2}}} \approx \sqrt{ \frac{s^2}{n_{1}} + \frac{s^2}{n_{2}}}\]

\(s^2 = \cfrac{(n_1 - 1)s^2_1 + (n_2 - 1)s^2_2}{n_1 + n_2 - 2}\) — это обобщенная дисперсия по двум выборкам.

Результирующая \(t\)-статистика подчиняется \(t\)-распределению с \(df = n_1 + n_2 - 2\).

Осторожно, равенство дисперсий в группах —
это часто нереалистичное предположение!

Если группы независимы и дисперсии в них неизвестны, то получается

\[SE_{\bar{x}_1 - \bar{x}_2} = \sqrt{ \frac{\sigma^2_{1}}{n_{1}} + \frac{\sigma^2_{2}}{n_{2}}} \approx \sqrt{ \frac{s^2_{1}}{n_{1}} + \frac{s^2_{2}}{n_{2}}}\]

Проблема в том, что эта величина лишь приблизительно следует t-распределению, если считать его степени свободы как обычно для двух групп \(df = n_1 + n_2 - 2\).

Это из-за того, что мы оцениваем две дисперсии
при помощи их стандартных отклонений.

Можно рассчитать по уравнению Уэлча-Саттеруэйта. Это решит проблему.

\[df_{ Welch-Satterthwaite} \approx \cfrac {\bigg(\cfrac{s^2_{1}}{n_{1}} + \cfrac{s^2_{2}}{n_{2}}\bigg)^2} {\cfrac{1}{n_{1} - 1}\bigg(\cfrac {s_{1}^2} {n_{1}}\bigg)^2 + \cfrac{1}{n_{2} - 1}\bigg(\cfrac {s_{2}^2} {n_{2}}\bigg)^2}\]

t-тестом Уэлча можно пользоваться, даже если дисперсии равны.

Он немного консервативнее, чем тест Стьюдента.

Почти такие же, как условия справедливости ЦПТ

• Наблюдения независимы друг от друга.
• Выборки независимы друг от друга (новое условие).
• Объем выборки достаточно велик или величины нормально распределены.

Синдром Кушинга — это нарушения уровня артериального давления, вызванные гиперсекрецией кортизола надпочечниками.

В датасете Cushings (пакет MASS) записаны данные о секреции двух метаболитов при разных типах синдрома (данные из кн. Aitchison, Dunsmore, 1975).

• Tetrahydrocortisone — секреция тетрагидрокортизона с мочой (мг/сут.)
• Pregnanetriol — секреция прегнантриола с мочой (мг/сут.)
• Type — тип синдрома:
  • a — аденома
  • b — двусторонняя гиперплазия
  • c — карцинома
  • u — не известно

Давайте сравним секрецию тетрагидрокортизона при аденома и двусторонней гиперплазии надпочечников. Различается ли она?

library(MASS)
data("Cushings")

Все ли правильно открылось?

head(Cushings)

##    Tetrahydrocortisone Pregnanetriol Type
## a1                 3.1         11.70    a
## a2                 3.0          1.30    a
## a3                 1.9          0.10    a
## a4                 3.8          0.04    a
## a5                 4.1          1.10    a
## a6                 1.9          0.40    a

str(Cushings)

## 'data.frame':    27 obs. of  3 variables:
##  $ Tetrahydrocortisone: num  3.1 3 1.9 3.8 4.1 1.9 8.3 3.8 3.9 7.8 ...
##  $ Pregnanetriol      : num  11.7 1.3 0.1 0.04 1.1 0.4 1 0.2 0.6 1.2 ...
##  $ Type               : Factor w/ 4 levels "a","b","c","u": 1 1 1 1 1 1 2 2 2 2 ...

Есть ли пропущенные значения?

colSums(is.na(Cushings))

## Tetrahydrocortisone       Pregnanetriol                Type 
##                   0                   0                   0

Каковы объемы выборки в каждой группе?

table(Cushings$Type)

## 
##  a  b  c  u 
##  6 10  5  6

Обратите внимание, объемы выборок маленькие.

1.Наблюдения независимы друг от друга?

• Да, независимы. Это случайная выборка.

2.Выборки независимы друг от друга?

• Да, независимы. В группах разные люди (естественно, т.к. тип синдрома у человека может быть только какой-то один).

3.Объем выборки достаточно велик или величины нормально распределены?

• Объем выборки мал

table(Cushings$Type)

## 
##  a  b  c  u 
##  6 10  5  6

Нужно проверить форму распределения в обеих группах.

library(car)
qqPlot(Cushings$Tetrahydrocortisone[Cushings$Type == 'a'])
qqPlot(Cushings$Tetrahydrocortisone[Cushings$Type == 'b'])

Удовлетворительно. При таких малых объемах выборки сложно ожидать лучшего.

Сравним секрецию тетрагидрокортизона при помощи двухвыборочного t-теста.

t.test(x = значения_в_гр.1, 
       y = значения_в_гр.2)

tt <- t.test(x = Cushings$Tetrahydrocortisone[Cushings$Type == 'a'],
             y = Cushings$Tetrahydrocortisone[Cushings$Type == 'b'])
tt

## 
##  Welch Two Sample t-test
## 
## data:  Cushings$Tetrahydrocortisone[Cushings$Type == "a"] and Cushings$Tetrahydrocortisone[Cushings$Type == "b"]
## t = -4.1499, df = 10.685, p-value = 0.001719
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -7.988307 -2.438360
## sample estimates:
## mean of x mean of y 
##  2.966667  8.180000

## 
##  Welch Two Sample t-test
## 
## data:  Cushings$Tetrahydrocortisone[Cushings$Type == "a"] and Cushings$Tetrahydrocortisone[Cushings$Type == "b"]
## t = -4.1499, df = 10.685, p-value = 0.001719
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -7.988307 -2.438360
## sample estimates:
## mean of x mean of y 
##  2.966667  8.180000

• Секреция тетрагидрокортизона значимо различается у пациентов с аденомой и двусторонней гиперплазией надпочечников (\(t_{10.69} = -4.15\), \(p = < 0.05\))

Можно указать в скобках не сравнение с \(\alpha\), а само значение \(p\):
(\(t_{10.69} = -4.15\), \(p = 0.0017\)).

Только не надо безумствовать и указывать слишком много знаков…

t.test(formula = что_зависит ~ группы, data = данные)

Второй вариант синтаксиса можно использовать только, если у вас есть только две группы. Если групп больше, то можно отфильтровать нужные с помощью логического вектора:

t.test(formula = что_зависит ~ группы, data = данные,
       subset = логический_вектор)

t.test(formula = Tetrahydrocortisone ~ Type, data = Cushings, 
       subset = Cushings$Type %in% c('a', 'b'))

Результаты точно такие же, естественно.

Как называются отдельные элементы результатов можно узнать посмотрев их структуру при помощи функции str()

str(tt)

## List of 10
##  $ statistic  : Named num -4.15
##   ..- attr(*, "names")= chr "t"
##  $ parameter  : Named num 10.7
##   ..- attr(*, "names")= chr "df"
##  $ p.value    : num 0.00172
##  $ conf.int   : num [1:2] -7.99 -2.44
##   ..- attr(*, "conf.level")= num 0.95
##  $ estimate   : Named num [1:2] 2.97 8.18
##   ..- attr(*, "names")= chr [1:2] "mean of x" "mean of y"
##  $ null.value : Named num 0
##   ..- attr(*, "names")= chr "difference in means"
##  $ stderr     : num 1.26
##  $ alternative: chr "two.sided"
##  $ method     : chr "Welch Two Sample t-test"
##  $ data.name  : chr "Cushings$Tetrahydrocortisone[Cushings$Type == \"a\"] and Cushings$Tetrahydrocortisone[Cushings$Type == \"b\"]"
##  - attr(*, "class")= chr "htest"

tt$parameter # степени свободы

##       df 
## 10.68512

tt$p.value # доверительная вероятность

## [1] 0.00171934

tt$statistic # значение t-критерия

##         t 
## -4.149899

ggplot(data = Cushings[Cushings$Type %in% c('a', 'b'), ],
       aes(x = Type, y = Tetrahydrocortisone)) +
  stat_summary(fun.data = mean_cl_normal)

Файл aml.csv содержит данные о влиянии регулярной химиотерапии на продолжительность ремиссии.

Прочитаем эти данные

rem <- read.csv(file = "data/aml.csv", header = TRUE)
str(rem)

## 'data.frame':    23 obs. of  2 variables:
##  $ time : int  9 13 13 18 23 28 31 34 45 48 ...
##  $ group: int  1 1 1 1 1 1 1 1 1 1 ...

• В переменной time представлена продолжительность ремиссии в днях.
• group указывает, к какой экспериментальной группе принадлежал пациент. В группе 1 проводилась регулярная химиотерапия, в группе 2 - нет.

Ваша задача сравнить эти группы с помощью t-теста.

t.test(formula = time ~ group, data = rem)

## 
##  Welch Two Sample t-test
## 
## data:  time by group
## t = 1.2741, df = 12.082, p-value = 0.2266
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -12.19413  46.60322
## sample estimates:
## mean in group 1 mean in group 2 
##        38.45455        21.25000

Или так:

t.test(x = rem$time[rem$group == 1],
       y = rem$time[rem$group == 2])

ggplot(data = rem, aes(x = factor(group), y = time)) + 
  stat_summary(fun.data = mean_cl_normal)

Другой вариант кода

rem$group <- factor(rem$group)
ggplot(data = rem, aes(x = group, y = time)) + 
  stat_summary(fun.data = mean_cl_normal)

В датасете sleep содержатся данные об увеличении продолжительности сна по сравнению с контролем после применения двух снотворных препаратов (Cushny, Peebles, 1905, Student, 1908)

Одинаково ли два снотворных влияют на увеличение продолжительности сна?

data(sleep)
head(sleep)

##   extra group ID
## 1   0.7     1  1
## 2  -1.6     1  2
## 3  -0.2     1  3
## 4  -1.2     1  4
## 5  -0.1     1  5
## 6   3.4     1  6

# str(sleep)

1.Наблюдения независимы друг от друга?

• Да, независимы. Это случайная выборка людей.

2.Выборки независимы друг от друга?

• Нет!. В группах одни и те же люди (10 человек, каждый пил оба снотворных).
  ОСТОРОЖНО. Эти данные НЕ ПОДХОДЯТ для двухвыборочного t-теста, т.к. группы не являются взаимно независимыми. Нужно использовать парный t-тест.

3.Объем выборки достаточно велик или величины нормально распределены?

• Объем выборки мал

table(sleep$group)

## 
##  1  2 
## 10 10

Нужно проверить форму распределения в обеих группах.

library(car)
qqPlot(sleep$extra[sleep$group == 1])
qqPlot(sleep$extra[sleep$group == 2])

Хорошо. Можно аппроксимировать нормальным распределением

Сравним увеличение продолжительности сна при помощи парного t-теста.

t.test(formula = extra ~ group, data = sleep, paired = TRUE)

## 
##  Paired t-test
## 
## data:  extra by group
## t = -4.0621, df = 9, p-value = 0.002833
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.4598858 -0.7001142
## sample estimates:
## mean of the differences 
##                   -1.58

## 
##  Paired t-test
## 
## data:  extra by group
## t = -4.0621, df = 9, p-value = 0.002833
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.4598858 -0.7001142
## sample estimates:
## mean of the differences 
##                   -1.58

• Различия изменения продолжительности сна при применении двух препаратов были достоверны (\(t_{9} = -4.06\), \(p = < 0.01\))

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

t.test(formula = extra ~ group, data = sleep)

ggplot(data = sleep, aes(x = factor(group), y = extra)) + 
  stat_summary(fun.data = mean_cl_normal)

• По центральной предельной теореме выборочные средние нормально распределены \[\bar X \sim N (\mu, \sigma / \sqrt{n})\].
• Доверительный интервал к среднему значению можно построить из нормального распределения, если известна \(\sigma\) в генеральной совокупности, или из \(t\)-распределения, если она не известна.
• При тестировании нулевой гипотезы оценивают вероятность получения данного или более экстремального значения тестовой статистики при условии, что \(H_0\) справедлива. Если эта вероятность меньше выбранного критического уровня значимости, то нулевую гипотезу отвергают.
• Для сравнения выборочных средних лучше использовать t-критерий в модификации Уэлча, который учитывает, что дисперсии в группах могут быть разными.
• Средние в независимых выборках нужно сравнивать двухвыборочным t-тестом.
• Если выборки зависимы — нужно это учесть, поэтому для сравнения средних используется парный t-тест.

• OpenIntro: Statistics
• Quinn, Keough, 2002
• Sokal, Rohlf, 1995
• Zar, 1999