我想比较观察到的双变量(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 研究中,人变量实际上在正态分布下只有很少的列可以将他们的项目案例分类,如下所示:
实际上,人员变量实际上是按等级排序的项目案例,因此 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 添加。
这是结果;正如人们所期望的那样完美匹配: