Significant differences between multiple variables in R

Question

I have a dataset of particle concentrations recorded at 5 different height. I want to find out whether the differences are significant. For each height, N=15.

What test would be appropriate to use?

I used pairwise.t.test, but am not sure if this is the right solution, as sampling size is quiet small. I also tried pairwise.wilcox.test which returns different p-values and errors "cannot compute exact p-value with ties". Is this due to the small sampling size and can I use it?

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")

Answer 1

If you only want to know if any group means differ significantly, you might want to use an analysis of variance (ANOVA).

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

results in

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

So the result is slightly non-significant at an alpha-level of 5%. The effect size, however, is large at a generalized eta-squared (ges) of 0.123.

The problem with pairwise tests (such as the t-test you mentioned) is that the alpha-error cumulates. In order to account for that alpha error inflation, you would need to reduce the alpha level of the individual tests, resulting in a drastically reduced power.

If the data comes from dependent measures (aka within data), i.e. you measured the same subject multiple times at these heights, you might use a within-subject analysis.

Addition: For a quick visualization, you might want to try

boxplot(df$values~df$height)

Answer 2

You could vectorize the wilcox.exact function of the exactRankTests package, which is able to deal with ties. With that you may apply it on column permutations using outer.

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

First we want to reshape the data into wide format to get columns.

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

Now apply the function on the matrix to get another matrix that gives the p values of the differences.

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

To see significance at a glance just do

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

Answer 3

I agree completely with the answer by @marvinschmitt however I'll show my approach regarding such data.

1. How does the data look like?

boxplot(df$values~df$height)

2. Don't forget about factors! Otherwise results will be wrong.

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

3. Let's bulid a model:

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

and check:

a) Normality:

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

b) Variance:

plot(model.lm, 1)

You can read about these diagnostics plots here

4. Analysis of variance:

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

5. Posthoc test:

(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

You can also take a look about non-parametric multiple testing:

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