I'm doing a 3-way ANOVA that shows me a significant interactive effect of my 3 factors, I anted to do a post-hoc test with emmeans_test to see if the pairs of my different groups are significantly different. The problem is that my data frame look like this (not the entire data frame, but to see how many different groups I have) :
| temperature | quality | quantity | growth |
|---|---|---|---|
| 20 | S | 0.1 | 150 |
| 20 | S | 0.3 | 50 |
| 20 | S | 0.6 | 160 |
| 20 | S | 0.9 | 80 |
| 20 | S | 1.5 | 95 |
| 28 | S | 0.1 | 120 |
| 28 | S | 0.3 | 54 |
| 28 | S | 0.6 | 70 |
| 28 | S | 0.9 | 94 |
| 28 | S | 1.5 | 74 |
| 20 | Y | 0.1 | 200 |
| 20 | Y | 0.3 | 64 |
| 20 | Y | 0.6 | 80 |
| 20 | Y | 0.9 | 87 |
| 20 | Y | 1.5 | 95 |
| 28 | Y | 0.1 | 87 |
| 28 | Y | 0.3 | 55 |
| 28 | Y | 0.6 | 92 |
| 28 | Y | 0.9 | 108 |
| 28 | Y | 1.5 | 121 |
I want to compare by pairs to see if 2 groups are different so I used emmeans_test, but when I run my script, the function group_by group my factor quantity as one unique quantity (0.830) whereas I have 5 differents :
library(rstatix)
> res1 <- donnees_tot_g_J4 %>%
+ group_by(quantity, quality) %>%
+ emmeans_test(growth_rate ~ temperature, p.adjust.method = "bonferroni")
> res1
# A tibble: 2 × 11
quantity quality term .y. group1 group2 df statistic p p.adj p.adj…¹
* <dbl> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
1 0.830 S temperature growth_rate 20 28 98 4.22 0.0000537 5.37e-5 ****
2 0.830 Y temperature growth_rate 20 28 98 3.05 0.00292 2.92e-3 **
# … with abbreviated variable name ¹p.adj.signif
And here is my code :
structure(list( quality = c("S","S", "S", "S", "S", "S", "S", "S",
"S", "S", "S", "Y", "Y", "Y",
"Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "S", "S", "S", "S", "S",
"S", "S", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "S", "S", "S", "S",
"S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S",
"S", "S", "S", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y",
"Y", "S", "S", "S", "S", "S", "S", "S", "Y", "Y", "Y", "Y", "Y",
"Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y",
"Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y",
"S"), quantity = c(0.3, 0.3, 0.3, 0.3, 0.6, 0.6, 0.6, 0.9, 0.9,
0.9, 0.9, 0.3, 0.3, 0.3, 0.6, 0.6, 0.6, 0.6, 0.9, 0.9, 0.9, 0.9,
0.1, 0.1, 0.1, 0.1, 1.5, 1.5, 1.5, 0.1, 0.1, 0.1, 0.1, 1.5, 1.5,
1.5, 0.1, 0.1, 0.1, 0.1, 0.3, 0.3, 0.3, 0.3, 0.3, 0.6, 0.6, 0.6,
0.6, 0.6, 0.6, 0.6, 0.9, 0.9, 0.9, 0.9, 0.1, 0.1, 0.1, 0.6, 0.6,
0.6, 0.6, 0.9, 0.9, 0.9, 0.9, 0.3, 0.3, 0.3, 1.5, 1.5, 1.5, 1.5,
0.3, 0.3, 0.3, 0.3, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5,
1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5,
1.5, 1.5, 1.5, 1.5, 1.5, 0.1), temperature = c(20, 20, 20, 20,
20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28,
28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28,
28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28,
28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28,
28, 28, 28, 28, 28, 28),growth_rate = c(0.303800071834841, 0.405269291837872, 0.316864747623708,
0.303800071834841, 0.352098906853808, 0.37289228981091, 0.382620123855391,
0.444775904433979, 0.409541745203879, 0.396477069415012,
0.472100291701688, 0.451919285802152, 0.444775904433979,
0.437402052060633, 0.52882661034269, 0.533795710392497, 0.541052497667245,
0.510531280665662, 0.533795710392497, 0.622339329512085,
0.627290489167077, 0.601386996896272, 0.148530350448613,
0.22732788723803, 0.191530869757777, 0.181562374726831, 0.354194957377592,
0.375858237686643, 0.34006122020639, 0.17110742318728, 0.160116090892034,
0.160116090892034, 0.210176015665816, 0.609895236463043,
0.591124373918946, 0.58187326824239, 0.0378568479840844,
0.126892202895244, 0.0152090656484838, 0.0152090656484838,
0.0585326533474545, 0.0585326533474545, 0.0775525505364441,
0.126892202895244, 0.0951622248242906, 0.0585326533474545,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 142.231, 0, 0, 142.231,
142.231, 142.231, 142.231, 142.231, 142.231, 142.231, 142.231,
142.231, 142.231, 142.231, 142.231, 142.231, 142.231, 142.231,
142.231, 142.231, 142.231, 142.231, 142.231, 142.231, 142.231,
142.231, 142.231, 0)), row.names = c(NA, -106L), class = c("tbl_df",
"tbl", "data.frame"))
Your quantity variable is numeric, so it gets reduced to its mean by default. If you want the values kept separate, add cov.reduce = FALSE to the emmeans_test call.
I had mistakenly assumed that emmeans_test was patterned after emmeans::emmeans.
Here is (perhaps) an equivalent analysis using the emmeans package itself. (I am using only the first 87 rows of the data, as the data in the OP is corrupt).
> mod <- lm(growth_rate ~ temperature * quantity * quality, data = dat)
> library(emmeans)
> EMM <- emmeans(mod, ~ temperature | quality * quantity, cov.reduce = FALSE)
> test(pairs(EMM), by = NULL, adjust = "bonferroni")
contrast quality quantity estimate SE df t.ratio p.value
temperature20 - temperature28 S 0.1 -0.877 19.5 79 -0.045 1.0000
temperature20 - temperature28 Y 0.1 -56.369 21.6 79 -2.610 0.1084
temperature20 - temperature28 S 0.3 -14.744 15.8 79 -0.935 1.0000
temperature20 - temperature28 Y 0.3 -69.216 17.8 79 -3.896 0.0020
temperature20 - temperature28 S 0.6 -35.543 12.9 79 -2.755 0.0728
temperature20 - temperature28 Y 0.6 -88.486 14.0 79 -6.340 <.0001
temperature20 - temperature28 S 0.9 -56.342 14.9 79 -3.776 0.0031
temperature20 - temperature28 Y 0.9 -107.756 14.2 79 -7.604 <.0001
temperature20 - temperature28 S 1.5 -97.941 27.4 79 -3.573 0.0060
temperature20 - temperature28 Y 1.5 -146.297 24.4 79 -6.007 <.0001
P value adjustment: bonferroni method for 10 tests
However, it is hard to know what model is actually being used by emmesns_test;
and since you seem to want to use both quantity and temperature as factors,
even though they are numeric predictors, perhaps you want the model
lm(growth_rate ~ factor(temperature) * factor(quantity) * quality, data = dat)