Блог

У меня есть набор данных о концентрациях частиц, зарегистрированных на 5 разных высотах. Я хочу выяснить, являются ли различия значительными. Для каждой высоты N = 15.

Какой тест было бы целесообразно использовать?

Я использовал pairwise.t.test, но не уверен, что это правильное решение, так как размер выборки очень мал. Я также пробовал попарно.wilcox.test, который возвращает разные значения p и ошибки «не удается вычислить точное значение p с помощью связей». Связано ли это с небольшим размером выборки и могу ли я его использовать?

mydata:

 structure(list(height = c(1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 
5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 
1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 
2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 
3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 
4L, 5L), values = c(1.67, 3.33, 6.67, 10, 15, 25, 20, 11.67, 
16.67, 18.33, 1.67, 0, 1.67, 5, 3.33, 5, 73.33, 8.33, 5, 5, 10, 
5, 6.67, 6.67, 3.33, 18.33, 18.33, 6.67, 38.33, 0, 23.33, 10, 
15, 11.67, 5, 11.67, 8.33, 1.67, 15, 3.33, 13.33, 10, 10, 3.33, 
10, 8.33, 21.67, 10, 41.67, 8.33, 3.33, 36.67, 15, 11.67, 8.33, 
8.33, 8.33, 5, 5, 0, 1.67, 8.33, 16.67, 3.33, 10, 16.67, 8.33, 
8.33, 25, 1.67, 6.67, 26.67, 3.33, 11.67, 1.67)), row.names = c(NA, 
-75L), class = "data.frame")
 

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

 library(afex)
df$id = 1:nrow(df)
aov_ez(data=df, id="id", between="height", dv="values")

 

приводит к

 Anova Table (Type 3 tests)

Response: values
  Effect    df    MSE      F  ges p.value
1 height 4, 70 118.38 2.45   .123    .054
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘ ’ 0.1 ‘ ’ 1
 

Таким образом, результат незначителен при альфа-уровне 5%. Однако размер эффекта велик при обобщенном eta-квадрате ( ges ) 0,123.

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

Если данные получены из зависимых измерений (иначе внутри данных), т. Е. Вы измеряли один и тот же объект несколько раз на этих высотах, вы можете использовать анализ внутри объекта.

Дополнение: Для быстрой визуализации вы можете попробовать

 boxplot(df$values~df$height)
 

Вы могли бы векторизовать wilcox.exact функцию exactRankTests пакета, которая способна работать со связями. При этом вы можете применить его к перестановкам столбцов, используя outer .

 wilcox.testv <- Vectorize(function(x, y) 
  exactRankTests::wilcox.exact(m[,x], m[,y])$p.value)
 

Сначала мы хотим преобразовать данные в широкий формат, чтобы получить столбцы.

 m <- as.matrix(reshape(transform(d, id=cumsum(height == 1)), timevar="height", 
                       direction="wide")[-1])
m
#    values.1 values.2 values.3 values.4 values.5
# 1      1.67     3.33     6.67    10.00    15.00
# 6     25.00    20.00    11.67    16.67    18.33
# 11     1.67     0.00     1.67     5.00     3.33
# 16     5.00    73.33     8.33     5.00     5.00
# 21    10.00     5.00     6.67     6.67     3.33
# 26    18.33    18.33     6.67    38.33     0.00
# 31    23.33    10.00    15.00    11.67     5.00
# 36    11.67     8.33     1.67    15.00     3.33
# 41    13.33    10.00    10.00     3.33    10.00
# 46     8.33    21.67    10.00    41.67     8.33
# 51     3.33    36.67    15.00    11.67     8.33
# 56     8.33     8.33     5.00     5.00     0.00
# 61     1.67     8.33    16.67     3.33    10.00
# 66    16.67     8.33     8.33    25.00     1.67
# 71     6.67    26.67     3.33    11.67     1.67
 

Теперь примените функцию к матрице, чтобы получить другую матрицу, которая дает p значений различий.

 cols <- 1:ncol(m)
res <- outer(cols, cols, wilcox.testv)
res
#           [,1]       [,2]      [,3]       [,4]       [,5]
# [1,] 1.0000000 0.32724202 0.6582911 0.47820691 0.14360144
# [2,] 0.3272420 1.00000000 0.1431578 0.81358101 0.01930055
# [3,] 0.6582911 0.14315777 1.0000000 0.29689457 0.18766290
# [4,] 0.4782069 0.81358101 0.2968946 1.00000000 0.02072233
# [5,] 0.1436014 0.01930055 0.1876629 0.02072233 1.00000000
 

Чтобы сразу увидеть значимость, просто выполните

 alpha <- .05
res < alpha
#       [,1]  [,2]  [,3]  [,4]  [,5]
# [1,] FALSE FALSE FALSE FALSE FALSE
# [2,] FALSE FALSE FALSE FALSE  TRUE
# [3,] FALSE FALSE FALSE FALSE FALSE
# [4,] FALSE FALSE FALSE FALSE  TRUE
# [5,] FALSE  TRUE FALSE  TRUE FALSE
 

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

1. Как выглядят данные?

 boxplot(df$values~df$height)
 

2. Не забывайте о факторах! В противном случае результаты будут неверными.

 str(df)
df$height <- as.factor(df$height)
 

3. Давайте построим модель:

 model.lm = lm(values ~ height, data=df)
 

и проверьте:

а) Нормальность:

 hist(resid(model.lm))
plot(model.lm, 2)
 

б) Дисперсия:

 plot(model.lm, 1)
 

Вы можете прочитать об этих диаграммах диагностики здесь

4. Дисперсионный анализ:

 a1 <- aov(model.lm)
summary(a1)
 

5. Постхок-тест:

 (TukeyHSD(a1, 'height', conf.level=0.95))
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = model.lm)

$height
          diff        lwr        upr     p adj
2-1   6.888000  -4.236727 18.0127273 0.4204011
3-1  -1.888000 -13.012727  9.2367273 0.9893557
4-1   3.667333  -7.457394 14.7920606 0.8870422
5-1  -4.112000 -15.236727  7.0127273 0.8382557
3-2  -8.776000 -19.900727  2.3487273 0.1885170
4-2  -3.220667 -14.345394  7.9040606 0.9265135
5-2 -11.000000 -22.124727  0.1247273 0.0540926
4-3   5.555333  -5.569394 16.6800606 0.6307915
5-3  -2.224000 -13.348727  8.9007273 0.9803501
5-4  -7.779333 -18.904061  3.3453940 0.2972209
 

Вы также можете взглянуть на непараметрическое множественное тестирование:

 kruskal.test(values ~ height, data=df)