2
$\begingroup$

can anyone help me to understand this. I have data to compare among 3 experiments with 2 treatments each.I want to compare a measured variable between the 2 conditions for each experiment. I first ran an ANOVA (significant) then a post hoc tukey test. One of comparison is not significant although the values between the conditions are different (please see the boxplot, exp1). I pasted a part of the post hoc output with the result confusing me Exp1:U-Exp1:B where pvalue>0.05.

enter image description here

can anyone help me understand this?

test <- structure(list(exp = c("Exp1", "Exp1", "Exp1", "Exp1", "Exp1", 
                   "Exp1", "Exp2", "Exp2", "Exp2", "Exp2", "Exp2", "Exp2", "Exp2", 
                   "Exp2", "Exp3", "Exp3", "Exp3", "Exp3", "Exp3"), cond = c("B", 
                                                                             "B", "B", "U", "U", "U", "B", "B", "B", "B", "U", "U", "U", "U", 
                                                                             "B", "B", "B", "U", "U"), variable = c(0.00838203, 0.0103495, 
                                                                                                                    0.00757493, 0.02157368, 0.0132083, 0.01336677, 0.03054078, 0.01570897, 
                                                                                                                    0.028895, 0.02730669, 0.05822746, 0.05476223, 0.05476223, 0.05814691, 
                                                                                                                    0.00358898, 0.00721144, 0.01070452, 0.00348329, 0.00613196)), class = c("tbl_df", 
                                                                                                                                                                                            "tbl", "data.frame"), row.names = c(NA, -19L))


ggboxplot(test , x="exp", y='variable',fill='cond')+facet_wrap(.~exp,scales = 'free')

#run an anova 
# Compute the analysis of variance
anova <- aov(variable ~ exp*cond, data = test )
# Summary of the analysis
summary(anova)

#apply a post hoc tukey test ukey HSD (Tukey Honest Significant Differences, R function: TukeyHSD()) 
TukeyHSD(anova)

>  Tukey multiple comparisons of means
    95% family-wise confidence level
> 
> Fit: aov(formula = variable ~ exp * cond, data = test)
> 
> $exp
>     >               diff         lwr           upr     p adj
>     > Exp2-Exp1  0.028634582  0.02268322  0.0345859455 0.0000000
>     > Exp3-Exp1 -0.006185164 -0.01285797  0.0004876468 0.0706136 
>     > Exp3-Exp2 -0.034819746 -0.04110199 -0.0285375024 0.0000000
>     > 
>     > $cond
>           diff        lwr        upr   p adj 
> U-B 0.01473809 0.01059542 0.01888076 3.5e-06
> 
> $`exp:cond`
>                       diff          lwr          upr     p adj 
> Exp2:B-Exp1:B  0.016844040  0.006273786  0.027414294 0.0016210
> Exp3:B-Exp1:B -0.001600507 -0.012900584  0.009699570 0.9964757
> Exp1:U-Exp1:B  0.007280763 -0.004019314  0.018580840 0.3277292

  #checking anova assumption validity 
# 2. Homogeneity of variances 
plot(anova, 1)

#Bartlett's test or Levene's test to check the homogeneity of variances. 
library(car)
leveneTest(variable ~ exp*cond, data = test )

# 2. Normality
plot(anova, 2)
#confirm with shapiro 
# Extract the residuals
anova_residuals <- residuals(object = anova)
# Run Shapiro-Wilk test (high P value= normalit is oK)
shapiro.test(x = anova_residuals )
Improve this question
$\endgroup$
4
  • $\begingroup$ what is your sample size? $\endgroup$
    – rep_ho
    Commented Jul 18, 2022 at 14:14
  • $\begingroup$ Given the low sample size, I suspect this might be due to the fact that in Exp1 the difference of means is much lower. You can try to follow this post to compute Tukey's HSD manually, and maybe that will clarify your doubts $\endgroup$
    – Gian
    Commented Jul 18, 2022 at 14:17
  • 1
    $\begingroup$ You have no more than three data points for each condition, sometimes even only two. Don't visualize such a small number of data points using boxplots! As you see, the boxplots are hiding a lot of information. Much better to simply plot the data points themselves - then you would see at a glance that you probably have far too little data and far too much variance to draw any valid conclusions. Boxplots are dangerous. $\endgroup$
    – Stephan Kolassa
    Commented Jul 18, 2022 at 15:08
  • $\begingroup$ Thank you all for you explanations and suggestions. Indeed the boxplot mislead me. $\endgroup$
    – user18309655
    Commented Jul 19, 2022 at 9:06

