在 R 中，如何编写对数据帧列表运行 t 检验的函数

Question

以iris数据集为例，我想写一个用户自定义函数

运行 pairwise t-test 所有 4 列免除每个数据拆分的 Species 列
将结果导出为一个 csv 文件的 3 个工作表

下面是我的尝试：

library(tidyr)
library(reshape) # for melting /stacking the data
library(multcomp) # for pairwise test 
library(xlsx) # export excel file with worksheet

options(scipen = 100)

# dataset
iris

data_stats <- function(data){
    # melt the dataframe
     df <- melt(data, id.vars=c('Species'),var='group')
    
# split the dataframe into three list of dataframe
    dfsplit<-split(df,df$column)
    
    #  pairwise t-test
    results <- pairwise.t.test(dfsplit$value, dfsplit$group,p.adjust.method = "BH")
    
    # export each result as a worksheet of an excel file
    write.xlsx(results, file="Results.xlsx", sheetName="versicolor_stats", row.names=FALSE)
    write.xlsx(results, file="Results.xlsx", sheetName="virginica_stats", append=TRUE, row.names=FALSE)
    write.xlsx(results, file="Results.xlsx", sheetName="setosa_stats", append=TRUE, row.names=FALSE)
    
}
# testing the code on iris data
data_stats(iris)

请发表评论并分享您的代码。谢谢

Answer 1

这是 tidyverse 的选项 - 使用 pivot_longer 重塑为 'long' 格式，然后使用 group_modify 执行 pairwise.t.test 、tidy 输出和 unnest list 输出

library(dplyr)
library(tidyr)
library(broom)
ttest_out <- iris %>% 
  pivot_longer(cols = -Species) %>% 
  group_by(Species) %>%
  group_modify(~ .x %>%
     summarise(out = list(pairwise.t.test(value, name) %>% 
       tidy))) %>% 
   ungroup %>% 
   unnest(out)

-输出

ttest_out
# A tibble: 18 × 4
   Species    group1       group2         p.value
   <fct>      <chr>        <chr>            <dbl>
 1 setosa     Petal.Width  Petal.Length 1.77e- 54
 2 setosa     Sepal.Length Petal.Length 2.77e-132
 3 setosa     Sepal.Length Petal.Width  1.95e-156
 4 setosa     Sepal.Width  Petal.Length 1.61e- 86
 5 setosa     Sepal.Width  Petal.Width  1.13e-123
 6 setosa     Sepal.Width  Sepal.Length 4.88e- 71
 7 versicolor Petal.Width  Petal.Length 5.35e- 90
 8 versicolor Sepal.Length Petal.Length 3.78e- 52
 9 versicolor Sepal.Length Petal.Width  5.02e-125
10 versicolor Sepal.Width  Petal.Length 1.36e- 45
11 versicolor Sepal.Width  Petal.Width  3.46e- 44
12 versicolor Sepal.Width  Sepal.Length 1.25e- 95
13 virginica  Petal.Width  Petal.Length 1.39e- 90
14 virginica  Sepal.Length Petal.Length 6.67e- 22
15 virginica  Sepal.Length Petal.Width  3.47e-110
16 virginica  Sepal.Width  Petal.Length 2.35e- 68
17 virginica  Sepal.Width  Petal.Width  1.87e- 19
18 virginica  Sepal.Width  Sepal.Length 2.47e- 92

Answer 2

更新：这道题的统计部分（将pairwise.t.test应用于鸢尾花数据集）已经answered previously on SO. And here另一种解法。

接受的解决方案运行一系列成对的 t-tests 并生成一列 p-values（正如问题所问），但它不是一组非常有意义的 t-tests。您可能会怀疑，从我们看到 p-values 像 2.77e-132 并且 group1 和 group2 是连续变量，而不是因子的水平这一事实来看。

这些测试评估的假设是，对于每个物种，萼片是否与花瓣相同，长度是否与宽度相同。成对 t-test 程序旨在比较所有级别的单个连续变量（例如萼片长度）和一个因子（例如物种）。

首先，让我们将 pairwise.t.test 应用到列 Sepal.Length，以便稍后检查我们是否得到正确的 p-values。

library("broom")
library("tidyverse")

pairwise.t.test(iris$Sepal.Width, iris$Species)
#> 
#>  Pairwise comparisons using t tests with pooled SD 
#> 
#> data:  iris$Sepal.Width and iris$Species 
#> 
#>            setosa  versicolor
#> versicolor < 2e-16 -         
#> virginica  9.1e-10 0.0031    
#> 
#> P value adjustment method: holm

如果您看过鸢尾花数据集，您就会知道这些 p-values“有道理”：Virginica 和 Versicolor 彼此之间的相似性比 Setosa 更相似。

现在让我们对四个数字列进行整齐的测试。

t_pvals <- iris %>%
  pivot_longer(
    -Species,
    names_to = "Variable",
    values_to = "x"
  ) %>%
  # The trick to performing the right tests is to group the tibble by Variable,
  # not by Species because Species is the grouping variable for the t-tests.
  group_by(
    Variable
  ) %>%
  group_modify(
    ~ tidy(pairwise.t.test(.x$x, .x$Species))
  ) %>%
  ungroup()
t_pvals
#> # A tibble: 12 × 4
#>    Variable     group1     group2      p.value
#>    <chr>        <chr>      <chr>         <dbl>
#>  1 Petal.Length versicolor setosa     1.05e-68
#>  2 Petal.Length virginica  setosa     1.23e-90
#>  3 Petal.Length virginica  versicolor 1.81e-31
#>  4 Petal.Width  versicolor setosa     2.51e-57
#>  5 Petal.Width  virginica  setosa     2.39e-85
#>  6 Petal.Width  virginica  versicolor 8.82e-37
#>  7 Sepal.Length versicolor setosa     1.75e-15
#>  8 Sepal.Length virginica  setosa     6.64e-32
#>  9 Sepal.Length virginica  versicolor 2.77e- 9
#> 10 Sepal.Width  versicolor setosa     5.50e-17
#> 11 Sepal.Width  virginica  setosa     9.08e-10
#> 12 Sepal.Width  virginica  versicolor 3.15e- 3

Sepal.Width 比较的 p-values 在底部。我们得到了预期的 p-values！

接下来我们格式化 p-values 以便它们看起来更舒适。

t_pvals <- t_pvals %>%
  mutate(
    across(
      p.value, rstatix::p_format,
      accuracy = 0.05
    )
  )
t_pvals
#> # A tibble: 12 × 4
#>    Variable     group1     group2     p.value
#>    <chr>        <chr>      <chr>      <chr>  
#>  1 Petal.Length versicolor setosa     <0.05  
#>  2 Petal.Length virginica  setosa     <0.05  
#>  3 Petal.Length virginica  versicolor <0.05  
#>  4 Petal.Width  versicolor setosa     <0.05  
#>  5 Petal.Width  virginica  setosa     <0.05  
#>  6 Petal.Width  virginica  versicolor <0.05  
#>  7 Sepal.Length versicolor setosa     <0.05  
#>  8 Sepal.Length virginica  setosa     <0.05  
#>  9 Sepal.Length virginica  versicolor <0.05  
#> 10 Sepal.Width  versicolor setosa     <0.05  
#> 11 Sepal.Width  virginica  setosa     <0.05  
#> 12 Sepal.Width  virginica  versicolor <0.05

最后我们将结果保存到文件中。

t_pvals %>%
  write_csv("pairwse-t-tests-on-iris-data.csv")

在 R 中，如何编写对数据帧列表运行 t 检验的函数

In R, how to write function that runs t-test on list of dataframes

excel

split

r

list

dataframe