从post hoc测试中获得整洁的输出

Question

考虑这个数据帧dat1

dat1 <- data.frame(Region = rep(c("r1","r2"), each = 100),
                   State = rep(c("NY","MA","FL","GA"), each = 10),
                   Loc = rep(c("a","b","c","d","e","f","g","h"),each = 5),
                   ID = rep(c(1:10), each=2),
                   var1 = rnorm(200),
                   var2 = rnorm(200),
                   var3 = rnorm(200),
                   var4 = rnorm(200),
                   var5 = rnorm(200))

我有类似于上面创建的dat1的数据帧。Region、State和Loc是每个观察值ID的分组变量，并且对每个观察值var1:var5进行5次测量。对于每个分组变量，我对每个var执行单变量方差分析。当发现显著差异时，我使用multcompView包中的TukeyHSD()函数和multcompLetters()函数在组上生成CLD。由于我想对每个分组变量执行此操作，因此我正在尝试编写一个函数来避免重复和输入错误。下面展示了我对这篇文章的看法：

library(tidyverse)
library(multcomp)
library(multcompView)

Tuk <- function(dat,groupvar,var){
  TUK <- TukeyHSD(aov(lm(get(var) ~ get(groupvar), data=dat)))
  names(TUK)[[1]] <- paste0(groupvar)
  lets<-multcompLetters(extract_p(TUK$groupvar))
  lets
}
#assuming all 5 vars were significant in the anovas, I would then run this for each grouping variable as follows:
vars <- paste0(names(dat1[,5:9]))
#by Region
lapply(vars, FUN=Tuk, dat=dat1, groupvar="Region")
#by State
lapply(vars, FUN=Tuk, dat=dat1, groupvar="State")
#by Loc
lapply(vars, FUN=Tuk, dat=dat1, groupvar="Loc")

代码在函数外部工作。该函数将创建模型，但我不知道如何格式化它，以便识别multcompLetters(extract_p())部件的groupvar是什么？我如何解决这个问题，以及如何让函数输出一个整洁的表格，该表格显示了我一次给出的每个组和每个变量的字母。例如，对于使用所有5个变量的State，它将如下所示

     NY   MA   FL   GA
var1  a    ab   c    a
var2  a    ab   b    c
var3  a    c    ab   bc  
var4  ab   c    ab   ab 
var5  a    b    c     b

另外，有没有一种合理的方法来让这个函数生成显示CLD字母的组的箱图(对于每个变量)？

Answer 1

假设这个情节真的是你想要的，这是否让你非常接近var1 ~ State的sing情节，Indrajeet在构建这个包方面做了很好的工作，我讨厌重新发明轮子。

dat1 <- data.frame(Region = rep(c("r1","r2"), each = 100),
                   State = rep(c("NY","MA","FL","GA"), each = 10),
                   Loc = rep(c("a","b","c","d","e","f","g","h"),each = 5),
                   ID = rep(c(1:10), each=2),
                   var1 = rnorm(200),
                   var2 = rnorm(200),
                   var3 = rnorm(200),
                   var4 = rnorm(200),
                   var5 = rnorm(200))

library(ggstatsplot)

ggbetweenstats(dat1, State, var1, 
               plot.type = "box", 
               pairwise.comparisons = TRUE, 
               pairwise.display = "everything")

#> Note: Shapiro-Wilk Normality Test for var1: p-value = 0.183
#> 
#> Note: Bartlett's test for homogeneity of variances for factor State: p-value = 0.373
#> 

​

​

Answer 2

你接受了答案，但只需做一点工作，你就可以得到你最初要求的东西……

# https://stackoverflow.com/questions/62050403/getting-tidy-output-from-post-hoc-tests/62052651#62052651

library(multcompView)
library(dplyr)
library(purrr)
library(tidyselect)

set.seed(1111)
dat1 <- data.frame(Region = rep(c("r1","r2"), each = 100),
                   State = rep(c("NY","MA","FL","GA"), each = 10),
                   Loc = rep(c("a","b","c","d","e","f","g","h"),each = 5),
                   ID = rep(c(1:10), each = 2),
                   var1 = rnorm(200),
                   var2 = rnorm(200),
                   var3 = rnorm(200),
                   var4 = rnorm(200),
                   var5 = rnorm(200))

