This vignette can be cited as:
To cite package 'statsExpressions' in publications use:
Patil, I., (2021). statsExpressions: R Package for Tidy Dataframes
and Expressions with Statistical Details. Journal of Open Source
Software, 6(61), 3236, https://doi.org/10.21105/joss.03236
A BibTeX entry for LaTeX users is
@Article{,
doi = {10.21105/joss.03236},
year = {2021},
publisher = {{The Open Journal}},
volume = {6},
number = {61},
pages = {3236},
author = {Indrajeet Patil},
title = {{statsExpressions: {R} Package for Tidy Dataframes and Expressions with Statistical Details}},
journal = {{Journal of Open Source Software}},
}One-sample tests
# for reproducibility
set.seed(123)
# ----------------------- parametric -----------------------
one_sample_test(mtcars, wt, test.value = 3)
#> # A tibble: 1 × 15
#> mu statistic df.error p.value method alternative effectsize
#> <dbl> <dbl> <dbl> <dbl> <chr> <chr> <chr>
#> 1 3 1.26 31 0.218 One Sample t-test two.sided Hedges' g
#> estimate conf.level conf.low conf.high conf.method conf.distribution n.obs
#> <dbl> <dbl> <dbl> <dbl> <chr> <chr> <int>
#> 1 0.217 0.95 -0.127 0.557 ncp t 32
#> expression
#> <list>
#> 1 <language>
# ----------------------- non-parametric -------------------
one_sample_test(mtcars, wt, test.value = 3, type = "nonparametric")
#> # A tibble: 1 × 12
#> statistic p.value method alternative effectsize
#> <dbl> <dbl> <chr> <chr> <chr>
#> 1 319 0.308 Wilcoxon signed rank test two.sided r (rank biserial)
#> estimate conf.level conf.low conf.high conf.method n.obs expression
#> <dbl> <dbl> <dbl> <dbl> <chr> <int> <list>
#> 1 0.208 0.95 -0.184 0.543 normal 32 <language>
# ----------------------- robust ---------------------------
one_sample_test(mtcars, wt, test.value = 3, type = "robust")
#> # A tibble: 1 × 10
#> statistic p.value n.obs method effectsize
#> <dbl> <dbl> <int> <chr> <chr>
#> 1 1.18 0.275 32 Bootstrap-t method for one-sample test Trimmed mean
#> estimate conf.level conf.low conf.high expression
#> <dbl> <dbl> <dbl> <dbl> <list>
#> 1 3.20 0.95 2.85 3.54 <language>
# ----------------------- Bayesian -------------------------
one_sample_test(mtcars, wt, test.value = 3, type = "bayes")
#> # A tibble: 1 × 16
#> term effectsize estimate conf.level conf.low conf.high pd
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 Difference Bayesian t-test 0.195 0.95 -0.165 0.555 0.86
#> prior.distribution prior.location prior.scale bf10 method
#> <chr> <dbl> <dbl> <dbl> <chr>
#> 1 cauchy 0 0.707 0.387 Bayesian t-test
#> conf.method log_e_bf10 n.obs expression
#> <chr> <dbl> <int> <list>
#> 1 ETI -0.950 32 <language>Two-sample tests
within-subjects design
# ----------------------- within-subjects -------------------------------------
# data
df <- dplyr::filter(bugs_long, condition %in% c("LDLF", "LDHF"))
# for reproducibility
set.seed(123)
# ----------------------- parametric ---------------------------------------
two_sample_test(df, condition, desire, subject.id = subject, paired = TRUE, type = "parametric")
#> # A tibble: 1 × 16
#> term group statistic df.error p.value method alternative
#> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr>
#> 1 desire condition 6.65 90 0.00000000222 Paired t-test two.sided
#> effectsize estimate conf.level conf.low conf.high conf.method
#> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
#> 1 Hedges' g 0.691 0.95 0.462 0.917 ncp
#> conf.distribution n.obs expression
#> <chr> <int> <list>
#> 1 t 91 <language>
# ----------------------- non-parametric -----------------------------------
two_sample_test(df, condition, desire, subject.id = subject, paired = TRUE, type = "nonparametric")
#> # A tibble: 1 × 14
#> parameter1 parameter2 statistic p.value method
#> <chr> <chr> <dbl> <dbl> <chr>
#> 1 desire condition 2250. 0.0000000241 Wilcoxon signed rank test
#> alternative effectsize estimate conf.level conf.low conf.high
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 two.sided r (rank biserial) 0.761 0.95 0.642 0.844
#> conf.method n.obs expression
#> <chr> <int> <list>
#> 1 normal 91 <language>
# ----------------------- robust --------------------------------------------
two_sample_test(df, condition, desire, subject.id = subject, paired = TRUE, type = "robust")
#> # A tibble: 1 × 15
#> statistic df.error p.value
#> <dbl> <dbl> <dbl>
#> 1 6.46 54 0.0000000313
#> method
#> <chr>
#> 1 Yuen's test on trimmed means for dependent samples
#> effectsize estimate conf.level
#> <chr> <dbl> <dbl>
#> 1 Algina-Keselman-Penfield robust standardized difference 0.533 0.95
#> conf.low conf.high mu small medium large n.obs expression
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <list>
#> 1 0.369 0.707 0 0.1 0.3 0.5 91 <language>
# ----------------------- Bayesian ---------------------------------------
two_sample_test(df, condition, desire, subject.id = subject, paired = TRUE, type = "bayes")
#> # A tibble: 1 × 16
#> term effectsize estimate conf.level conf.low conf.high pd
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 Difference Bayesian t-test 1.63 0.95 1.13 2.11 1
#> prior.distribution prior.location prior.scale bf10 method
#> <chr> <dbl> <dbl> <dbl> <chr>
#> 1 cauchy 0 0.707 4762370. Bayesian t-test
#> conf.method log_e_bf10 n.obs expression
#> <chr> <dbl> <int> <list>
#> 1 ETI 15.4 91 <language>between-subjects design
# ----------------------- between-subjects -------------------------------------
# for reproducibility
set.seed(123)
# ----------------------- parametric ---------------------------------------
# unequal variance
two_sample_test(ToothGrowth, supp, len, type = "parametric")
#> # A tibble: 1 × 18
#> parameter1 parameter2 mean.parameter1 mean.parameter2 statistic df.error
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 len supp 20.7 17.0 1.92 55.3
#> p.value method alternative effectsize estimate conf.level
#> <dbl> <chr> <chr> <chr> <dbl> <dbl>
#> 1 0.0606 Welch Two Sample t-test two.sided Hedges' g 0.488 0.95
#> conf.low conf.high conf.method conf.distribution n.obs expression
#> <dbl> <dbl> <chr> <chr> <int> <list>
#> 1 -0.0217 0.993 ncp t 60 <language>
# equal variance
two_sample_test(ToothGrowth, supp, len, type = "parametric", var.equal = TRUE)
#> # A tibble: 1 × 18
#> parameter1 parameter2 mean.parameter1 mean.parameter2 statistic df.error
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 len supp 20.7 17.0 1.92 58
#> p.value method alternative effectsize estimate conf.level conf.low
#> <dbl> <chr> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 0.0604 Two Sample t-test two.sided Hedges' g 0.488 0.95 -0.0217
#> conf.high conf.method conf.distribution n.obs expression
#> <dbl> <chr> <chr> <int> <list>
#> 1 0.993 ncp t 60 <language>
# ----------------------- non-parametric -----------------------------------
two_sample_test(ToothGrowth, supp, len, type = "nonparametric")
#> # A tibble: 1 × 14
#> parameter1 parameter2 statistic p.value method alternative
#> <chr> <chr> <dbl> <dbl> <chr> <chr>
#> 1 len supp 576. 0.0645 Wilcoxon rank sum test two.sided
#> effectsize estimate conf.level conf.low conf.high conf.method n.obs
#> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <int>
#> 1 r (rank biserial) 0.279 0.95 -0.00812 0.523 normal 60
#> expression
#> <list>
#> 1 <language>
# ----------------------- robust --------------------------------------------
two_sample_test(ToothGrowth, supp, len, type = "robust")
#> # A tibble: 1 × 11
#> statistic df.error p.value
#> <dbl> <dbl> <dbl>
#> 1 2.29 33.5 0.0286
#> method
#> <chr>
#> 1 Yuen's test on trimmed means for independent samples
#> effectsize estimate conf.level
#> <chr> <dbl> <dbl>
#> 1 Algina-Keselman-Penfield robust standardized difference 0.683 0.95
#> conf.low conf.high n.obs expression
#> <dbl> <dbl> <int> <list>
#> 1 -0.00736 2.36 60 <language>
# ----------------------- Bayesian ---------------------------------------
two_sample_test(ToothGrowth, supp, len, type = "bayes")
#> # A tibble: 1 × 16
#> term effectsize estimate conf.level conf.low conf.high pd
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 Difference Bayesian t-test 3.16 0.95 -0.338 6.78 0.961
#> prior.distribution prior.location prior.scale bf10 method
#> <chr> <dbl> <dbl> <dbl> <chr>
#> 1 cauchy 0 0.707 1.20 Bayesian t-test
#> conf.method log_e_bf10 n.obs expression
#> <chr> <dbl> <int> <list>
#> 1 ETI 0.181 60 <language>One-way ANOVAs
within-subjects design
suppressPackageStartupMessages(library(afex))
# ----------------------- parametric ---------------------------------------
set.seed(123)
oneway_anova(
data = bugs_long,
x = condition,
y = desire,
paired = TRUE,
subject.id = subject,
type = "p"
)
#> # A tibble: 1 × 18
#> term sumsq sum.squares.error df df.error meansq statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 condition 233. 984. 2.63 229. 4.30 20.6 8.27e-11
#> method effectsize estimate
#> <chr> <chr> <dbl>
#> 1 ANOVA estimation for factorial designs using 'afex' Omega2 (partial) 0.0783
#> conf.level conf.low conf.high conf.method conf.distribution n.obs expression
#> <dbl> <dbl> <dbl> <chr> <chr> <int> <list>
#> 1 0.95 0.0280 1 ncp F 88 <language>
# ----------------------- non-parametric -----------------------------------
set.seed(123)
oneway_anova(
data = bugs_long,
x = condition,
y = desire,
paired = TRUE,
subject.id = subject,
type = "np"
)
#> # A tibble: 1 × 15
#> parameter1 parameter2 statistic df.error p.value method
#> <chr> <chr> <dbl> <dbl> <dbl> <chr>
#> 1 desire condition 55.8 3 4.56e-12 Friedman rank sum test
#> effectsize estimate conf.level conf.low conf.high conf.method
#> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
#> 1 Kendall's W 0.211 0.95 0.148 1 percentile bootstrap
#> conf.iterations n.obs expression
#> <int> <int> <list>
#> 1 100 88 <language>
# ----------------------- robust --------------------------------------------
set.seed(123)
oneway_anova(
data = bugs_long,
x = condition,
y = desire,
paired = TRUE,
subject.id = subject,
type = "r"
)
#> # A tibble: 1 × 12
#> statistic df df.error p.value
#> <dbl> <dbl> <dbl> <dbl>
#> 1 21.0 2.73 145. 1.15e-10
#> method
#> <chr>
#> 1 A heteroscedastic one-way repeated measures ANOVA for trimmed means
#> effectsize estimate
#> <chr> <dbl>
#> 1 Algina-Keselman-Penfield robust standardized difference average 0.664
#> conf.level conf.low conf.high n.obs expression
#> <dbl> <dbl> <dbl> <int> <list>
#> 1 0.95 0.466 0.971 88 <language>
# ----------------------- Bayesian ---------------------------------------
set.seed(123)
oneway_anova(
data = bugs_long,
x = condition,
y = desire,
paired = TRUE,
subject.id = subject,
type = "bayes"
)
#> # A tibble: 8 × 19
#> term pd prior.distribution prior.location prior.scale effect
#> <chr> <dbl> <chr> <dbl> <dbl> <chr>
#> 1 mu 1 cauchy 0 0.707 fixed
#> 2 condition-HDHF 1 cauchy 0 0.707 fixed
#> 3 condition-HDLF 0.861 cauchy 0 0.707 fixed
#> 4 condition-LDHF 0.995 cauchy 0 0.707 fixed
#> 5 condition-LDLF 1 cauchy 0 0.707 fixed
#> 6 sig2 1 cauchy 0 1 fixed
#> 7 g_condition 1 cauchy 0 1 fixed
#> 8 g_.rowid 1 cauchy 0 1 fixed
#> bf10 method log_e_bf10 effectsize
#> <dbl> <chr> <dbl> <chr>
#> 1 1372773377. Bayes factors for linear models 21.0 Bayesian R-squared
#> 2 1372773377. Bayes factors for linear models 21.0 Bayesian R-squared
#> 3 1372773377. Bayes factors for linear models 21.0 Bayesian R-squared
#> 4 1372773377. Bayes factors for linear models 21.0 Bayesian R-squared
#> 5 1372773377. Bayes factors for linear models 21.0 Bayesian R-squared
#> 6 1372773377. Bayes factors for linear models 21.0 Bayesian R-squared
#> 7 1372773377. Bayes factors for linear models 21.0 Bayesian R-squared
#> 8 1372773377. Bayes factors for linear models 21.0 Bayesian R-squared
#> estimate std.dev conf.level conf.low conf.high conf.method component n.obs
#> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr> <int>
#> 1 0.529 0.0332 0.95 0.461 0.587 HDI conditional 88
#> 2 0.529 0.0332 0.95 0.461 0.587 HDI conditional 88
#> 3 0.529 0.0332 0.95 0.461 0.587 HDI conditional 88
#> 4 0.529 0.0332 0.95 0.461 0.587 HDI conditional 88
#> 5 0.529 0.0332 0.95 0.461 0.587 HDI conditional 88
#> 6 0.529 0.0332 0.95 0.461 0.587 HDI conditional 88
#> 7 0.529 0.0332 0.95 0.461 0.587 HDI conditional 88
#> 8 0.529 0.0332 0.95 0.461 0.587 HDI conditional 88
#> expression
#> <list>
#> 1 <language>
#> 2 <language>
#> 3 <language>
#> 4 <language>
#> 5 <language>
#> 6 <language>
#> 7 <language>
#> 8 <language>between-subjects design
# ----------------------- parametric ---------------------------------------
# unequal variance
set.seed(123)
oneway_anova(
data = iris,
x = Species,
y = Sepal.Length,
type = "p"
)
#> # A tibble: 1 × 14
#> statistic df df.error p.value
#> <dbl> <dbl> <dbl> <dbl>
#> 1 139. 2 92.2 1.51e-28
#> method effectsize estimate
#> <chr> <chr> <dbl>
#> 1 One-way analysis of means (not assuming equal variances) Omega2 0.743
#> conf.level conf.low conf.high conf.method conf.distribution n.obs expression
#> <dbl> <dbl> <dbl> <chr> <chr> <int> <list>
#> 1 0.95 0.671 1 ncp F 150 <language>
# equal variance
set.seed(123)
oneway_anova(
data = iris,
x = Species,
y = Sepal.Length,
var.equal = TRUE,
type = "p"
)
#> # A tibble: 1 × 14
#> statistic df df.error p.value method effectsize
#> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
#> 1 119. 2 147 1.67e-31 One-way analysis of means Omega2
#> estimate conf.level conf.low conf.high conf.method conf.distribution n.obs
#> <dbl> <dbl> <dbl> <dbl> <chr> <chr> <int>
#> 1 0.612 0.95 0.534 1 ncp F 150
#> expression
#> <list>
#> 1 <language>
# ----------------------- non-parametric -----------------------------------
set.seed(123)
oneway_anova(
data = iris,
x = Species,
y = Sepal.Length,
type = "np"
)
#> # A tibble: 1 × 15
#> parameter1 parameter2 statistic df.error p.value
#> <chr> <chr> <dbl> <int> <dbl>
#> 1 Sepal.Length Species 96.9 2 8.92e-22
#> method effectsize estimate conf.level conf.low
#> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 Kruskal-Wallis rank sum test Epsilon2 (rank) 0.651 0.95 0.595
#> conf.high conf.method conf.iterations n.obs expression
#> <dbl> <chr> <int> <int> <list>
#> 1 1 percentile bootstrap 100 150 <language>
# ----------------------- robust --------------------------------------------
set.seed(123)
oneway_anova(
data = iris,
x = Species,
y = Sepal.Length,
type = "r"
)
#> # A tibble: 1 × 12
#> statistic df df.error p.value
#> <dbl> <dbl> <dbl> <dbl>
#> 1 112. 2 53.8 0
#> method
#> <chr>
#> 1 A heteroscedastic one-way ANOVA for trimmed means
#> effectsize estimate conf.level conf.low conf.high
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 Explanatory measure of effect size 0.846 0.95 0.759 0.933
#> n.obs expression
#> <int> <list>
#> 1 150 <language>
# ----------------------- Bayesian ---------------------------------------
set.seed(123)
oneway_anova(
data = iris,
x = Species,
y = Sepal.Length,
type = "bayes"
)
#> # A tibble: 6 × 17
#> term pd prior.distribution prior.location prior.scale bf10
#> <chr> <dbl> <chr> <dbl> <dbl> <dbl>
#> 1 mu 1 cauchy 0 0.707 1.87e28
#> 2 Species-setosa 1 cauchy 0 0.707 1.87e28
#> 3 Species-versicolor 0.936 cauchy 0 0.707 1.87e28
#> 4 Species-virginica 1 cauchy 0 0.707 1.87e28
#> 5 sig2 1 cauchy 0 0.707 1.87e28
#> 6 g_Species 1 cauchy 0 0.707 1.87e28
#> method log_e_bf10 effectsize estimate std.dev
#> <chr> <dbl> <chr> <dbl> <dbl>
#> 1 Bayes factors for linear models 65.1 Bayesian R-squared 0.612 0.0311
#> 2 Bayes factors for linear models 65.1 Bayesian R-squared 0.612 0.0311
#> 3 Bayes factors for linear models 65.1 Bayesian R-squared 0.612 0.0311
#> 4 Bayes factors for linear models 65.1 Bayesian R-squared 0.612 0.0311
#> 5 Bayes factors for linear models 65.1 Bayesian R-squared 0.612 0.0311
#> 6 Bayes factors for linear models 65.1 Bayesian R-squared 0.612 0.0311
#> conf.level conf.low conf.high conf.method n.obs expression
#> <dbl> <dbl> <dbl> <chr> <int> <list>
#> 1 0.95 0.544 0.667 HDI 150 <language>
#> 2 0.95 0.544 0.667 HDI 150 <language>
#> 3 0.95 0.544 0.667 HDI 150 <language>
#> 4 0.95 0.544 0.667 HDI 150 <language>
#> 5 0.95 0.544 0.667 HDI 150 <language>
#> 6 0.95 0.544 0.667 HDI 150 <language>Contingency table analyses
if (identical(Sys.getenv("NOT_CRAN"), "true")) {
#### -------------------- association test ------------------------ ####
# ------------------------ frequentist ---------------------------------
# unpaired
set.seed(123)
contingency_table(
data = mtcars,
x = am,
y = vs,
paired = FALSE
)
# paired
paired_data <- tibble(
response_before = structure(c(1L, 2L, 1L, 2L), levels = c("no", "yes"), class = "factor"),
response_after = structure(c(1L, 1L, 2L, 2L), levels = c("no", "yes"), class = "factor"),
Freq = c(65L, 25L, 5L, 5L)
)
set.seed(123)
contingency_table(
data = paired_data,
x = response_before,
y = response_after,
paired = TRUE,
counts = Freq
)
# ------------------------ Bayesian -------------------------------------
# unpaired
set.seed(123)
contingency_table(
data = mtcars,
x = am,
y = vs,
paired = FALSE,
type = "bayes"
)
# paired
set.seed(123)
contingency_table(
data = paired_data,
x = response_before,
y = response_after,
paired = TRUE,
counts = Freq,
type = "bayes"
)
#### -------------------- goodness-of-fit test -------------------- ####
# ------------------------ frequentist ---------------------------------
set.seed(123)
contingency_table(
data = as.data.frame(HairEyeColor),
x = Eye,
counts = Freq
)
# ------------------------ Bayesian -------------------------------------
set.seed(123)
contingency_table(
data = as.data.frame(HairEyeColor),
x = Eye,
counts = Freq,
ratio = c(0.2, 0.2, 0.3, 0.3),
type = "bayes"
)
}
#> # A tibble: 1 × 4
#> bf10 prior.scale method expression
#> <dbl> <dbl> <chr> <list>
#> 1 4.17e55 1 Bayesian one-way contingency table analysis <language>Correlation analyses
# for reproducibility
set.seed(123)
# ----------------------- parametric -----------------------
corr_test(mtcars, wt, mpg, type = "parametric")
#> # A tibble: 1 × 14
#> parameter1 parameter2 effectsize estimate conf.level conf.low
#> <chr> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 wt mpg Pearson correlation -0.868 0.95 -0.934
#> conf.high statistic df.error p.value method n.obs conf.method
#> <dbl> <dbl> <int> <dbl> <chr> <int> <chr>
#> 1 -0.744 -9.56 30 1.29e-10 Pearson correlation 32 normal
#> expression
#> <list>
#> 1 <language>
# ----------------------- non-parametric -------------------
corr_test(mtcars, wt, mpg, type = "nonparametric")
#> # A tibble: 1 × 13
#> parameter1 parameter2 effectsize estimate conf.level conf.low
#> <chr> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 wt mpg Spearman correlation -0.886 0.95 -0.945
#> conf.high statistic p.value method n.obs conf.method expression
#> <dbl> <dbl> <dbl> <chr> <int> <chr> <list>
#> 1 -0.774 10292. 1.49e-11 Spearman correlation 32 normal <language>
# ----------------------- robust ---------------------------
corr_test(mtcars, wt, mpg, type = "robust")
#> # A tibble: 1 × 14
#> parameter1 parameter2 effectsize estimate conf.level
#> <chr> <chr> <chr> <dbl> <dbl>
#> 1 wt mpg Winsorized Pearson correlation -0.864 0.95
#> conf.low conf.high statistic df.error p.value method
#> <dbl> <dbl> <dbl> <int> <dbl> <chr>
#> 1 -0.932 -0.738 -9.41 30 1.84e-10 Winsorized Pearson correlation
#> n.obs conf.method expression
#> <int> <chr> <list>
#> 1 32 normal <language>
# ----------------------- Bayesian -------------------------
corr_test(mtcars, wt, mpg, type = "bayes")
#> # A tibble: 1 × 17
#> parameter1 parameter2 effectsize estimate conf.level
#> <chr> <chr> <chr> <dbl> <dbl>
#> 1 wt mpg Bayesian Pearson correlation -0.843 0.95
#> conf.low conf.high pd rope.percentage prior.distribution prior.location
#> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
#> 1 -0.934 -0.734 1 0 beta 1.41
#> prior.scale bf10 method n.obs conf.method
#> <dbl> <dbl> <chr> <int> <chr>
#> 1 1.41 56223033. Bayesian Pearson correlation 32 HDI
#> expression
#> <list>
#> 1 <language>Meta-analysis
library(metaplus)
# renaming columns to `{statsExpressions}` conventions
df <- dplyr::rename(mag, estimate = yi, std.error = sei)
# ----------------------- parametric ---------------------------------------
set.seed(123)
meta_analysis(df, type = "parametric")
#> # A tibble: 1 × 14
#> term effectsize estimate std.error conf.level conf.low
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 Overall meta-analytic summary estimate -0.767 0.212 0.95 -1.18
#> conf.high statistic p.value weight method conf.method
#> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
#> 1 -0.351 -3.62 0.000295 NA Meta-analysis using 'metafor' Wald
#> n.obs expression
#> <int> <list>
#> 1 16 <language>
# ----------------------- robust --------------------------------------------
set.seed(123)
meta_analysis(df, type = "robust")
#> # A tibble: 1 × 14
#> term effectsize estimate std.error conf.low conf.high
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 Overall meta-analytic summary estimate -0.746 0.233 -1.26 -0.344
#> statistic p.value weight conf.level method
#> <dbl> <dbl> <dbl> <dbl> <chr>
#> 1 -3.20 0.000501 NA 0.95 Robust meta-analysis using 'metaplus'
#> conf.method n.obs expression
#> <chr> <int> <list>
#> 1 Wald 16 <language>
# ----------------------- Bayesian ---------------------------------------
# suppress warnings about divergent transitions after warmup
set.seed(123)
suppressWarnings(meta_analysis(df, type = "bayes"))
#> # A tibble: 2 × 20
#> term effectsize estimate std.error conf.level
#> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 Overall meta-analytic posterior estimate -0.643 0.220 0.95
#> 2 tau meta-analytic posterior estimate 0.484 0.182 0.95
#> conf.low conf.high weight bf10 rhat ess component prior.distribution
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
#> 1 -1.11 -0.242 NA 53.0 1 3507 meta Student's t
#> 2 0.205 0.909 NA 53.0 1 3460. meta Inverse gamma
#> prior.location prior.scale method conf.method
#> <dbl> <dbl> <chr> <chr>
#> 1 0 0.707 Bayesian meta-analysis using 'metaBMA' ETI
#> 2 1 0.15 Bayesian meta-analysis using 'metaBMA' ETI
#> log_e_bf10 n.obs expression
#> <dbl> <int> <list>
#> 1 3.97 16 <language>
#> 2 3.97 16 <language>Centrality description
# for reproducibility
set.seed(123)
# ----------------------- parametric -----------------------
centrality_description(iris, Species, Sepal.Length, type = "parametric")
#> # A tibble: 3 × 12
#> Species Sepal.Length std.dev iqr min max skewness kurtosis n.obs
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
#> 1 setosa 5.01 0.352 0.400 4.3 5.8 0.120 -0.253 50
#> 2 versicolor 5.94 0.516 0.7 4.9 7 0.105 -0.533 50
#> 3 virginica 6.59 0.636 0.750 4.9 7.9 0.118 0.0329 50
#> missing.obs expression n.expression
#> <int> <list> <chr>
#> 1 0 <language> "setosa\n(n = 50)"
#> 2 0 <language> "versicolor\n(n = 50)"
#> 3 0 <language> "virginica\n(n = 50)"
# ----------------------- non-parametric -------------------
centrality_description(mtcars, am, wt, type = "nonparametric")
#> # A tibble: 2 × 12
#> am wt mad iqr min max skewness kurtosis n.obs missing.obs
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <int>
#> 1 0 3.52 0.452 0.41 2.46 5.42 1.15 1.06 19 0
#> 2 1 2.32 0.682 0.942 1.51 3.57 0.269 -0.654 13 0
#> expression n.expression
#> <list> <chr>
#> 1 <language> "0\n(n = 19)"
#> 2 <language> "1\n(n = 13)"
# ----------------------- robust ---------------------------
centrality_description(ToothGrowth, supp, len, type = "robust")
#> # A tibble: 2 × 12
#> supp len std.dev iqr min max skewness kurtosis n.obs missing.obs
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <int>
#> 1 OJ 21.7 6.61 10.9 8.2 30.9 -0.580 -0.831 30 0
#> 2 VC 16.6 8.27 12.5 4.2 33.9 0.306 -0.700 30 0
#> expression n.expression
#> <list> <chr>
#> 1 <language> "OJ\n(n = 30)"
#> 2 <language> "VC\n(n = 30)"
# ----------------------- Bayesian -------------------------
centrality_description(sleep, group, extra, type = "bayes")
#> # A tibble: 2 × 11
#> group extra iqr min max skewness kurtosis n.obs missing.obs expression
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <int> <list>
#> 1 1 0.0579 2.8 -1.6 3.7 0.581 -0.630 10 0 <language>
#> 2 2 0.973 3.82 -0.1 5.5 0.386 -1.42 10 0 <language>
#> n.expression
#> <chr>
#> 1 "1\n(n = 10)"
#> 2 "2\n(n = 10)"Pairwise comparisons for one-way design
Between-subjects design
# ----------------------- parametric -----------------------
# if `var.equal = TRUE`, then Student's *t*-test will be run
pairwise_comparisons(
data = ggplot2::msleep,
x = vore,
y = brainwt,
type = "parametric",
var.equal = TRUE,
paired = FALSE,
p.adjust.method = "bonferroni"
)
#> # A tibble: 6 × 6
#> group1 group2 p.value p.adjust.method test expression
#> <chr> <chr> <dbl> <chr> <chr> <list>
#> 1 carni herbi 1 Bonferroni Student's t <language>
#> 2 carni insecti 1 Bonferroni Student's t <language>
#> 3 carni omni 1 Bonferroni Student's t <language>
#> 4 herbi insecti 1 Bonferroni Student's t <language>
#> 5 herbi omni 0.979 Bonferroni Student's t <language>
#> 6 insecti omni 1 Bonferroni Student's t <language>
# if `var.equal = FALSE`, then Games-Howell test will be run
pairwise_comparisons(
data = ggplot2::msleep,
x = vore,
y = brainwt,
type = "parametric",
var.equal = FALSE,
paired = FALSE,
p.adjust.method = "bonferroni"
)
#> # A tibble: 6 × 9
#> group1 group2 statistic p.value alternative distribution p.adjust.method
#> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
#> 1 carni herbi 2.17 1 two.sided q Bonferroni
#> 2 carni insecti -2.17 1 two.sided q Bonferroni
#> 3 carni omni 1.10 1 two.sided q Bonferroni
#> 4 herbi insecti -2.41 1 two.sided q Bonferroni
#> 5 herbi omni -1.87 1 two.sided q Bonferroni
#> 6 insecti omni 2.19 1 two.sided q Bonferroni
#> test expression
#> <chr> <list>
#> 1 Games-Howell <language>
#> 2 Games-Howell <language>
#> 3 Games-Howell <language>
#> 4 Games-Howell <language>
#> 5 Games-Howell <language>
#> 6 Games-Howell <language>
# ----------------------- non-parametric -------------------
pairwise_comparisons(
data = ggplot2::msleep,
x = vore,
y = brainwt,
type = "nonparametric",
paired = FALSE,
p.adjust.method = "none"
)
#> # A tibble: 6 × 9
#> group1 group2 statistic p.value alternative distribution p.adjust.method
#> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
#> 1 carni herbi 0.582 0.561 two.sided z None
#> 2 carni insecti 1.88 0.0595 two.sided z None
#> 3 carni omni 1.14 0.254 two.sided z None
#> 4 herbi insecti 1.63 0.102 two.sided z None
#> 5 herbi omni 0.717 0.474 two.sided z None
#> 6 insecti omni 1.14 0.254 two.sided z None
#> test expression
#> <chr> <list>
#> 1 Dunn <language>
#> 2 Dunn <language>
#> 3 Dunn <language>
#> 4 Dunn <language>
#> 5 Dunn <language>
#> 6 Dunn <language>
# ----------------------- robust ---------------------------
pairwise_comparisons(
data = ggplot2::msleep,
x = vore,
y = brainwt,
type = "robust",
paired = FALSE,
p.adjust.method = "fdr"
)
#> # A tibble: 6 × 10
#> group1 group2 estimate conf.level conf.low conf.high p.value p.adjust.method
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
#> 1 carni herbi -0.0323 0.95 -0.248 0.184 0.790 FDR
#> 2 carni insecti 0.0451 0.95 -0.0484 0.139 0.552 FDR
#> 3 carni omni 0.00520 0.95 -0.114 0.124 0.898 FDR
#> 4 herbi insecti 0.0774 0.95 -0.133 0.288 0.552 FDR
#> 5 herbi omni 0.0375 0.95 -0.182 0.257 0.790 FDR
#> 6 insecti omni -0.0399 0.95 -0.142 0.0625 0.552 FDR
#> test expression
#> <chr> <list>
#> 1 Yuen's trimmed means <language>
#> 2 Yuen's trimmed means <language>
#> 3 Yuen's trimmed means <language>
#> 4 Yuen's trimmed means <language>
#> 5 Yuen's trimmed means <language>
#> 6 Yuen's trimmed means <language>
# ----------------------- Bayesian -------------------------
pairwise_comparisons(
data = ggplot2::msleep,
x = vore,
y = brainwt,
type = "bayes",
paired = FALSE
)
#> # A tibble: 6 × 18
#> group1 group2 term effectsize estimate conf.level conf.low
#> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 carni herbi Difference Bayesian t-test -0.376 0.95 -1.36
#> 2 carni insecti Difference Bayesian t-test 0.0339 0.95 -0.0440
#> 3 carni omni Difference Bayesian t-test -0.0447 0.95 -0.241
#> 4 herbi insecti Difference Bayesian t-test 0.353 0.95 -0.733
#> 5 herbi omni Difference Bayesian t-test 0.364 0.95 -0.313
#> 6 insecti omni Difference Bayesian t-test -0.0766 0.95 -0.337
#> conf.high pd prior.distribution prior.location prior.scale bf10
#> <dbl> <dbl> <chr> <dbl> <dbl> <dbl>
#> 1 0.508 0.800 cauchy 0 0.707 0.540
#> 2 0.127 0.812 cauchy 0 0.707 0.718
#> 3 0.145 0.688 cauchy 0 0.707 0.427
#> 4 1.62 0.742 cauchy 0 0.707 0.540
#> 5 1.10 0.848 cauchy 0 0.707 0.571
#> 6 0.164 0.743 cauchy 0 0.707 0.545
#> conf.method log_e_bf10 n.obs expression test
#> <chr> <dbl> <int> <list> <chr>
#> 1 ETI -0.617 29 <language> Student's t
#> 2 ETI -0.332 14 <language> Student's t
#> 3 ETI -0.851 26 <language> Student's t
#> 4 ETI -0.616 25 <language> Student's t
#> 5 ETI -0.560 37 <language> Student's t
#> 6 ETI -0.606 22 <language> Student's tWithin-subjects design
# ----------------------- parametric -----------------------
pairwise_comparisons(
data = bugs_long,
x = condition,
y = desire,
subject.id = subject,
type = "parametric",
paired = TRUE,
p.adjust.method = "BH"
)
#> # A tibble: 6 × 6
#> group1 group2 p.value p.adjust.method test expression
#> <chr> <chr> <dbl> <chr> <chr> <list>
#> 1 HDHF HDLF 1.06e- 3 FDR Student's t <language>
#> 2 HDHF LDHF 7.02e- 2 FDR Student's t <language>
#> 3 HDHF LDLF 3.95e-12 FDR Student's t <language>
#> 4 HDLF LDHF 6.74e- 2 FDR Student's t <language>
#> 5 HDLF LDLF 1.99e- 3 FDR Student's t <language>
#> 6 LDHF LDLF 6.66e- 9 FDR Student's t <language>
# ----------------------- non-parametric -------------------
pairwise_comparisons(
data = bugs_long,
x = condition,
y = desire,
subject.id = subject,
type = "nonparametric",
paired = TRUE,
p.adjust.method = "BY"
)
#> # A tibble: 6 × 9
#> group1 group2 statistic p.value alternative distribution p.adjust.method
#> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
#> 1 HDHF HDLF 4.78 1.44e- 5 two.sided t BY
#> 2 HDHF LDHF 2.44 4.47e- 2 two.sided t BY
#> 3 HDHF LDLF 8.01 5.45e-13 two.sided t BY
#> 4 HDLF LDHF 2.34 4.96e- 2 two.sided t BY
#> 5 HDLF LDLF 3.23 5.05e- 3 two.sided t BY
#> 6 LDHF LDLF 5.57 4.64e- 7 two.sided t BY
#> test expression
#> <chr> <list>
#> 1 Durbin-Conover <language>
#> 2 Durbin-Conover <language>
#> 3 Durbin-Conover <language>
#> 4 Durbin-Conover <language>
#> 5 Durbin-Conover <language>
#> 6 Durbin-Conover <language>
# ----------------------- robust ---------------------------
pairwise_comparisons(
data = bugs_long,
x = condition,
y = desire,
subject.id = subject,
type = "robust",
paired = TRUE,
p.adjust.method = "hommel"
)
#> # A tibble: 6 × 11
#> group1 group2 estimate conf.level conf.low conf.high p.value p.crit
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 HDHF HDLF 1.03 0.95 0.140 1.92 0.00999 0.0127
#> 2 HDHF LDHF 0.454 0.95 -0.104 1.01 0.0520 0.025
#> 3 HDHF LDLF 1.95 0.95 1.09 2.82 0.000000564 0.00851
#> 4 HDLF LDHF -0.676 0.95 -1.61 0.256 0.0520 0.05
#> 5 HDLF LDLF 0.889 0.95 0.0244 1.75 0.0203 0.0169
#> 6 LDHF LDLF 1.35 0.95 0.560 2.14 0.000102 0.0102
#> p.adjust.method test expression
#> <chr> <chr> <list>
#> 1 Hommel Yuen's trimmed means <language>
#> 2 Hommel Yuen's trimmed means <language>
#> 3 Hommel Yuen's trimmed means <language>
#> 4 Hommel Yuen's trimmed means <language>
#> 5 Hommel Yuen's trimmed means <language>
#> 6 Hommel Yuen's trimmed means <language>
# ----------------------- Bayesian -------------------------
pairwise_comparisons(
data = bugs_long,
x = condition,
y = desire,
subject.id = subject,
type = "bayes",
paired = TRUE,
bf.prior = 0.77
)
#> # A tibble: 6 × 18
#> group1 group2 term effectsize estimate conf.level conf.low
#> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 HDHF HDLF Difference Bayesian t-test 1.13 0.95 0.476
#> 2 HDHF LDHF Difference Bayesian t-test 0.458 0.95 -0.0405
#> 3 HDHF LDLF Difference Bayesian t-test 2.14 0.95 1.63
#> 4 HDLF LDHF Difference Bayesian t-test -0.656 0.95 -1.31
#> 5 HDLF LDLF Difference Bayesian t-test 0.980 0.95 0.385
#> 6 LDHF LDLF Difference Bayesian t-test 1.66 0.95 1.15
#> conf.high pd prior.distribution prior.location prior.scale bf10
#> <dbl> <dbl> <chr> <dbl> <dbl> <dbl>
#> 1 1.73 1.00 cauchy 0 0.77 3.95e+ 1
#> 2 0.953 0.965 cauchy 0 0.77 5.42e- 1
#> 3 2.65 1 cauchy 0 0.77 1.22e+10
#> 4 0.0590 0.966 cauchy 0 0.77 6.50e- 1
#> 5 1.57 1.00 cauchy 0 0.77 1.72e+ 1
#> 6 2.17 1 cauchy 0 0.77 4.78e+ 6
#> conf.method log_e_bf10 n.obs expression test
#> <chr> <dbl> <int> <list> <chr>
#> 1 ETI 3.68 88 <language> Student's t
#> 2 ETI -0.612 88 <language> Student's t
#> 3 ETI 23.2 88 <language> Student's t
#> 4 ETI -0.430 88 <language> Student's t
#> 5 ETI 2.84 88 <language> Student's t
#> 6 ETI 15.4 88 <language> Student's tSuggestions
If you find any bugs or have any suggestions/remarks, please file an issue on GitHub: https://github.com/IndrajeetPatil/statsExpressions/issues