For 循环使用 t-stat 函数创建列表答案

【问题标题】：For loop using a t-stat function to create a listFor 循环使用 t-stat 函数创建列表
【发布时间】：2011-11-11 22:31:40
【问题描述】：

我正在使用以下函数来计算数据帧 (x) 中数据的 t-stat：

    wilcox.test.all.genes<-function(x,s1,s2) {
     x1<-x[s1]
     x2<-x[s2]
     x1<-as.numeric(x1)
     x2<-as.numeric(x2)
     wilcox.out<-wilcox.test(x1,x2,exact=F,alternative="two.sided",correct=T)
     out<-as.numeric(wilcox.out$statistic)
     return(out)
    }

我需要编写一个 for 循环来迭代特定次数。对于每次迭代，需要对列进行洗牌，执行上述函数并将最大 t-stat 值保存到列表中。

我知道我可以使用sample() 函数来打乱数据框的列，并使用max() 函数来识别最大t-stat 值，但我不知道如何将它们放在一起实现一个可行的代码。

【问题讨论】：

下面的答案有效，但我仍然不清楚你实际上想要做什么。将您的真实测试统计数据与这些排列的分布输出进行比较将为您提供经验 p 值。但是，您提到取最大值，如果您正在进行许多类似的测试，通常会这样做以纠正多重比较。我的困惑是你的 x 是一个向量，你只执行一个测试。无论如何，如果你澄清一下，我很乐意相应地编辑我的答案。
我使用的数据框有来自两类个人的数据，我正在根据类拆分数据。该函数是对来自两个类的数据的 Wilcox 检验。初始测试完成后，我想打乱数据框的列，再次执行该功能，从新数据中获取最大测试统计量并将此数字保存到列表中。我想执行 500 次改组，从而从改组数据中生成 500 个最大测试统计信息的列表。我将使用此测试统计数据列表来确定 95% 的测试统计数据，以便与原始数据集进行比较。谢谢
我明白了。那正是我所想。如果您提供示例数据会有所帮助，因为我们无法知道您有数据框。此外，还有一些混乱，因为要进行这样的置换测试，您需要置换类标签（即 rows），而不是列。请参阅下面我编辑的答案以获取您的解决方案。祝你好运！

标签： list r

【解决方案1】：

您正在尝试生成经验 p 值，并针对您正在进行的多重比较进行校正，因为您的数据中有多个列。首先，让我们模拟一个示例数据集：

# Simulate data
n.row = 100
n.col = 10

set.seed(12345)
group = factor(sample(2, n.row, replace=T))
data  = data.frame(matrix(rnorm(n.row*n.col), nrow=n.row))

计算每一列的 Wilcoxon 检验，但我们将重复多次，同时置换观察值的类成员资格。这为我们提供了此检验统计量的经验零分布。

# Re-calculate columnwise test statisitics many times while permuting class labels
perms = replicate(500, apply(data[sample(nrow(data)), ], 2, function(x) wilcox.test(x[group==1], x[group==2], exact=F, alternative="two.sided", correct=T)$stat))

通过折叠多重比较来计算最大检验统计量的零分布。

# For each permuted replication, calculate the max test statistic across the multiple comparisons
perms.max = apply(perms, 2, max)

通过简单地对结果进行排序，我们现在可以确定 p=0.05 的临界值。

# Identify critical value 
crit = sort(perms.max)[round((1-0.05)*length(perms.max))]

我们还可以将分布与临界值一起绘制。

# Plot 
dev.new(width=4, height=4)
hist(perms.max)
abline(v=crit, col='red')

最后，将真实检验统计量与此分布进行比较将为您提供经验 p 值，并通过将家庭误差控制到 p

> length(which(perms.max>1600))/length(perms.max)
[1] 0.074

【讨论】：