I'm puzzled about some of the results that I got, after plotting the data.

I have a dataset where 38 participants are tested on three conditions "time" and are thus tested three times. I figured this is a classic repeated measures ANOVA problem.

Here's a reproducible example of my data:

data_ex <- data.frame( pnum = rep(1:38, times=3),
                   time = rep(c("t1", "t2", "t3"),  each=38),
                   score = c(0.6, 0.8, 1.0, 1.0, 0.3, 0.7, 0.9, 0.3, 0.8,
                             0.3, 0.7, 0.8, 0.8, 1.0, 0.9, 0.7, 0.8, 0.7,
                             1.0, 0.9, 0.4, 0.8, 0.7, 1.0, 1.0, 0.4, 0.9,
                             0.6, 0.3, 0.7, 1.0, 0.9, 0.4, 0.6, 0.3, 0.6,
                             0.8, 0.9, 0.7, 0.9, 0.8, 0.8, 0.6, 1.0, 0.9,
                             0.8, 1.0, 0.4, 0.9, 0.5, 0.6, 0.8, 0.7, 0.6,
                             0.8, 0.7, 0.8, 0.8, 0.4, 0.9, 0.9, 0.9, 1.0,
                             0.6, 0.6, 0.7, 0.3, 0.8, 1.0, 0.9, 0.9, 0.6,
                             0.3, 0.3, 0.9, 0.8, 0.6, 0.8, 1.0, 1.0, 0.6,
                             0.9, 0.9, 0.7, 1.0, 0.4, 0.7, 0.9, 0.4, 0.9,
                             0.9, 0.8, 0.8, 0.9, 0.8, 0.8, 0.7, 1.0, 0.7,
                             1.0, 1.0, 0.9, 1.0, 0.7, 0.5, 0.7, 1.0, 1.0, 
                             0.4, 0.8, 0.4, 0.8, 1.0, 0.8)) 

To perform a repeated measures ANOVA in the tidyverse syntax, I used the anova_test() function from the rstatix package, which apparently is a wrapper around car::Anova().

data_ex_anova  <- anova_test(data = data_ex, dv = score, wid = pnum, within = time)

ANOVA Table (type III tests)

  Effect DFn DFd     F     p p<.05   ges
1   time   2  74 3.593 0.032     * 0.025

Here I find a significant effect. I further check where the difference comes from with the pairwise comparisons

data_ex_anova_pw <- data_ex %>% pairwise_t_test(score ~ time, paired = TRUE, p.adjust.method = "bonferroni")

# A tibble: 3 x 10
  .y.   group1 group2    n1    n2 statistic    df     p p.adj p.adj.signif
* <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl> <dbl> <chr>       
1 score t1      t2         38    38    -0.492    37 0.626 1     ns          
2 score t1      t3         38    38    -2.88     37 0.007 0.02  *           
3 score t2      t3         38    38    -1.94     37 0.06  0.181 ns  

Seems like time 1 and 3 differ.

However, when plotting the results, I get suspicious:

ggplot(data=data_ex,aes(x = as.factor(time), y = score )) +  
  geom_boxplot(size = 1.25, color = "black", outlier.size = NA) + ylim(0,1.1) +
  geom_jitter(size = 3, alpha = 0.5, stroke = 1.5, width = 0.01, height = 0.02) +
  stat_summary(fun.y=mean, geom="point", shape =24, fill = "white", size = 5)

Boxplot of my data

The boxplots are completely overlapping! The medians are basically the same and all the boxes are within each others ranges. Now why this is surprising to me, is when I read about boxplots, this should be an indication of non-significance. Only with large samples, this is not clearly the case. Now since I only have 38 observations per time point, I would not consider this a large sample. Could anyone enlighten me what's going on? Are the results valid or am I indeed performing the wrong test?

PS: I also dabbled with the lmer package, to hopefully account for the repeated participant in all time points, which led to a similar finding:

