I'm getting not significance differences in ANOVA lsmeans, but there are some really significance data. What's wrong with my scripts?

df <- structure(list(value = c(1.693086732, 0.25167691, 1.100272527, 1.60428654, 0.908237338, 1.449864567, 1.06604818, 0.596785144, 0.652925021, 0.453697295, 0.544252785, 1.464221767, 1.043720641, 0.735035158, 0.938875327, 0.712832947, 1.701854524, 1.021094251, 0.564349482, 2.326316679, 1.10170484, 1.075217638, 1.397442796, 0.501086703, 0.675502908, 0.846651623, 1.578086856, 1.857360967, 1.194232629, 1.875837087, 1.106882408, 1.112407609, 1.30479321, 0.637491754, 1.281566883, 1.103115742, 1.895286629, 1.623933836, 0.941989812, 1.30636425, 0.69977606, 1.937055334, 0.666069131, 0.829396619, 0.892844633, 0.573255443, 1.27370148, 0.531593222, 2.782899244, 0.972928201, 0.729463812, 1.121965821, 2.55117084, 0.999302442, 1.138902544, 1.656807007, 0.545349299, 0.550315908, 2.346074577, 0.637551271), gene = c("WT", "WT", "WT", "WT", "WT", "WT", "WT", "WT", "WT", "WT", "WT", "WT", "aox2", "aox2", "aox2", "aox2", "aox2", "aox2", "aox2", "aox2", "aox2", "aox2", "aox2", "aox2", "aox5", "aox5", "aox5", "aox5", "aox5", "aox5", "aox5", "aox5", "aox5", "aox5", "aox5", "aox5", "aox7", "aox7", "aox7", "aox7", "aox7", "aox7", "aox7", "aox7", "aox7", "aox7", "aox7", "aox7", "aox9", "aox9", "aox9", "aox9", "aox9", "aox9", "aox9", "aox9", "aox9", "aox9", "aox9", "aox9"), time = c("0t", "0t", "0t", "1t", "1t", "1t", "3t", "3t", "3t", "5t", "5t", "5t", "0t", "0t", "0t", "1t", "1t", "1t", "3t", "3t", "3t", "5t", "5t", "5t", "0t", "0t", "0t", "1t", "1t", "1t", "3t", "3t", "3t", "5t", "5t", "5t", "0t", "0t", "0t", "1t", "1t", "1t", "3t", "3t", "3t", "5t", "5t", "5t", "0t", "0t", "0t", "1t", "1t", "1t", "3t", "3t", "3t", "5t", "5t", "5t")), row.names = c(NA, -60L), class = c("data.table", "data.frame"))

Scripts I have used:

library(emmeans)
fit <- aov(value ~ gene*time, df)
summary(fit)
em <- emmeans(fit, ~ gene | time)
pairs(em)
pairs(em, adjust=NULL)

When I'm plotting the data in bar chart with Std Err, I see there are really significance in some samples, especially at 3 day. I've already tested this in pairwise t-test and those samples are significantly different at p 0.05.

This is with the pairwise-test, even excel gives the sig results at time point 3.

Could you please try to fix my ANOVA and ls means.

Solution

Not sure why you're expecting a Tukey HSD pairwise comparison to be significant; none of the main effects in the Anova are, nor is there anything obvious in the plot.

For completeness, you can have Tukey HSD significant even when the main effect isn't, but it's rare, and it's important to decide ahead of time which test is more relevant to your question, not pick the one that gives "better" results.

summary(fit)
##             Df Sum Sq Mean Sq F value Pr(>F)
## gene         4  0.873  0.2184   0.625  0.648
## time         3  1.670  0.5568   1.593  0.206
## gene:time   12  1.577  0.1314   0.376  0.965
## Residuals   40 13.984  0.3496     

library(ggplot2)
sem <- summary(em)
sem$gene <- relevel(factor(sem$gene), 'WT')
ggplot(sem) + aes(gene, emmean, ymin=lower.CL, ymax=upper.CL) + geom_pointrange() +
  facet_wrap(~time, nrow=1) + ggtitle("emmeans with 95% CIs")

The plot shows raw data in red (that'll be important).

By the way, none of the pairwise t-tests are significant at the .05 level either, not even without correction for multiple comparisons. As in your code, I test pairwise for genes at each time.

library(broom)
library(tidyr)
library(dplyr)
options(digits=2)
lapply(split(df, df$time), function(x) {
  data.frame(time=x$time[1], tidy(pairwise.t.test(x$value, x$gene, p.adjust.method="none")))
}) %>% bind_rows() %>% spread(time, p.value)
##    group1 group2   0t   1t   3t   5t
## 1    aox5   aox2 0.82 0.32 0.71 0.97
## 2    aox7   aox2 0.32 0.73 0.21 0.70
## 3    aox7   aox5 0.43 0.50 0.37 0.67
## 4    aox9   aox2 0.31 0.40 0.60 0.71
## 5    aox9   aox5 0.42 0.86 0.88 0.74
## 6    aox9   aox7 0.99 0.62 0.45 0.45
## 7      WT   aox2 0.85 0.72 0.20 0.74
## 8      WT   aox5 0.97 0.51 0.34 0.71
## 9      WT   aox7 0.41 0.99 0.95 0.96
## 10     WT   aox9 0.40 0.63 0.42 0.49

However, if you run this without pooling the variance (or again, without correcting for multiple comparisons), you do get the result you expected, with the p-value for aox5 and aox7 at 3t = 0.016.

However, this is the WRONG thing to do for those two reasons; first, you really should correct for multiple comparisons, and second, you don't have nearly enough data to estimate the variance well enough, so pooling is important. You can see that in the raw data; those two happen to have tighter data than the others, but that's almost certainly due to random chance.

##   group1 group2   0t   1t    3t   5t
## 1    aox5   aox2 0.70 0.25 0.794 0.96
## 2    aox7   aox2 0.17 0.73 0.413 0.61
## 3    aox7   aox5 0.32 0.49 0.016 0.53
## 4    aox9   aox2 0.46 0.52 0.744 0.79
## 5    aox9   aox5 0.56 0.89 0.868 0.80
## 6    aox9   aox7 0.99 0.71 0.428 0.59
## 7      WT   aox2 0.82 0.65 0.398 0.70
## 8      WT   aox5 0.97 0.36 0.096 0.65
## 9      WT   aox7 0.41 0.99 0.892 0.95
## 10     WT   aox9 0.57 0.69 0.409 0.63

So the short answer is mostly that the ANOVA with pairwise tests doesn't match up with your pairwise t-tests because ANOVA pools the SD but your t-tests didn't, but the bigger answer is that you should both pool and correct for multiple comparisons, so the code you had at the beginning, namely:

pairs(em)

is the right thing to do.