blog/how-to-do-a-t-test-or-anova-for-many-variables-at-once-in-r-and-communicate-the-results-in-a-better-way/

[utterances-bot]

How to do a t-test or ANOVA for more than one variable at once in R - Stats and R

Learn how to compare samples for multiple variables at once in R thanks to a Student t-test or ANOVA and communicate the results in a better way

https://statsandr.com/blog/how-to-do-a-t-test-or-anova-for-many-variables-at-once-in-r-and-communicate-the-results-in-a-better-way/

[AntoineSoetewey]

Comment written by CARLOS ORTEGA on March 19, 2020 19:59:08:

Thanks for your entry, very elaborated and clear.
For further reference, and I used it very recently, I strongly recommend you to take a look to "compareGroups" package.

It will make the first analysis (not the graphical one) automatically even if the dataset includes categoricals.
https://cloud.r-project.org/web/packages/compareGroups/index.html

Thanks again,
Carlos.

[AntoineSoetewey]

Comment written by Duleep Samuel on March 20, 2020 08:02:54:

thanks for sharing a very useful site

[AntoineSoetewey]

Comment written by Duleep Samuel on March 20, 2020 08:02:54:

thanks for sharing a very useful site

Comment written by Antoine Soetewey on March 20, 2020 08:05:25:

Glad you like it!

[AntoineSoetewey]

Comment written by CARLOS ORTEGA on March 19, 2020 19:59:08:

Thanks for your entry, very elaborated and clear.
For further reference, and I used it very recently, I strongly recommend you to take a look to "compareGroups" package.

It will make the first analysis (not the graphical one) automatically even if the dataset includes categoricals.
https://cloud.r-project.org/web/packages/compareGroups/index.html

Thanks again,
Carlos.

Comment written by Antoine Soetewey on March 20, 2020 09:38:54:

Thanks Carlos for your comment.

I added a reference to the package you mentioned, but I'll definitely discuss it in more details in an article about descriptive statistics and/or about presenting odds ratio.

Regards,
Antoine

[AntoineSoetewey]

Comment written by Chuck Powell on March 20, 2020 12:18:19:

