两个不相关的正态变量之间的样本相关系数分布如何?

机器算法验证 r 分布 相关性 皮尔逊-r 斯皮尔曼罗
2022-03-04 02:51:00

我想比较观察到的双变量(Pearson'sρ和斯皮尔曼的ρ) 相关系数与随机数据的预期值。

假设我们测量了 36 个案例,涉及非常多的变量(1000 个)。(我知道这很奇怪,它被称为Q 方法论。进一步假设每个变量(严格)正态分布 在案例中。(同样,非常奇怪,但确实如此,因为作为人员变量的人将排序项目案例排序在正态分布。)

所以,如果人们随机排序,我们应该得到:

m <- sapply(X = 1:1000, FUN = function(x) rnorm(36))

现在——因为这是 Q 方法——我们将所有人员变量关联起来:

cors <- cor(x = m, method = "pearson")

然后我们尝试绘制它,并将 Pearson 相关系数在随机数据中的分布叠加起来,这实际上应该非常接近于我们的假数据中观察到的相关性:

library(ggplot2)
cor.data <- cors[upper.tri(cors, diag = FALSE)]  # we're only interested in one of the off-diagonals, otherwise there'd be duplicates
cor.data <- as.data.frame(cor.data)  # that's how ggplot likes it
colnames(cor.data) <- "pearson"
g <- ggplot(data = cor.data, mapping = aes(x = pearson))
g <- g + xlim(-1,1)  # actual limits of pearsons r
g <- g + geom_histogram(mapping = aes(y = ..density..))
g <- g + stat_function(fun = dt, colour = "red", args = list(df = 36-1))
g

这给出了:

密度图

叠加曲线显然是错误的。(另请注意,虽然很奇怪,但 y 轴密度实际上是正确的:因为 x 值非常小,这就是面积总和为 1 的方式)。

我(模糊地)记得 t 分布在这种情况下是相关的,但我不知道如何正确地对其进行参数化。特别是,自由度是由相关数(1000^2/2-500)还是这些相关性所基于的观察数(36)给出的?

无论哪种方式,上面的叠加曲线显然是错误的。

我也很困惑,因为皮尔逊 r 的概率分布需要有界(没有超出 (-) 1 的值)——但 t 分布是没有界的。

哪个分布描述了 Pearson 的ρ在这种情况下?


奖金:

上面的数据实际上是理想化的:在我真正的 Q 研究中,人变量实际上在正态分布下只有很少的列可以将他们的项目案例分类,如下所示:

q-排序

实际上,人员变量实际上是按等级排序的项目案例,因此 Pearson's 不适用。作为一个粗略和肮脏的修复,我选择了斯皮尔曼的ρ, 反而。Spearman 的概率分布是否相同ρ?


更新:如果有人感兴趣,下面是实现@amoeba 精彩响应的 R 代码:

library(ggplot2)
cor.data <- cors[upper.tri(cors, diag = FALSE)]  # we're only interested in one of the off-diagonals, otherwise there'd be duplicates
cor.data <- as.data.frame(cor.data)  # that's how ggplot likes it
summary(cor.data)
colnames(cor.data) <- "pearson"
pearson.p <- function(r, n) {
  pofr <- ((1-r^2)^((n-4)/2))/beta(a = 1/2, b = (n-2)/2)
  return(pofr)
}
g <- NULL
g <- ggplot(data = cor.data, mapping = aes(x = pearson))
g <- g + xlim(-1,1)  # actual limits of pearsons r
g <- g + geom_histogram(mapping = aes(y = ..density..))
g <- g + stat_function(fun = pearson.p, colour = "red", args = list(n = nrow(m)))
g

至关重要的是pearson.p函数和最后一个 ggplot2 添加。

这是结果;正如人们所期望的那样完美匹配:

在此处输入图像描述

1个回答

作为一般评论,您的问题通常非常清晰且说明性很好,但往往过于解释您的主题(“Q 方法论”或其他任何内容),可能会在此过程中失去一些读者。

在这种情况下,您似乎在问:

样本的概率分布是多少(n=36) 两个不相关的高斯变量之间的皮尔逊相关系数?

答案很容易找到,例如 Wikipedia 关于 Pearson 相关系数的文章确切的分布可以写成任何n以及人口相关性的任何值ρ就超几何函数而言。这个公式很吓人,我不想在这里复制它。在你的情况下ρ=0它大大简化如下(参见同一篇 Wiki 文章):

p(r)=(1r2)(n4)/2Beta(1/2,(n2)/2).

在你的随机情况下36×1000矩阵n=36. 我们可以检查一下公式:

相关系数分布

这里蓝线显示了随机生成的相关矩阵的非对角元素的直方图,红线显示了上面的分布。合身是完美的。

请注意,分布可能呈现高斯分布,但它不能完全是高斯分布,因为它仅定义在[1,1]而正态分布有无限的支持。我用黑色虚线绘制了具有相同方差的正态分布;您可以看到它与红线非常相似,但在峰值处略高。


Matlab代码

n = 36;
p = 1000;

X = randn(n,p);
C = corr(X);
offDiagElements = C(logical(triu(C,1)));

figure
step = 0.01;
x = -1:step:1;
h = histc(offDiagElements, x);
stairs(x,h/sum(h)/step)
hold on

r = -1:0.01:1;
plot(r, 1/beta(1/2,(n-2)/2)*(1-r.^2).^((n-4)/2), 'r')

sigma2 = var(offDiagElements);
plot(r, 1/sqrt(sigma2*2*pi)*exp(-r.^2/(2*sigma2)), 'k--')

斯皮尔曼相关系数

我不知道有关样本 Spearman 相关性分布的理论结果。但在上面的模拟中,很容易将 Pearson 的相关性替换为 Spearman 的相关性:

C = corr(X, 'type', 'Spearman');

这似乎根本没有改变分布。

更新: @Glen_b 在聊天中指出“分布不能相同,因为 Spearman 的分布是离散的,而 Pearson 的分布是连续的”。这是真的,可以通过我的代码清楚地看到较小的值n. 奇怪的是,如果使用足够大的直方图箱以使离散性消失,则直方图开始与 Pearson 的直方图完美重叠。我不确定如何在数学上精确地表达这种关系。