1 Answer 1

Reset to default
6
$\begingroup$

Update: Interesting discussion (don't skip the comments) about box-and-whiskers for small samples: How should we do boxplots with small samples?

Your plot is misleading rather than revealing.

  • The box and whiskers plot is a 5-number summary: it shows the minimum, lower quartile, median, upper quartile, and maximum. You use a box-plot to summarize as few as 2 values. It's not possible to create a meaningful 5-number summary of 2 numbers. As @StephanKolassa explains, boxplots can be dangerous.
  • Each panel has its own scale for the y-axis. This highlights differences within each experiment but keep in mind that "statistical differences" are measured on the scale of the residual error which is the same for all observations.

Here is a different plot of the same data. If you had made this plot, you might not have written the question.

enter image description here

You wrote the question, however, so let's point a couple more things.

  • Three pairs of condition U measurements are almost exactly the same. (I added jitter to avoid overlap.) Suggests that variance is not the same under the two conditions.
  • It may be easier to make comparisons with the emmeans package.
mod <- lm(variable ~ exp * cond, data = test)
pairs(emmeans(mod, ~ cond | exp))

#> exp = Exp1:
#>  contrast estimate      SE df t.ratio p.value
#>  B - U    -0.00728 0.00341 13  -2.137  0.0522
#> 
#> exp = Exp2:
#>  contrast estimate      SE df t.ratio p.value
#>  B - U    -0.03086 0.00295 13 -10.458  <.0001
#> 
#> exp = Exp3:
#>  contrast estimate      SE df t.ratio p.value
#>  B - U     0.00236 0.00381 13   0.620  0.5462
Improve this answer
$\endgroup$
5
  • $\begingroup$ Hello @dipetkov, thank you for your detailed answer and your suggestions. In the beginning, what made me expect significant result was looking at the mean ± sd of the values : B=0.88a ± 0.14 vs. U=1.60a ± 0.48. But actually you are right choice of plot type was not appropriate and it is clear that I have 2 data points in U that are much closer to the values in B. Thanks again for the suggestions and explanations. $\endgroup$
    – user18309655
    Commented Jul 19, 2022 at 9:18
  • $\begingroup$ How did you compute the standard deviations in your first computation? Was it separately for each condition? $\endgroup$
    – dipetkov
    Commented Jul 19, 2022 at 10:01
  • $\begingroup$ I computed the sd for each condition per experiment separately using the function std in r. $\endgroup$
    – user18309655
    Commented Jul 19, 2022 at 12:02
  • $\begingroup$ This is one more thing to clarify: In the model every observation has the same error variance $\sigma^2$; for each combination of experiment $i$ & condition $j$, the standard error of the sample mean is a function of $\sigma^2$ and the group size $n_{ij}$: $\sigma^2/n_{ij}$. Since your design is almost balanced and there are 10 observations under B and 9 observations under U, the estimated means for the two conditions have similar standard error as well. So in a way the difference between $\text{SE}_B = 0.14$ and $\text{SE}_U = 0.48$ is too big under the equal variance assumption of ANOVA. $\endgroup$
    – dipetkov
    Commented Jul 19, 2022 at 18:21
  • $\begingroup$ Thank you for link for the discussion ! Glad to see that my unintended wrong choice of plot triggered this very interesting discussion! $\endgroup$
    – user18309655
    Commented Jul 21, 2022 at 15:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.