使用组均值跨多个组进行一次样本t.test

Question

我感兴趣的是，当使用两个样本的t检验与来自相同数据的单样本t检验进行均值比较时，结果可能会有什么不同。比较了使用两种不同的技术(“源”)从多个总体("Group_ID")中随机抽样“计数”的两种评估。p>0.05的结果表明，均值在统计上没有差异(即，两种评估技术给出的结果相似)。source=A的原始数据是可靠的，并且始终可用。source=B的原始数据可能不可用，但将始终提供平均值。使用测试数据集，我想通过使用一个样本t.test而不是两个样本t.test来检查t.test结果如何不同。

使用dplyr和broom函数，我已经确定了如何在许多情况下("Group_ID")进行多个双样本t测试，在这些情况下进行了两次评估。来自两个源的数据被合并，以创建包含由两个源(“源”)之一标识的原始计数(“计数”)的三列数据帧。

glimpse(Data)

Observations: 2,552
Variables: 3
$ Group_ID    <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
$ count       <dbl> 7, 8, 5, 5, 7, 3, 4, 2, 8, 11, 12, 1, 3, 5, 5, 12, 1, 5...
$ source      <chr> "B", "B", "B", "B", "B", "B", "B", "B",...

双样本测试比较了两个来源的平均值，并使用以下方法考虑了两个数据集中的差异：

Data_Stats <- Data %>% 
    group_by(Group_ID) %>%
    do(tidy(t.test(count ~ source, alt="two.sided", conf=0.95, var.eq=FALSE, paired=FALSE, data = .)))

结果：

Observations: 38
Variables: 11
Groups: Group_ID [38]
$ Group_ID <fct> 1, 1029, 1032, 1033, 1041, 1044, 1064, 1065, 1067, 1080, 1081, 1083, 1084, 117, 127, 180, 2...
$ estimate    <dbl> -0.4250000, -6.5000000, -1.1944444, 0.3437500, -5.2250000, -1.4375000, -1.6250000, -1.48387...
$ estimate1   <dbl> 5.250000, 9.166667, 5.833333, 6.156250, 5.375000, 3.937500, 2.075000, 6.000000, 4.108108, 9...
$ estimate2   <dbl> 5.675000, 15.666667, 7.027778, 5.812500, 10.600000, 5.375000, 3.700000, 7.483871, 6.540541,...
$ statistic   <dbl> -0.42469044, -3.42643903, -1.19922603, 0.32509809, -3.36599817, -1.94947775, -2.47005992, -...
$ p.value     <dbl> 6.723526e-01, 1.480949e-03, 2.355386e-01, 7.463614e-01, 1.509467e-03, 5.649284e-02, 1.57612...
$ parameter   <dbl> 70.66653, 38.05593, 55.46167, 54.07284, 47.94418, 53.45682, 75.67387, 44.38453, 49.87894, 6...
$ conf.low    <dbl> -2.420560, -10.340117, -3.190125, -1.776090, -8.346179, -2.916196, -2.935370, -3.969856, -5...
$ conf.high   <dbl> 1.570559944, -2.659882919, 0.801236360, 2.463590202, -2.103821060, 0.041195928, -0.31462953...
$ method      <chr> "Welch Two Sample t-test", "Welch Two Sample t-test", "Welch Two Sample t-test", "Welch Two...
$ alternative <chr> "two.sided", "two.sided", "two.sided", "two.sided", "two.sided", "two.sided", "two.sided", ...

我知道我可以使用以下命令为每个案例获取方法：

Data_means <- Data %>% 
    group_by(Group_ID) %>% 
    summarize(count_mean = mean(count))

我正在寻找一个建议，如何最好地使用来自source=2的每个Group_ID的平均值作为t.test函数中"mu=“调用的值，以便与对38个不同的source=1使用单样本t检验的平均值进行比较？

Answer 1

这是一个不适合与dplyr进行数据争论的问题，使用split.data.frame和lapply创建测试可能更容易，但这里介绍了如何在不拼接不同数据帧的情况下解决这个问题。

首先，我需要一些与您问题中的数据相同的可重现数据：

library(tidyverse)

set.seed(69)
df <- data.frame(Group_ID = factor(rep(1:4, each = 40)),
                 count = sample(10, 160, T) + rep(1:2, 80) + rep(1:4, each = 40),
                 source = factor(rep(c("A", "B"), 80)))

现在，我们可以使用与您使用的方法类似的方法来获得两个样本的p值。然后，我们使用的技巧是在每个ID中获取每个源的平均计数，然后对数据帧进行解组并复制平均值，但将其移位，使"B“平均值位于"A”行。接下来，我们可以使用它作为"A“的单样本t检验的平均值。当我们总结时，我们同时有1个样本和2个样本数据。