Anova(lmer(score~time + (1|pnum), data = data_ex))

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: score
      Chisq Df Pr(>Chisq)  
time 6.4703  1    0.01097 *
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

EDIT: I changed the variable "time" to a factor, results are still the same.

  • 1
    $\begingroup$ It looks like you made a programming mistake by entering time as a number rather than a factor. The biggest clue is the unexpected value of DFd in the output. $\endgroup$
    – whuber
    May 11, 2020 at 11:51
  • $\begingroup$ Changing time into a factor did not change the results..... what value would you expect to find for the DFd? $\endgroup$
    – Inkling
    May 11, 2020 at 12:02
  • $\begingroup$ "when I read about boxplots, this should be an indication of non-significance" - can you specify where you read this (sounds strange to me, especially for within-subject tests)? Also, be aware that ANOVA tests for difference in means, not medians, so identical medians should not necessarily be a problem $\endgroup$
    – mesolimbic
    May 11, 2020 at 12:10

3 Answers 3


The boxplots ignore the repeated nature of the data.

If you want a plot of these data (and you should want one!) you can make a plot where the x axis is time, the y axis is score and each participant gets a line. With n = 38 this ought to be readable, but if not, you can separate the data into two parts, either based on some relevant independent variable or on starting value.

In your particular case, it is possible that every single person went up by a small amount between 1 and 3.

Another point: If you know when the measurements happened (e.g. time 1 = day 1, time 2 = day 5, time 3 = day 21 or whatever) it may be better to use that rather than a factor the way you have it.

I also would be wary of RM-ANOVA. It makes strong assumptions (in particular, it assumes sphericity) that are often not reasonable with repeated data. I'd go with either generalized estimating equations (GEE) or a multilevel model.

  • 1
    $\begingroup$ Thank you for the explanation. Would the same be true if the factor time in this case was a variable called condition, in which I have three measurements by the same participant, but not per se split over time (e.g. score on positive, negative and neutral trials of the same task). I think a RM-ANOVA (leaving out the sphericity issue for now) would fit, but plotting the scores with lines between them would seem odd, right? $\endgroup$
    – Inkling
    May 11, 2020 at 14:03
  • $\begingroup$ You would still want a MLM or GEE. The plotting issue would be interesting. I'm not sure what would be best. $\endgroup$
    – Peter Flom
    May 12, 2020 at 10:29

The correct way to make a plot of this sort of data as Peter Flom said would be to make a line for each individual, like this:

ggplot(data_ex, aes(x = time, y = score, color = as.factor(pnum), group = pnum)) + 
  geom_line() +
  facet_wrap(.~pnum) +

which gets this image:


if you take a closer look, you can study how each individual behaves, doing it "fast" and bad, you get something like this, where I draw a red square over the individuals that somehow obtained a bigger score over time (around 19), yellow in those who stayed the same (around 14) and green in the ones that go down (5):

plot with tendency

Even with a lot of doubts about some individuals you can see some tendency.

That visual tendency gets corroborated by a model:

data.table::setorder(data_ex, pnum)
m1 <- geeglm(score ~ time,data = data_ex, id = pnum, corstr = "unstructured")

geeglm(formula = score ~ time, data = data_ex, id = pnum, corstr = "unstructured")

            Estimate Std.err    Wald Pr(>|W|)    
(Intercept)  0.71842 0.03756 365.857  < 2e-16 ***
timet2       0.01579 0.03169   0.248  0.61836    
timet3       0.07632 0.02617   8.506  0.00354 ** 
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation structure = unstructured 
Estimated Scale Parameters:

            Estimate  Std.err
(Intercept)    0.043 0.006599
    Link = identity 

Estimated Correlation Parameters:
          Estimate Std.err
alpha.1:2   0.6464 0.13019
alpha.1:3   0.7306 0.09916
alpha.2:3   0.4571 0.12413
Number of clusters:   38  Maximum cluster size: 3