Take a look at ggstatsplot (https://github.com/IndrajeetPatil/ggstatsplot) [full disclosure I contributed a small amount to the package in the past]

[AntoineSoetewey]

Dear Chuck,

I updated the article with a detailed section about the {ggstatsplot} package. See this section. As a contributor to the package, feel free to let me know if my feature request is reasonable or not.

Regards,
Antoine

[AntoineSoetewey]

Comment written by Chuck Powell on March 20, 2020 12:18:19:

Take a look at ggstatsplot (https://github.com/IndrajeetPatil/ggstatsplot) [full disclosure I contributed a small amount to the package in the past]

Comment written by Antoine Soetewey on March 20, 2020 14:58:17:

The package is really complete, thanks for your input Chuck!

I have mentioned it in the article and I'll definitely discuss about it again in a future article.

[AntoineSoetewey]

Comment written by Shodan on April 02, 2020 17:16:00:

Thank you for this I have found it very helpful. However, I have some questions.

Questions

• In the first graph, is the top value shown "1.9e-07" a p-value?
• If so, why isn't it written out as such?
• If son, is it considered significant?
• If so, why aren't asterisks shown?
• The second and third numbers are p-values 
• Are they significant? 
• If so, why aren't asterisks shown? 
• What level of significance do the 4 asterisks shown  represent?
• Could a graph have comparisons with p-values at different levels of significance (e.g., 0.05, 0.01, 0.001 etc.)?

Any help interpreting these graphs would be very helpful.

Thank you. Jason, the #rstatsnewbie

[AntoineSoetewey]

Comment written by Shodan on April 02, 2020 17:16:00:

Thank you for this I have found it very helpful. However, I have some questions.

Questions

• In the first graph, is the top value shown "1.9e-07" a p-value?
• If so, why isn't it written out as such?
• If son, is it considered significant?
• If so, why aren't asterisks shown?
• The second and third numbers are p-values
• Are they significant?
• If so, why aren't asterisks shown?
• What level of significance do the 4 asterisks shown  represent?
• Could a graph have comparisons with p-values at different levels of significance (e.g., 0.05, 0.01, 0.001 etc.)?

Any help interpreting these graphs would be very helpful.

Thank you. Jason, the #rstatsnewbie

Comment written by Antoine Soetewey on April 03, 2020 08:00:14:

Dear Shodan,

• In the first graph, is the top value shown "1.9e-07" a p-value? Yes it is a p-value. 1.9e-07 is the scientific notation for 0.00000019, so the p-value is very small. 
• If so, why isn't it written out as such? Simply because writing the scientific notation is shorter than to write 0.00000019 (so it is more concise, especially if the p-value is even smaller). 
• If so, is it considered significant? The p-value is so small, we reject the null hypothesis (even for a small alpha) so the 2 species versicolor and virginica are significantly different from each other in terms of Sepal.Length. 
• If so, why aren't asterisks shown? The package doesn't allow to write both the p-values AND the asterisks, it is either the p-values or the asterisks. As I find p-values more informative, I chose to display them and not the asterisks. Furthermore, due to multiple testing, the p-values should be adjusted accordingly so the asterisks may not be correct if you have multiple tests. For instance, usually a p-value of 0.04 is significant at 5% level, so an asterisk is usually displayed. However, if you performed 2 t-tests, the significance level is alpha = 0.05 / 2 = 0.025, so your p-value of 0.04 is not significant anymore, so displaying an asterisk is not correct in that case.
• The second and third numbers are p-values? All numeric values in the plots are p-values. 
• Are they significant? It depends on the significance level you choose. Usually, we choose a significance level of 5%. So if p-values are below 0.05, they are significant and we reject the null hypothesis (H0), if p-values >= 0.05, they are not significant and we fail to reject H0. 
• If so, why aren't asterisks shown? Same reasons than above. 
• What level of significance do the 4 asterisks shown represent? I don't see any asterisk in this article. In R, if a p-value is less than 0.05, it is flagged with one star (*), if a p-value is less than 0.01, it is flagged with two stars (**), and if a p-value is less than 0.001, it is flagged with three stars (***), so I assume a p-value with 4 stars means it is smaller than 0.0001.
• Could a graph have comparisons with p-values at different levels of significance (e.g., 0.05, 0.01, 0.001 etc.)? Not sure I understand your question here but I'll try. Of course a p-value shown in a plot allows you to compare it to any significance levels of your choice. It's the point of the p-value, so you can conclude (reject or do not reject H0) at different significance levels.

Hope this helps.

Regards,
Antoine

[AntoineSoetewey]

Comment written by Pat on May 29, 2020 21:49:01:

I have a question on your ANOVA section for this article. I was under the impression that you couldn't use a t-test multiple times without making a correction to your p value, but the values shown in the graphs are non-adjusted is that right? - EDIT - It looks like you have that addressed in the introductory note on 'multiplicity'.  I'm assuming that since your t-test values are so small, they qualify as fitting the Bonferroni correction for multiple t-tests.

[AntoineSoetewey]

Dear Pat,

Thanks for your question.

You’re right the p-values for the multiple comparisons in the ANOVA section are not adjusted.

But as you mention, I usually use the Bonferroni correction by adjusting the alpha. To do so, you can simply divide your desired alpha level by the number of comparisons, and compare it to the p-values obtained.

Hope this helps.

Regards,
Antoine

[Reineany]

Dear Pat, Antoine, my question is connected to the problem of multiple testing in T-tests.
My question is: what alternative analyses you recommend to be done in a scenario with one IV having two levels (gender) on multiple dependent variables? I am so hesitant to perform a simple Anova?
I adopted your code on performing multiple t-tests in one and corrected with bonferroni using your code for p.adjust. Thank you so much for providing this routine.
I did wonder however if I should have performed the Hotelling's T²? In the internet, I could not find a R-implementation and wondered, if you had heard of the Hotelling's T² and have any experience with it?

[AntoineSoetewey]

Hello @Reineany,

If you have multiple dependent variables in function of one IV:

• having two levels, you could use Hotelling's $T^2$
• having more than two levels, you could use a multivariate ANOVA (MANOVA)

Sources: 1 (see page 1 "The multivariate equivalent of the t-test is Hotelling’s $T^2$. [...] MANOVA is the multivariate analogue to Hotelling's $T^2$), and 2.

In the internet, I could not find a R-implementation and wondered, if you had heard of the Hotelling's T² and have any experience with it?

I have heard about it, unfortunately I don't have any good resources about its implementation in R (I have not looked very far though). If you happen to find some, feel free to let me know, I would be interested to learn more about it.

Regards,
Antoine

[Reineany]

As always - clear, concise, incredibly helpful. THANK YOU!

[grahamreidpsychology]

Thank you so much for sharing your code. It is excellent! I have 6 ANOVAS that I want to run on my depedent variable and plot as boxplots. I have written a loop using lapply(), which allows me to run the 6 ANOVAS and plot the boxplots no problem. However, I have no idea how to loop so that I can do the pairwise comparisons and add the significant ones to the boxplots. Would you have any idea how to do this or how to ammend your code (which is similar to mine) to add these comparisons in one loop ?

[AntoineSoetewey]

Thank you so much for sharing your code. It is excellent! I have 6 ANOVAS that I want to run on my depedent variable and plot as boxplots. I have written a loop using lapply(), which allows me to run the 6 ANOVAS and plot the boxplots no problem. However, I have no idea how to loop so that I can do the pairwise comparisons and add the significant ones to the boxplots. Would you have any idea how to do this or how to ammend your code (which is similar to mine) to add these comparisons in one loop ?

Hello,

Hard for me to reply without seeing your code. My code can be used to run ANOVAs and do multiple comparisons, but it will display all comparisons and not just the significant ones.

I believe editing my code to display only the significant comparisons would be too much of a struggle (and not very efficient). However, here are two resources that may be of interest to you:

• Grouped analysis with grouped_ggbetweenstats
• Repeating function execution across multiple columns in a dataframe

Hope this helps.

Regards,
Antoine

[Hanagit-art]

Dear Antoine,
I tried to perform a Kruskal-Wallis test with your code but it did not work. I am totally new in R, so that's maybe that bug :) But I would really appreciate your help to solve this problem. I have cancer patient data. I need to divide continuous variables into 3 groups according to treatment response.

Thank you in advance!
Hana

[AntoineSoetewey]

Dear Hana,

Where are you stuck exactly?

1. You managed to divide your continuous variable into 3 groups but you are unable to do the Kruskal-Wallis test?
2. Or you are stuck at dividing your continuous variable into 3 groups?

If you are stuck at 2., see this section. If you are stuck at 1., perhaps this article about ANOVA will help you (because ANOVA is quite similar to Kruskal-Wallis in R).

If nothing helps, you can send me your data, I'll try to have a look when I find the time.

Regards,
Antoine

[Hanagit-art]

Hi Antoine, thank you a lot for your prompt answer. I will try it and let you know.   Kind regards Hana

[Hanagit-art]

Dear Antoine, I reviewed my code and realized that I forget to set categories. Now everything is working.   Thank you for this blog! It is really helpful! test4$Response <- as.factor(test4$Response) test4$Response = factor(test4$Response, labels = c(A", "B", "C"))

[AntoineSoetewey]

Glad you fixed it Hana!

[zahradm]

Dear Antoine,

Thanks for this useful post.

I have some questions, is it possible for you to help me?

There is a problem in computer science called "graph coloring". Our goal in graph coloring is to color a graph with a number of colors, provided that each connected vertex has a different color.
For solving this problem, we can use logic. So we need to convert graph coloring to CNF form in logic. There are two approaches for encoding, vertex encoding, and edge encoding.
We want to compare the performance of two encoding methods. The main problem is, solvers can affect the performance of encoding. Actually, maybe one solver has better performance with vertex encoding, and or another solver has better performance with edge encoding.

Our goal is to compare two methods of encoding, but I don't know how to delete the solver effect and just compare encoding methods.
In this problem, we have some benchmarks (random graphs) and we have the average solving time of these benchmarks that is obtained from running one solver for 10 iterations on the benchmark.
If we want to compare solvers, we can use t-test, but is it valid to use this method for comparing encoding methods?
And is it possible for statistics to choose the best encoding methods? Or just it compares these and says they are different or not?

Observe that I don't know if the time of running is the best variable to compare, or if I can identify a better metric!

Best regards,
Zahra

[AntoineSoetewey]

Dear Zahra,

Thanks for your question.

I’m not familiar at all with graph coloring.

To make sure I understood well, let me summarize your problem in my words:

• your variable of interest, the solving time, is a quantitative continuous variable
• you have a dependent variable, the encoding method (vertex and edge encoding), which forms the two samples

Now you want to compare the two methods to see whether one is better than the other (in that case, if one method has a significantly smaller time than the other).

If my understanding is correct, yes it’s possible.

Since you have two groups to compare, it would be a Student’s t-test (provided that the assumptions are met, otherwise use the non-parametric version).

Now you have to choose between independent and paired samples:

• If you have two independent samples to compare (because each subject performs only one method), you can use the Student’s t-test for independent samples.
• If your samples are dependent (because each subject performs the two methods, so you measured two times the solving time for each subject), you should use the Student’s t-test for paired samples (to take into account the dependency between the two samples).

See how to perform these two tests in R and by hand in this article. If you want to use the non-parametric version, see this article.

Hope this helps.

Regards,
Antoine

[zahradm]

Hi Antoine,

Many thanks for your response to the question.

Yes, you understood well, and run time has been considered as a quantitative continuous variable.

Student's t-test is a good idea, but as mentioned, different solvers can affect encoding performance.
Is it possible to compare encoding performance without considering solver effects? For instance, I design an experiment that tests different solvers by using different approaches to encoding, and then I remove the effect of the solvers by using "Block Analysis of Variance".

Several researchers have researched comparing solvers, for example in this article, and claim that comparing only the running time is not very reliable for choosing the best solver.
What are your thoughts? How can I find the best features for comparison encoding methods? Comparing solvers and encoding methods differs in terms of problem assumptions, but we can get an idea from solver comparisons.
My problem is comparing encoding methods, not comparing solvers.

And for the last question, how to determine the number of random graphs (samples)? In this problem, I generated 50 graphs. Is it sufficient? Because I don't have to make assumptions about the population of random graphs, and I have chosen the number according to my system's computational limitations!

Kind regards,
Zahra

[AntoineSoetewey]

Again, my knowledge in graph coloring and the issue of solvers is limited so I may not fully grasp your research question, but why don’t you run a multiple linear regression or a two-way ANOVA (if the variable solver is also a categorical variable)?

With the solving time as dependent variable, and the encoding method and the solver as explanatory variables, you would be able to check whether the encoding method is significantly associated with solving time, while controlling for the effect of the solver. You could even add an interaction effect to see if the relationship between solving time and encoding method is dependent on the level of the variable solver.

Regarding the sample size, usually the higher the better of course, but 50 for the types of analysis I mentioned above is already not bad and you could definitely find some good results with such a sample size.

Regards,
Antoine

[zahradm]

Thanks so much.
Your explanations were very useful. I will check them.

Kind regards,
Zahra

[Reineany]

  Dear Antoine, thank you so much - if I understand correctly, using your version of t-tests with bonferroni adjustment would suffice? I have no experience with Hotelling's either which is why I'd prefer to perform your version.   Next I have an embarassing question which I would prefer to ask in private.... :/ If I want to know if the number of people in one group are significantly different from the number of people in another group, then I use the chi-square, right? It basically is, I guess, a frequency of occurences test, not computed on means but on sums but I really --- cannot for the life of me remember how to check for whether group sizes are significantly different :/ My example: female adults with epilepsy, female adults with physical disability, female adults with psychological disability, female adults with chronic pain. Each group contains a number of females but how do I check if there are significant more female adults in one group compared to the other - per school that they attend (4 different schools). :/ I am so embarrassed.....

[AntoineSoetewey]

Dear @Reineany,

First of all, note that your question is not private since it is posted on GitHub and displayed in the comments of the article.
That being said, there is no stupid question so I continue in this thread rather than in private (so everyone can benefit from the discussion).

If I want to know if the number of people in one group are significantly different from the number of people in another group, then I use the chi-square, right?

Correct. Here you have one qualitative variable (the group); see this section for more information about the chi-square goodness of fit test.

Each group contains a number of females but how do I check if there are significant more female adults in one group compared to the other - per school that they attend (4 different schools).

If I understand correctly, now you have two qualitative variables: 1) the school they attend and 2) the group. In this case, you can use the Chi-square test of independence. See more about this test in this article.

Hope this helps.

Regards,
Antoine

[Bellis1511]

At first, thank you for thevery nice explanations. :) I tried the comparison of 15 groups, with the last code here at the side. I get the plots, but there are just in half of the boxplots jitter points. Any idea why? Maybe, it is limited?

[AntoineSoetewey]

Yes it is possible to upload pictures via GitHub: #74. You can simply drag and drop your picture into your message.

Anyway, I have seen your email with your code and plot. Only half of the boxplots have jittered points because by default the color palette accepts only up to 7 groups.

If you want to draw points for more groups, you should use another color palette. For instance, see the last two arguments of this code:

ggbetweenstats(dat,
               x = group,
               y = value,
               centrality.plotting = FALSE,
               plot.type = "box",
               type = "nonparametric",
               package = "ggsci",
               palette = "category20_d3")

This gives:

You'll find many other packages and palettes by running View(paletteer::palettes_d_names).

Hope this helps.

Regards,
Antoine

[Bellis1511]

Oh thanks, that helps a lot :)

But I got another question. Is it also possible to change the code so that the 95% confidence interval is also calculated and displayed in the plot?

[AntoineSoetewey]

At the moment it’s displayed in the title. To display it on the plot, you’ll have to edit the code of the ggbetweenstats() function.

[Bellis1511]

Im so sorry, but I need your help for the last time. I run the code with the new lines from you and got an error massage:

Error in `ensym()`:
! Can't convert to a symbol.

Here is my code, I think I did somethring wrong adjusting y:

plotlist <-
  purrr::pmap(
    .l = list(
      data = list(as_tibble(indiv_char)),
      x = x,
      y = as.list(colnames(indiv_char)[cols]),
      plot.type = "box", # for boxplot
      type = type, # parametric or nonparametric
      pairwise.comparisons = TRUE, # to run post-hoc tests if more than 2 groups
      pairwise.display = "significant", # show only significant differences
      bf.message = FALSE, # remove message about Bayes Factor
      centrality.plotting = FALSE # remove central measure
    ),
    .f = ifelse(paired, # automatically use ggwithinstats if paired samples, ggbetweenstats otherwise
                ggstatsplot::ggwithinstats,
                ggstatsplot::ggbetweenstats(indiv_char,
                                            x = group,
                                            y = as.list(colnames(indiv_char)[cols]),
                                            centrality.plotting = FALSE,
                                            plot.type = "box",
                                            type = "nonparametric",
                                            package = "ggsci",
                                            palette = "category20_d3")))

[AntoineSoetewey]

Try with adding the two arguments package = "ggsci" and palette = "category20_d3" inside the pmap() function instead of inside the ifelse() function, like this:

plotlist <-
  purrr::pmap(
    .l = list(
      data = list(as_tibble(indiv_char)),
      x = x,
      y = as.list(colnames(indiv_char)[cols]),
      plot.type = "box", # for boxplot
      type = type, # parametric or nonparametric
      pairwise.comparisons = TRUE, # to run post-hoc tests if more than 2 groups
      pairwise.display = "significant", # show only significant differences
      bf.message = FALSE, # remove message about Bayes Factor
      centrality.plotting = FALSE, # remove central measure
      package = "ggsci",
      palette = "category20_d3"
    ),
    .f = ifelse(paired, # automatically use ggwithinstats if paired samples, ggbetweenstats otherwise
                ggstatsplot::ggwithinstats,
                ggstatsplot::ggbetweenstats
    )
  )

# print all plots together with statistical results
for (i in 1:length(plotlist)) {
  print(plotlist[[i]])
}

Also make sure that all your dependent variables are quantitative.

[Linophy]

Dear Antoine Soetewey,

I was wondering if you have the answers to the following questions of mine (I guess I should also ask the author of ggstatsplot):

1. Do you know if it is possible / find a way to specify "Dunnett’s test" in ggstatsplot? as I understand it, the default setting in ggstatsplot will do all the possible comparisons which might not be wanted.
2. Do you know if it is possible/ find a way to specify specific comparisons in ggstatsplot? I am thinking of Bonferroni correction-like approach with the advantages of having stats information included in the plot.
   My apologies if the answers to these questions might have been explained in one of your other post. Oh, lastly, this link does not seem to work (at least from my PC): https://antoinesoetewey.shinyapps.io/ggstatsplotShiny/.

Best wishes

[AntoineSoetewey]

Dear @Linophy,

Thanks for your questions.

1. It is not currently possible, unfortunately. See implementation of Dunnett's test IndrajeetPatil/ggstatsplot#499
2. Yes it is possible, with the p.adjust.method argument. For example:

ggbetweenstats(
  data             = movies_long,
  x                = mpaa,
  y                = length,
  p.adjust.method  = "bonferroni")

Regarding the broken link, thanks for spotting it, I removed the link from the post.
I have indeed unpublished this Shiny app because I reached the maximum number of Shiny apps allowed in my current shinyapps.io plan. However, you'll find the code of the app on GitHub: https://github.com/AntoineSoetewey/ggstatsplotShiny. Feel free to copy it if you want to run the app on your PC.

Hope this helps.

Regards,
Antoine

[Linophy]

Cool. Thank you very much ( :