rstatix provides a simple, pipe-friendly framework for common statistical tests, coherent with the tidyverse design philosophy. Each function takes a data frame and a formula, and returns the result as a tidy tibble — so a test flows straight into dplyr pipelines, into reporting tables (gtsummary, gt, broom), and into ggpubr plots, without reshaping.
This vignette walks one analysis from start to finish: describe the data, check the assumptions, run the test, quantify the effect with a confidence interval, and report it. The same handful of verbs covers t-tests, Wilcoxon tests, ANOVA, Kruskal-Wallis, correlation, and their effect sizes.
We use ToothGrowth, a base R dataset: the length of odontoblasts (len) in 60 guinea pigs, each given one of three vitamin C doses (dose) by one of two delivery methods (supp).
df <- ToothGrowth
df$dose <- factor(df$dose)
head(df)
#> len supp dose
#> 1 4.2 VC 0.5
#> 2 11.5 VC 0.5
#> 3 7.3 VC 0.5
#> 4 5.8 VC 0.5
#> 5 6.4 VC 0.5
#> 6 10.0 VC 0.5get_summary_stats() returns the summaries you actually report, for one or many variables, grouped or not.
df %>%
group_by(dose) %>%
get_summary_stats(len, type = "mean_sd")
#> # A tibble: 3 × 5
#> dose variable n mean sd
#> <fct> <fct> <dbl> <dbl> <dbl>
#> 1 0.5 len 20 10.6 4.5
#> 2 1 len 20 19.7 4.42
#> 3 2 len 20 26.1 3.77Before comparing the three doses, check_test_assumptions() tests normality (Shapiro-Wilk, per group) and homogeneity of variance (Levene), and returns a one-row recommendation: which omnibus test and which post-hoc test the data call for.
df %>% check_test_assumptions(len ~ dose)
#> # A tibble: 1 × 8
#> .y. normality.p homogeneity.p normal equal.variance significance omnibus
#> <chr> <dbl> <dbl> <lgl> <lgl> <dbl> <chr>
#> 1 len 0.164 0.528 TRUE TRUE 0.05 anova_test
#> # ℹ 1 more variable: posthoc <chr>Here both assumptions hold, so the recommendation is a one-way ANOVA followed by Tukey HSD. When normality fails it points to Kruskal-Wallis with Dunn’s test; when only the variances differ, to Welch ANOVA with Games-Howell. (Choosing a test by first testing its assumptions on the same data is a known trade-off, noted in ?check_test_assumptions — pre-specifying the test is always defensible.)
Follow the recommendation. anova_test() returns the ANOVA table; asking for the partial eta-squared effect size with a confidence interval adds pes and its bounds.
df %>% anova_test(len ~ dose, effect.size = "pes", ci = 0.95)
#> ANOVA Table (type II tests)
#>
#> Effect DFn DFd F p p<.05 pes conf.low conf.high
#> 1 dose 2 57 67.4 9.53e-16 * 0.703 0.567 0.785Dose has a large, highly significant effect on tooth length.
posthoc_test() runs the post-hoc appropriate to the same decision tree and shows which test it picked and why, so the routing is transparent.
df %>% posthoc_test(len ~ dose)
#> Post-hoc test chosen: Tukey HSD
#> Normality (Shapiro-Wilk, min across groups): p = 0.164 -> normal
#> Homogeneity of variance (Levene): p = 0.528 -> equal variances
#>
#> # A tibble: 3 × 9
#> term group1 group2 null.value estimate conf.low conf.high p.adj
#> * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 dose 0.5 1 0 9.13 5.90 12.4 2.00e- 8
#> 2 dose 0.5 2 0 15.5 12.3 18.7 1.12e-11
#> 3 dose 1 2 0 6.37 3.14 9.59 4.25e- 5
#> # ℹ 1 more variable: p.adj.signif <chr>Every pairwise difference is significant after adjustment. If you prefer to drive the pairwise step yourself, the individual verbs compose the same way — tukey_hsd(), or dunn_test() / games_howell_test() — and any pairwise test can be piped through adjust_pvalue() and add_significance().
A p-value tells you an effect exists; an effect size tells you how large. For a pairwise mean comparison that is Cohen’s d:
df %>% cohens_d(len ~ dose)
#> # A tibble: 3 × 7
#> .y. group1 group2 effsize n1 n2 magnitude
#> * <chr> <chr> <chr> <dbl> <int> <int> <ord>
#> 1 len 0.5 1 -2.05 20 20 large
#> 2 len 0.5 2 -3.73 20 20 large
#> 3 len 1 2 -1.55 20 20 largeEach metric has a matching verb: eta_squared() for ANOVA (on the fitted model), cramer_v() for association, cliff_delta() and wilcox_effsize() for the non-parametric case, kruskal_effsize() for Kruskal-Wallis. Each can also return a confidence interval — ci = TRUE for the formula-based verbs, a confidence level such as ci = 0.95 for eta_squared().
For a manuscript, get_test_label(style = "apa") writes the APA-7 in-text sentence, including the effect-size confidence interval when the result carries one:
model <- df %>% anova_test(len ~ dose, effect.size = "pes", ci = 0.95)
get_test_label(model, type = "text", style = "apa")
#> [1] "F(2, 57) = 67.42, p < .001, eta2[p] = .70, 95% CI [.57, .79]"The test results are rstatix_test objects with tidy() and glance() methods, which hand them as plain tibbles to the reporting-table ecosystem (broom, gtsummary, gt):
df %>% anova_test(len ~ dose) %>% tidy()
#> # A tibble: 1 × 7
#> Effect DFn DFd F p `p<.05` ges
#> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
#> 1 dose 2 57 67.4 9.53e-16 * 0.703
df %>% anova_test(len ~ dose) %>% glance()
#> # A tibble: 1 × 2
#> method n
#> <chr> <int>
#> 1 anova_test 1The workflow is group_by()-aware: run the same test within every level of a grouping variable, then adjust and flag, in a single pipe. Here we compare the two delivery methods within each dose.
df %>%
group_by(dose) %>%
t_test(len ~ supp) %>%
adjust_pvalue(method = "holm") %>%
add_significance()
#> # A tibble: 3 × 11
#> dose .y. group1 group2 n1 n2 statistic df p p.adj
#> <fct> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
#> 1 0.5 len OJ VC 10 10 3.17 15.0 0.00636 0.0127
#> 2 1 len OJ VC 10 10 4.03 15.4 0.00104 0.00312
#> 3 2 len OJ VC 10 10 -0.0461 14.0 0.964 0.964
#> # ℹ 1 more variable: p.adj.signif <chr>The delivery method matters at the two lower doses but not at the highest.
rstatix returns the data; ggpubr draws it. The tidy result of any test pipes straight into a plot: add_xy_position() places the significance brackets and stat_pvalue_manual() adds them, so the p-values on the figure are the ones you computed. See the ggpubr documentation.
Extended tutorials — comparing means, ANOVA designs, correlation analysis, effect sizes — are on the package website: https://rpkgs.datanovia.com/rstatix/.
sessionInfo()
#> R version 4.5.1 (2025-06-13)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sequoia 15.6.1
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
#>
#> locale:
#> [1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8
#>
#> time zone: Europe/Paris
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] dplyr_1.2.1 rstatix_1.1.0
#>
#> loaded via a namespace (and not attached):
#> [1] jsonlite_2.0.0 compiler_4.5.1 tidyselect_1.2.1 tidyr_1.3.2
#> [5] jquerylib_0.1.4 systemfonts_1.3.1 textshaping_1.0.4 yaml_2.3.12
#> [9] fastmap_1.2.0 R6_2.6.1 generics_0.1.4 Formula_1.2-5
#> [13] knitr_1.51 htmlwidgets_1.6.4 backports_1.5.0 tibble_3.3.1
#> [17] car_3.1-3 desc_1.4.3 bslib_0.10.0 pillar_1.11.1
#> [21] rlang_1.2.0 utf8_1.2.6 cachem_1.1.0 broom_1.0.13
#> [25] xfun_0.57 fs_2.0.1 sass_0.4.10 otel_0.2.0
#> [29] cli_3.6.6 withr_3.0.2 pkgdown_2.2.0 magrittr_2.0.4
#> [33] digest_0.6.39 lifecycle_1.0.5 vctrs_0.7.2 evaluate_1.0.5
#> [37] glue_1.8.0 ragg_1.5.0 abind_1.4-8 carData_3.0-5
#> [41] rmarkdown_2.31 purrr_1.2.2 tools_4.5.1 pkgconfig_2.0.3
#> [45] htmltools_0.5.9