从事后测试中获得整洁的输出答案

【问题标题】：Getting tidy output from post hoc tests从事后测试中获得整洁的输出
【发布时间】：2020-09-14 22:11:03
【问题描述】：

考虑这个数据框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 进行单变量方差分析。当发现显着差异时，我使用 TukeyHSD() 函数和 multcompView 包中的 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")

代码在函数之外工作。该函数将创建模型，但我不知道如何格式化它，以便它识别 groupvar 是什么 multcompLetters(extract_p()) 部分？我该如何解决这个问题，以及如何让它输出一个整洁的表格，显示每个组和我一次给它的每个变量的字母。例如，使用所有 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 字母的组的箱线图（针对每个变量）？

【问题讨论】：

要使函数正常工作，请将 TUK$groupvar 替换为 TUK[[groupvar]]。
您的问题是否相关？目标是生成箱线图并标记 p.values 的差异吗？
@Chuck P 是的，我希望该函数最终生成箱线图并标记 p.values 的差异。 @stefan 对 groupvar 问题的解决方案有效，我仍然希望在使用它来生成绘图方面获得帮助
嗯，让我向您展示一个现有的包装，在我们投入时间推出您自己的包装之前，它可以满足您的需求。我对multcompLetters 不是很熟悉，但对单向 AOV 的绘图和图表非常熟悉。

标签： r function functional-programming tidyeval tukey

【解决方案1】：

您接受了答案，但只需记录一下您就可以通过一些工作得到您最初要求的内容...

# 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

我切断了结果

【讨论】：

这看起来很有希望。不幸的是，当我使用我的数据尝试此操作时，我收到错误消息“请引用所有变量”。我不确定为什么会发生这种情况，因为我不理解这部分代码。您能否解释一下何时出现错误消息？
如果您传递给函数的变量列表具有裸变量名而不是作为字符串传递，则会显示错误消息。因此，请确保您使用的任何 outcomes2 和 groupings2 都是字符串

【解决方案2】：

假设情节确实是您正在寻找的内容，这是否会让您非常接近 var1 ~ State 的歌唱情节，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
#>

【讨论】：