Package {rstatix}

Pipe operator

Description

See magrittr::%>% for details.

Usage

lhs %>% rhs

Manova exported from car package

Description

See car::Manova for details. See the Datanovia tutorial MANOVA in R for a worked walkthrough.

See Also

The Datanovia tutorial: MANOVA in R.

Compact Letter Display of All-Pairwise Comparisons

Description

Adds the compact letter display (CLD) to a pairwise comparison result. Groups that do not share a letter are significantly different. This is a convenient way to annotate plots (e.g. one letter per box/bar) after an all-pairwise post-hoc test such as tukey_hsd(), dunn_test(), games_howell_test(), conover_test(), wilcox_test() or t_test().

The letters are computed with the insert-and-absorb algorithm (Piepho, 2004) using base R only, so no additional package is required (the results match multcompView::multcompLetters()).

Usage

add_cld(test, p.col = NULL, threshold = 0.05, reversed = FALSE, ...)

Arguments

Value

a tibble with one row per group and the following columns: any grouping variables (for a grouped test), .y. (the outcome variable, when present), group (the group level) and cld (the compact letter display). Groups sharing a letter are not significantly different.

References

Piepho, H.-P. (2004) An Algorithm for a Letter-Based Representation of All-Pairwise Comparisons. Journal of Computational and Graphical Statistics, 13(2), 456-466.

See Also

tukey_hsd, dunn_test, games_howell_test, add_significance

Examples

# Tukey HSD post-hoc, then compact letter display
res <- ToothGrowth %>%
  mutate(dose = factor(dose)) %>%
  tukey_hsd(len ~ dose)
res %>% add_cld()

# Works on rank-based post-hocs too
ToothGrowth %>% dunn_test(len ~ dose) %>% add_cld()

# Grouped pairwise test -> one CLD per group
ToothGrowth %>%
  mutate(dose = factor(dose)) %>%
  group_by(supp) %>%
  tukey_hsd(len ~ dose) %>%
  add_cld()

Add P-value Significance Symbols

Description

Add p-value significance symbols into a data frame.

Usage

add_significance(
  data,
  p.col = NULL,
  output.col = NULL,
  cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
  symbols = c("****", "***", "**", "*", "ns")
)

Arguments

Value

a data frame

Examples

# Perform pairwise comparisons and adjust p-values
ToothGrowth %>%
 t_test(len ~ dose) %>%
 adjust_pvalue() %>%
 add_significance("p.adj")

Adjust P-values for Multiple Comparisons

Description

A pipe-friendly function to add an adjusted p-value column into a data frame. Supports grouped data.

Usage

adjust_pvalue(data, p.col = NULL, output.col = NULL, method = "holm")

Arguments

Details

For grouped data (and, equivalently, when a test is run on data grouped with dplyr::group_by() using an in-test p.adjust.method), the p-value adjustment is computed within each group separately, not across all groups. If you instead want a single family of comparisons adjusted across all groups, run the test without adjustment (p.adjust.method = "none") and then call adjust_pvalue() on the combined result (see the grouped example below).

Value

a data frame

Examples

# Perform pairwise comparisons and adjust p-values
ToothGrowth %>%
 t_test(len ~ dose) %>%
 adjust_pvalue()

# Grouped data: adjustment within vs across groups
# Per-group adjustment (within each supp level):
ToothGrowth %>%
  group_by(supp) %>%
  t_test(len ~ dose)                    # in-test holm, adjusted within each group

# One family across ALL comparisons (all groups together):
ToothGrowth %>%
  group_by(supp) %>%
  t_test(len ~ dose, p.adjust.method = "none") %>%
  adjust_pvalue(method = "holm")

Create Nice Summary Tables of ANOVA Results

Description

Create beautiful summary tables of ANOVA test results obtained from either Anova() or aov().

The results include ANOVA table, generalized effect size and some assumption checks.

See the Datanovia tutorial One-Way ANOVA in R for a worked walkthrough.

Usage

anova_summary(
  object,
  effect.size = "ges",
  detailed = FALSE,
  observed = NULL,
  ci = NULL
)

Arguments

Value

return an object of class anova_test a data frame containing the ANOVA table for independent measures ANOVA. However, for repeated/mixed measures ANOVA, it is a list containing the following components are returned:

• ANOVA: a data frame containing ANOVA results

• Mauchly's Test for Sphericity: If any within-Ss variables with more than 2 levels are present, a data frame containing the results of Mauchly's test for Sphericity. Only reported for effects that have more than 2 levels because sphericity necessarily holds for effects with only 2 levels.

• Sphericity Corrections: If any within-Ss variables are present, a data frame containing the Greenhouse-Geisser and Huynh-Feldt epsilon values, and corresponding corrected p-values.

The returned object might have an attribute called args if you compute ANOVA using the function anova_test(). The attribute args is a list holding the arguments used to fit the ANOVA model, including: data, dv, within, between, type, model, etc.

The following abbreviations are used in the different results tables:

• DFn Degrees of Freedom in the numerator (i.e. DF effect).

• DFd Degrees of Freedom in the denominator (i.e., DF error).

• SSn Sum of Squares in the numerator (i.e., SS effect).

• SSd Sum of Squares in the denominator (i.e.,SS error).

• F F-value.

• p p-value (probability of the data given the null hypothesis).

• p<.05 Highlights p-values less than the traditional alpha level of .05.

• ges Generalized Eta-Squared measure of effect size.

• GGe Greenhouse-Geisser epsilon.

• p[GGe] p-value after correction using Greenhouse-Geisser epsilon.

• p[GGe]<.05 Highlights p-values (after correction using Greenhouse-Geisser epsilon) less than the traditional alpha level of .05.

• HFe Huynh-Feldt epsilon.

• p[HFe] p-value after correction using Huynh-Feldt epsilon.

• p[HFe]<.05 Highlights p-values (after correction using Huynh-Feldt epsilon) less than the traditional alpha level of .05.

• W Mauchly's W statistic

Author(s)

Alboukadel Kassambara, alboukadel.kassambara@gmail.com

See Also

anova_test(), factorial_design() The Datanovia tutorial: One-Way ANOVA in R.

Examples

# Load data
#:::::::::::::::::::::::::::::::::::::::
data("ToothGrowth")
df <- ToothGrowth
df$dose <- as.factor(df$dose)

# Independent measures ANOVA
#:::::::::::::::::::::::::::::::::::::::::
# Compute ANOVA and display the summary
res.anova <- Anova(lm(len ~ dose*supp, data = df))
anova_summary(res.anova)

# Display both SSn and SSd using detailed = TRUE
# Show generalized eta squared using effect.size = "ges"
anova_summary(res.anova, detailed = TRUE, effect.size = "ges")

# Show partial eta squared using effect.size = "pes"
anova_summary(res.anova, detailed = TRUE, effect.size = "pes")

# Repeated measures designs using car::Anova()
#:::::::::::::::::::::::::::::::::::::::::
# Prepare the data
df$id <- as.factor(rep(1:10, 6)) # Add individuals ids
head(df)

# Easily perform repeated measures ANOVA using the car package
design <- factorial_design(df, dv = len, wid = id, within = c(supp, dose))
res.anova <- Anova(design$model, idata = design$idata, idesign = design$idesign, type = 3)
anova_summary(res.anova)

# Repeated measures designs using stats::Aov()
#:::::::::::::::::::::::::::::::::::::::::
res.anova <- aov(len ~ dose*supp + Error(id/(supp*dose)), data = df)
anova_summary(res.anova)

Anova Test

Description

Provides a pipe-friendly framework to perform different types of ANOVA tests, including:

• Independent measures ANOVA: between-Subjects designs,

• Repeated measures ANOVA: within-Subjects designs

• Mixed ANOVA: Mixed within within- and between-Subjects designs, also known as split-plot ANOVA and

• ANCOVA: Analysis of Covariance.

The function is an easy to use wrapper around Anova() and aov(). It makes ANOVA computation handy in R and It's highly flexible: can support model and formula as input. Variables can be also specified as character vector using the arguments dv, wid, between, within, covariate.

The results include ANOVA table, generalized effect size and some assumption checks.

Usage

anova_test(
  data,
  formula,
  dv,
  wid,
  between,
  within,
  covariate,
  type = NULL,
  effect.size = "ges",
  error = NULL,
  white.adjust = FALSE,
  observed = NULL,
  detailed = FALSE,
  ci = NULL
)

get_anova_table(x, correction = c("auto", "GG", "HF", "none"))

## S3 method for class 'anova_test'
print(x, ...)

## S3 method for class 'anova_test'
plot(x, ...)

Arguments

Details

Contrasts. By default, R uses treatment contrasts (contr.treatment), where each factor level is compared to the first level used as baseline; the current setting can be checked with options('contrasts').

How anova_test() handles contrasts depends on the interface you use:

• When you use the formula interface (or the dv + between/within arguments), anova_test() fits the model internally with options(contrasts = c('contr.sum', 'contr.poly')), restoring your global option afterwards. This gives orthogonal contrasts, where every level is compared to the overall mean, and type-III results that match the most commonly used commercial softwares, like SPSS.

