Get started with rstatix

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.

library(rstatix)
library(dplyr)

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.5

Describe the data

get_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.77

Check the assumptions

Before 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.)

Run the omnibus test

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.785

Dose has a large, highly significant effect on tooth length.

Post-hoc comparisons

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().

Quantify the effect

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 large

Each 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().

Report it

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     1

Grouped analysis in one call

The 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.

Plotting and further reading

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.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] dplyr_1.2.1       rstatix_1.1.0.999 rmarkdown_2.31   
#> 
#> loaded via a namespace (and not attached):
#>  [1] jsonlite_2.0.0   compiler_4.6.1   tidyselect_1.2.1 tidyr_1.3.2     
#>  [5] jquerylib_0.1.4  yaml_2.3.12      fastmap_1.2.0    R6_2.6.1        
#>  [9] generics_0.1.4   Formula_1.2-5    knitr_1.51       backports_1.5.1 
#> [13] tibble_3.3.1     car_3.1-5        maketools_1.3.2  bslib_0.11.0    
#> [17] pillar_1.11.1    rlang_1.3.0      utf8_1.2.6       cachem_1.1.0    
#> [21] broom_1.0.13     xfun_0.60        sass_0.4.10      sys_3.4.3       
#> [25] otel_0.2.0       cli_3.6.6        withr_3.0.3      magrittr_2.0.5  
#> [29] digest_0.6.39    lifecycle_1.0.5  vctrs_0.7.3      evaluate_1.0.5  
#> [33] glue_1.8.1       buildtools_1.0.0 abind_1.4-8      carData_3.0-6   
#> [37] purrr_1.2.2      tools_4.6.1      pkgconfig_2.0.3  htmltools_0.5.9