About the overlapping boxplots and the signification, is it possible that there is some confusion with overlapping boxplots and overlapping confidence intervals of a mean? It is a widespread myth that overlapping mean CI implies no statistical differences (although it is not true)

  • $\begingroup$ In principle this is an important graph to look at. In practice, there is a slight preponderance of increases over decreases, as expected from the means, and so it doesn't seem to add much information over comparison of distributions, $\endgroup$
    – Nick Cox
    May 11, 2020 at 14:00
  • $\begingroup$ Thank you @Avl, as I already commented under Peter Flom's reply: Would you think it's valid to plot lines in the case the factor time was a variable called condition, in which I have three measurements by the same participant, but not per se split over time (e.g. score on positive, negative and neutral trials of the same task). I think a RM-ANOVA would fit, but plotting the scores with lines between them would seem odd, right? $\endgroup$
    – Inkling
    May 11, 2020 at 14:17
  • $\begingroup$ @Inkling Yes, I think even if you change the time for a factor, what would make a different concept of problem, the resolution and the plot without the lines would still be OK. I think you could do RM-ANOVA or gee indistintinctly as the dependent variable is continuous. $\endgroup$
    – Avl
    May 11, 2020 at 16:37
  • $\begingroup$ @NickCox I was aiming at seeing how the individuals evolved over time. If I understand you say that it doesn't add information over the distribution of the scores in each time. I didn't think that the confusion was about that, but about how the repeated measures were behaving. $\endgroup$
    – Avl
    May 11, 2020 at 16:42
  • 1
    $\begingroup$ It's not the same information: the distributions are implicit in the time trajectories (but hard for most people to reduce mentally) but not conversely. It is not inevitable that increasing mean scores from one time to another imply that most individuals increase their scores over time, as a few individuals who improved dramatically might compensate for a larger number who got worse -- yet the distribution graphs also don't support that either, as variability of scores decreases over time. The positive point here is that this kind of graph is an essential check on how individuals change. $\endgroup$
    – Nick Cox
    May 12, 2020 at 18:48

Your plot -- although not actively misleading -- nevertheless doesn't do justice to all the fine structure of your data.

A box plot can't always work well when data are granular with lots of ties: here values are reported only as multiples of 0.1. My bias is that the jittering you use can't be as clear as the stacking of identical values that is possible in a sample this small.

Further, there is a hint in the values, accidentally obscured a little by showing 0.3, 0.6 and 0.9 only as labelled points on the score axis, together with the negative skewness, that scores just possibly have an upper limit at 1. Whether that is so and what that implies about the analysis of variance are open questions.

The display here uses a slightly unorthodox box plot with whiskers just to the extremes. The common convention of drawing whiskers only to points within 1.5 IQR of the nearer quartile I find a little oversold. In any case it can't work especially well with so many ties. But which whisker rule is used is immaterial whenever any display, like that below, shows all the data any way.

With so many ties, the medians can't work well at showing the levels of the distributions. To that extent the box plots are misleading. Naturally showing means as well is a good idea. Here they are diamonds.

This display leaves aside the important detail of which individual is which, explored in another answer.

enter image description here

(Graphic prepared using Stata)

  • $\begingroup$ Thank you for your insights. I was wondering about your comment on the upper limit being 1. I know for a fact this is the case, scores are able to range from 0 to 1 as they represent a False Alarm rate (#false positive answers / #trials), what would be the implicitations on the analysis, knowing this? $\endgroup$
    – Inkling
    May 11, 2020 at 14:10
  • $\begingroup$ The limit of 1 undermines any assumption of normal distributions conditional on predictors, does it not? How much that bites depends what you want to do. P-values are likely to be off, but they are usually wrong any way. $\endgroup$
    – Nick Cox
    May 11, 2020 at 14:13
  • $\begingroup$ If your data are essentially binomial, then a different model may make more sense. $\endgroup$
    – Nick Cox
    May 11, 2020 at 20:34

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.