# You want just the letters which you can get by ...
multcompLetters(TukeyHSD(aov(var1  ~ State, data = dat1))$State[,4])$Letters
#>  GA  MA  NY  FL 
#> "a" "a" "a" "a"

# Your function redone...
Tuk3 <- function(data,
                 groupvar,
                 var) {
  lst <- as.list(match.call())
  if (is.symbol(lst$groupvar) || is.symbol(lst$var)) {
    stop("Please quote all variables")
  }
  if (!is.call(groupvar)) {
    grouplabel <- rlang::as_name(rlang::enquo(groupvar))
  }

  data <-
    dplyr::select(
      .data = data,
      var = {{ var }},
      groupvar = {{ groupvar }}
    )

  aov_object <- aov(var ~ groupvar, data = data)
  aov_results <- broom::tidy(aov_object) %>%
    mutate(term = if_else(term != "Residuals", grouplabel, term))
  tukey_results <- broom::tidy(TukeyHSD(aov_object)) %>%
    mutate(term = grouplabel)

  # multcompLetters is annoying and wants named vectors
  p_values <- tukey_results %>% pull(adj.p.value)
  names(p_values) <- tukey_results %>% pull(comparison)
  letters_results <-  data.frame(as.list(multcompLetters(p_values)$Letters))
  return(letters_results)
}

# works for one
Tuk3(data = dat1, groupvar = "State", var = "var1")
#>   GA MA NY FL
#> 1  a  a  a  a

# you could do this manually but I do it a lot so I have a function
variables_list <- CGPfunctions::cross2_var_vectors(dat1, 1:3, 5:9)

# Make the names nice
outcomes2 <- variables_list$lista
groupings2 <- variables_list$listb
names(groupings2) <- unlist(groupings2)
names(outcomes2) <- paste(unlist(outcomes2), "~", unlist(groupings2))

# get all 15 results and a final map_dfr to make one tibble
map2(.x = outcomes2,
     .y = groupings2,
     .f = ~ Tuk3(dat = dat1,
                 var = tidyselect::all_of(.x),
                 groupvar = tidyselect::all_of(.y))) %>% 
       map_dfr(~ rbind(.), .id = "Which_ANOVA")
#>      Which_ANOVA   r2   r1   GA   MA   NY   FL    b    c    d    e    f    g
#> 1  var1 ~ Region    a    a <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
#> 2  var2 ~ Region    a    a <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
#> 3  var3 ~ Region    a    a <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
#> 4  var4 ~ Region    a    a <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
#> 5  var5 ~ Region    a    a <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
#> 6   var1 ~ State <NA> <NA>    a    a    a    a <NA> <NA> <NA> <NA> <NA> <NA>
#> 7   var2 ~ State <NA> <NA>    a    a    a    a <NA> <NA> <NA> <NA> <NA> <NA>
#> 8   var3 ~ State <NA> <NA>    a    a    a    a <NA> <NA> <NA> <NA> <NA> <NA>
#> 9   var4 ~ State <NA> <NA>    a    a    a    a <NA> <NA> <NA> <NA> <NA> <NA>
#> 10  var5 ~ State <NA> <NA>    a    b   ab   ab <NA> <NA> <NA> <NA> <NA> <NA>
#> 11    var1 ~ Loc <NA> <NA> <NA> <NA> <NA> <NA>    a    a    a    a    a    a
#> 12    var2 ~ Loc <NA> <NA> <NA> <NA> <NA> <NA>    a    a    a    a    a    a
#> 13    var3 ~ Loc <NA> <NA> <NA> <NA> <NA> <NA>    a    a    a    a    a    a
#> 14    var4 ~ Loc <NA> <NA> <NA> <NA> <NA> <NA>    a    a    a    a    a    a
#> 15    var5 ~ Loc <NA> <NA> <NA> <NA> <NA> <NA>   ab   ab   ab    a   ab    b

我切断了搜索结果