r - 循环 t.tests 以获取 r 中的数据帧子集

Question

我有一个包含 32 个变量的数据框“math.numeric”。每行代表一个学生，每个变量是一个属性。学生们根据他们的期末成绩被分成了 5 组。

数据如下:

head(math.numeric)
school sex age address famsize Pstatus Medu Fedu Mjob Fjob reason ... group
1      1   18  2       1       1       4    4    1    5    1          2
1      1   17  2       1       2       1    1    1    3    1          2
1      1   15  2       2       2       1    1    1    3    3          3
1      1   15  2       1       2       4    2    2    4    2          4
1      1   16  2       1       2       3    3    3    3    2          3
1      2   16  2       2       2       4    3    4    3    4          4

我正在对第 1 组与所有其他组的每个变量进行 t 检验，以识别该组显着不同的属性。我希望提取每个测试的 p 值，例如:

t.test(subset(math.numeric$school, math.numeric$group == 1),
      subset(math.numeric$school, math.numeric$group != 1))$p.value
t.test(subset(math.numeric$sex, math.numeric$group == 1), 
        subset(math.numeric$sex, math.numeric$group != 1))$p.value
t.test(subset(math.numeric$age, math.numeric$group == 1), 
        subset(math.numeric$age, math.numeric$group != 1))$p.value

我一直试图弄清楚如何创建一个循环来执行此操作，而不是一次写出每个测试。我尝试了一个 for 循环和 lapply，但到目前为止我还没有运气。

我对此很陌生，因此我们将不胜感激。

考特尼

Answer 1

您的示例数据不足以对所有子组实际执行 t 检验。出于这个原因，我采用 iris 数据集，其中包含 3 种植物:Setosa、Versicolor 和 Virginica。这些是我的小组。您将不得不相应地调整您的代码。下面我将展示如何测试一组与所有其他组、一组与另一组以及各个组的所有组合。

一个组与所有其他组的组合:

首先，假设我想将 Versicolor 和 Virginica 与 Setosa 进行比较，即 Setosa 是我的 group 1，所有其他组都应与之进行比较。以下是实现您想要的简单方法:

sapply(names(iris)[-ncol(iris)], function(x){
             t.test(iris[iris$Species=="setosa", x], 
                    iris[iris$Species!="setosa", x])$p.value 
                    })
Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
7.709331e-32 1.035396e-13 1.746188e-69 1.347804e-60 

在这里，我提供了数据集中不同变量的名称names(iris) - 不包括指示分组变量的列[-ncol(iris)] (因为它是最后一列) - 作为 sapply 的向量，它将相应的名称作为参数传递给我定义的函数。

一组与其他组:

如果您想对所有组进行分组比较，以下操作可能会有所帮助:首先，创建一个包含您将要做的所有组 x 变量组合的数据框，不包括分组变量本身和引用组，类(class)。这可以通过以下方式实现:

comps <- expand.grid(unique(iris$Species)[-1], # excluding Setosa as reference group
                     names(iris)[-ncol(iris)] # excluding group column
                     )
head(comps)
        Var1         Var2
1 versicolor Sepal.Length
2  virginica Sepal.Length
3 versicolor  Sepal.Width
4  virginica  Sepal.Width
5 versicolor Petal.Length
6  virginica Petal.Length

这里，Var1 是不同的物种，Var2 是要进​​行比较的不同变量。在这种情况下，引用 group 1 或 Setosa 是隐含的。现在，我可以使用 apply 来创建测试。我通过使用 comps 的每一行作为带有两个元素的参数来做到这一点，其中第一个指示轮到哪个组，第二个参数指示应该比较哪个变量。这些将用于对原始数据帧进行子集化。

comps$pval <- apply(comps, 1, function(x) {
    t.test(iris[iris$Species=="setosa", x[2]], iris[iris$Species==x[1], x[2]])$p.value 
    } )

其中第 1 组又名 Setosa 在函数中被硬编码。这给了我一个数据框，其中包含所有组合的 p 值(以 Setosa 作为引用组)，以便它们易于查找:

head(comps)
        Var1         Var2         pval
1 versicolor Sepal.Length 3.746743e-17
2  virginica Sepal.Length 3.966867e-25
3 versicolor  Sepal.Width 2.484228e-15
4  virginica  Sepal.Width 4.570771e-09
5 versicolor Petal.Length 9.934433e-46
6  virginica Petal.Length 9.269628e-50

所有组的组合:

您可以轻松扩展上述内容以生成一个数据框，其中包含每个组组合的 t 检验的 p 值。一种方法是:

comps <- expand.grid(unique(iris$Species), unique(iris$Species), names(iris)[-ncol(iris)])

现在有三列。前两个是组，第三个是变量:

head(comps)
        Var1       Var2         Var3
1     setosa     setosa Sepal.Length
2 versicolor     setosa Sepal.Length
3  virginica     setosa Sepal.Length
4     setosa versicolor Sepal.Length
5 versicolor versicolor Sepal.Length
6  virginica versicolor Sepal.Length

您可以使用它来进行测试:

comps$pval <- apply(comps, 1, function(x) {
  t.test(iris[iris$Species==x[1], x[3]], iris[iris$Species==x[2], x[3]])$p.value 
} )

我收到一条错误消息:我该怎么办？

t.test 如果样本量太小或一组的值是恒定的，则可能会抛出错误消息。这是有问题的，因为它可能只发生在特定的群体中，而且您可能事先不知道它是哪一个群体。然而，该错误会中断对 apply 的整个函数调用，您将看不到任何结果。

避免这种情况并识别有问题的组的一种方法是将函数 t.test 包裹在 dplyr::failwith 周围(另请参见 ?tryCatch)。为了说明这是如何工作的，请考虑以下几点:

smalln <- data.frame(a=1, b=2)
t.test(smalln$a, smalln$b)
> Error in t.test.default(smalln$a, smalln$b) : not enough 'x' observations

failproof.t <- failwith(default="Some default of your liking", t.test, quiet = T)
failproof.t(smalln$a, smalln$b)
[1] "Some default of your liking"

这样，每当 t.test 抛出错误时，您会得到一个字符作为结果，并且计算会继续与其他组一起进行。不用说，您还可以将 default 设置为数字或其他任何值。它不必是一个字符。

统计免责声明:说了这么多，请注意进行几次 t 检验不一定是好的统计实践。您可能希望调整 p 值以考虑多次测试，或者您可能希望使用进行联合测试的替代测试程序。