使用dplyr创建t.test表格？

Question

假设我有这样的数据：

set.seed(031915)
myDF <- data.frame(
  Name= rep(c("A", "B"), times = c(10,10)),
  Group = rep(c("treatment", "control", "treatment", "control"), times = c(5,5,5,5)),
  X = c(rnorm(n=5,mean = .05, sd = .001), rnorm(n=5,mean = .02, sd = .001),
        rnorm(n=5,mean = .08, sd = .02), rnorm(n=5,mean = .03, sd = .02))
)

我希望创建一个t.test表格，其中有一个"A"的行和一个"B"的行。

我可以编写自己的函数来实现这一点：

ttestbyName <- function(Name) {
  b <- t.test(myDF$X[myDF$Group == "treatment" & myDF$Name==Name], 
              myDF$X[myDF$Group == "control" & myDF$Name==Name], 
              conf.level = 0.90)
  dataNameX <- data.frame(Name = Name,
                          treatment = round(b$estimate[[1]], digits = 4),
                          control = round(b$estimate[[2]], digits = 4),
                          CI = paste('(',round(b$conf.int[[1]], 
                                                  digits = 4),', ',
                                        round(b$conf.int[[2]], 
                                              digits = 4), ')',
                                        sep=""),
                          pvalue = round(b$p.value, digits = 4),
                          ntreatment = nrow(myDF[myDF$Group == "treatment" & myDF$Name==Name,]),
                          ncontrol = nrow(myDF[myDF$Group == "control" & myDF$Name==Name,]))
}
library(parallel)
Test_by_Name <- mclapply(unique(myDF$Name), ttestbyName)
Test_by_Name <- do.call("rbind", Test_by_Name)

并且输出如下：

 Name treatment control               CI pvalue ntreatment ncontrol
1    A    0.0500  0.0195 (0.0296, 0.0314) 0.0000          5        5
2    B    0.0654  0.0212  (0.0174, 0.071) 0.0161          5        5

我想知道是否有更简洁的方法使用dplyr完成这个任务。我考虑过使用groupby，但我有点迷失。

谢谢！

- Ignacio

Answer 1

虽然没有太大的改进，但是这里有一个提高:

library(dplyr)

ttestbyName <- function(myName) {
  bt <- filter(myDF, Group=="treatment", Name==myName)
  bc <- filter(myDF, Group=="control", Name==myName)

  b <- t.test(bt$X, bc$X, conf.level=0.90)

  dataNameX <- data.frame(Name = myName,
                      treatment = round(b$estimate[[1]], digits = 4),
                      control = round(b$estimate[[2]], digits = 4),
                      CI = paste('(',round(b$conf.int[[1]], 
                                           digits = 4),', ',
                                 round(b$conf.int[[2]], 
                                       digits = 4), ')',
                                 sep=""),
                      pvalue = round(b$p.value, digits = 4),
                      ntreatment = nrow(bt),  # changes only in
                      ncontrol = nrow(bc))    # these 2 nrow() args
}

你应该使用 data.table 中的 rbindlist 函数替代 do.call 函数：

library(data.table)
Test_by_Name <- lapply(unique(myDF$Name), ttestbyName)
Test_by_Name <- rbindlist(Test_by_Name)

或者更好的方法是使用%>% 管道：

Test_by_Name <- myDF$Name %>% 
                unique %>% 
                lapply(., ttestbyName) %>% 
                rbindlist

> Test_by_Name
 Name treatment control               CI pvalue ntreatment ncontrol
1:    A    0.0500  0.0195 (0.0296, 0.0314) 0.0000          5        5
2:    B    0.0654  0.0212  (0.0174, 0.071) 0.0161          5        5

- desertnaut

Answer 2

这是一个老问题，但broom包已经为此提供了确切的解决方案（以及其他统计测试）：

library(broom)
library(dplyr)
myDF %>% group_by(Name) %>%
         do(tidy(t.test(X~Group, data = .)))
Source: local data frame [2 x 9]
Groups: Name [2]

    Name    estimate  estimate1  estimate2  statistic      p.value
  (fctr)       (dbl)      (dbl)      (dbl)      (dbl)        (dbl)
1      A -0.03050475 0.01950384 0.05000860 -63.838440 1.195226e-09
2      B -0.04423181 0.02117864 0.06541046  -3.104927 1.613625e-02
Variables not shown: parameter (dbl), conf.low (dbl), conf.high (dbl)

- jeremycg

---

如果我想对多个组进行单样本单尾检验怎么办？ - Polar Bear

你好，我正在寻找一种更完整的方法，即能够打印p.adjust、置信区间和显著性水平...有任何线索吗？ - 12666727b9

Answer 3

library(tidyr)
library(dplyr)
myDF %>% group_by(Group) %>% mutate(rowname=1:n())%>% 
  spread(Group, X) %>% 
  group_by(Name) %>%
  do(b = t.test(.$control, .$treatment)) %>%  
  mutate(
         treatment = round(b[['estimate']][[2]], digits = 4),
         control = round(b[['estimate']][[1]], digits = 4),
         CI = paste0("(", paste(b[['conf.int']], collapse=", "), ")"),
         pvalue = b[['p.value']]
         )
#  Name treatment control                                        CI       pvalue
#1    A    0.0500  0.0195 (-0.031677109707283, -0.0293323994902097) 1.195226e-09
#2    B    0.0654  0.0212 (-0.0775829100729602, -0.010880719830447) 1.613625e-02

您可以手动添加ncontrol和ntreatment。

- Khashaa

Answer 4

您可以使用自定义的t.test函数和do完成：

my.t.test <- function(data, formula, ...)
{
    tt <- t.test(formula=formula, data=data, ...)
    ests <- tt$estimate
    names(ests) <- sub("mean in group ()", "\\1",names(ests))
    counts <- xtabs(formula[c(1,3)],data)
    names(counts) <- paste0("n",names(counts))
    cbind(
          as.list(ests),
          data.frame(
            CI = paste0("(", paste(tt$conf.int, collapse=", "), ")"),
            pvalue = tt$p.value,
            stringsAsFactors=FALSE
            ),
          as.list(counts)
    )
}

myDF %>% group_by(Name) %>% do(my.t.test(.,X~Group))
Source: local data frame [2 x 7]
Groups: Name

  Name    control  treatment                                        CI       pvalue ncontrol ntreatment
1    A 0.01950384 0.05000860 (-0.031677109707283, -0.0293323994902097) 1.195226e-09        5          5
2    B 0.02117864 0.06541046 (-0.0775829100729602, -0.010880719830447) 1.613625e-02        5          5

- James

---

这对我很有帮助，但我发现在调用data.frame()时执行ci_lb = tt$conf.int[[1]], ci_ub = tt$conf.int[[2]]更快、更方便，因为(i)结果将是数字而不是字符，(ii)拆分塞入单个列中的多个值有点困难。 - Curt F.