2

考虑这个数据框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进行 5 次测量var1:var5。对于每个分组变量，我在每个var. 当发现显着差异时，我使用该TukeyHSD()函数和包中的multcompLetters()函数在multcompView组上生成 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())？我该如何解决这个问题，以及如何让它输出一个整洁的表格，该表格显示每个组和我一次给它的每个变量的字母。例如，State使用所有 5 个变量看起来像这样

     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 字母的组的箱线图（针对每个变量）？

4

2 回答 2

3

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

于 2020-05-27T20:51:41.353 回答
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

我切断了结果

于 2020-05-29T15:52:42.097 回答