logo标签列表

在R中循环遍历数据框子集进行t检验

rloopssubsetp-value
3

我有一个名为“math.numeric”的数据框，其中包含32个变量。每一行代表一个学生，每个变量都是一个属性。根据他们的最终成绩，学生已被分为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

我正在对每个变量执行t检验，以比较第1组和其他所有组之间的重要差异。我希望获取每个检验的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，但到目前为止都没有成功。

我对此还比较新，所以非常感谢任何帮助。

Courtney

你应该包含一些样本数据以使你的问题具有可重现性，并且更容易被帮助解决。 - MrFlick
3个回答
5
您提供的数据示例不足以对所有子组进行T检验。因此，我使用包含 3 种植物（Setosa、Versicolor 和 Virginica）的iris数据集作为我的分组。您需要相应地调整代码。下面我将展示如何测试一个组与所有其他组，一个组与每个其他组以及所有单独组的组合之间的差异。
一个组与所有其他组的比较：
首先，假设我想比较Versicolor和Virginica与Setosa，即Setosa是我的“第一组”，应该与所有其他组进行比较。实现您想要的简单方法如下：
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变量组合的数据框，当然不包括分组变量本身和参考组。可以通过以下方式实现：
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 
    } )

函数中硬编码了第一组，也就是Setosa。这使得我得到了一个数据框，其中包含了所有组合（以Setosa作为参考组）的p值，方便查阅。

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值以考虑多次测试，或者使用进行联合测试的替代测试程序。

3
这个怎么样？
pvals <- numeric() #the vector of p values
k <- 1 #in case you choose to use a subset not continuing from 1

# "for(i in seq(1,5))" is just doing the pvalues for the first 5 columns. You could do a 
# list, like "c(1,2,4)" (in place of "seq(1,5)"), to do tests for columns 1, 2, and 4. 
# To do all of the columns, try "for(i in seq(1,(ncol(math.numeric)-1)))".

for(i in seq(1,5)){

  # using your code to grab the p-values and store them in the kth element of "pvals"
  pvals[k] <- t.test(subset(math.numeric[,i], math.numeric$group == 1),
         subset(math.numeric[,i], math.numeric$group != 1))$p.value    

  #iterating the "pvals" vector entry counter
  k=k+1
}
pvals #printing the p values for each test
1
考虑按组将数据框分割，并在列之间使用mapply()。输出变为测试统计量的编译列表：统计量、参数、p值、置信区间等。
# FILTER ROWS AND SUBSET NUMERIC COLS
group1df <- df[df$group==1, 1:ncol(df)-1]
othgroupdf <- df[df$group!=1, 1:ncol(df)-1]

# T-TEST FCT
tfct <- function(v1, v2){
      t.test(v1, v2) 
}

# RUN T-TESTS BY COL, SAVE RESULTS TO LIST
ttests <- mapply(tfct, group1df, othgroupdf)
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接