Feature suggestion for Bayesian muliple comparisons

[ibecav]

Hi Indrajeet,

Was working on multiple comparisons for bayes analysis see https://ibecav.netlify.com/post/pairwise-bayesian-comparisons-even-faster/ when it struck me that a nice future feature for ggstatsplot::ggbetweenstats might be to display BF10 similar to the way you do p values for parametric.

I hacked at your code a little bit signif_column and pairwise_p and the reprex kind of shows what's possible...

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 (richarddmorey@gmail.com).
#> 
#> 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

^{Created on 2019-07-11 by the reprex package (v0.3.0)}

[ibecav]

library(tidyverse)
library(ggsignif)
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 (richarddmorey@gmail.com).
#> 
#> Type BFManual() to open the manual.
#> ************
xxx <- msleep %>% filter(!is.na(vore), !is.na(brainwt))
magic1 <- function(x = NULL, y = NULL) {
  results <- ttestBF(x = x, y = y) %>% 
    extractBF() %>% 
    mutate(p.value = bf)
  return(results)
}
ggplot(xxx, aes(x=vore, y=brainwt)) +
  geom_boxplot() +
  geom_signif(comparisons = list(c("carni", "herbi"),
                                 c("carni", "insecti"),
                                 c("carni", "omni")
                                ),
              test = "magic1",
              step_increase = .1)

^{Created on 2019-07-12 by the reprex package (v0.3.0)}

[IndrajeetPatil]

closed in favor of https://github.com/IndrajeetPatil/pairwiseComparisons/issues/10