df                                                                      %>%
group_by(Group_ID)                                                      %>%
mutate(two_group_pval = t.test(count ~ source)$p.value)                 %>%
group_by(Group_ID, source)                                              %>%
mutate(mean_A = mean(count))                                            %>%
arrange(source, .by_group = T)                                          %>%
group_by(Group_ID)                                                      %>%
mutate(mean_B = lead(mean_A, length(which(source == "B"))))             %>%
filter(source == "A")                                                   %>%
group_by(Group_ID)                                                      %>%
mutate(one_group = t.test(count, mu = mean(mean_B, na.rm = T))$p.value) %>%
summarise(observations = length(count),
          mean_A = mean(mean_A, na.rm = T),
          mean_B = mean(mean_B, na.rm = T),
          one_sample_p_value = mean(one_group),
          two_sample_p_value = mean(two_group_pval))

#> # A tibble: 4 x 6
#>   Group_ID observations mean_A mean_B one_sample_p_value two_sample_p_value
#>   <fct>           <int>  <dbl>  <dbl>              <dbl>              <dbl>
#> 1 1                  20    6.6   8.65            0.00341             0.0201
#> 2 2                  20    8.3   9.85            0.0364              0.0999
#> 3 3                  20    9.9  11.5             0.0103              0.0600
#> 4 4                  20   10.4  11.4             0.122               0.290 

您会注意到，在我的数据中，两个样本测试的p值更高。这是因为样本是从均匀分布而不是钟形曲线中提取的，所以t检验的假设不成立。在依赖于单样本或双样本t检验之前，您应该检查您自己的数据是否近似正态分布。如果它们不正常，你应该切换到wilcox.test。

编辑

基于OP的需求变得更加清晰，上面的代码不适用于给定的Group_ID的两个级别的示例。以下是如何使用自解释代码在任意级别的base R中解决该问题：

multi_ss_t_test <- function(x, y) as.numeric(t.test(x$count, mu = y)$p.value)
multi_ts_t_test <- function(x, y) as.numeric(t.test(x$count, y$count)$p.value)

source_dfs <- split.data.frame(df, df$source)

A_groups   <- split.data.frame(source_dfs$A, source_dfs$A$Group_ID)
B_groups   <- split.data.frame(source_dfs$B, source_dfs$B$Group_ID)

B_means    <- tapply(source_dfs$B$count, source_dfs$B$Group_ID, mean)
A_means    <- tapply(source_dfs$A$count, source_dfs$A$Group_ID, mean)

ss_pvals   <- mapply(multi_ss_t_test, A_groups, B_means)
ts_pvals   <- mapply(multi_ts_t_test, A_groups, B_groups)

result     <- data.frame(group = levels(df$Group_ID),
                         source_A_mean = A_means,
                         source_B_means = B_means,
                         one_sample_pval = ss_pvals,
                         two_sample_pval = ts_pvals)

现在，如果我提供一个具有38个级别的数据框，您可以看到它将为每个Group_ID输出一个样本和两个样本的p值

set.seed(69)
df <- data.frame(Group_ID = factor(rep(rep(1:38, each = 20), 2)),
                 count = c(sample(10:40, 760, T), sample(12:42, 760, T)),
                 source = rep(c("A", "B"), each = 760))

运行上面的程序，就会得到以下结果：

as_tibble(result)
#> # A tibble: 38 x 5
#>    group source_A_mean source_B_means one_sample_pval two_sample_pval
#>    <fct>         <dbl>          <dbl>           <dbl>           <dbl>
#>  1 1              24.5           25.6          0.620           0.724 
#>  2 2              23.6           25.0          0.407           0.608 
#>  3 3              24.2           29.7          0.0154          0.0438
#>  4 4              26             22.6          0.123           0.264 
#>  5 5              26.6           25.6          0.531           0.683 
#>  6 6              26.8           24.7          0.303           0.414 
#>  7 7              25.5           26            0.807           0.852 
#>  8 8              24.3           28.9          0.0167          0.0887
#>  9 9              23.6           26.6          0.137           0.255 
#> 10 10             25.8           27.1          0.533           0.651 
#>#... with 28 more rows

如果您只需要p值，一种更有效的方法是：

get_pvals <- function(x)
{
  c(one_sample_p_value = t.test(x$count ~ x$source)$p.value,
    two_sample_p_value = t.test(x$count[x$source == "A"], 
                                mu = mean(x$count[x$source == "B"]))$p.value) 
}

split(df, list(group = df$Group_ID)) %>% sapply(get_pvals) %>% t %>% as.data.frame()