library(tidyverse)
library(ggstatsplot)
#> Registered S3 methods overwritten by 'broom.mixed':
#> method from
#> augment.lme broom
#> augment.merMod broom
#> glance.lme broom
#> glance.merMod broom
#> glance.stanreg broom
#> tidy.brmsfit broom
#> tidy.gamlss broom
#> tidy.lme broom
#> tidy.merMod broom
#> tidy.rjags broom
#> tidy.stanfit broom
#> tidy.stanreg broom
#> Registered S3 methods overwritten by 'lme4':
#> method from
#> cooks.distance.influence.merMod car
#> influence.merMod car
#> dfbeta.influence.merMod car
#> dfbetas.influence.merMod car
library(BayesFactor)
#> Loading required package: coda
#> Loading required package: Matrix
#>
#> Attaching package: 'Matrix'
#> The following object is masked from 'package:tidyr':
#>
#> expand
#> ************
#> Welcome to BayesFactor 0.9.12-4.2. If you have questions, please contact Richard Morey ([email protected]).
#>
#> Type BFManual() to open the manual.
#> ************
# function body
xpairwise_p <- function(data,
x,
y,
type = "parametric",
tr = 0.1,
paired = FALSE,
var.equal = FALSE,
p.adjust.method = "holm",
k = 2,
messages = TRUE,
...) {
ellipsis::check_dots_used()
# ---------------------------- data cleanup -------------------------------
# creating a dataframe
data <-
dplyr::select(
.data = data,
x = !!rlang::enquo(x),
y = !!rlang::enquo(y)
) %>%
dplyr::mutate(.data = ., x = droplevels(as.factor(x))) %>%
tibble::as_tibble(x = .)
# return(data)
# ---------------------------- parametric ---------------------------------
if (type %in% c("parametric", "p")) {
if (isTRUE(var.equal) || isTRUE(paired)) {
# anova model
aovmodel <- stats::aov(formula = y ~ x, data = data)
# safeguarding against edge cases
aovmodel$model %<>%
dplyr::mutate(
.data = .,
x = forcats::fct_relabel(
.f = x,
.fun = ~ stringr::str_replace(
string = .x,
pattern = "-",
replacement = "_"
)
)
)
# extracting and cleaning up Tukey's HSD output
df_tukey <- stats::TukeyHSD(x = aovmodel, conf.level = 0.95) %>%
broomExtra::tidy(x = .) %>%
dplyr::select(.data = ., comparison, estimate) %>%
tidyr::separate(
data = .,
col = comparison,
into = c("group1", "group2"),
sep = "-"
) %>%
dplyr::rename(.data = ., mean.difference = estimate) %>%
dplyr::mutate_at(
.tbl = .,
.vars = dplyr::vars(dplyr::matches("^group[0-9]$")),
.funs = ~ stringr::str_replace(
string = .,
pattern = "_",
replacement = "-"
)
)
# tidy dataframe with results from pairwise tests
df_tidy <- broomExtra::tidy(
stats::pairwise.t.test(
x = data$y,
g = data$x,
p.adjust.method = p.adjust.method,
paired = paired,
alternative = "two.sided",
na.action = na.omit
)
) %>%
signif_column(data = ., p = p.value)
# combining mean difference and results from pairwise t-test
df <-
dplyr::full_join(
x = df_tukey,
y = df_tidy,
by = c("group1", "group2")
) %>% # the group columns need to be swapped to be consistent
dplyr::rename(.data = ., group2 = group1, group1 = group2) %>%
dplyr::select(.data = ., group1, group2, dplyr::everything())
# display message about the post hoc tests run
if (isTRUE(messages)) {
message(cat(
crayon::green("Note: "),
crayon::blue(
"The parametric pairwise multiple comparisons test used-\n",
"Student's t-test.\n",
"Adjustment method for p-values: "
),
crayon::yellow(p.adjust.method),
sep = ""
))
}
} else {
# dataframe with Games-Howell test results
df <-
games_howell(data = data, x = x, y = y) %>%
p_adjust_column_adder(df = ., p.adjust.method = p.adjust.method) %>%
dplyr::select(.data = ., -conf.low, -conf.high)
# display message about the post hoc tests run
if (isTRUE(messages)) {
message(cat(
crayon::green("Note: "),
crayon::blue(
"The parametric pairwise multiple comparisons test used-\n",
"Games-Howell test.\n",
"Adjustment method for p-values: "
),
crayon::yellow(p.adjust.method),
sep = ""
))
}
}
}
# ---------------------------- nonparametric ----------------------------
if (type %in% c("nonparametric", "np")) {
if (!isTRUE(paired)) {
# running Dwass-Steel-Crichtlow-Fligner test using `jmv` package
jmv_pairs <-
jmv::anovaNP(
data = data,
deps = "y",
group = "x",
pairs = TRUE
)
# extracting the pairwise tests and formatting the output
df <-
as.data.frame(x = jmv_pairs$comparisons[[1]]) %>%
tibble::as_tibble(x = .) %>%
dplyr::rename(
.data = .,
group1 = p1,
group2 = p2,
p.value = p
) %>%
ggstatsplot:::p_adjust_column_adder(df = ., p.adjust.method = p.adjust.method)
# letting the user know which test was run
if (isTRUE(messages)) {
message(cat(
crayon::green("Note: "),
crayon::blue(
"The nonparametric pairwise multiple comparisons test used-\n",
"Dwass-Steel-Crichtlow-Fligner test.\n",
"Adjustment method for p-values: "
),
crayon::yellow(p.adjust.method),
sep = ""
))
}
} else {
# converting the entered long format data to wide format
data_wide <- long_to_wide_converter(
data = data,
x = x,
y = y
)
# running Durbin-Conover test using `jmv` package
jmv_pairs <-
jmv::anovaRMNP(
data = data_wide,
measures = names(data_wide[, -1]),
pairs = TRUE
)
# extracting the pairwise tests and formatting the output
df <-
as.data.frame(x = jmv_pairs$comp) %>%
tibble::as_tibble(x = .) %>%
dplyr::select(.data = ., -sep) %>%
dplyr::rename(
.data = .,
group1 = i1,
group2 = i2,
statistic = stat,
p.value = p
) %>%
ggstatsplot:::p_adjust_column_adder(df = ., p.adjust.method = p.adjust.method)
# letting the user know which test was run
if (isTRUE(messages)) {
message(cat(
crayon::green("Note: "),
crayon::blue(
"The nonparametric pairwise multiple comparisons test used-\n",
"Durbin-Conover test.\n",
"Adjustment method for p-values: "
),
crayon::yellow(p.adjust.method),
sep = ""
))
}
}
}
# ---------------------------- robust ----------------------------------
if (type %in% c("robust", "r")) {
if (!isTRUE(paired)) {
# object with all details about pairwise comparisons
rob_pairwise_df <-
WRS2::lincon(
formula = y ~ x,
data = data,
tr = tr
)
} else {
# converting to long format and then getting it back in wide so that the
# rowid variable can be used as the block variable for WRS2 functions
data_within <-
long_to_wide_converter(
data = data,
x = x,
y = y
) %>%
tidyr::gather(data = ., key, value, -rowid) %>%
dplyr::arrange(.data = ., rowid)
# running pairwise multiple comparison tests
rob_pairwise_df <-
with(
data = data_within,
expr = WRS2::rmmcp(
y = value,
groups = key,
blocks = rowid,
tr = tr
)
)
}
# extracting the robust pairwise comparisons and tidying up names
rob_df_tidy <-
suppressMessages(tibble::as_tibble(
x = rob_pairwise_df$comp,
.name_repair = "unique"
)) %>%
dplyr::rename(
.data = .,
group1 = Group...1,
group2 = Group...2
)
# cleaning the raw object and getting it in the right format
df <-
dplyr::full_join(
# dataframe comparing comparison details
x = rob_df_tidy %>%
ggstatsplot:::p_adjust_column_adder(df = ., p.adjust.method = p.adjust.method) %>%
tidyr::gather(
data = .,
key = "key",
value = "rowid",
group1:group2
),
# dataframe with factor level codings
y = rob_pairwise_df$fnames %>%
tibble::enframe(x = ., name = "rowid"),
by = "rowid"
) %>%
dplyr::select(.data = ., -rowid) %>%
tidyr::spread(data = ., key = "key", value = "value") %>%
dplyr::select(.data = ., group1, group2, dplyr::everything())
# for paired designs, there will be an unnecessary column to remove
if (("p.crit") %in% names(df)) {
df %<>% dplyr::select(.data = ., -p.crit)
}
# renaming confidence interval names
df %<>% dplyr::rename(.data = ., conf.low = ci.lower, conf.high = ci.upper)
# message about which test was run
if (isTRUE(messages)) {
message(cat(
crayon::green("Note: "),
crayon::blue(
"The robust pairwise multiple comparisons test used-\n",
"Yuen's trimmed means comparisons test.\n",
"Adjustment method for p-values: "
),
crayon::yellow(p.adjust.method),
sep = ""
))
}
}
# ---------------------------- bayes factor --------------------------------
# print a message telling the user that this is currently not supported
if (type %in% c("bf", "bayes")) {
# anova model
aovmodel <- stats::aov(formula = y ~ x, data = data)
# safeguarding against edge cases
aovmodel$model %<>%
dplyr::mutate(
.data = .,
x = forcats::fct_relabel(
.f = x,
.fun = ~ stringr::str_replace(
string = .x,
pattern = "-",
replacement = "_"
)
)
)
# extracting and cleaning up Tukey's HSD output
df_tukey <- stats::TukeyHSD(x = aovmodel, conf.level = 0.95) %>%
broomExtra::tidy(x = .) %>%
dplyr::select(.data = ., comparison, estimate) %>%
tidyr::separate(
data = .,
col = comparison,
into = c("group1", "group2"),
sep = "-"
) %>%
dplyr::rename(.data = ., mean.difference = estimate) %>%
dplyr::mutate_at(
.tbl = .,
.vars = dplyr::vars(dplyr::matches("^group[0-9]$")),
.funs = ~ stringr::str_replace(
string = .,
pattern = "_",
replacement = "-"
)
)
# return(df_tukey)
# combining mean difference and results from pairwise t-test
df <- df_tukey %>% # the group columns need to be swapped to be consistent
dplyr::rename(.data = ., group2 = group1, group1 = group2) %>%
dplyr::select(.data = ., group1, group2, dplyr::everything())
g1_list <- df_tukey %>% pull(group2) %>% as.character()
g2_list <- df_tukey %>% pull(group1) %>% as.character()
# return(g1_list)
bfresults <- map2(
g1_list,
g2_list,
function(a, b)
data %>%
filter(!is.na(x)) %>%
filter(!is.na(y)) %>%
filter(x %in% c(a, b)) %>%
droplevels() %>%
as.data.frame()
) %>%
map(.x = ., ~ ttestBF(
formula = y ~ x,
data = .
)) %>%
map(.x = ., ~ extractBF(x = .)) %>%
map_dbl(.x = ., ~ .[, "bf"])
# return(bfresults)
df$bfactor <- bfresults
# display message about the post hoc tests run
if (isTRUE(messages)) {
message(cat(
crayon::green("Note: "),
crayon::blue(
"The pairwise multiple comparisons test used - \n",
"BayesFactor::ttestBF.\n"
),
sep = ""
))
}
} # bayes
# ---------------------------- cleanup ----------------------------------
# if there are factors, covert them to character to make life easy
df %<>%
dplyr::mutate_if(
.tbl = .,
.predicate = is.factor,
.funs = ~ as.character(.)
)
# return(df)
if (type %in% c("parametric", "nonparametric", "robust", "p", "np", "r")) {
df %<>%
purrrlyr::by_row(
.d = .,
..f = ~ specify_decimal_p(
x = .$p.value,
k = k,
p.value = TRUE
),
.collate = "rows",
.to = "label",
.labels = TRUE
) %>%
dplyr::mutate(
.data = .,
p.value.label = dplyr::case_when(
label == "< 0.001" ~ "p <= 0.001",
TRUE ~ paste("p = ", label, sep = "")
)
) %>%
dplyr::select(.data = ., -label)
}
# return
return(tibble::as_tibble(df))
}
bf_column <- function(data = NULL, bf) {
# if dataframe is provided
if (!is.null(data)) {
# storing variable name to be assigned later
p_lab <- colnames(dplyr::select(
.data = data,
!!rlang::enquo(bf)
))
# preparing dataframe
df <-
dplyr::select(
.data = data,
# column corresponding to bf-values
bf = !!rlang::enquo(bf),
dplyr::everything()
)
} else {
# if only vector is provided
df <-
base::cbind.data.frame(bf = bf)
}
# make sure the bf-value column is numeric; if not, convert it to numeric
if (!is.numeric(df$bf)) {
# display message about conversion
base::message(cat(
crayon::green("Note:"),
crayon::blue(
"The entered vector is of class",
crayon::yellow(class(df$bf)[[1]]),
"; attempting to convert it to numeric."
)
))
# conversion
df$bf <- as.numeric(as.character(df$bf))
}
# add new support column based on
# Wagenmakers, Wetzels, Borsboom, & Van Der Maas, 2011
df %<>%
dplyr::mutate(
.data = .,
support = dplyr::case_when(
# first condition
bf < .01 ~ "extreme BF01",
bf < .03 & bf >= .01 ~ "very strong BF01",
bf < .1 & bf >= .03 ~ "strong BF01",
bf < 1/3 & bf >= .1 ~ "moderate BF01",
bf < 1 & bf >= 1/3 ~ "anectdotal BF01",
bf >= 1 & bf < 3 ~ "anectdotal BF10",
bf >= 3 & bf < 10 ~ "moderate BF10",
bf >= 10 & bf < 30 ~ "strong BF10",
bf >= 30 & bf < 100 ~ "very strong BF10",
bf >= 100 ~ "extreme BF10",
# fourth condition
bf < 0.001 ~ "***"
)
) %>%
tibble::as_data_frame(x = .) # convert to tibble dataframe
# change back from the generic bf-value to the original name that was provided by the user for the bf-value
if (!is.null(data)) {
# reordering the dataframe
df %<>%
dplyr::select(.data = ., -bf, -support, dplyr::everything())
# renaming the bf-value variable with the name provided by the user
colnames(df)[which(names(df) == "bf")] <- p_lab
}
# return the final tibble dataframe
return(df)
}
testme <- xpairwise_p(data = ggplot2::msleep,
x = vore,
y = brainwt,
type = "bf")
#> Note: The pairwise multiple comparisons test used -
#> BayesFactor::ttestBF.
#>
bf_column(testme, bfactor)
#> Warning: `as_data_frame()` is deprecated, use `as_tibble()` (but mind the new semantics).
#> This warning is displayed once per session.
#> # A tibble: 6 x 5
#> group1 group2 mean.difference bfactor support
#> <chr> <chr> <dbl> <dbl> <chr>
#> 1 carni herbi 0.542 0.540 anectdotal BF01
#> 2 carni insecti -0.0577 0.718 anectdotal BF01
#> 3 carni omni 0.0665 0.427 anectdotal BF01
#> 4 herbi insecti -0.600 0.540 anectdotal BF01
#> 5 herbi omni -0.476 0.571 anectdotal BF01
#> 6 insecti omni 0.124 0.545 anectdotal BF01
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent