R

兩個不相關的正態變量之間的樣本相關係數分佈如何?

  • January 22, 2016

我想比較觀察到的雙變量(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 添加。

這是結果;正如人們所期望的那樣完美匹配:

在此處輸入圖像描述

作為一般評論,您的問題通常非常清晰且說明良好,但往往過於解釋您的主題(“Q 方法論”或其他任何內容),可能會在此過程中失去一些讀者。

在這種情況下,您似乎在問:

樣本的概率分佈是多少( $ n=36 $ ) 兩個不相關的高斯變量之間的皮爾遜相關係數?

答案很容易找到,例如 Wikipedia 關於 Pearson 相關係數的文章。確切的分佈可以寫成任何 $ n $ 以及人口相關性的任何值 $ \rho $ 就超幾何函數而言。這個公式很嚇人,我不想在這裡複製它。在你的情況下 $ \rho=0 $ 它大大簡化如下(參見同一篇 Wiki 文章):

$$ p(r) = \frac{(1-r^2)^{(n-4)/2}}{\operatorname{Beta}(1/2, (n-2)/2)}. $$

在你的隨機情況下 $ 36\times 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 的直方圖完美重疊。我不確定如何在數學上精確地表達這種關係。

引用自:https://stats.stackexchange.com/questions/191937

comments powered by Disqus