• When you instead pass a pre-fitted model (lm() or aov()), anova_test() does not change its contrasts: the model keeps whatever contrasts were in effect when it was fitted (R's default contr.treatment unless you set otherwise). Fitting with the default treatment contrasts and then requesting type = 3 can therefore give different results from the formula interface.

To reproduce the formula-interface (SPSS) result from a pre-fitted model, or to obtain the same result with car::Anova() directly, set options(contrasts = c('contr.sum', 'contr.poly')) before fitting the model and use type = 3.

Value

return an object of class anova_test a data frame containing the ANOVA table for independent measures ANOVA.

However, for repeated/mixed measures ANOVA, a list containing the following components are returned: ANOVA table, Mauchly's Test for Sphericity, Sphericity Corrections. These table are described more in the documentation of the function anova_summary().

The returned object has an attribute called args, which is a list holding the arguments used to fit the ANOVA model, including: data, dv, within, between, type, model, etc.

Functions

• anova_test(): perform anova test

• get_anova_table(): extract anova table from an object of class anova_test. When within-subject factors are present, either sphericity corrected or uncorrected degrees of freedom can be reported.

Note

Default sums-of-squares type differs between the two interfaces for unbalanced designs. When type is not supplied:

• the formula interface (anova_test(data, y ~ a*b)) uses type II for between-subjects designs, regardless of balance;

• the dv=/between=/within= interface uses type II for balanced between-subjects designs but switches to type III for unbalanced between-subjects designs with more than one factor (both interfaces use type III for repeated-measures designs).

For a balanced design the SS types coincide, so the two interfaces agree. For an unbalanced factorial design they can give different main-effect F/p values unless you pass type explicitly — which is recommended for reproducibility. Example:

 d <- mtcars %>% dplyr::mutate(cyl = factor(cyl), am = factor(am))
 # differ (unbalanced): formula -> type II, dv/between -> type III
 d %>% anova_test(mpg ~ cyl * am)
 d %>% anova_test(dv = mpg, between = c(cyl, am))
 # agree once type is explicit:
 d %>% anova_test(mpg ~ cyl * am, type = 3)
 d %>% anova_test(dv = mpg, between = c(cyl, am), type = 3)
 

Author(s)

Alboukadel Kassambara, alboukadel.kassambara@gmail.com

See Also

anova_summary(), factorial_design() The Datanovia tutorials: One-Way ANOVA in R, Repeated Measures ANOVA in R, Mixed ANOVA in R, ANCOVA in R.

Examples

# Load data
#:::::::::::::::::::::::::::::::::::::::
data("ToothGrowth")
df <- ToothGrowth

# One-way ANOVA test
#:::::::::::::::::::::::::::::::::::::::::
df %>% anova_test(len ~ dose)

# Grouped One-way ANOVA test
#:::::::::::::::::::::::::::::::::::::::::
df %>%
  group_by(supp) %>%
  anova_test(len ~ dose)

# Two-way ANOVA test
#:::::::::::::::::::::::::::::::::::::::::
df %>% anova_test(len ~ supp*dose)

# Two-way repeated measures ANOVA
#:::::::::::::::::::::::::::::::::::::::::
df$id <- rep(1:10, 6) # Add individuals id
# Use formula

df %>% anova_test(len ~ supp*dose + Error(id/(supp*dose)))



# or use character vector
df %>% anova_test(dv = len, wid = id, within = c(supp, dose))

# Extract ANOVA table and apply correction
#:::::::::::::::::::::::::::::::::::::::::
res.aov <- df %>% anova_test(dv = len, wid = id, within = c(supp, dose))
get_anova_table(res.aov, correction = "GG")


# Use model as arguments
#:::::::::::::::::::::::::::::::::::::::::
.my.model <- lm(yield ~ block + N*P*K, npk)
anova_test(.my.model)

Convert a Correlation Test Data Frame into a Correlation Matrix

Description

Convert a correlation test data frame, returned by the cor_test(), into a correlation matrix format.

See the Datanovia tutorial Correlation Matrix in R for a worked walkthrough.

Usage

as_cor_mat(x)

Arguments

Value

Returns a data frame containing the matrix of the correlation coefficients. The output has an attribute named "pvalue", which contains the matrix of the correlation test p-values.

See Also

cor_mat(), cor_test() The Datanovia tutorial: Correlation Matrix in R.

Examples

# Pairwise correlation tests between variables
#:::::::::::::::::::::::::::::::::::::::::::::::
res.cor.test <- mtcars %>%
  select(mpg, disp, hp, drat, wt, qsec) %>%
  cor_test()
res.cor.test

# Convert the correlation test into a correlation matrix
#:::::::::::::::::::::::::::::::::::::::::::::::
res.cor.test %>% as_cor_mat()

Exact Binomial Test

Description

Performs exact binomial test and pairwise comparisons following a significant exact multinomial test. Wrapper around the R base function link[stats]{binom.test}() that returns a data frame as a result.

Usage

binom_test(
  x,
  n,
  p = 0.5,
  alternative = "two.sided",
  conf.level = 0.95,
  detailed = FALSE
)

pairwise_binom_test(
  x,
  p.adjust.method = "holm",
  alternative = "two.sided",
  conf.level = 0.95
)

pairwise_binom_test_against_p(
  x,
  p = rep(1/length(x), length(x)),
  p.adjust.method = "holm",
  alternative = "two.sided",
  conf.level = 0.95
)

Arguments

Value

return a data frame containing the p-value and its significance. with some the following columns:

• group, group1, group2: the categories or groups being compared.

• statistic: the number of successes.

• parameter: the number of trials.

• p: p-value of the test.

• p.adj: the adjusted p-value.

• method: the used statistical test.

• p.signif, p.adj.signif: the significance level of p-values and adjusted p-values, respectively.

• estimate: the estimated probability of success.

• alternative: a character string describing the alternative hypothesis.

• conf.low,conf.high: Lower and upper bound on a confidence interval for the probability of success.

The returned object has an attribute called args, which is a list holding the test arguments.

Functions

• binom_test(): performs exact binomial test. Wrapper around the R base function binom.test that returns a dataframe as a result.

• pairwise_binom_test(): performs pairwise comparisons (binomial test) following a significant exact multinomial test.

• pairwise_binom_test_against_p(): performs pairwise comparisons (binomial test) following a significant exact multinomial test for given probabilities.

See Also

multinom_test

Examples

# Exact binomial test
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Data: 160 mice with cancer including 95 male and 65 female
# Q1: Does cancer affect more males than females?
binom_test(x = 95, n = 160)
# => yes, there are a significant difference


# Q2: compare the observed proportion of males
# to an expected proportion (p = 3/5)
binom_test(x = 95, n = 160, p = 3/5)
# => there are no significant difference

# Multinomial test
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Data
tulip <- c(red = 81, yellow = 50, white = 27)
# Question 1: are the color equally common ?
# this is a test of homogeneity
res <- multinom_test(tulip)
res
attr(res, "descriptives")

# Pairwise comparisons between groups
pairwise_binom_test(tulip, p.adjust.method = "bonferroni")


# Question 2: comparing observed to expected proportions
# this is a goodness-of-fit test
expected.p <- c(red = 0.5, yellow = 0.33, white = 0.17)
res <- multinom_test(tulip, expected.p)
res
attr(res, "descriptives")

# Pairwise comparisons against a given probabilities
pairwise_binom_test_against_p(tulip, expected.p)

Box's M-test for Homogeneity of Covariance Matrices

Description

Performs the Box's M-test for homogeneity of covariance matrices obtained from multivariate normal data according to one grouping variable. The test is based on the chi-square approximation.

See the Datanovia tutorial Homogeneity of Variance Test in R for a worked walkthrough.

Usage

box_m(data, group)

Arguments

Value

A data frame containing the following components:

See Also

The Datanovia tutorial: Homogeneity of Variance Test in R.

Examples

data(iris)
box_m(iris[, -5], iris[, 5])

Check One-Way Assumptions and Recommend the Test

Description

For a one-way, independent-groups design (outcome ~ group), check the two assumptions that decide which family of tests is appropriate — normality (Shapiro-Wilk, per group) and homogeneity of variance (Levene) — and return the verdicts together with the recommended omnibus and post-hoc test:

• each group normal and equal variances: anova_test() + tukey_hsd();

• each group normal but unequal variances: welch_anova_test() + games_howell_test();

• at least one group not normal: kruskal_test() + dunn_test().

The result is a tidy one-row tibble, so the same single assumption check can drive both the omnibus and the post-hoc coherently — run the recommended omnibus, then pass the result to posthoc_test() via its .assumptions argument to avoid re-checking.

A note on choosing a test from the data. Selecting the test by first testing its assumptions on the same data is convenient but has a known cost: the assumption gate is least reliable exactly when it matters (Shapiro-Wilk has little power at small n and rejects trivial departures at large n), and conditioning the choice on it makes the p-value of the test finally run no longer the exact nominal quantity. A common alternative is to skip the gate and use a robust method unconditionally — Welch ANOVA with Games-Howell (which reduce to the classic result when variances are equal) or a rank-based test. Treat this recommendation as guidance, not a substitute for judgement.

See the Datanovia tutorial Statistical Tests and Assumptions in R for a worked walkthrough.

Usage

check_test_assumptions(data, formula, significance = 0.05)

Arguments

Value

a one-row tibble with the columns .y. (the outcome), normality.p (the smallest Shapiro-Wilk p across groups), homogeneity.p (Levene's p), the logical verdicts normal and equal.variance, the significance used, and the recommended omnibus and posthoc test names.

See Also

posthoc_test(), shapiro_test(), levene_test(). The Datanovia tutorial: Statistical Tests and Assumptions in R.

Examples

df <- ToothGrowth
df$dose <- as.factor(df$dose)
df %>% check_test_assumptions(len ~ dose)

Chi-squared Test for Count Data

Description

Performs chi-squared tests, including goodness-of-fit, homogeneity and independence tests.

chisq_test() also accepts a pipe-friendly data-frame interface for the test of independence between two categorical variables: pass a data frame as x and the two columns either positionally (data %>% chisq_test(var1, var2)) or via vars (data %>% chisq_test(vars = c("var1", "var2"))). The contingency table is built internally. Note that in the positional form the second column occupies the correct argument slot, so use the vars form (or the table interface) if you need to set correct/simulate.p.value.

See the Datanovia tutorial Chi-Square Test of Independence in R for a worked walkthrough.

Usage

chisq_test(
  x,
  y = NULL,
  correct = TRUE,
  p = rep(1/length(x), length(x)),
  rescale.p = FALSE,
  simulate.p.value = FALSE,
  B = 2000,
  vars = NULL
)

pairwise_chisq_gof_test(x, p.adjust.method = "holm", ...)

pairwise_chisq_test_against_p(
  x,
  p = rep(1/length(x), length(x)),
  p.adjust.method = "holm",
  ...
)

chisq_descriptives(res.chisq)

expected_freq(res.chisq)

observed_freq(res.chisq)

pearson_residuals(res.chisq)

std_residuals(res.chisq)

Arguments

Value

return a data frame with some the following columns:

• n: the number of participants.

• group, group1, group2: the categories or groups being compared.

• statistic: the value of Pearson's chi-squared test statistic.

• df: the degrees of freedom of the approximate chi-squared distribution of the test statistic. NA if the p-value is computed by Monte Carlo simulation.

• p: p-value.

• p.adj: the adjusted p-value.

• method: the used statistical test.

• p.signif, p.adj.signif: the significance level of p-values and adjusted p-values, respectively.

• observed: observed counts.

• expected: the expected counts under the null hypothesis.

The returned object has an attribute called args, which is a list holding the test arguments.

Functions

• chisq_test(): performs chi-square tests including goodness-of-fit, homogeneity and independence tests.

• pairwise_chisq_gof_test(): perform pairwise comparisons between groups following a global chi-square goodness of fit test.

• pairwise_chisq_test_against_p(): perform pairwise comparisons after a global chi-squared test for given probabilities. For each group, the observed and the expected proportions are shown. Each group is compared to the sum of all others.

• chisq_descriptives(): returns the descriptive statistics of the chi-square test. These include, observed and expected frequencies, proportions, residuals and standardized residuals. Only available for a single (ungrouped) chisq_test() result.

• expected_freq(): returns the expected counts from the chi-square test result.

• observed_freq(): returns the observed counts from the chi-square test result.

• pearson_residuals(): returns the Pearson residuals, (observed - expected) / sqrt(expected).

• std_residuals(): returns the standardized residuals

See Also

The Datanovia tutorial: Chi-Square Test of Independence in R, Chi-Square Goodness-of-Fit Test in R.

Examples

# Chi-square goodness of fit test
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tulip <- c(red = 81, yellow = 50, white = 27)
# Q1: Are the colors equally common?
chisq_test(tulip)
pairwise_chisq_gof_test(tulip)
# Q2: comparing observed to expected proportions
chisq_test(tulip, p = c(1/2, 1/3, 1/6))
pairwise_chisq_test_against_p(tulip, p = c(0.5, 0.33, 0.17))

# Homogeneity of proportions between groups
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Data: Titanic
xtab <- as.table(rbind(
  c(203, 118, 178, 212),
  c(122, 167, 528, 673)
))
dimnames(xtab) <- list(
  Survived = c("Yes", "No"),
  Class = c("1st", "2nd", "3rd", "Crew")
)
xtab
# Chi-square test
chisq_test(xtab)
# Compare the proportion of survived between groups
pairwise_prop_test(xtab)

# Test of independence using the data-frame interface
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
df <- data.frame(
  gender = rep(c("M", "F"), each = 100),
  smoker = rep(c("yes", "no", "yes", "no"), times = c(30, 70, 60, 40))
)
# Positional columns
df %>% chisq_test(gender, smoker)
# Equivalent, using vars (keeps `correct` settable)
df %>% chisq_test(vars = c("gender", "smoker"), correct = FALSE)

Cliff's Delta Effect Size for Ordinal / Non-parametric Comparisons

Description

Compute Cliff's delta, a non-parametric effect size for the difference between two groups. It is the standardized version of the Mann-Whitney statistic and estimates the probability that a randomly drawn value from one group exceeds a randomly drawn value from the other, minus the reverse probability: \delta = (\#\{x > y\} - \#\{x < y\}) / (n_1 n_2). It ranges from -1 to 1 and, unlike the rank-biserial r, is unaffected by ties beyond their contribution to the counts.

See the Datanovia tutorial Wilcoxon Test in R for a worked walkthrough.

Usage

cliff_delta(
  data,
  formula,
  comparisons = NULL,
  ref.group = NULL,
  ci = FALSE,
  conf.level = 0.95,
  ci.type = "perc",
  nboot = 1000,
  ...,
  boot.parallel = getOption("boot.parallel", "no"),
  boot.ncpus = getOption("boot.ncpus", 1L)
)

Arguments

Details

The magnitude thresholds are those of Romano et al. (2006): |delta| < 0.147 "negligible", < 0.33 "small", < 0.474 "medium", otherwise "large". Cliff's delta is algebraically identical to the rank-biserial correlation, so the point estimate equals effectsize::rank_biserial().

Value

a tibble with one row per comparison and the columns .y., group1, group2, effsize (Cliff's delta), n1, n2 and magnitude; conf.low / conf.high are added when ci = TRUE.

References

Cliff, N. (1993). Dominance statistics: Ordinal analyses to answer ordinal questions. Psychological Bulletin, 114(3), 494-509.

Romano, J., Kromrey, J. D., Coraggio, J., & Skowronek, J. (2006). Appropriate statistics for ordinal level data. Annual meeting of the Florida Association of Institutional Research.

See Also

The Datanovia tutorial: Wilcoxon Test in R.

Examples

# Two-samples Cliff's delta
ToothGrowth %>% cliff_delta(len ~ supp)

# Pairwise comparisons
ToothGrowth %>% cliff_delta(len ~ dose)

# Grouped data
ToothGrowth %>%
  dplyr::group_by(supp) %>%
  cliff_delta(len ~ dose)

Cochran's Q Test

Description

Performs the Cochran's Q test for unreplicated randomized block design experiments with a binary response variable and paired data. This test is analogue to the friedman.test() with 0,1 coded response. It's an extension of the McNemar Chi-squared test for comparing more than two paired proportions.

See the Datanovia tutorial Cochran’s Q Test in R for a worked walkthrough.

Usage

cochran_qtest(data, formula)

Arguments

See Also

The Datanovia tutorial: Cochran’s Q Test in R.

Examples

# Generate a demo data
mydata <- data.frame(
  outcome = c(0,1,1,0,0,1,0,1,1,1,1,1,0,0,1,1,0,1,0,1,1,0,0,1,0,1,1,0,0,1),
  treatment = gl(3,1,30,labels=LETTERS[1:3]),
  participant = gl(10,3,labels=letters[1:10])
)
mydata$outcome <- factor(
  mydata$outcome, levels = c(1, 0),
  labels = c("success", "failure")
  )
# Cross-tabulation
xtabs(~outcome + treatment, mydata)

# Compare the proportion of success between treatments
cochran_qtest(mydata, outcome ~ treatment|participant)

# pairwise comparisons between groups
pairwise_mcnemar_test(mydata, outcome ~ treatment|participant)

Compute Cohen's d Measure of Effect Size

Description

Compute the effect size for t-test. T-test conventional effect sizes, proposed by Cohen, are: 0.2 (small effect), 0.5 (moderate effect) and 0.8 (large effect).

Cohen's d is calculated as the difference between means or mean minus mu divided by the estimated standardized deviation.

For independent samples t-test, there are two possibilities implemented. If the t-test did not make a homogeneity of variance assumption, (the Welch test), the variance term will mirror the Welch test, otherwise a pooled estimate is used.

If a paired samples t-test was requested, then effect size desired is based on the standard deviation of the differences.

It can also return confidence intervals for the effect size, either by bootstrap (the default) or by an analytic, deterministic method (see ci.method).

See the Datanovia tutorial Cohen’s d Effect Size in R for a worked walkthrough.

Usage

cohens_d(
  data,
  formula,
  comparisons = NULL,
  ref.group = NULL,
  paired = FALSE,
  mu = 0,
  var.equal = FALSE,
  hedges.correction = FALSE,
  ci = FALSE,
  conf.level = 0.95,
  ci.type = "perc",
  nboot = 1000,
  boot.parallel = getOption("boot.parallel", "no"),
  boot.ncpus = getOption("boot.ncpus", 1L),
  id = NULL,
  ci.method = c("boot", "analytic")
)

Arguments

Details

Quantification of the effect size magnitude is performed using the thresholds defined in Cohen (1992). The magnitude is assessed using the thresholds provided in (Cohen 1992), i.e. |d| < 0.2 "negligible", |d| < 0.5 "small", |d| < 0.8 "medium", otherwise "large".

Value

return a data frame with some of the following columns:

• .y.: the y variable used in the test.

• group1,group2: the compared groups in the pairwise tests.

• n,n1,n2: Sample counts.

• effsize: estimate of the effect size (d value).

• magnitude: magnitude of effect size.

• conf.low,conf.high: lower and upper bound of the effect size confidence interval.

References

• Cohen, J. (1988). Statistical power analysis for the behavioral sciences (2nd ed.). New York:Academic Press.

• Cohen, J. (1992). A power primer. Psychological Bulletin, 112, 155-159.

• Hedges, Larry & Olkin, Ingram. (1985). Statistical Methods in Meta-Analysis. 10.2307/1164953.

• Navarro, Daniel. 2015. Learning Statistics with R: A Tutorial for Psychology Students and Other Beginners (Version 0.5).

• Steiger, J. H. (2004). Beyond the F test: Effect size confidence intervals and tests of close fit in the analysis of variance and contrast analysis. Psychological Methods, 9(2), 164-182.

• Cumming, G., & Finch, S. (2001). A primer on the understanding, use, and calculation of confidence intervals that are based on central and noncentral distributions. Educational and Psychological Measurement, 61(4), 532-574.

See Also

The Datanovia tutorial: Cohen’s d Effect Size in R.

Examples

# One-sample t test effect size
ToothGrowth %>% cohens_d(len ~ 1, mu = 0)

# Two indepedent samples t-test effect size
ToothGrowth %>% cohens_d(len ~ supp, var.equal = TRUE)

# Paired samples effect size
df <- data.frame(
  id = 1:5,
  pre  = c(110, 122, 101, 120, 140),
  post = c(150, 160, 110, 140, 155)
)
df <- df %>% gather(key = "treatment", value = "value", -id)
head(df)

df %>% cohens_d(value ~ treatment, paired = TRUE)

Conover's All-Pairs Rank Comparison Test

Description

Performs Conover's test (also known as the Conover-Iman test) for pairwise multiple comparisons of the ranked data, following a significant Kruskal-Wallis test. It is closely related to dunn_test(), but uses the pooled within-group rank variance and refers the test statistic to a t-distribution (with N - k degrees of freedom) instead of the standard normal distribution. The Conover-Iman test is generally more powerful than Dunn's test, but should only be used as a post-hoc procedure when the Kruskal-Wallis test is itself significant (Conover, 1999).

If a reference group is specified (via ref.group), then each of the remaining group levels is compared only to the reference (control) group, and the p-value adjustment for multiple comparisons is computed over only these k - 1 comparisons (instead of all k(k - 1)/2 pairwise comparisons), exactly as for dunn_test().

See the Datanovia tutorial Kruskal-Wallis Test in R for a worked walkthrough.

Usage

conover_test(
  data,
  formula,
  p.adjust.method = "holm",
  ref.group = NULL,
  detailed = FALSE
)

Arguments

Details

The Conover-Iman pairwise statistic for comparing groups i and j is

t_{ij} = \frac{\bar{R}_i - \bar{R}_j}{\sqrt{S^2 \, \frac{N - 1 - H}{N - k} \left(\frac{1}{n_i} + \frac{1}{n_j}\right)}}

where \bar{R} are the mean ranks, H is the (tie-corrected) Kruskal-Wallis statistic, N is the total sample size, k is the number of groups, and S^2 is the variance of the ranks (S^2 = N(N+1)/12 when there are no ties; otherwise S^2 = \frac{1}{N - 1}\left(\sum r^2 - \frac{N(N+1)^2}{4}\right)). The statistic is referred to a t-distribution with N - k degrees of freedom.

In the returned table each row is oriented with i = group2 and j = group1: estimate is \bar{R}_{group2} - \bar{R}_{group1} and statistic carries its sign, the same convention as dunn_test().

The results match PMCMRplus::kwAllPairsConoverTest().

Value

return a data frame with some of the following columns:

• .y.: the y (outcome) variable used in the test.

• group1,group2: the compared groups in the pairwise tests.

• n1,n2: Sample counts.

• estimate: mean ranks difference.

• estimate1, estimate2: show the mean rank values of the two groups, respectively.

• statistic: Test statistic (t-value) used to compute the p-value.

• df: degrees of freedom (N - k, the same for every comparison).

• p: p-value.

• p.adj: the adjusted p-value.

• method: the statistical test used to compare groups.

• p.adj.signif: the significance level of the adjusted p-values.

The returned object has an attribute called args, which is a list holding the test arguments.

References

Conover, W. J. (1999) Practical Nonparametric Statistics, 3rd edition. Wiley.

Conover, W. J. and Iman, R. L. (1979) On multiple-comparisons procedures. Technical Report LA-7677-MS, Los Alamos Scientific Laboratory.

See Also

dunn_test, kruskal_test The Datanovia tutorial: Kruskal-Wallis Test in R.

Examples

# Simple test
ToothGrowth %>% conover_test(len ~ dose)

# Comparison against a reference (control) group
# each group is compared to the reference; the p-value
# adjustment corrects for only these k - 1 comparisons
ToothGrowth %>% conover_test(len ~ dose, ref.group = "0.5")

# Grouped data
ToothGrowth %>%
  group_by(supp) %>%
  conover_test(len ~ dose)

Factors

Description

Provides pipe-friendly functions to convert simultaneously multiple variables into a factor variable.

Helper functions are also available to set the reference level and the levels order.

Usage

convert_as_factor(data, ..., vars = NULL, make.valid.levels = FALSE)

set_ref_level(data, name, ref)

reorder_levels(data, name, order)

Arguments

Functions

• convert_as_factor(): Convert one or multiple variables into factor.

• set_ref_level(): Change a factor reference level or group.

• reorder_levels(): Change the order of a factor levels

Examples

# Create a demo data
df <- tibble(
  group = c("a", "a", "b", "b", "c", "c"),
  time = c("t1", "t2", "t1", "t2", "t1", "t2"),
  value = c(5, 6, 1, 3, 4, 5)
)
df
# Convert group and time into factor variable
result <- df %>% convert_as_factor(group, time)
result
# Show group levels
levels(result$group)

# Set c as the reference level (the first one)
result <- result %>%
  set_ref_level("group", ref = "c")
levels(result$group)

# Set the order of levels
result <- result %>%
  reorder_levels("group", order = c("b", "c", "a"))
levels(result$group)

Replace Correlation Coefficients by Symbols

Description

Take a correlation matrix and replace the correlation coefficients by symbols according to the level of the correlation.

See the Datanovia tutorial Correlation Matrix in R for a worked walkthrough.

Usage

cor_as_symbols(
  x,
  cutpoints = c(0, 0.25, 0.5, 0.75, 1),
  symbols = c(" ", ".", "+", "*")
)

Arguments

See Also

cor_mat() The Datanovia tutorial: Correlation Matrix in R.

Examples

# Compute correlation matrix
#::::::::::::::::::::::::::::::::::::::::::
cor.mat <- mtcars %>%
  select(mpg, disp, hp, drat, wt, qsec) %>%
  cor_mat()

# Replace correlation coefficient by symbols
#::::::::::::::::::::::::::::::::::::::::::
cor.mat %>%
  cor_as_symbols() %>%
  pull_lower_triangle()

Reshape Correlation Data

Description

Reshape correlation analysis results. Key functions:

• cor_gather(): takes a correlation matrix and collapses (i.e. melt) it into a paired list (long format).

• cor_spread(): spread a long correlation data format across multiple columns. Particularly, it takes the results of cor_test and transforms it into a correlation matrix.

See the Datanovia tutorial Correlation Matrix in R for a worked walkthrough.

Usage

cor_gather(data, drop.na = TRUE)

cor_spread(data, value = "cor")

Arguments

Functions

• cor_gather(): takes a correlation matrix and collapses (or melt) it into long format data frame (paired list)

• cor_spread(): spread a long correlation data frame into wide format. Expects the columns "var1", "var2" and "cor" in the data. (correlation matrix).

See Also

cor_mat(), cor_reorder() The Datanovia tutorial: Correlation Matrix in R.

Examples

# Data preparation
#::::::::::::::::::::::::::::::::::::::::::
mydata <- mtcars %>%
  select(mpg, disp, hp, drat, wt, qsec)
head(mydata, 3)

# Reshape a correlation matrix
#::::::::::::::::::::::::::::::::::::::::::
# Compute a correlation matrix
cor.mat <- mydata %>% cor_mat()
cor.mat

# Collapse the correlation matrix into long format
# paired list data frame
long.format <- cor.mat %>% cor_gather()
long.format

# Spread a correlation data format
#::::::::::::::::::::::::::::::::::::::::::
# Spread the correlation coefficient value
long.format %>% cor_spread(value = "cor")
# Spread the p-value
long.format %>% cor_spread(value = "p")

Add Significance Levels To a Correlation Matrix

Description

Combines correlation coefficients and significance levels in a correlation matrix data.

See the Datanovia tutorial Correlation Matrix in R for a worked walkthrough.

Usage

cor_mark_significant(
  x,
  cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
  symbols = c("****", "***", "**", "*", "")
)

Arguments

Value

a data frame containing the lower triangular part of the correlation matrix marked by significance symbols.

See Also

The Datanovia tutorial: Correlation Matrix in R.

Examples

mtcars %>%
  select(mpg, disp, hp, drat, wt, qsec) %>%
  cor_mat() %>%
  cor_mark_significant()

Compute Correlation Matrix with P-values

Description

Compute correlation matrix with p-values. Numeric columns in the data are detected and automatically selected for the analysis. You can also specify variables of interest to be used in the correlation analysis.

See the Datanovia tutorial Correlation Matrix in R for a worked walkthrough.

Usage

cor_mat(
  data,
  ...,
  vars = NULL,
  method = "pearson",
  alternative = "two.sided",
  conf.level = 0.95
)

cor_pmat(
  data,
  ...,
  vars = NULL,
  method = "pearson",
  alternative = "two.sided",
  conf.level = 0.95
)

cor_get_pval(x)

Arguments

Value

a data frame

Functions

• cor_mat(): compute correlation matrix with p-values. Returns a data frame containing the matrix of the correlation coefficients. The output has an attribute named "pvalue", which contains the matrix of the correlation test p-values.

• cor_pmat(): compute the correlation matrix but returns only the p-values of the tests.

• cor_get_pval(): extract a correlation matrix p-values from an object of class cor_mat(). P-values are not adjusted.

See Also

cor_test(), cor_reorder(), cor_gather(), cor_select(), cor_as_symbols(), pull_triangle(), replace_triangle() The Datanovia tutorial: Correlation Matrix in R.

Examples

# Data preparation
#:::::::::::::::::::::::::::::::::::::::::::
mydata <- mtcars %>%
  select(mpg, disp, hp, drat, wt, qsec)
head(mydata, 3)

# Compute correlation matrix
#::::::::::::::::::::::::::::::::::::::::::
# Correlation matrix between all variables
cor.mat <- mydata %>% cor_mat()
cor.mat

# Specify some variables of interest
mydata %>% cor_mat(mpg, hp, wt)

# Or remove some variables in the data
# before the analysis
mydata %>% cor_mat(-mpg, -hp)

# Significance levels
#::::::::::::::::::::::::::::::::::::::::::
cor.mat %>% cor_get_pval()


# Visualize
#::::::::::::::::::::::::::::::::::::::::::
# Insignificant correlations are marked by crosses
cor.mat %>%
  cor_reorder() %>%
  pull_lower_triangle() %>%
  cor_plot(label = TRUE)

# Gather/collapse correlation matrix into long format
#::::::::::::::::::::::::::::::::::::::::::
cor.mat %>% cor_gather()

Visualize Correlation Matrix Using Base Plot

Description

Provide a tibble-friendly framework to visualize a correlation matrix. Wrapper around the R base function corrplot(). Compared to corrplot(), it can handle directly the output of the functions cor_mat() (in rstatix), rcorr() (in Hmisc), correlate() (in corrr) and cor() (in stats).

The p-values contained in the outputs of the functions cor_mat() and rcorr() are automatically detected and used in the visualization.

See the Datanovia tutorial Correlation Matrix in R for a worked walkthrough.

Usage

cor_plot(
  cor.mat,
  method = "circle",
  type = "full",
  palette = NULL,
  p.mat = NULL,
  significant.level = 0.05,
  insignificant = c("cross", "blank"),
  label = FALSE,
  font.label = list(),
  ...
)

Arguments

See Also

cor_as_symbols() The Datanovia tutorial: Correlation Matrix in R.

Examples

# Compute correlation matrix
#::::::::::::::::::::::::::::::::::::::::::
cor.mat <- mtcars %>%
  select(mpg, disp, hp, drat, wt, qsec) %>%
  cor_mat()

# Visualize correlation matrix
#::::::::::::::::::::::::::::::::::::::::::
# Full correlation matrix,
# insignificant correlations are marked by crosses
cor.mat %>% cor_plot()

# Reorder by correlation coefficient
# pull lower triangle and visualize
cor.lower.tri <- cor.mat %>%
  cor_reorder() %>%
  pull_lower_triangle()
cor.lower.tri %>% cor_plot()

# Change visualization methods
#::::::::::::::::::::::::::::::::::::::::::
cor.lower.tri %>%
  cor_plot(method = "pie")

cor.lower.tri %>%
  cor_plot(method = "color")

cor.lower.tri %>%
  cor_plot(method = "number")

# Show the correlation coefficient: label = TRUE
# Blank the insignificant correlation
#::::::::::::::::::::::::::::::::::::::::::
cor.lower.tri %>%
  cor_plot(
    method = "color",
    label = TRUE,
    insignificant = "blank"
  )

# Change the color palettes
#::::::::::::::::::::::::::::::::::::::::::

# Using custom color palette
# Require ggpubr: install.packages("ggpubr")
if(require("ggpubr")){
  my.palette <- get_palette(c("red", "white", "blue"), 200)
  cor.lower.tri %>%
   cor_plot(palette = my.palette)
}

# Using RcolorBrewer color palette
if(require("ggpubr")){
  my.palette <- get_palette("PuOr", 200)
  cor.lower.tri %>%
   cor_plot(palette = my.palette)
}

Reorder Correlation Matrix

Description

reorder correlation matrix, according to the coefficients, using the hierarchical clustering method.

See the Datanovia tutorial Correlation Matrix in R for a worked walkthrough.

Usage

cor_reorder(x)

Arguments

Value

a data frame

See Also

cor_mat(), cor_gather(), cor_spread() The Datanovia tutorial: Correlation Matrix in R.

Examples

# Compute correlation matrix
#::::::::::::::::::::::::::::::::::::::::::
cor.mat <- mtcars %>%
  select(mpg, disp, hp, drat, wt, qsec) %>%
  cor_mat()

# Reorder by correlation and get p-values
#::::::::::::::::::::::::::::::::::::::::::
# Reorder
cor.mat %>%
  cor_reorder()
# Get p-values of the reordered cor_mat
cor.mat %>%
  cor_reorder() %>%
  cor_get_pval()

Subset Correlation Matrix

Description

See the Datanovia tutorial Correlation Matrix in R for a worked walkthrough.

Usage

cor_select(x, ..., vars = NULL)

Arguments

Value

a data frame

See Also

cor_mat(), pull_triangle(), replace_triangle(). The Datanovia tutorial: Correlation Matrix in R.

Examples

# Compute correlation matrix
#::::::::::::::::::::::::::::::::::::::::::
cor.mat <- mtcars %>%
  select(mpg, disp, hp, drat, wt, qsec) %>%
  cor_mat()

# Subsetting correlation matrix
#::::::::::::::::::::::::::::::::::::::::::

# Select some variables of interest
cor.mat %>%
  cor_select(mpg, drat, wt)

# Remove variables
cor.mat %>%
  cor_select(-mpg, -wt)

Correlation Test

Description

Provides a pipe-friendly framework to perform correlation test between paired samples, using Pearson, Kendall or Spearman method. Wrapper around the function cor.test().

Can also performs multiple pairwise correlation analyses between more than two variables or between two different vectors of variables. Using this function, you can also compute, for example, the correlation between one variable vs many.

See the Datanovia tutorial Correlation Test in R for a worked walkthrough.

Usage

cor_test(
  data,
  ...,
  vars = NULL,
  vars2 = NULL,
  alternative = "two.sided",
  method = "pearson",
  conf.level = 0.95,
  use = "pairwise.complete.obs"
)

Arguments

Value

return a data frame with the following columns:

• var1, var2: the variables used in the correlation test.

• cor: the correlation coefficient.

• statistic: Test statistic used to compute the p-value.

• df: the degrees of freedom (Pearson method only).

• p: p-value.

• conf.low,conf.high: Lower and upper bounds on a confidence interval.

• method: the method used to compute the statistic.

Functions

• cor_test(): correlation test between two or more variables.

Note

cor_test() does not support weighted correlations: it wraps cor.test(), which has no weights argument. Passing weights = therefore raises an error rather than silently returning the unweighted result. For a weighted Pearson correlation use base R, e.g. stats::cov.wt(data[, c("x", "y")], wt = data$w, cor = TRUE)$cor.

See Also

cor_mat(), as_cor_mat(), rstatix-programming (using variable names held in strings) The Datanovia tutorial: Correlation Test in R.

Examples

# Correlation between the specified variable vs
# the remaining numeric variables in the data
#:::::::::::::::::::::::::::::::::::::::::
mtcars %>% cor_test(mpg)

# Correlation test between two variables
#:::::::::::::::::::::::::::::::::::::::::
mtcars %>% cor_test(wt, mpg)

# Pairwise correlation between multiple variables
#:::::::::::::::::::::::::::::::::::::::::
mtcars %>% cor_test(wt, mpg, disp)

# Grouped data
#:::::::::::::::::::::::::::::::::::::::::
iris %>%
  group_by(Species) %>%
 cor_test(Sepal.Width, Sepal.Length)

# Multiple correlation test
#:::::::::::::::::::::::::::::::::::::::::
# Correlation between one variable vs many
mtcars %>% cor_test(
  vars = "mpg",
  vars2 = c("disp", "hp", "drat")
 )

# Correlation between two vectors of variables
# Each element in vars is tested against all elements in vars2
mtcars %>% cor_test(
  vars = c("mpg", "wt"),
  vars2 = c("disp", "hp", "drat")
 )

Convert a Table of Counts into a Data Frame of cases

Description

converts a contingency table or a data frame of counts into a data frame of individual observations.

Usage

counts_to_cases(x, count.col = "Freq")

Arguments

Value

a data frame of cases

Examples

# Create a cross-tabulation demo data
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xtab <- as.table(
  rbind(c(20, 5), c(16,9))
)
dimnames(xtab) <- list(
  before = c("non.smoker", "smoker"),
  after = c("non.smoker", "smoker")
)
xtab

# Convert into a data frame of cases
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
df <- counts_to_cases(xtab)
head(df)

Compute Cramer's V

Description

Compute Cramer's V, which measures the strength of the association between categorical variables.

See the Datanovia tutorial Chi-Square Test of Independence in R for a worked walkthrough.

Usage

cramer_v(x, y = NULL, correct = FALSE, ..., ci = FALSE, conf.level = 0.95)

Arguments

Details

Cramer's V is V = \sqrt{\chi^2 / (N (k - 1))} (Cramer, 1946), where \chi^2 is the Pearson chi-square statistic, N the total count and k the smaller of the two table dimensions.

The confidence interval is obtained by inverting the noncentral chi-square distribution (Smithson, 2003; Steiger, 2004): the noncentrality parameters \lambda whose distributions place the observed chi-square statistic at the 1 - \alpha/2 and \alpha/2 quantiles are found by root finding, and each is converted with V = \sqrt{\lambda / (N (k - 1))}. The interval is computed from the same chi-square statistic as the point estimate, so correct = TRUE shifts both.

The bounds are clipped to [0, 1], the range of Cramer's V. They are NA, with a warning, when the statistic or its degrees of freedom are undefined – for instance when the table has an empty row or column, or when simulate.p.value = TRUE is passed on to chisq.test(), which then reports no degrees of freedom.

The interval usually brackets the reported effsize, but it does not in the near-independence corner: when the observed chi-square falls below the \alpha/2 quantile of its central distribution, no noncentrality is consistent with the data at that quantile, both bounds collapse to [0, 0], and the (necessarily positive) point estimate lies above them. This is a property of the noncentral inversion rather than of this implementation, and it only arises for effect sizes indistinguishable from zero.

At the default correct = FALSE, the results match effectsize::cramers_v(adjust = FALSE, ci = , alternative = "two.sided") away from the near-independence corner above (where numerical inversions differ), and DescTools::CramerV(conf.level = ) to about four decimals (its inversion uses a looser tolerance), except at a chi-square of exactly zero, where DescTools returns NA bounds and this function returns the collapsed [0, 0] interval. Neither package applies Yates' continuity correction, so correct = TRUE values have no counterpart there (DescTools's own correct argument selects the Bergsma bias correction, a different adjustment).

Value

By default, a single numeric value: Cramer's V.

When ci = TRUE, a one-row tibble with the columns effsize (Cramer's V), conf.low and conf.high – the same confidence interval columns that anova_test(ci = ) returns.

References

Cramer, H. (1946). Mathematical Methods of Statistics. Princeton University Press.

Smithson, M. (2003). Confidence Intervals. Sage Publications.

Steiger, J. H. (2004). Beyond the F test: Effect size confidence intervals and tests of close fit in the analysis of variance and contrast analysis. Psychological Methods, 9, 164-182.

See Also

The Datanovia tutorial: Chi-Square Test of Independence in R.

Examples

# Data preparation
df <- as.table(rbind(c(762, 327, 468), c(484, 239, 477)))
dimnames(df) <- list(
  gender = c("F", "M"),
  party = c("Democrat","Independent", "Republican")
)
df
# Compute cramer's V
cramer_v(df)

# Add a confidence interval
cramer_v(df, ci = TRUE)

# Yates' continuity correction only affects 2x2 tables. It belongs to the
# test, not to the effect size, so it is off by default.
tab <- as.table(rbind(c(20, 30), c(35, 15)))
cramer_v(tab)
cramer_v(tab, correct = TRUE)

Arrange Rows by Column Values

Description

Order the rows of a data frame by values of specified columns. Wrapper arround the arrange() function. Supports standard and non standard evaluation.

Usage

df_arrange(data, ..., vars = NULL, .by_group = FALSE)

Arguments

Value

a data frame

Examples

df <- head(ToothGrowth)
df

# Select column using standard evaluation
df %>% df_arrange(vars = c("dose", "len"))

# Select column using non-standard evaluation
df %>% df_arrange(dose, desc(len))

Get User Specified Variable Names

Description

Returns user specified variable names. Supports standard and non standard evaluation.

Usage

df_get_var_names(data, ..., vars = NULL)

Arguments

Value

a character vector

Examples

# Non standard evaluation
ToothGrowth %>%
 df_get_var_names(dose, len)

# Standard evaluation
ToothGrowth %>%
 df_get_var_names(vars = c("len", "dose"))

Group a Data Frame by One or more Variables

Description

Group a data frame by one or more variables. Supports standard and non standard evaluation.

Usage

df_group_by(data, ..., vars = NULL)

Arguments

Examples

# Non standard evaluation
by_dose <- head(ToothGrowth) %>%
   df_group_by(dose)
by_dose

# Standard evaluation
head(ToothGrowth) %>%
   df_group_by(vars = c("dose", "supp"))

Functions to Label Data Frames by Grouping Variables

Description

Functions to label data frame rows by one or multiple grouping variables.

Usage

df_label_both(data, ..., vars = NULL, label_col = "label", sep = c(", ", ":"))

df_label_value(data, ..., vars = NULL, label_col = "label", sep = ", ")

Arguments

Value

a modified data frame with a column containing row labels.

Functions

• df_label_both(): Displays both the variable name and the factor value.

• df_label_value(): Displays only the value of a factor.

Examples

# Data preparation
df <- head(ToothGrowth)

# Labelling: Non standard evaluation
df %>%
  df_label_both(dose, supp)

# Standard evaluation
df %>%
  df_label_both(dose, supp)

# Nesting the data then label each subset by groups
ToothGrowth %>%
  df_nest_by(dose, supp) %>%
  df_label_both(supp, dose)

Nest a Tibble By Groups

Description

Nest a tibble data frame using grouping specification. Supports standard and non standard evaluation.

Usage

df_nest_by(data, ..., vars = NULL)

Arguments

Value

A tbl with one row per unique combination of the grouping variables. The first columns are the grouping variables, followed by a list column of tibbles with matching rows of the remaining columns.

Examples

# Non standard evaluation
ToothGrowth %>%
 df_nest_by(dose, supp)

# Standard evaluation
ToothGrowth %>%
 df_nest_by(vars = c("dose", "supp"))

Select Columns in a Data Frame

Description

A wrapper around the select() function for selection data frame columns. Supports standard and non standard evaluations. Usefull to easily program with dplyr

Usage

df_select(data, ..., vars = NULL)

Arguments

Value

a data frame

Examples

df <- head(ToothGrowth)
df

# Select column using standard evaluation
df %>% df_select(vars = c("dose", "len"))

# Select column using non-standard evaluation
df %>% df_select(dose, len)

Split a Data Frame into Subset

Description

Split a data frame by groups into subsets or data panel. Very similar to the function df_nest_by(). The only difference is that, it adds label to each data subset. Labels are the combination of the grouping variable levels. The column holding labels are named "label".

Usage

df_split_by(
  data,
  ...,
  vars = NULL,
  label_col = "label",
  labeller = df_label_both,
  sep = c(", ", ":")
)

Arguments

Value

A tbl with one row per unique combination of the grouping variables. The first columns are the grouping variables, followed by a list column of tibbles with matching rows of the remaining columns, and a column named label, containing labels.

Examples

# Split a data frame
# :::::::::::::::::::::::::::::::::::::::::::::::::
# Create a grouped data
res <- ToothGrowth %>%
  df_split_by(dose, supp)
res

# Show subsets
res$data

# Add panel/subset labels
res <- ToothGrowth %>%
  df_split_by(dose, supp)
res

Unite Multiple Columns into One

Description

Paste together multiple columns into one. Wrapper arround unite() that supports standard and non standard evaluation.

Usage

df_unite(data, col, ..., vars = NULL, sep = "_", remove = TRUE, na.rm = FALSE)

df_unite_factors(
  data,
  col,
  ...,
  vars = NULL,
  sep = "_",
  remove = TRUE,
  na.rm = FALSE
)

Arguments

Functions

• df_unite(): Unite multiple columns into one.

• df_unite_factors(): Unite factor columns. First, order factors levels then merge them into one column. The output column is a factor.

Examples

# Non standard evaluation
head(ToothGrowth) %>%
 df_unite(col = "dose_supp", dose, supp)

# Standard evaluation
head(ToothGrowth) %>%
 df_unite(col = "dose_supp", vars = c("dose", "supp"))

Alternative to dplyr::do for Doing Anything

Description

Provides a flexible alternative to the dplyr:do() function. Technically it uses nest() + mutate() + map() to apply arbitrary computation to a grouped data frame.

The output is a data frame. If the applied function returns a data frame, then the output will be automatically unnested. Otherwise, the output includes the grouping variables and a column named ".results." (by default), which is a "list-columns" containing the results for group combinations.

Usage

doo(data, .f, ..., result = ".results.")

Arguments

Value

a data frame

Examples

# Custom function
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
stat_test <- function(data, formula){
  t.test(formula, data) %>%
    tidy()
}
# Example 1: pipe-friendly stat_test().
# Two possibilities of usage are available
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Use this
ToothGrowth %>%
  group_by(dose) %>%
  doo(~stat_test(data =., len ~ supp))

# Or this
ToothGrowth %>%
  group_by(dose) %>%
  doo(stat_test, len ~ supp)

# Example 2: R base function t.test() (not pipe friendly)
# One possibility of usage
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
comparisons <- ToothGrowth %>%
   group_by(dose) %>%
   doo(~t.test(len ~ supp, data =.))
comparisons
comparisons$.results.

# Example 3: R base function combined with tidy()
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ToothGrowth %>%
   group_by(dose) %>%
   doo(~t.test(len ~ supp, data =.) %>% tidy())

Dunn's Test of Multiple Comparisons

Description

Performs Dunn's test for pairwise multiple comparisons of the ranked data. The mean rank of the different groups is compared. Used for post-hoc test following Kruskal-Wallis test.

The default of the rstatix::dunn_test() function is to perform a two-sided Dunn test like the well known commercial softwares, such as SPSS and GraphPad. This is not the case for some other R packages (dunn.test and jamovi), where the default is to perform one-sided test. This discrepancy is documented at https://github.com/kassambara/rstatix/issues/50.

If a reference group is specified (via ref.group), then each of the remaining group levels is compared only to the reference (control) group, and the p-value adjustment for multiple comparisons is computed over only these k - 1 comparisons (instead of all k(k - 1)/2 pairwise comparisons). Note that this affects the adjusted p-values: it is not equivalent to filtering the full pairwise result afterwards, which would still adjust over all pairs.

See the Datanovia tutorial Kruskal-Wallis Test in R for a worked walkthrough.

Usage

dunn_test(
  data,
  formula,
  p.adjust.method = "holm",
  ref.group = NULL,
  detailed = FALSE,
  effect.size = FALSE
)

Arguments

Details

DunnTest performs the post hoc pairwise multiple comparisons procedure appropriate to follow up a Kruskal-Wallis test, which is a non-parametric analog of the one-way ANOVA. The Wilcoxon rank sum test, itself a non-parametric analog of the unpaired t-test, is possibly intuitive, but inappropriate as a post hoc pairwise test, because (1) it fails to retain the dependent ranking that produced the Kruskal-Wallis test statistic, and (2) it does not incorporate the pooled variance estimate implied by the null hypothesis of the Kruskal-Wallis test.

Value

return a data frame with some of the following columns:

• .y.: the y (outcome) variable used in the test.

• group1,group2: the compared groups in the pairwise tests.

• n1,n2: Sample counts.

• estimate: mean ranks difference.

• estimate1, estimate2: show the mean rank values of the two groups, respectively.

• statistic: Test statistic (z-value) used to compute the p-value.

• p: p-value.

• p.adj: the adjusted p-value.

• method: the statistical test used to compare groups.

• p.signif, p.adj.signif: the significance level of p-values and adjusted p-values, respectively.

The returned object has an attribute called args, which is a list holding the test arguments.

References

Dunn, O. J. (1964) Multiple comparisons using rank sums Technometrics, 6(3):241-252.

See Also

The Datanovia tutorial: Kruskal-Wallis Test in R.

Examples

# Simple test
ToothGrowth %>% dunn_test(len ~ dose)

# Comparison against a reference (control) group
# each group is compared to the reference; the p-value
# adjustment corrects for only these k - 1 comparisons
ToothGrowth %>% dunn_test(len ~ dose, ref.group = "0.5")

# Grouped data
ToothGrowth %>%
  group_by(supp) %>%
  dunn_test(len ~ dose)

Dunnett's Many-to-One Comparisons Test

Description

Performs Dunnett's test for comparing each of several treatment groups against a single control (reference) group. Unlike all-pairwise post-hoc tests, Dunnett's procedure controls the family-wise error rate over only the k - 1 treatment-vs-control comparisons, using the exact multivariate-t distribution (which accounts for the correlation between the comparisons that share the control group).

This is a pipe-friendly wrapper around emmeans::emmeans() + emmeans::contrast() (with adjust = "mvt"), so the emmeans package must be installed. The results match DescTools::DunnettTest() and multcomp::glht(..., mcp(... = "Dunnett")).

See the Datanovia tutorial One-Way ANOVA in R for a worked walkthrough.

Usage

dunnett_test(
  data,
  formula,
  ref.group = NULL,
  conf.level = 0.95,
  detailed = FALSE
)

Arguments

Value

a data frame with some of the following columns:

• .y.: the outcome variable used in the test.

• group1,group2: the compared groups; group1 is the control (reference) and group2 is the treatment, consistent with the ref.group convention of t_test()/wilcox_test()/ dunn_test()/emmeans_test().

• n1,n2: sample sizes of the control and treatment groups.

• estimate: the estimated mean difference group1 - group2 (control minus treatment).

• conf.low,conf.high: simultaneous confidence interval for the difference.

• statistic: the t-statistic.

• df: degrees of freedom.

• p.adj: the Dunnett-adjusted p-value.

• method: the statistical test used.

• p.adj.signif: the significance level of the adjusted p-value.

The estimate, confidence interval, se and method columns are returned only when detailed = TRUE.

The returned object has an attribute called args, which is a list holding the test arguments.

References

Dunnett, C. W. (1955) A multiple comparison procedure for comparing several treatments with a control. Journal of the American Statistical Association, 50, 1096-1121.

See Also

tukey_hsd(), games_howell_test(), emmeans_test() The Datanovia tutorial: One-Way ANOVA in R.

Examples

if (requireNamespace("emmeans", quietly = TRUE)) {

# Compare each dose to the control dose ("0.5")
ToothGrowth %>% dunnett_test(len ~ dose)

# Detailed output (estimate + simultaneous confidence interval)
ToothGrowth %>% dunnett_test(len ~ dose, detailed = TRUE)

# Grouped data
ToothGrowth %>%
  group_by(supp) %>%
  dunnett_test(len ~ dose)

}

Pairwise Comparisons of Estimated Marginal Means

Description

Performs pairwise comparisons between groups using the estimated marginal means. Pipe-friendly wrapper arround the functions emmans() + contrast() from the emmeans package, which need to be installed before using this function. This function is useful for performing post-hoc analyses following ANOVA/ANCOVA tests.

See the Datanovia tutorial One-Way ANOVA in R for a worked walkthrough.

Usage

emmeans_test(
  data,
  formula,
  covariate = NULL,
  ref.group = NULL,
  comparisons = NULL,
  p.adjust.method = "bonferroni",
  conf.level = 0.95,
  model = NULL,
  detailed = FALSE
)

get_emmeans(emmeans.test)

Arguments

Value

return a data frame with some the following columns:

• .y.: the y variable used in the test.

• group1,group2: the compared groups in the pairwise tests.

• statistic: Test statistic (t.ratio) used to compute the p-value.

• df: degrees of freedom.

• p: p-value.

• p.adj: the adjusted p-value.

• method: the statistical test used to compare groups.

• p.signif, p.adj.signif: the significance level of p-values and adjusted p-values, respectively.

• estimate: estimate of the effect size, that is the difference between the two emmeans (estimated marginal means).

• conf.low,conf.high: Lower and upper bound on a confidence interval of the estimate.

The returned object has an attribute called args, which is a list holding the test arguments. It has also an attribute named "emmeans", a data frame containing the groups emmeans.

Functions

• get_emmeans(): returns the estimated marginal means from an object of class emmeans_test

See Also

The Datanovia tutorial: One-Way ANOVA in R.

Examples

if (requireNamespace("emmeans", quietly = TRUE)) {

# Data preparation
df <- ToothGrowth
df$dose <- as.factor(df$dose)

# Pairwise comparisons
res <- df %>%
 group_by(supp) %>%
 emmeans_test(len ~ dose, p.adjust.method = "bonferroni")
res

# Display estimated marginal means
attr(res, "emmeans")

# Show details
df %>%
 group_by(supp) %>%
 emmeans_test(len ~ dose, p.adjust.method = "bonferroni", detailed = TRUE)

# Marginal means averaged over another factor (e.g. a 2x3 design).
# Fit the full model and pass it with `model =` so that the estimated
# marginal means for `dose` are averaged over `supp` (instead of fitting
# `len ~ dose` alone, which would ignore `supp`):
model <- lm(len ~ supp * dose, data = df)
df %>% emmeans_test(len ~ dose, model = model)

# Repeated-measures / mixed designs: pass a fitted within-subject model
# (e.g. stats::aov() with an Error() term, or nlme::lme()) with `model =`:

set.seed(123)
d <- data.frame(
  id    = factor(rep(1:10, 3)),
  time  = factor(rep(c("t1", "t2", "t3"), each = 10)),
  score = rnorm(30)
)
rm_model <- stats::aov(score ~ time + Error(id / time), data = d)
d %>% emmeans_test(score ~ time, model = rm_model)


}

Effect Size for ANOVA

Description

Compute eta-squared and partial eta-squared for all terms in an ANOVA model.

See the Datanovia tutorial One-Way ANOVA in R for a worked walkthrough.

Usage

eta_squared(model, ci = NULL)

partial_eta_squared(model, ci = NULL)

Arguments

Value

a named numeric vector of effect sizes, one per model term; or, when ci is a confidence level, a tibble with the columns Effect, effsize, conf.low and conf.high.

Functions

• eta_squared(): compute eta squared

• partial_eta_squared(): compute partial eta squared.

References

Steiger, J. H. (2004). Beyond the F test: Effect size confidence intervals and tests of close fit in the analysis of variance and contrast analysis. Psychological Methods, 9, 164-182.

See Also

The Datanovia tutorial: One-Way ANOVA in R.

Examples

# Data preparation
df <- ToothGrowth
df$dose <- as.factor(df$dose)

# Compute ANOVA
res.aov <- aov(len ~ supp*dose, data = df)
summary(res.aov)

# Effect size
eta_squared(res.aov)
partial_eta_squared(res.aov)

# Effect size with confidence interval
eta_squared(res.aov, ci = 0.95)
partial_eta_squared(res.aov, ci = 0.95)

Build Factorial Designs for ANOVA

Description

Provides helper functions to build factorial design for easily computing ANOVA using the Anova() function. This might be very useful for repeated measures ANOVA, which is hard to set up with the car package.

Usage

factorial_design(data, dv, wid, between, within, covariate)

Arguments

Value

a list with the following components:

• the specified arguments: dv, wid, between, within

• data: the original data (long format) or independent ANOVA. The wide format is returned for repeated measures ANOVA.

• idata: an optional data frame giving the levels of factors defining the intra-subject model for multivariate repeated-measures data.

• idesign: a one-sided model formula using the “data” in idata and specifying the intra-subject design.

• repeated: logical. Value is TRUE when the data is a repeated design.

• lm_formula: the formula used to build the lm model.

• lm_data: the data used to build the lm model. Can be either in a long format (i.e., the original data for independent measures ANOVA) or in a wide format (case of repeated measures ANOVA).

• model: the lm model

Author(s)

Alboukadel Kassambara, alboukadel.kassambara@gmail.com

See Also

anova_test(), anova_summary()

Examples

# Load data
#:::::::::::::::::::::::::::::::::::::::
data("ToothGrowth")
df <- ToothGrowth
head(df)

# Repeated measures designs
#:::::::::::::::::::::::::::::::::::::::::
# Prepare the data
df$id <- rep(1:10, 6) # Add individuals id
head(df)
# Build factorial designs
design <- factorial_design(df, dv = len, wid = id, within = c(supp, dose))
design
# Easily perform repeated measures ANOVA using the car package
res.anova <- Anova(design$model, idata = design$idata, idesign = design$idesign, type = 3)
summary(res.anova, multivariate = FALSE)

# Independent measures designs
#:::::::::::::::::::::::::::::::::::::::::
# Build factorial designs
df$id <- 1:nrow(df)
design <- factorial_design(df, dv = len, wid = id, between = c(supp, dose))
design
# Perform ANOVA
Anova(design$model, type = 3)

Fisher's Exact Test for Count Data

Description

Performs Fisher's exact test for testing the null of independence of rows and columns in a contingency table.

Wrappers around the R base function fisher.test() but have the advantage of performing pairwise and row-wise fisher tests, the post-hoc tests following a significant chi-square test of homogeneity for 2xc and rx2 contingency tables.

See the Datanovia tutorial Fisher’s Exact Test in R for a worked walkthrough.

Usage

fisher_test(
  xtab,
  workspace = 2e+05,
  alternative = "two.sided",
  conf.int = TRUE,
  conf.level = 0.95,
  simulate.p.value = FALSE,
  B = 2000,
  detailed = FALSE,
  ...
)

pairwise_fisher_test(xtab, p.adjust.method = "holm", detailed = FALSE, ...)

row_wise_fisher_test(xtab, p.adjust.method = "holm", detailed = FALSE, ...)

Arguments

Value

return a data frame with some the following columns:

• group: the categories in the row-wise proportion tests.

• p: p-value.

• p.adj: the adjusted p-value.

• method: the used statistical test.

• p.signif, p.adj.signif: the significance level of p-values and adjusted p-values, respectively.

• estimate: an estimate of the odds ratio. Only present in the 2 by 2 case.

• alternative: a character string describing the alternative hypothesis.

• conf.low,conf.high: a confidence interval for the odds ratio. Only present in the 2 by 2 case and if argument conf.int = TRUE.

The returned object has an attribute called args, which is a list holding the test arguments.

Functions

• fisher_test(): performs Fisher's exact test for testing the null of independence of rows and columns in a contingency table with fixed marginals. Wrapper around the function fisher.test().

• pairwise_fisher_test(): pairwise comparisons between proportions, a post-hoc tests following a significant Fisher's exact test of homogeneity for 2xc design.

• row_wise_fisher_test(): performs row-wise Fisher's exact test of count data, a post-hoc tests following a significant chi-square test of homogeneity for rx2 contingency table. The test is conducted for each category (row).

See Also

The Datanovia tutorial: Fisher’s Exact Test in R.

Examples

# Comparing two proportions
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Data: frequencies of smokers between two groups
xtab <- as.table(rbind(c(490, 10), c(400, 100)))
dimnames(xtab) <- list(
  group = c("grp1", "grp2"),
  smoker = c("yes", "no")
)
xtab
# compare the proportion of smokers
fisher_test(xtab, detailed = TRUE)

# Homogeneity of proportions between groups
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# H0: the proportion of smokers is similar in the four groups
# Ha:  this proportion is different in at least one of the populations.
#
# Data preparation
grp.size <- c( 106, 113, 156, 102 )
smokers  <- c( 50, 100, 139, 80 )
no.smokers <- grp.size - smokers
xtab <- as.table(rbind(
  smokers,
  no.smokers
))
dimnames(xtab) <- list(
  Smokers = c("Yes", "No"),
  Groups = c("grp1", "grp2", "grp3", "grp4")
)
xtab

# Compare the proportions of smokers between groups
fisher_test(xtab, detailed = TRUE)

# Pairwise comparison between groups
pairwise_fisher_test(xtab)


# Pairwise proportion tests
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Data: Titanic
xtab <- as.table(rbind(
  c(122, 167, 528, 673),
  c(203, 118, 178, 212)
))
dimnames(xtab) <- list(
  Survived = c("No", "Yes"),
  Class = c("1st", "2nd", "3rd", "Crew")
)
xtab
# Compare the proportion of survived between groups
pairwise_fisher_test(xtab)

# Row-wise proportion tests
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Data: Titanic
xtab <- as.table(rbind(
  c(180, 145), c(179, 106),
  c(510, 196), c(862, 23)
))
dimnames(xtab) <- list(
  Class = c("1st", "2nd", "3rd", "Crew"),
  Gender = c("Male", "Female")
)
xtab
# Compare the proportion of males and females in each category
row_wise_fisher_test(xtab)

# A r x c table  Agresti (2002, p. 57) Job Satisfaction
Job <- matrix(c(1,2,1,0, 3,3,6,1, 10,10,14,9, 6,7,12,11), 4, 4,
              dimnames = list(income = c("< 15k", "15-25k", "25-40k", "> 40k"),
                             satisfaction = c("VeryD", "LittleD", "ModerateS", "VeryS")))
fisher_test(Job)
fisher_test(Job, simulate.p.value = TRUE, B = 1e5)

Fligner-Killeen Test

Description

Provides a pipe-friendly framework to perform the Fligner-Killeen test, a non-parametric (rank-based) test of the homogeneity of group variances. It is robust against departures from normality and is a useful alternative to levene_test(). Wrapper around the function fligner.test().

See the Datanovia tutorial Homogeneity of Variance Test in R for a worked walkthrough.

Usage

fligner_test(data, formula, ...)

Arguments

Value

return a data frame with the following columns:

• .y.: the y variable used in the test.

• n: sample count.

• statistic: the Fligner-Killeen test statistic (a chi-squared statistic) used to compute the p-value.

• df: the degrees of freedom.

• p: p-value.

• method: the statistical test used to compare groups.

See Also

levene_test The Datanovia tutorial: Homogeneity of Variance Test in R.

Examples

# Load data
#:::::::::::::::::::::::::::::::::::::::
data("ToothGrowth")
df <- ToothGrowth

# Fligner-Killeen test
#:::::::::::::::::::::::::::::::::::::::::
df %>% fligner_test(len ~ dose)

# Grouped data
df %>%
  group_by(supp) %>%
  fligner_test(len ~ dose)

Compute Frequency Table

Description

compute frequency table.

Usage

freq_table(data, ..., vars = NULL, na.rm = TRUE)

Arguments

Value

a data frame

Examples

data("ToothGrowth")
ToothGrowth %>% freq_table(supp, dose)

# Grouped data: frequencies are computed within each group
ToothGrowth %>% group_by(supp) %>% freq_table(dose)

Conover's All-Pairs Comparisons Test for Friedman Rank Sums

Description

Performs Conover's all-pairs comparison test (also known as the Durbin-Conover test) for a two-way balanced complete block design, following a significant Friedman rank sum test. It is the repeated-measures analogue of Conover's test for the Kruskal-Wallis design: the within-block ranks are compared pairwise using the pooled rank variance and the statistic is referred to a t-distribution with (b - 1)(k - 1) degrees of freedom (b blocks, k treatments). It should only be used as a post-hoc procedure when the Friedman test is itself significant (Conover, 1999).

If a reference group is specified (via ref.group), then each of the remaining treatments is compared only to the reference (control) treatment, and the p-value adjustment for multiple comparisons is computed over only these k - 1 comparisons (as for dunn_test()).

See the Datanovia tutorial Friedman Test in R for a worked walkthrough.

Usage

friedman_conover_test(
  data,
  formula,
  p.adjust.method = "holm",
  ref.group = NULL,
  detailed = FALSE
)

Arguments

Details

For a balanced complete block design with b blocks and k treatments, the observations within each block are ranked. Let R_j be the sum of the within-block ranks for treatment j and let A = \sum r^2 be the sum of the squared within-block ranks. The pairwise statistic for treatments i and j is

t_{ij} = \frac{R_i - R_j}{\sqrt{\dfrac{2\,(b\,A - \sum_j R_j^2)}{(b - 1)(k - 1)}}}

referred to a t-distribution with (b - 1)(k - 1) degrees of freedom. This is the Conover (1999) post-hoc, also known as the Durbin-Conover test. In the returned table each row is oriented with i = group2 and j = group1: estimate is R_{group2} - R_{group1} and statistic carries its sign, the same convention as conover_test() and dunn_test().

The p-values match PMCMRplus::frdAllPairsConoverTest(). That function reports the t statistic for the reversed comparison, so its sign is the opposite of the one returned here; the magnitude is the same.

Value

return a data frame with some of the following columns:

• .y.: the y (outcome) variable used in the test.

• group1,group2: the compared treatments in the pairwise tests.

• n1,n2: the number of blocks (subjects) contributing to each treatment.

• estimate: the rank-sum difference.

• estimate1, estimate2: the rank sums of the two treatments, respectively.

• statistic: Test statistic (t-value) used to compute the p-value.

• df: degrees of freedom ((b - 1)(k - 1)).

• p: p-value.

• p.adj: the adjusted p-value.

• method: the statistical test used to compare groups.

• p.adj.signif: the significance level of the adjusted p-values.

The returned object has an attribute called args, which is a list holding the test arguments.

References

Conover, W. J. (1999) Practical Nonparametric Statistics, 3rd edition. Wiley.

See Also

friedman_test, friedman_nemenyi_test, friedman_effsize The Datanovia tutorial: Friedman Test in R.

Examples

# A balanced complete block design: 3 treatments measured on 6 subjects
df <- data.frame(
  id        = factor(rep(1:6, 3)),
  treatment = factor(rep(c("A", "B", "C"), each = 6)),
  score     = c(4, 6, 3, 5, 4, 5,    7, 8, 6, 7, 9, 6,    6, 9, 7, 8, 8, 9)
)

# Omnibus Friedman test
df %>% friedman_test(score ~ treatment | id)

# Conover (Durbin-Conover) all-pairs post-hoc
df %>% friedman_conover_test(score ~ treatment | id)

# Comparison against a reference (control) treatment
df %>% friedman_conover_test(score ~ treatment | id, ref.group = "A")

Friedman Test Effect Size (Kendall's W Value)

Description

Compute the effect size estimate (referred to as w) for Friedman test: W = X2/N(K-1); where W is the Kendall's W value; X2 is the Friedman test statistic value; N is the sample size. k is the number of measurements per subject.

The Kendall’s W coefficient assumes the value from 0 (indicating no relationship) to 1 (indicating a perfect relationship).

Kendalls uses the Cohen’s interpretation guidelines of 0.1 - < 0.3 (small effect), 0.3 - < 0.5 (moderate effect) and >= 0.5 (large effect)

Confidence intervals are calculated by bootstap.

See the Datanovia tutorial Friedman Test in R for a worked walkthrough.

Usage

friedman_effsize(
  data,
  formula,
  ci = FALSE,
  conf.level = 0.95,
  ci.type = "perc",
  nboot = 1000,
  ...,
  boot.parallel = getOption("boot.parallel", "no"),
  boot.ncpus = getOption("boot.ncpus", 1L)
)

Arguments

Value

return a data frame with some of the following columns:

• .y.: the y variable used in the test.

• n: Sample counts.

• effsize: estimate of the effect size.

• magnitude: magnitude of effect size.

• conf.low,conf.high: lower and upper bound of the effect size confidence interval.

References

Maciej Tomczak and Ewa Tomczak. The need to report effect size estimates revisited. An overview of some recommended measures of effect size. Trends in Sport Sciences. 2014; 1(21):19-25.

See Also

The Datanovia tutorial: Friedman Test in R.

Examples

# Load data
#:::::::::::::::::::::::::::::::::::::::
data("ToothGrowth")
df <- ToothGrowth %>%
    filter(supp == "VC") %>%
    mutate(id = rep(1:10, 3))
head(df)

# Friedman test effect size
#:::::::::::::::::::::::::::::::::::::::::
df %>% friedman_effsize(len ~ dose | id)

Nemenyi Post-Hoc Test for Friedman Rank Sums

Description

Performs the Nemenyi (Wilcoxon-Nemenyi-McDonald-Thompson) all-pairs post-hoc test for a two-way balanced complete block design, following a significant Friedman rank sum test. The treatment rank sums are compared pairwise and the test statistic is referred to the studentized range distribution, which already accounts for the multiplicity of the all-pairs comparisons (so, as for tukey_hsd(), there is no separate p-value adjustment step). It should only be used as a post-hoc procedure when the Friedman test is itself significant.

The Nemenyi test is the rank-based, repeated-measures analogue of Tukey's HSD. Unlike friedman_conover_test() (the Durbin-Conover test), it does not borrow the residual rank variance and is therefore more conservative.

See the Datanovia tutorial Friedman Test in R for a worked walkthrough.

Usage

friedman_nemenyi_test(data, formula, detailed = FALSE)

Arguments

Details

For a balanced complete block design with b blocks and k treatments, the observations within each block are ranked. Let R_j be the sum of the within-block ranks for treatment j. The pairwise statistic for treatments i and j is

q_{ij} = \frac{|R_i - R_j|}{\sqrt{b\,k\,(k+1)/12}}

and the p-value is obtained from the studentized range distribution with k groups and infinite degrees of freedom.

The p-values match PMCMRplus::frdAllPairsNemenyiTest(). That function reports the magnitude of the statistic; the value returned here carries the sign of the rank-sum difference between the two groups, so the two agree only where that difference is positive.

Value

return a data frame with some of the following columns:

• .y.: the y (outcome) variable used in the test.

• group1,group2: the compared treatments in the pairwise tests.

• n1,n2: the number of blocks (subjects) contributing to each treatment.

• estimate: the rank-sum difference.

• estimate1, estimate2: the rank sums of the two treatments, respectively.

• statistic: the studentized-range test statistic.

• p.adj: the p-value (already adjusted for multiple comparisons via the studentized range distribution).

• method: the statistical test used to compare groups.

• p.adj.signif: the significance level of the adjusted p-values.

The returned object has an attribute called args, which is a list holding the test arguments.

References

Nemenyi, P. (1963) Distribution-free Multiple Comparisons. PhD Thesis, Princeton University.

Hollander, M., Wolfe, D. A. (1973) Nonparametric Statistical Methods. Wiley.

See Also

friedman_test, friedman_conover_test, friedman_effsize The Datanovia tutorial: Friedman Test in R.

Examples

# A balanced complete block design: 3 treatments measured on 6 subjects
df <- data.frame(
  id        = factor(rep(1:6, 3)),
  treatment = factor(rep(c("A", "B", "C"), each = 6)),
  score     = c(4, 6, 3, 5, 4, 5,    7, 8, 6, 7, 9, 6,    6, 9, 7, 8, 8, 9)
)

# Omnibus Friedman test
df %>% friedman_test(score ~ treatment | id)

# Nemenyi all-pairs post-hoc
df %>% friedman_nemenyi_test(score ~ treatment | id)

Friedman Rank Sum Test

Description

Provides a pipe-friendly framework to perform a Friedman rank sum test, which is the non-parametric alternative to the one-way repeated measures ANOVA test. Wrapper around the function friedman.test().

See the Datanovia tutorial Friedman Test in R for a worked walkthrough.

Usage

friedman_test(data, formula, ...)

Arguments

Value

return a data frame with the following columns:

• .y.: the y (dependent) variable used in the test.

• n: sample count.

• statistic: the value of Friedman's chi-squared statistic, used to compute the p-value.

• p: p-value.

• method: the statistical test used to compare groups.

See Also

The Datanovia tutorial: Friedman Test in R.

Examples

# Load data
#:::::::::::::::::::::::::::::::::::::::
data("ToothGrowth")
df <- ToothGrowth %>%
    filter(supp == "VC") %>%
    mutate(id = rep(1:10, 3))
head(df)

# Friedman rank sum test
#:::::::::::::::::::::::::::::::::::::::::
df %>% friedman_test(len ~ dose | id)

Games Howell Post-hoc Tests

Description

Performs Games-Howell test, which is used to compare all possible combinations of group differences when the assumption of homogeneity of variances is violated. This post hoc test provides confidence intervals for the differences between group means and shows whether the differences are statistically significant.

The test is based on Welch’s degrees of freedom correction and uses Tukey’s studentized range distribution for computing the p-values. The test compares the difference between each pair of means with appropriate adjustment for the multiple testing. So there is no need to apply additional p-value corrections.

See the Datanovia tutorial One-Way ANOVA in R for a worked walkthrough.

Usage

games_howell_test(
  data,
  formula,
  conf.level = 0.95,
  detailed = FALSE,
  effect.size = FALSE
)

Arguments

Details

The Games-Howell method is an improved version of the Tukey-Kramer method and is applicable in cases where the equivalence of variance assumption is violated. It is a t-test using Welch’s degree of freedom. This method uses a strategy for controlling the type I error for the entire comparison and is known to maintain the preset significance level even when the size of the sample is different. However, the smaller the number of samples in each group, the it is more tolerant the type I error control. Thus, this method can be applied when the number of samples is six or more.

Because the test relies on Welch's variance correction, comparisons involving a group with zero variance (constant values) or undefined variance (a single observation) can be undefined; such pairs are returned as NA (with a warning) while all other comparisons are computed as usual.

Value

return a data frame with some of the following columns:

• .y.: the y (outcome) variable used in the test.

• group1,group2: the compared groups in the pairwise tests.

• n1,n2: Sample counts.

• estimate, conf.low, conf.high: mean difference and its confidence intervals.

• statistic: Test statistic (t-value) used to compute the p-value.

• df: degrees of freedom calculated using Welch’s correction.

• p.adj: adjusted p-value using Tukey's method.

• method: the statistical test used to compare groups.

• p.adj.signif: the significance level of p-values.

The returned object has an attribute called args, which is a list holding the test arguments.

References

• Aaron Schlege, https://rpubs.com/aaronsc32/games-howell-test.

• Sangseok Lee, Dong Kyu Lee. What is the proper way to apply the multiple comparison test?. Korean J Anesthesiol. 2018;71(5):353-360.

See Also

The Datanovia tutorial: One-Way ANOVA in R.

Examples

# Simple test
ToothGrowth %>% games_howell_test(len ~ dose)

# Grouped data
ToothGrowth %>%
  group_by(supp) %>%
  games_howell_test(len ~ dose)

Create a List of Possible Comparisons Between Groups

Description

Create a list of possible pairwise comparisons between groups. If a reference group is specified, only comparisons against reference will be kept.

Usage

get_comparisons(data, variable, ref.group = NULL)

Arguments

Value

a list of all possible pairwise comparisons.

Examples

# All possible pairwise comparisons
ToothGrowth %>%
  get_comparisons("dose")

# Comparisons against reference groups
ToothGrowth %>%
  get_comparisons("dose", ref.group = "0.5")

# Comparisons against  all (basemean)
ToothGrowth %>%
  get_comparisons("dose", ref.group = "all")

Compute Mode

Description

Compute the mode in a given vector. Mode is the most frequent value.

Usage

get_mode(x)

Arguments

Examples

# Mode of numeric vector
x <- c(1:5, 6, 6, 7:10)
get_mode(x)

# Bimodal
x <- c(1:5, 6, 6, 7, 8, 9, 9, 10)
get_mode(x)

# No mode
x <- c(1, 2, 3, 4, 5)
get_mode(x)

# Nominal vector
fruits <-  c(rep("orange", 10), rep("apple", 5), rep("lemon", 2))
get_mode(fruits)

Extract Label Information from Statistical Tests

Description

Extracts label information from statistical tests. Useful for labelling plots with test outputs.

See the Datanovia tutorial P-values from Tests on ggplots in R (rstatix) for a worked walkthrough.

Usage

get_pwc_label(stat.test, type = c("expression", "text"))

get_test_label(
  stat.test,
  description = NULL,
  p.col = "p",
  type = c("expression", "text"),
  correction = c("auto", "GG", "HF", "none"),
  row = NULL,
  detailed = FALSE,
  style = c("classic", "apa")
)

create_test_label(
  statistic.text,
  statistic,
  p,
  parameter = NA,
  description = NULL,
  n = NA,
  effect.size = NA,
  effect.size.text = NA,
  type = c("expression", "text"),
  detailed = FALSE,
  style = c("classic", "apa"),
  effect.size.ci = NA,
  effect.size.bounded = TRUE,
  effect.size.ci.level = 0.95
)

get_n(stat.test)

get_description(stat.test)

Arguments

Value

a text label or an expression to pass to a plotting function.

Functions

• get_pwc_label(): Extract label from pairwise comparisons.

• get_test_label(): Extract labels for statistical tests.

• create_test_label(): Create labels from user specified test results.

• get_n(): Extracts sample counts (n) from an rstatix test outputs. Returns a numeric vector.

• get_description(): Extracts the description of an rstatix test outputs. Returns a character vector.

References

American Psychological Association (2020). Publication Manual of the American Psychological Association (7th ed.).

See Also

The Datanovia tutorial: P-values from Tests on ggplots in R (rstatix).

Examples

# Load data
#:::::::::::::::::::::::::::::::::::::::
data("ToothGrowth")
df <- ToothGrowth

# One-way ANOVA test
#:::::::::::::::::::::::::::::::::::::::::
anov <- df %>% anova_test(len ~ dose)
get_test_label(anov, detailed = TRUE, type = "text")

# Two-way ANOVA test
#:::::::::::::::::::::::::::::::::::::::::
anov <- df %>% anova_test(len ~ supp*dose)
get_test_label(anov, detailed = TRUE, type = "text",
   description = "Two Way ANOVA")


# Kruskal-Wallis test
#:::::::::::::::::::::::::::::::::::::::::
kruskal<- df %>% kruskal_test(len ~ dose)
get_test_label(kruskal, detailed = TRUE, type = "text")

# Wilcoxon test
#:::::::::::::::::::::::::::::::::::::::::
# Unpaired test
wilcox <- df %>% wilcox_test(len ~ supp)
get_test_label(wilcox, detailed = TRUE, type = "text")
# Paired test
wilcox <- df %>% wilcox_test(len ~ supp, paired = TRUE)
get_test_label(wilcox, detailed = TRUE, type = "text")

# T test
#:::::::::::::::::::::::::::::::::::::::::
ttest <- df %>% t_test(len ~ dose)
get_test_label(ttest, detailed = TRUE, type = "text")


# Pairwise comparisons labels
#:::::::::::::::::::::::::::::::::::::::::
get_pwc_label(ttest, type = "text")


# Create test labels
#:::::::::::::::::::::::::::::::::::::::::
create_test_label(
  statistic.text = "F", statistic = 71.82,
  parameter = "4, 294",
  p = "<0.0001",
  description = "ANOVA",
  type = "text"
)


# Extract infos
#:::::::::::::::::::::::::::::::::::::::::
stat.test <- df %>% t_test(len ~ dose)
get_n(stat.test)
get_description(stat.test)

Compute Summary Statistics

Description

Compute summary statistics for one or multiple numeric variables.

See the Datanovia tutorial Descriptive Statistics in R for a worked walkthrough.

Usage

get_summary_stats(
  data,
  ...,
  type = c("full", "common", "robust", "five_number", "mean_sd", "mean_se", "mean_ci",
    "median_iqr", "median_mad", "quantile", "mean", "median", "min", "max"),
  show = NULL,
  probs = seq(0, 1, 0.25),
  digits = 3
)

Arguments

Value

A data frame containing descriptive statistics, such as:

• n: the number of individuals

• min: minimum

• max: maximum

• median: median

• mean: mean

• q1, q3: the first and the third quartile, respectively.

• iqr: interquartile range

• mad: median absolute deviation (see ?MAD)

• sd: standard deviation of the mean

• se: standard error of the mean

• ci: 95 percent confidence interval of the mean

When requested through show, the output can also contain:

• skewness: bias-corrected sample skewness

• kurtosis: bias-corrected sample excess kurtosis (0 for a normal distribution).

Both use the type-2 (bias-corrected) estimator, matching e1071 with type = 2: skewness = g_1\sqrt{n(n-1)}/(n-2) and kurtosis = [(n+1)g_2 + 6] (n-1)/[(n-2)(n-3)], where g_1 = m_3/m_2^{1.5} and g_2 = m_4/m_2^2 - 3. Skewness is NA for n < 3 and kurtosis for n < 4.

See Also

rstatix-programming for selecting columns by names held in strings (!!, {{ }}, vars=, all_of()). The Datanovia tutorial: Descriptive Statistics in R.

Examples

# Full summary statistics
data("ToothGrowth")
ToothGrowth %>% get_summary_stats(len)

# Summary statistics of grouped data
# Show only common summary
ToothGrowth %>%
  group_by(dose, supp) %>%
  get_summary_stats(len, type = "common")

# Robust summary statistics
ToothGrowth %>% get_summary_stats(len, type = "robust")

# Five number summary statistics
ToothGrowth %>% get_summary_stats(len, type = "five_number")

# Compute only mean and sd
ToothGrowth %>% get_summary_stats(len, type = "mean_sd")

# Compute full summary statistics but show only mean, sd, median, iqr
ToothGrowth %>%
    get_summary_stats(len, show = c("mean", "sd", "median", "iqr"))

# Include skewness and kurtosis (computed on demand via show)
ToothGrowth %>%
    get_summary_stats(len, show = c("mean", "sd", "skewness", "kurtosis"))

Autocompute P-value Positions For Plotting Significance

Description

Compute p-value x and y positions for plotting significance levels.

See the Datanovia tutorials P-values from Tests on ggplots in R (rstatix) and Auto P-values in ggplot with geom_pwc (ggpubr) for worked walkthroughs.

Usage

get_y_position(
  data,
  formula,
  fun = "max",
  ref.group = NULL,
  comparisons = NULL,
  step.increase = 0.12,
  y.trans = NULL,
  stack = FALSE,
  scales = c("fixed", "free", "free_y")
)

add_y_position(
  test,
  fun = "max",
  step.increase = 0.12,
  data = NULL,
  formula = NULL,
  ref.group = NULL,
  comparisons = NULL,
  y.trans = NULL,
  stack = FALSE,
  scales = c("fixed", "free", "free_y")
)

add_x_position(
  test,
  x = NULL,
  group = NULL,
  dodge = 0.8,
  scales = c("fixed", "free", "free_y")
)

add_xy_position(
  test,
  x = NULL,
  group = NULL,
  dodge = 0.8,
  stack = FALSE,
  fun = "max",
  step.increase = 0.12,
  scales = c("fixed", "free", "free_y"),
  ...
)

Arguments

Functions

• get_y_position(): compute the p-value y positions

• add_y_position(): add p-value y positions to an object of class rstatix_test

• add_x_position(): compute and add p-value x positions.

• add_xy_position(): compute and add both x and y positions.

Added columns

add_y_position() adds y.position; add_x_position() adds x, xmin, xmax; and add_xy_position() adds all of these. A column named groups (the pair of compared groups, used internally to position dodged grouped comparisons) is also added. groups is a reserved internal column: do not rely on it and avoid using groups as one of your own column names in the test object, since it is overwritten here. (It may be renamed to a dotted, less collision-prone name in a future major version.)

See Also

The Datanovia tutorials: P-values from Tests on ggplots in R (rstatix), Auto P-values in ggplot with geom_pwc (ggpubr).

Examples

# Data preparation
#::::::::::::::::::::::::::::::::::::
df <- ToothGrowth
df$dose <- as.factor(df$dose)
df$group <- factor(rep(c(1, 2), 30))
head(df)

# Stat tests
#::::::::::::::::::::::::::::::::::::
stat.test <- df %>%
  t_test(len ~ dose)
stat.test

# Add the test into box plots
#::::::::::::::::::::::::::::::::::::
stat.test <- stat.test %>%
  add_y_position()

 if(require("ggpubr")){
   ggboxplot(df, x = "dose", y = "len") +
     stat_pvalue_manual(stat.test, label = "p.adj.signif", tip.length = 0.01)
  }

Identify Univariate Outliers Using Boxplot Methods

Description

Detect outliers using boxplot methods. Boxplots are a popular and an easy method for identifying outliers. There are two categories of outlier: (1) outliers and (2) extreme points.

Values above Q3 + 1.5xIQR or below Q1 - 1.5xIQR are considered as outliers. Values above Q3 + 3xIQR or below Q1 - 3xIQR are considered as extreme points (or extreme outliers).

Q1 and Q3 are the first and third quartile, respectively. IQR is the interquartile range (IQR = Q3 - Q1).

Generally speaking, data points that are labelled outliers in boxplots are not considered as troublesome as those considered extreme points and might even be ignored. Note that, any NA and NaN are automatically removed before the quantiles are computed.

See the Datanovia tutorial Descriptive Statistics in R for a worked walkthrough.

Usage

identify_outliers(data, ..., variable = NULL)

is_outlier(x, coef = 1.5)

is_extreme(x)

Arguments

Value

• identify_outliers(). Returns the input data frame with two additional columns: "is.outlier" and "is.extreme", which hold logical values.

• is_outlier() and is_extreme(). Returns logical vectors.

Functions

• identify_outliers(): takes a data frame and extract rows suspected as outliers according to a numeric column. The following columns are added "is.outlier" and "is.extreme".

• is_outlier(): detect outliers in a numeric vector. Returns logical vector.

• is_extreme(): detect extreme points in a numeric vector. An alias of is_outlier(), where coef = 3. Returns logical vector.

See Also

The Datanovia tutorial: Descriptive Statistics in R.

Examples

# Generate a demo data
set.seed(123)
demo.data <- data.frame(
  sample = 1:20,
  score = c(rnorm(19, mean = 5, sd = 2), 50),
  gender = rep(c("Male", "Female"), each = 10)
)

# Identify outliers according to the variable score
demo.data %>%
  identify_outliers(score)

# Identify outliers by groups
demo.data %>%
  group_by(gender) %>%
  identify_outliers("score")

Kruskal-Wallis Effect Size

Description

Compute the effect size for Kruskal-Wallis test as the eta squared based on the H-statistic: eta2[H] = (H - k + 1)/(n - k); where H is the value obtained in the Kruskal-Wallis test; k is the number of groups; n is the total number of observations.

The eta-squared estimate assumes values from 0 to 1 and multiplied by 100 indicates the percentage of variance in the dependent variable explained by the independent variable. The interpretation values commonly in published litterature are: 0.01- < 0.06 (small effect), 0.06 - < 0.14 (moderate effect) and >= 0.14 (large effect).

Note that eta2[H] is a bias-corrected estimator, so the raw formula can return a small negative value for a near-null effect (very small H). In that case the estimate is floored to 0, keeping the reported effect size within its valid [0, 1] range.

Confidence intervals are calculated by bootstap.

See the Datanovia tutorial Kruskal-Wallis Test in R for a worked walkthrough.

Usage

kruskal_effsize(
  data,
  formula,
  ci = FALSE,
  conf.level = 0.95,
  ci.type = "perc",
  nboot = 1000,
  boot.parallel = getOption("boot.parallel", "no"),
  boot.ncpus = getOption("boot.ncpus", 1L),
  method = c("eta2", "epsilon2")
)

Arguments

Value

return a data frame with some of the following columns:

• .y.: the y variable used in the test.

• n: Sample counts.

• effsize: estimate of the effect size.

• magnitude: magnitude of effect size.

• conf.low,conf.high: lower and upper bound of the effect size confidence interval.

References

Maciej Tomczak and Ewa Tomczak. The need to report effect size estimates revisited. An overview of some recommended measures of effect size. Trends in Sport Sciences. 2014; 1(21):19-25.

http://imaging.mrc-cbu.cam.ac.uk/statswiki/FAQ/effectSize

http://www.psy.gla.ac.uk/~steve/best/effect.html

See Also

The Datanovia tutorial: Kruskal-Wallis Test in R.

Examples

# Load data
#:::::::::::::::::::::::::::::::::::::::
data("ToothGrowth")
df <- ToothGrowth

# Kruskal-wallis rank sum test
#:::::::::::::::::::::::::::::::::::::::::
df %>% kruskal_effsize(len ~ dose)

# Grouped data
df %>%
  group_by(supp) %>%
  kruskal_effsize(len ~ dose)

Kruskal-Wallis Test

Description

Provides a pipe-friendly framework to perform Kruskal-Wallis rank sum test. Wrapper around the function kruskal.test().

See the Datanovia tutorial Kruskal-Wallis Test in R for a worked walkthrough.

Usage

kruskal_test(data, formula, ...)

Arguments

Value

return a data frame with the following columns:

• .y.: the y variable used in the test.

• n: sample count.

• statistic: the kruskal-wallis rank sum statistic used to compute the p-value.

• p: p-value.

• method: the statistical test used to compare groups.

See Also

The Datanovia tutorial: Kruskal-Wallis Test in R.

Examples

# Load data
#:::::::::::::::::::::::::::::::::::::::
data("ToothGrowth")
df <- ToothGrowth

# Kruskal-wallis rank sum test
#:::::::::::::::::::::::::::::::::::::::::
df %>% kruskal_test(len ~ dose)

# Grouped data
df %>%
  group_by(supp) %>%
  kruskal_test(len ~ dose)

Two-Sample Kolmogorov-Smirnov Test

Description

Provides a pipe-friendly framework to perform the two-sample Kolmogorov-Smirnov test, comparing the (empirical) distributions of a numeric variable between two groups. Wrapper around the R base function ks.test().

When the grouping factor contains more than two levels, pairwise Kolmogorov-Smirnov tests are automatically performed, with p-value adjustment.

See the Datanovia tutorial Normality Test in R for a worked walkthrough.

Usage

ks_test(
  data,
  formula,
  comparisons = NULL,
  ref.group = NULL,
  p.adjust.method = "holm",
  alternative = "two.sided",
  exact = NULL,
  detailed = FALSE
)

Arguments

Value

return a data frame with some of the following columns:

• .y.: the y variable used in the test.

• group1,group2: the compared groups in the pairwise tests.

• n1,n2: sample counts.

• statistic: the value of the test statistic D (the maximum difference between the two empirical cumulative distribution functions).

• p: p-value.

• p.adj: the adjusted p-value.

• method: the statistical test used to compare groups.

• p.signif, p.adj.signif: the significance level of p-values and adjusted p-values, respectively.

• alternative: the alternative hypothesis.

The returned object has an attribute called args, which is a list holding the test arguments.

See Also

wilcox_test(), t_test() The Datanovia tutorial: Normality Test in R.

Examples

# Two-samples test
#:::::::::::::::::::::::::::::::::::::::::
ToothGrowth %>% ks_test(len ~ supp)

# Pairwise comparisons (more than two groups)
#:::::::::::::::::::::::::::::::::::::::::
ToothGrowth %>% ks_test(len ~ dose)

# Comparison against a reference group
#:::::::::::::::::::::::::::::::::::::::::
ToothGrowth %>% ks_test(len ~ dose, ref.group = "0.5")

Levene's Test

Description

Provide a pipe-friendly framework to easily compute Levene's test for homogeneity of variance across groups.

Wrapper around the function leveneTest(), which can additionally handles a grouped data.

See the Datanovia tutorial Homogeneity of Variance Test in R for a worked walkthrough.

Usage

levene_test(data, formula, center = median)

Arguments

Value

a data frame with the following columns: df1, df2 (df.residual), statistic and p.

See Also

The Datanovia tutorial: Homogeneity of Variance Test in R.

Examples

# Prepare the data
data("ToothGrowth")
df <- ToothGrowth
df$dose <- as.factor(df$dose)
# Compute Levene's Test
df %>% levene_test(len ~ dose)

# Grouped data
df %>%
  group_by(supp) %>%
  levene_test(len ~ dose)

Compute Mahalanobis Distance and Flag Multivariate Outliers

Description

Pipe-friendly wrapper around to the function mahalanobis(), which returns the squared Mahalanobis distance of all rows in x. Compared to the base function, it automatically flags multivariate outliers.

Mahalanobis distance is a common metric used to identify multivariate outliers. The larger the value of Mahalanobis distance, the more unusual the data point (i.e., the more likely it is to be a multivariate outlier).

The distance tells us how far an observation is from the center of the cloud, taking into account the shape (covariance) of the cloud as well.

To detect outliers, the calculated Mahalanobis distance is compared against a chi-square (X^2) distribution with degrees of freedom equal to the number of dependent (outcome) variables and an alpha level of 0.001.

The threshold to declare a multivariate outlier is determined using the function qchisq(0.999, df) , where df is the degree of freedom (i.e., the number of dependent variable used in the computation).

See the Datanovia tutorial Normality Test in R for a worked walkthrough.

Usage

mahalanobis_distance(data, ...)

Arguments

Value

Returns the input data frame with two additional columns: 1) "mahal.dist": Mahalanobis distance values; and 2) "is.outlier": logical values specifying whether a given observation is a multivariate outlier

See Also

The Datanovia tutorial: Normality Test in R.

Examples

# Compute mahalonobis distance and flag outliers if any
iris %>%
  doo(~mahalanobis_distance(.))

# Compute distance by groups and filter outliers
iris %>%
 group_by(Species) %>%
 doo(~mahalanobis_distance(.)) %>%
 filter(is.outlier == TRUE)

Make Clean Names

Description

Pipe-friendly function to make syntactically valid names out of character vectors.

Usage

make_clean_names(data)

Arguments

Value

a data frame or a vector depending on the input data

Examples

# Vector
make_clean_names(c("a and b", "a-and-b"))
make_clean_names(1:10)

# data frame
df <- data.frame(
`a and b` = 1:4,
`c and d` = 5:8,
 check.names = FALSE
)
df
make_clean_names(df)

McNemar's Chi-squared Test for Count Data

Description

Performs McNemar chi-squared test to compare paired proportions.

Wrappers around the R base function mcnemar.test(), but provide pairwise comparisons between multiple groups

See the Datanovia tutorial McNemar’s Test in R for a worked walkthrough.

Usage

mcnemar_test(x, y = NULL, correct = TRUE)

pairwise_mcnemar_test(
  data,
  formula,
  type = c("mcnemar", "exact"),
  correct = TRUE,
  p.adjust.method = "bonferroni"
)

Arguments

Value

return a data frame with the following columns:

• n: the number of participants.

• statistic: the value of McNemar's statistic.

• df the degrees of freedom of the approximate chi-squared distribution of the test statistic.

• p: p-value.

• p.adj: the adjusted p-value.

• method: the used statistical test.

• p.signif: the significance level of p-values.

For pairwise_mcnemar_test(), the statistic and df columns are returned for type = "mcnemar" (the default); they are omitted for type = "exact", which is based on an exact binomial test and has no chi-squared statistic or degrees of freedom.

The returned object has an attribute called args, which is a list holding the test arguments.

Functions

• mcnemar_test(): performs McNemar's chi-squared test for comparing two paired proportions

• pairwise_mcnemar_test(): performs pairwise McNemar's chi-squared test between multiple groups. Could be used for post-hoc tests following a significant Cochran's Q test.

See Also

The Datanovia tutorial: McNemar’s Test in R.

Examples

# Comparing two paired proportions
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Data: frequencies of smokers before and after interventions
xtab <- as.table(
  rbind(c(25, 6), c(21,10))
)
dimnames(xtab) <- list(
  before = c("non.smoker", "smoker"),
  after = c("non.smoker", "smoker")
)
xtab

# Compare the proportion of smokers
mcnemar_test(xtab)

# Comparing multiple related proportions
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Generate a demo data
mydata <- data.frame(
  outcome = c(0,1,1,0,0,1,0,1,1,1,1,1,0,0,1,1,0,1,0,1,1,0,0,1,0,1,1,0,0,1),
  treatment = gl(3,1,30,labels=LETTERS[1:3]),
  participant = gl(10,3,labels=letters[1:10])
)
mydata$outcome <- factor(
  mydata$outcome, levels = c(1, 0),
  labels = c("success", "failure")
  )
# Cross-tabulation
xtabs(~outcome + treatment, mydata)

# Compare the proportion of success between treatments
cochran_qtest(mydata, outcome ~ treatment|participant)

# pairwise comparisons between groups
pairwise_mcnemar_test(mydata, outcome ~ treatment|participant)

Exact Multinomial Test

Description

Performs an exact multinomial test. Alternative to the chi-square test of goodness-of-fit-test when the sample size is small.

Usage

multinom_test(x, p = rep(1/length(x), length(x)), detailed = FALSE)

Arguments

Value

return a data frame containing the p-value and its significance.

The returned object has an attribute called args, which is a list holding the test arguments.

See Also

binom_test

Examples

# Data
tulip <- c(red = 81, yellow = 50, white = 27)

# Question 1: are the color equally common ?
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# this is a test of homogeneity
res <- multinom_test(tulip)
res

attr(res, "descriptives")

# Pairwise comparisons between groups
pairwise_binom_test(tulip, p.adjust.method = "bonferroni")


# Question 2: comparing observed to expected proportions
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# this is a goodness-of-fit test
expected.p <- c(red = 0.5, yellow = 0.33, white = 0.17)
res <- multinom_test(tulip, expected.p)
res
attr(res, "descriptives")

# Pairwise comparisons against a given probabilities
pairwise_binom_test_against_p(tulip, expected.p)

Omega Squared for ANOVA

Description

Compute (classic, full) omega-squared and partial omega-squared for all terms in a between-subjects ANOVA model. Omega squared is a less-biased alternative to eta squared, estimating the population effect size rather than the sample one.

See the Datanovia tutorial One-Way ANOVA in R for a worked walkthrough.

Usage

omega_squared(model)

partial_omega_squared(model)

Arguments

Value

a named numeric vector of effect sizes, one per model term. A negative point estimate (which can arise for a term with F < 1) is floored at 0, since omega squared estimates a non-negative proportion of variance.

Functions

• omega_squared(): compute the classic (full) omega squared.

• partial_omega_squared(): compute partial omega squared.

References

Olejnik, S., & Algina, J. (2003). Generalized eta and omega squared statistics: Measures of effect size for some common research designs. Psychological Methods, 8(4), 434-447.

See Also

eta_squared(), anova_test(). The Datanovia tutorial: One-Way ANOVA in R.

Examples

# Data preparation
df <- ToothGrowth
df$dose <- as.factor(df$dose)

# Fit the model
res.aov <- aov(len ~ supp * dose, data = df)

# Effect size
omega_squared(res.aov)
partial_omega_squared(res.aov)

Rounding and Formatting p-values

Description

Round and format p-values. Can also mark significant p-values with stars.

Usage

p_round(x, ..., digits = 3)

p_format(
  x,
  ...,
  new.col = FALSE,
  digits = 2,
  accuracy = 1e-04,
  decimal.mark = ".",
  leading.zero = TRUE,
  trailing.zero = FALSE,
  add.p = FALSE,
  space = FALSE
)

p_mark_significant(
  x,
  ...,
  new.col = FALSE,
  cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 1),
  symbols = c("****", "***", "**", "*", "")
)

p_detect(data, type = c("all", "p", "p.adj"))

p_names()

p_adj_names()

Arguments

Value

a vector or a data frame containing the rounded/formatted p-values.

Functions

• p_round(): round p-values

• p_format(): format p-values. Add a symbol "<" for small p-values.

• p_mark_significant(): mark p-values with significance levels

• p_detect(): detects and returns p-value column names in a data frame.

• p_names(): returns known p-value column names

• p_adj_names(): returns known adjust p-value column names

Examples

# Round and format a vector of p-values
# ::::::::::::::::::::::::::::::::::::::::::::
# Format
p <- c(0.5678, 0.127, 0.045, 0.011, 0.009, 0.00002, NA)
p_format(p)

# Specify the accuracy
p_format(p, accuracy = 0.01)

# Add p and remove the leading zero
p_format(p, add.p = TRUE, leading.zero = FALSE)

# Remove space before and after "=" or "<".
p_format(p, add.p = TRUE, leading.zero = FALSE, space = FALSE)

# Mark significant p-values
# ::::::::::::::::::::::::::::::::::::::::::::
p_mark_significant(p)

# Round, the mark significant
p %>% p_round(digits = 2) %>% p_mark_significant()

# Format, then mark significant
p %>% p_format(digits = 2) %>% p_mark_significant()

# Perform stat test, format p and mark significant
# ::::::::::::::::::::::::::::::::::::::::::::
ToothGrowth %>%
  group_by(dose) %>%
  t_test(len ~ supp) %>%
  p_format(digits = 2, leading.zero = FALSE) %>%
  p_mark_significant()

Choose and Run the Appropriate Post-Hoc Test

Description

Given a one-way, independent-groups design (outcome ~ group), check the ANOVA assumptions and run the post-hoc test they imply, following the standard decision tree:

• each group normal and variances equal: Tukey HSD (tukey_hsd());

• each group normal but variances unequal: Games-Howell (games_howell_test());

• at least one group not normal: Dunn's test (dunn_test()).

Normality is assessed per group with the Shapiro-Wilk test applied to each group's values, routing on the smallest p-value across groups (a single non-normal group sends the data to the non-parametric test). This is deliberately not the pooled model residuals, which unequal variances would make non-normal and so hide the Games-Howell case. Homogeneity of variance is assessed with Levene's test (levene_test()). Both are judged at the significance level. The function returns the chosen test's usual pairwise result, with the selected method and the assumption verdicts attached (and shown when the result is printed), so the routing is transparent rather than hidden.

To check the same assumptions before the omnibus test and read off the recommended omnibus/post-hoc pair, use check_test_assumptions(); its result can be passed back here through .assumptions so the checks are not repeated. Choosing a test by first testing its assumptions on the same data has a known cost — see the note in ?check_test_assumptions.

See the Datanovia tutorial Statistical Tests and Assumptions in R for a worked walkthrough.

Usage

posthoc_test(
  data,
  formula,
  significance = 0.05,
  ...,
  .assumptions = NULL,
  omnibus = NULL
)

## S3 method for class 'posthoc_test'
print(x, ...)

Arguments

Value

the pairwise comparison table returned by the selected post-hoc test (a tukey_hsd, games_howell_test or dunn_test object), additionally classed posthoc_test. The selected method and the assumption verdicts are stored in the attributes "posthoc.method" and "assumptions", and printed above the table.

See Also

check_test_assumptions(), tukey_hsd(), games_howell_test(), dunn_test(), levene_test(), shapiro_test(). The Datanovia tutorial: Statistical Tests and Assumptions in R.

Examples

df <- ToothGrowth
df$dose <- as.factor(df$dose)

# Assumptions hold here, so Tukey HSD is chosen
df %>% posthoc_test(len ~ dose)

Proportion Test

Description

Performs proportion tests to either evaluate the homogeneity of proportions (probabilities of success) in several groups or to test that the proportions are equal to certain given values.

Wrappers around the R base function prop.test() but have the advantage of performing pairwise and row-wise z-test of two proportions, the post-hoc tests following a significant chi-square test of homogeneity for 2xc and rx2 contingency tables.

See the Datanovia tutorial Proportion Z-Test in R for a worked walkthrough.

Usage

prop_test(
  x,
  n,
  p = NULL,
  alternative = c("two.sided", "less", "greater"),
  correct = TRUE,
  conf.level = 0.95,
  detailed = FALSE
)

pairwise_prop_test(xtab, p.adjust.method = "holm", ...)

row_wise_prop_test(xtab, p.adjust.method = "holm", detailed = FALSE, ...)

Arguments

Value

return a data frame with some the following columns:

• n: the number of participants.

• group: the categories in the row-wise proportion tests.

• statistic: the value of Pearson's chi-squared test statistic.

• df: the degrees of freedom of the approximate chi-squared distribution of the test statistic.

• p: p-value.

• p.adj: the adjusted p-value.

• method: the used statistical test.

• p.signif, p.adj.signif: the significance level of p-values and adjusted p-values, respectively.

• estimate: a vector with the sample proportions x/n.

• estimate1, estimate2: the proportion in each of the two populations.

• alternative: a character string describing the alternative hypothesis.

• conf.low,conf.high: Lower and upper bound on a confidence interval. a confidence interval for the true proportion if there is one group, or for the difference in proportions if there are 2 groups and p is not given, or NULL otherwise. In the cases where it is not NULL, the returned confidence interval has an asymptotic confidence level as specified by conf.level, and is appropriate to the specified alternative hypothesis.

The returned object has an attribute called args, which is a list holding the test arguments.

Functions

• prop_test(): performs one-sample and two-samples z-test of proportions. Wrapper around the function prop.test().

• pairwise_prop_test(): pairwise comparisons between proportions, a post-hoc tests following a significant chi-square test of homogeneity for 2xc design. Wrapper around pairwise.prop.test()

• row_wise_prop_test(): performs row-wise z-test of two proportions, a post-hoc tests following a significant chi-square test of homogeneity for rx2 contingency table. The z-test of two proportions is calculated for each category (row).

See Also

The Datanovia tutorial: Proportion Z-Test in R.

Examples

# Comparing an observed proportion to an expected proportion
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
prop_test(x = 95, n = 160, p = 0.5, detailed = TRUE)

# Comparing two proportions
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Data: frequencies of smokers between two groups
xtab <- as.table(rbind(c(490, 10), c(400, 100)))
dimnames(xtab) <- list(
  group = c("grp1", "grp2"),
  smoker = c("yes", "no")
)
xtab
# compare the proportion of smokers
prop_test(xtab, detailed = TRUE)

# Homogeneity of proportions between groups
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# H0: the proportion of smokers is similar in the four groups
# Ha:  this proportion is different in at least one of the populations.
#
# Data preparation
grp.size <- c( 106, 113, 156, 102 )
smokers  <- c( 50, 100, 139, 80 )
no.smokers <- grp.size - smokers
xtab <- as.table(rbind(
  smokers,
  no.smokers
))
dimnames(xtab) <- list(
  Smokers = c("Yes", "No"),
  Groups = c("grp1", "grp2", "grp3", "grp4")
)
xtab

# Compare the proportions of smokers between groups
prop_test(xtab, detailed = TRUE)

# Pairwise comparison between groups
pairwise_prop_test(xtab)


# Pairwise proportion tests
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Data: Titanic
xtab <- as.table(rbind(
  c(122, 167, 528, 673),
  c(203, 118, 178, 212)
))
dimnames(xtab) <- list(
  Survived = c("No", "Yes"),
  Class = c("1st", "2nd", "3rd", "Crew")
)
xtab
# Compare the proportion of survived between groups
pairwise_prop_test(xtab)

# Row-wise proportion tests
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Data: Titanic
xtab <- as.table(rbind(
  c(180, 145), c(179, 106),
  c(510, 196), c(862, 23)
))
dimnames(xtab) <- list(
  Class = c("1st", "2nd", "3rd", "Crew"),
  Gender = c("Male", "Female")
)
xtab
# Compare the proportion of males and females in each category
row_wise_prop_test(xtab)

Test for Trend in Proportions

Description

Performs chi-squared test for trend in proportion. This test is also known as Cochran-Armitage trend test.

Wrappers around the R base function prop.trend.test() but returns a data frame for easy data visualization.

See the Datanovia tutorial Cochran-Armitage Trend Test in R for a worked walkthrough.

Usage

prop_trend_test(xtab, score = NULL)

Arguments

Value

return a data frame with some the following columns:

• n: the number of participants.

• statistic: the value of Chi-squared trend test statistic.

• df: the degrees of freedom.

• p: p-value.

• method: the used statistical test.

• p.signif: the significance level of p-values and adjusted p-values, respectively.

The returned object has an attribute called args, which is a list holding the test arguments.

See Also

The Datanovia tutorial: Cochran-Armitage Trend Test in R.

Examples

# Proportion of renal stone (calculi) across age
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Data
xtab <- as.table(rbind(
  c(384, 536, 335),
  c(951, 869, 438)
))
dimnames(xtab) <- list(
  stone = c("yes", "no"),
  age = c("30-39", "40-49", "50-59")
)
xtab
# Compare the proportion of survived between groups
prop_trend_test(xtab)

Pull Lower and Upper Triangular Part of a Matrix

Description

Returns the lower or the upper triangular part of a (correlation) matrix.

Usage

pull_triangle(x, triangle = c("lower", "upper"), diagonal = FALSE)

pull_upper_triangle(x, diagonal = FALSE)

pull_lower_triangle(x, diagonal = FALSE)

Arguments

Value

an object of class cor_mat_tri, which is a data frame

Functions

• pull_triangle(): returns either the lower or upper triangular part of a matrix.

• pull_upper_triangle(): returns an object of class upper_tri, which is a data frame containing the upper triangular part of a matrix.

• pull_lower_triangle(): returns an object of class lower_tri, which is a data frame containing the lower triangular part of a matrix.

See Also

replace_triangle()

Examples

# Data preparation
#::::::::::::::::::::::::::::::::::::::::::
mydata <- mtcars %>%
  select(mpg, disp, hp, drat, wt, qsec)
head(mydata, 3)

# Compute correlation matrix and pull triangles
#::::::::::::::::::::::::::::::::::::::::::
# Correlation matrix
cor.mat <- cor_mat(mydata)
cor.mat

# Pull lower triangular part
cor.mat %>% pull_lower_triangle()

# Pull upper triangular part
cor.mat %>% pull_upper_triangle()

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

car

Anova()

dplyr

desc(), filter(), group_by(), mutate(), select()

generics

augment(), glance(), tidy()

tibble

tibble()

tidyr

drop_na(), gather(), spread()

Remove Non-Significant from Statistical Tests

Description

Filter out non-significant (NS) p-values from a statistical test. Can detect automatically p-value columns

Usage

remove_ns(stat.test, col = NULL, signif.cutoff = 0.05)

Arguments

Value

a data frame

Examples

# Statistical test
stat.test <- PlantGrowth %>% wilcox_test(weight ~ group)
# Remove ns: automatic detection of p-value columns
stat.test %>% remove_ns()
# Remove ns by the column p
stat.test %>% remove_ns(col ="p")

Replace Lower and Upper Triangular Part of a Matrix

Description

Replace the lower or the upper triangular part of a (correlation) matrix.

Usage

replace_triangle(x, triangle = c("lower", "upper"), by = "", diagonal = FALSE)

replace_upper_triangle(x, by = "", diagonal = FALSE)

replace_lower_triangle(x, by = "", diagonal = FALSE)

Arguments

Value

an object of class cor_mat_tri, which is a data frame

Functions

• replace_triangle(): replaces the specified triangle by empty or NA.

• replace_upper_triangle(): replaces the upper triangular part of a matrix. Returns an object of class lower_tri.

• replace_lower_triangle(): replaces the lower triangular part of a matrix. Returns an object of class lower_tri

See Also

pull_triangle()

Examples

# Compute correlation matrix and pull triangles
#::::::::::::::::::::::::::::::::::::::::::
# Correlation matrix
 cor.mat <- mtcars %>%
   select(mpg, disp, hp, drat, wt, qsec) %>%
   cor_mat()
 cor.mat

 # Replace upper triangle by NA
#::::::::::::::::::::::::::::::::::::::::::
cor.mat %>% replace_upper_triangle(by = NA)


# Replace upper triangle by NA and reshape the
# correlation matrix to have unique combinations of variables
#::::::::::::::::::::::::::::::::::::::::::
cor.mat %>%
  replace_upper_triangle(by = NA) %>%
  cor_gather()

Programming with rstatix (tidy evaluation)

Description

How to use rstatix functions programmatically — i.e. when the variable names are held in character strings or passed into your own wrapper functions, as is common in package development, loops, and Shiny apps.

rstatix has two kinds of interfaces, and each supports a standard way of "programming over variables":

• Selection interface (functions that select columns through ... or vars=, e.g. cor_test(), get_summary_stats(), cor_mat(), freq_table()). These support full tidy-evaluation: bare names, the injection operators !!/!!!, the embracing operator {{ }} inside your own functions, a character vector via vars=, and tidyselect helpers such as all_of()/any_of().

• Formula interface (tests that take a formula, e.g. t_test(), wilcox_test(), kruskal_test(), anova_test(), friedman_test()). A formula is an ordinary R object, so build it from strings with reformulate() or stats::as.formula(paste(...)) and pass it in.

Details

Injecting a string directly into a raw formula (e.g. t_test(df, y ~ {{var}})) is not supported: a formula is captured as a syntax tree, not a quosure, so the embracing/injection operators do not apply there. Build the formula instead with reformulate(rhs, lhs) — see the examples.

In examples below, helpers that rstatix does not re-export are namespaced (rlang::sym, dplyr::all_of, dplyr::across); attach rlang/dplyr and you can drop the prefixes.

See Also

cor_test(), get_summary_stats(), t_test().

Examples

# Selection interface -----------------------------------------------------

# A variable name held in a string, injected with !!
x <- "mpg"; y <- "wt"
mtcars %>% cor_test(!!rlang::sym(x), !!rlang::sym(y))

# Several names at once, spliced with !!!
vars <- c("mpg", "disp", "hp")
mtcars %>% cor_test(!!!rlang::syms(vars))

# A character vector via `vars =` (no rlang needed)
mtcars %>% cor_test(vars = c("mpg", "disp"))

# tidyselect helpers
iris %>% get_summary_stats(dplyr::all_of(c("Sepal.Length", "Sepal.Width")))

# Your own wrapper function: embrace the argument with {{ }}
my_summary <- function(data, var) {
  data %>% get_summary_stats({{ var }}, type = "mean_sd")
}
my_summary(iris, Sepal.Length)

# Formula interface -------------------------------------------------------

# Build the formula from strings with reformulate(rhs, lhs)
outcome <- "len"; group <- "supp"
ToothGrowth %>% t_test(reformulate(group, outcome))

# The same inside a wrapper function
my_test <- function(data, outcome, group) {
  data %>% t_test(reformulate(group, outcome))
}
my_test(ToothGrowth, "len", "supp")

# Programmatic grouping + a built formula (a common end-to-end pattern)
gv <- "supp"
ToothGrowth %>%
  group_by(dplyr::across(dplyr::all_of(gv))) %>%
  t_test(reformulate("dose", "len"))

References and related packages

Description

Where the methods in rstatix come from, and which other packages implement them.

Most of the statistical methods in rstatix are implemented in base R from the published formulas. A few call another package to compute the result, and one is adapted from another package's code; Adapted code below says which. Where a function records the source of its method, it does so in its own References section; cite those authors, not this package, when you report a result.

Method sources

• Cramer's V — Cramer, H. (1946) Mathematical Methods of Statistics. See cramer_v().

• Effect-size confidence intervals (partial eta squared, Cramer's V) — obtained by inverting a noncentral distribution: Smithson, M. (2003) Confidence Intervals; Steiger, J. H. (2004) Beyond the F test. Psychological Methods, 9, 164-182. See anova_test() and cramer_v().

• Conover's all-pairs test — Conover, W. J. (1999) Practical Nonparametric Statistics, 3rd edition. See conover_test() and friedman_conover_test().

• Nemenyi's all-pairs test — Nemenyi, P. (1963) Distribution-free Multiple Comparisons. See friedman_nemenyi_test().

• Compact letter display — Piepho, H.-P. (2004) An algorithm for a letter-based representation of all-pairwise comparisons. Journal of Computational and Graphical Statistics, 13, 456-466. See add_cld().

• Dunnett's many-to-one comparisons — Dunnett, C. W. (1955) A multiple comparison procedure for comparing several treatments with a control. Journal of the American Statistical Association, 50, 1096-1121. See dunnett_test().

Related packages

The following packages implement some of the same methods. rstatix compares its results against them while developing, and they offer functionality that rstatix does not:

• effectsize — a broad effect-size toolkit; effectsize::cramers_v() and effectsize::eta_squared() produce the same intervals as cramer_v(ci = TRUE), anova_test(ci = ), and eta_squared(ci = ) / partial_eta_squared(ci = ) when called with alternative = "two.sided", which is not their default.

• DescTools — DescTools::CramerV() and DescTools::DunnettTest().

• multcomp — multcomp::glht() for general linear hypotheses, including Dunnett contrasts.

• multcompView — multcompView::multcompLetters() for compact letter displays.

• PMCMRplus — a large collection of all-pairs and many-to-one nonparametric procedures, including the Conover and Nemenyi tests.

These packages are not dependencies of rstatix. The values they return are recorded in the rstatix test suite as fixed numbers, checked against the package they came from at the time they were written; a recorded number cannot detect a later change in the package that produced it.

Adapted code

Of the functions documented in this package, sign_test() is the exception to the description above. Its one- and two-sample test code, and the confidence interval it reports for the median, are adapted with modifications from DescTools::SignTest() and DescTools::MedianCI(), written by Andri Signorell. DescTools is distributed under GPL (>= 2); rstatix is distributed under GPL-2.

No other source file in rstatix carries a statistical method adapted from a contributed package. Argument-validation idioms are shared with base R — check_two_samples_test_args() follows the preamble of stats::wilcox.test(), as DescTools and PMCMRplus do — and those are not methods anyone cites. The test suite scans the sources for the phrases that declare an adaptation and fails when one appears in a file this section does not name; code copied without a word about it would not be caught that way.

Copying code is not the same as calling it. Several functions compute their result by calling another package, as each of them states in its own documentation: anova_test() uses car::Anova() for type II and type III sums of squares, and stats::aov() for type I; dunnett_test() and emmeans_test() use emmeans::emmeans(); and wilcox_effsize() takes its default statistic (r = Z / sqrt(N)) from coin::wilcoxsign_test() or coin::wilcox_test(), while its method = "rank_biserial" alternative is computed in base R.

See Also

cramer_v(), anova_test(), conover_test(), add_cld(), dunnett_test().

Sample n Rows By Group From a Table

Description

sample n rows by group from a table using the sample_n() function.

Usage

sample_n_by(data, ..., size = 1, replace = FALSE)

Arguments

Examples

ToothGrowth %>% sample_n_by(dose, supp, size = 2)

Shapiro-Wilk Normality Test

Description

Provides a pipe-friendly framework to performs Shapiro-Wilk test of normality. Support grouped data and multiple variables for multivariate normality tests. Wrapper around the R base function shapiro.test(). Can handle grouped data.

See the Datanovia tutorial Normality Test in R for a worked walkthrough.

Usage

shapiro_test(data, ..., vars = NULL)

mshapiro_test(data)

Arguments

Value

a data frame containing the value of the Shapiro-Wilk statistic and the corresponding p.value.

Functions

• shapiro_test(): univariate Shapiro-Wilk normality test

• mshapiro_test(): multivariate Shapiro-Wilk normality test. This is a modified copy of the mshapiro.test() function of the package mvnormtest, for internal convenience.

See Also

The Datanovia tutorial: Normality Test in R.

Examples

# Shapiro Wilk normality test for one variable
iris %>% shapiro_test(Sepal.Length)

# Shapiro Wilk normality test for two variables
iris %>% shapiro_test(Sepal.Length, Petal.Width)

# Multivariate normality test
mshapiro_test(iris[, 1:3])

Sign Test

Description

Performs one-sample and two-sample sign tests.

See the Datanovia tutorial Sign Test in R for a worked walkthrough.

Usage

sign_test(
  data,
  formula,
  comparisons = NULL,
  ref.group = NULL,
  p.adjust.method = "holm",
  alternative = "two.sided",
  mu = 0,
  conf.level = 0.95,
  detailed = FALSE
)

pairwise_sign_test(
  data,
  formula,
  comparisons = NULL,
  ref.group = NULL,
  p.adjust.method = "holm",
  detailed = FALSE,
  ...
)

Arguments

Value

return a data frame with some the following columns:

• .y.: the y variable used in the test.

• group1,group2: the compared groups in the pairwise tests.

• n,n1,n2: Sample counts.

• statistic: Test statistic used to compute the p-value. That is the S-statistic (the number of positive differences between the data and the hypothesized median), with names attribute "S".

• df, parameter: degrees of freedom. Here, the total number of valid differences.

• p: p-value.

• method: the statistical test used to compare groups.

• p.signif, p.adj.signif: the significance level of p-values and adjusted p-values, respectively.

• estimate: estimate of the effect size. It corresponds to the median of the differences.

• alternative: a character string describing the alternative hypothesis.

• conf.low,conf.high: Lower and upper bound on a confidence interval of the estimate.

The returned object has an attribute called args, which is a list holding the test arguments.

Functions

• sign_test(): Sign test

• pairwise_sign_test(): performs pairwise two sample Wilcoxon test.

Note

The sign test and the confidence interval for the median are adapted, with modifications, from DescTools::SignTest() and DescTools::MedianCI(), written by Andri Signorell. The results match DescTools::SignTest().

Source

Adapted from DescTools::SignTest() and DescTools::MedianCI() (Andri Signorell), distributed under GPL (>= 2). See rstatix-references.

See Also

The Datanovia tutorial: Sign Test in R.

Examples

# Load data
#:::::::::::::::::::::::::::::::::::::::
data("ToothGrowth")
df <- ToothGrowth

# One-sample test
#:::::::::::::::::::::::::::::::::::::::::
df %>% sign_test(len ~ 1, mu = 0)


# Two-samples paired test
#:::::::::::::::::::::::::::::::::::::::::
df %>% sign_test(len ~ supp)


# Compare supp levels after grouping the data by "dose"
#::::::::::::::::::::::::::::::::::::::::
df %>%
  group_by(dose) %>%
  sign_test(data =., len ~ supp) %>%
  adjust_pvalue(method = "bonferroni") %>%
  add_significance("p.adj")

# pairwise comparisons
#::::::::::::::::::::::::::::::::::::::::
# As dose contains more than two levels ==>
# pairwise test is automatically performed.
df %>% sign_test(len ~ dose)

# Comparison against reference group
#::::::::::::::::::::::::::::::::::::::::
# each level is compared to the ref group
df %>% sign_test(len ~ dose, ref.group = "0.5")

T-test

Description

Provides a pipe-friendly framework to performs one and two sample t-tests.

See the Datanovia tutorial T-Test in R for a worked walkthrough.

Usage

t_test(
  data,
  formula,
  comparisons = NULL,
  ref.group = NULL,
  p.adjust.method = "holm",
  paired = FALSE,
  var.equal = FALSE,
  alternative = "two.sided",
  mu = 0,
  conf.level = 0.95,
  detailed = FALSE,
  id = NULL,
  error.as.na = FALSE,
  effect.size = FALSE
)

pairwise_t_test(
  data,
  formula,
  comparisons = NULL,
  ref.group = NULL,
  p.adjust.method = "holm",
  paired = FALSE,
  pool.sd = !paired,
  detailed = FALSE,
  ...,
  effect.size = FALSE
)

Arguments

Details

- If a list of comparisons is specified, the result of the pairwise tests is filtered to keep only the comparisons of interest. The p-value is adjusted after filtering.

- For a grouped data, if pairwise test is performed, then the p-values are adjusted for each group level independently.

Value

return a data frame with some the following columns:

• .y.: the y variable used in the test.

• group1,group2: the compared groups in the pairwise tests.

• n,n1,n2: Sample counts.

• statistic: Test statistic used to compute the p-value.

• df: degrees of freedom.

• p: p-value.

• p.adj: the adjusted p-value.

• method: the statistical test used to compare groups.

• p.signif, p.adj.signif: the significance level of p-values and adjusted p-values, respectively.

• estimate: estimate of the effect size. It corresponds to the estimated mean or difference in means depending on whether it was a one-sample test or a two-sample test. For a two-sample test the difference is taken as estimate1 - estimate2, i.e. mean(group1) - mean(group2) (following t.test).

• estimate1, estimate2: show the mean values of the two groups, respectively, for independent samples t-tests.

• alternative: a character string describing the alternative hypothesis.

• conf.low,conf.high: Lower and upper bound on a confidence interval.

The returned object has an attribute called args, which is a list holding the test arguments.

Functions

• t_test(): t test

• pairwise_t_test(): performs pairwise two sample t-test. Wrapper around the R base function pairwise.t.test.

Note

On the sign of statistic and estimate when a ref.group is specified: the reference group is taken as group1 and the other group as group2, and the difference is computed as estimate = mean(group1) - mean(group2) = mean(ref.group) - mean(other) (the t.test convention). A positive statistic/estimate therefore means the value is higher in the reference group. To orient results so that a positive sign means "higher in the non-reference group", flip the sign yourself, e.g. mutate(statistic = -statistic, estimate = -estimate).

See Also

rstatix-programming for building the formula from variable names held in strings (e.g. reformulate()). The Datanovia tutorial: T-Test in R.

Examples

# Load data
#:::::::::::::::::::::::::::::::::::::::
data("ToothGrowth")
df <- ToothGrowth

# One-sample test
#:::::::::::::::::::::::::::::::::::::::::
df %>% t_test(len ~ 1, mu = 0)


# Two-samples unpaired test
#:::::::::::::::::::::::::::::::::::::::::
df %>% t_test(len ~ supp)

# Two-samples paired test
#:::::::::::::::::::::::::::::::::::::::::
df %>% t_test (len ~ supp, paired = TRUE)

# Compare supp levels after grouping the data by "dose"
#::::::::::::::::::::::::::::::::::::::::
df %>%
  group_by(dose) %>%
  t_test(data =., len ~ supp) %>%
  adjust_pvalue(method = "bonferroni") %>%
  add_significance("p.adj")

# pairwise comparisons
#::::::::::::::::::::::::::::::::::::::::
# As dose contains more than two levels ==>
# pairwise test is automatically performed.
df %>% t_test(len ~ dose)

# Comparison against reference group
#::::::::::::::::::::::::::::::::::::::::
# each level is compared to the ref group
df %>% t_test(len ~ dose, ref.group = "0.5")

# Comparison against all
#::::::::::::::::::::::::::::::::::::::::
df %>% t_test(len ~ dose, ref.group = "all")

Tidy an rstatix Test Result

Description

tidy() and glance() methods for objects of class rstatix_test — the result of a test function such as t_test(), wilcox_test(), anova_test() or kruskal_test(). The results are already tidy tibbles; these methods drop the internal rstatix classes and the stashed test arguments so the object passes cleanly to tools that dispatch on tidy / glance, such as broom, gtsummary and gt. Correlation results (cor_test(), cor_mat()) carry a different class and are not covered by these methods.

Usage

## S3 method for class 'rstatix_test'
tidy(x, ...)

## S3 method for class 'rstatix_test'
glance(x, ...)

Arguments

Value

tidy() returns the same result as a plain tibble, one row per comparison or model term, with the internal classes and the args attribute removed. glance() returns a one-row tibble with the test method and n, the number of rows in the result (the number of comparisons or model terms).

Examples

res <- ToothGrowth %>% t_test(len ~ dose)

# A plain tibble, ready for broom / gtsummary / gt
tidy(res)

# One-row summary
glance(res)

Tukey Honest Significant Differences

Description

Provides a pipe-friendly framework to performs Tukey post-hoc tests. Wrapper around the function TukeyHSD(). It is essentially a t-test that corrects for multiple testing.

Can handle different inputs formats: aov, lm, formula.

See the Datanovia tutorial One-Way ANOVA in R for a worked walkthrough.

Usage

tukey_hsd(x, ...)

## Default S3 method:
tukey_hsd(x, ...)

## S3 method for class 'lm'
tukey_hsd(x, ...)

## S3 method for class 'data.frame'
tukey_hsd(x, formula, ...)

Arguments

Value

a tibble data frame containing the results of the different comparisons.

Methods (by class)

• tukey_hsd(default): performs tukey post-hoc test from aov() results.

• tukey_hsd(lm): performs tukey post-hoc test from lm() model.

• tukey_hsd(data.frame): performs tukey post-hoc tests using data and formula as inputs. ANOVA will be automatically performed using the function aov()

See Also

The Datanovia tutorial: One-Way ANOVA in R.

Examples

# Data preparation
df <- ToothGrowth
df$dose <- as.factor(df$dose)
# Tukey HSD from ANOVA results
aov(len ~ dose, data = df) %>% tukey_hsd()

# two-way anova with interaction
aov(len ~ dose*supp, data = df) %>% tukey_hsd()

# Tukey HSD from lm() results
lm(len ~ dose, data = df) %>% tukey_hsd()

# Tukey HSD from data frame and formula
tukey_hsd(df, len ~ dose)

# Tukey HSD using grouped data
df %>%
  group_by(supp) %>%
  tukey_hsd(len ~ dose)

Welch One-Way ANOVA Test

Description

Tests for equal means in a one-way design (not assuming equal variance). A wrapper around the base function oneway.test(). This is is an alternative to the standard one-way ANOVA in the situation where the homogeneity of variance assumption is violated.

See the Datanovia tutorial One-Way ANOVA in R for a worked walkthrough.

Usage

welch_anova_test(data, formula)

Arguments

Value

return a data frame with the following columns:

• .y.: the y variable used in the test.

• n: sample count.

• statistic: the value of the test statistic.

• p: p-value.

• method: the statistical test used to compare groups.

See Also

The Datanovia tutorial: One-Way ANOVA in R.

Examples

# Load data
#:::::::::::::::::::::::::::::::::::::::
data("ToothGrowth")
df <- ToothGrowth
df$dose <- as.factor(df$dose)

# Welch one-way ANOVA test (not assuming equal variance)
#:::::::::::::::::::::::::::::::::::::::::
df %>% welch_anova_test(len ~ dose)

# Grouped data
#:::::::::::::::::::::::::::::::::::::::::
df %>%
  group_by(supp) %>%
  welch_anova_test(len ~ dose)

Wilcoxon Effect Size

Description

Compute Wilcoxon effect size (r) for:

• one-sample test (Wilcoxon one-sample signed-rank test);

• paired two-samples test (Wilcoxon two-sample paired signed-rank test) and

• independent two-samples test ( Mann-Whitney, two-sample rank-sum test).

It can also returns confidence intervals by bootstap.

The effect size r is calculated as Z statistic divided by square root of the sample size (N) (Z/\sqrt{N}). The Z value is extracted from either coin::wilcoxsign_test() (case of one- or paired-samples test) or coin::wilcox_test() (case of independent two-samples test).

Here, N is the number of independent observations contributing to the test: the total sample size for the independent two-samples test, and the number of pairs (equivalently, the number of difference scores) for the one-sample and paired tests. This is because the paired test reduces to a one-sample signed-rank test on the pairwise differences, so each pair counts once. This convention matches the default of rcompanion::wilcoxonPairedR() (its cases = TRUE setting).

Some references instead define N as the total number of observations, i.e. twice the number of pairs (Field, 2012; Tomczak & Tomczak, 2014), which yields a smaller r. If you need that convention for a paired test, divide the reported r (or the Z) by \sqrt 2; it is also available via rcompanion::wilcoxonPairedR(..., cases = FALSE).

The r value varies from 0 to close to 1. The interpretation values for r commonly in published litterature and on the internet are: 0.10 - < 0.3 (small effect), 0.30 - < 0.5 (moderate effect) and >= 0.5 (large effect).

See the Datanovia tutorial Wilcoxon Test in R for a worked walkthrough.

Usage

wilcox_effsize(
  data,
  formula,
  comparisons = NULL,
  ref.group = NULL,
  paired = FALSE,
  alternative = "two.sided",
  mu = 0,
  ci = FALSE,
  conf.level = 0.95,
  ci.type = "perc",
  nboot = 1000,
  detailed = FALSE,
  ...,
  boot.parallel = getOption("boot.parallel", "no"),
  boot.ncpus = getOption("boot.ncpus", 1L),
  method = c("r", "rank_biserial")
)

Arguments

Value

return a data frame with some of the following columns:

• .y.: the y variable used in the test.

• group1,group2: the compared groups in the pairwise tests.

• n,n1,n2: Sample counts.

• effsize: estimate of the effect size (r value).

• magnitude: magnitude of effect size.

• conf.low,conf.high: lower and upper bound of the effect size confidence interval.

• statistic: the Z statistic and p: the p-value (only when detailed = TRUE).

References

Maciej Tomczak and Ewa Tomczak. The need to report effect size estimates revisited. An overview of some recommended measures of effect size. Trends in Sport Sciences. 2014; 1(21):19-25.

See Also

The Datanovia tutorial: Wilcoxon Test in R.

Examples

if(require("coin")){

# One-sample Wilcoxon test effect size
ToothGrowth %>% wilcox_effsize(len ~ 1, mu = 0)

# Independent two-samples wilcoxon effect size
ToothGrowth %>% wilcox_effsize(len ~ supp)


# Paired-samples wilcoxon effect size
ToothGrowth %>% wilcox_effsize(len ~ supp, paired = TRUE)

# Pairwise comparisons
ToothGrowth %>% wilcox_effsize(len ~ dose)

# Grouped data
ToothGrowth %>%
  group_by(supp) %>%
  wilcox_effsize(len ~ dose)

}

Wilcoxon Tests

Description

Provides a pipe-friendly framework to performs one and two sample Wilcoxon tests.

See the Datanovia tutorial Wilcoxon Test in R for a worked walkthrough.

Usage

wilcox_test(
  data,
  formula,
  comparisons = NULL,
  ref.group = NULL,
  p.adjust.method = "holm",
  paired = FALSE,
  exact = NULL,
  alternative = "two.sided",
  mu = 0,
  conf.level = 0.95,
  detailed = FALSE,
  id = NULL,
  error.as.na = FALSE,
  effect.size = FALSE
)

pairwise_wilcox_test(
  data,
  formula,
  comparisons = NULL,
  ref.group = NULL,
  p.adjust.method = "holm",
  detailed = FALSE,
  ...,
  effect.size = FALSE
)

Arguments

Details

- pairwise_wilcox_test() applies the standard two sample Wilcoxon test to all possible pairs of groups. This method calls the wilcox.test(), so extra arguments are accepted.

- If a list of comparisons is specified, the result of the pairwise tests is filtered to keep only the comparisons of interest.The p-value is adjusted after filtering.

- For a grouped data, if pairwise test is performed, then the p-values are adjusted for each group level independently.

- a nonparametric confidence interval and an estimator for the pseudomedian (one-sample case) or for the difference of the location parameters x-y is computed, where x and y are the compared samples or groups. The column estimate and the confidence intervals are displayed in the test result when the option detailed = TRUE is specified in the wilcox_test() and pairwise_wilcox_test() functions. Read more about the calculation of the estimate in the details section of the R base function wilcox.test() documentation by typing ?wilcox.test in the R console.

- With effect.size = TRUE, an independent-samples test is annotated with Cliff's delta and a paired test with the matched-pairs rank-biserial correlation, (R^+ - R^-)/(R^+ + R^-) over the signed ranks of the paired differences (Kerby, 2014).

Value

return a data frame with some of the following columns:

• .y.: the y variable used in the test.

• group1,group2: the compared groups in the pairwise tests.

• n,n1,n2: Sample counts.

• statistic: Test statistic used to compute the p-value.

• p: p-value.

• p.adj: the adjusted p-value.

• method: the statistical test used to compare groups.

• p.signif, p.adj.signif: the significance level of p-values and adjusted p-values, respectively.

• estimate: an estimate of the location parameter (Only present if argument detailed = TRUE). This corresponds to the pseudomedian (for one-sample case) or to the difference of the location parameter (for two-samples case).

  • The pseudomedian of a distribution F is the median of the distribution of (u+v)/2, where u and v are independent, each with distribution F. If F is symmetric, then the pseudomedian and median coincide.

  • Note that in the two-sample case the estimator for the difference in location parameters does not estimate the difference in medians (a common misconception) but rather the median of the difference between a sample from x and a sample from y.

• conf.low, conf.high: a confidence interval for the location parameter. (Only present if argument conf.int = TRUE.)

The returned object has an attribute called args, which is a list holding the test arguments.

Functions

• wilcox_test(): Wilcoxon test

• pairwise_wilcox_test(): performs pairwise two sample Wilcoxon test.

Note

When a ref.group is specified, the reference group is taken as group1 and the other group as group2, and the comparison is computed as group1 versus group2 (i.e. ref.group versus the other group), following the wilcox.test convention. With detailed = TRUE, the estimate is the Hodges-Lehmann location shift of group1 relative to group2, so a positive estimate means the reference group is shifted higher; flip its sign (mutate(estimate = -estimate)) if you want a positive sign to mean "higher in the non-reference group". (The statistic is the rank-sum/signed-rank W, which is not a signed difference.)

References

Kerby, D. S. (2014). The simple difference formula: An approach to teaching nonparametric correlation. Comprehensive Psychology, 3, 11.IT.3.1.

See Also

The Datanovia tutorial: Wilcoxon Test in R.

Examples

# Load data
#:::::::::::::::::::::::::::::::::::::::
data("ToothGrowth")
df <- ToothGrowth

# One-sample test
#:::::::::::::::::::::::::::::::::::::::::
df %>% wilcox_test(len ~ 1, mu = 0)


# Two-samples unpaired test
#:::::::::::::::::::::::::::::::::::::::::
df %>% wilcox_test(len ~ supp)

# Two-samples paired test
#:::::::::::::::::::::::::::::::::::::::::
df %>% wilcox_test (len ~ supp, paired = TRUE)

# Compare supp levels after grouping the data by "dose"
#::::::::::::::::::::::::::::::::::::::::
df %>%
  group_by(dose) %>%
  wilcox_test(data =., len ~ supp) %>%
  adjust_pvalue(method = "bonferroni") %>%
  add_significance("p.adj")

# pairwise comparisons
#::::::::::::::::::::::::::::::::::::::::
# As dose contains more than two levels ==>
# pairwise test is automatically performed.
df %>% wilcox_test(len ~ dose)

# Comparison against reference group
#::::::::::::::::::::::::::::::::::::::::
# each level is compared to the ref group
df %>% wilcox_test(len ~ dose, ref.group = "0.5")

# Comparison against all
#::::::::::::::::::::::::::::::::::::::::
df %>% wilcox_test(len ~ dose, ref.group = "all")