Wie ist die Verteilung der Stichprobenkorrelationskoeffizienten zwischen zwei nicht korrelierten Normalvariablen?

8

Ich möchte beobachtete bivariate (Pearson's und Spearman's ) Korrelationskoeffizienten mit dem vergleichen, was von zufälligen Daten erwartet wird.ρρ

Angenommen, wir messen beispielsweise 36 Fälle über sehr viele Variablen (1000). (Ich weiß, dass dies seltsam ist, es wird Q-Methodik genannt . Nehmen wir weiter an, dass jede der Variablen (streng) normal über die Fälle verteilt ist . (Wieder sehr seltsam, aber wahr, weil Personen als Personenvariablen die Reihenfolge der Artikelfälle unter a ordnen Normalverteilung.)

Also, wenn die Menschen zufällig sortiert , sollten wir erhalten:

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

Nun - da dies eine Q-Methode ist - korrelieren wir alle Personenvariablen :

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

Dann versuchen wir, dies aufzuzeichnen und die Verteilung des Pearson-Korrelationskoeffizienten in Zufallsdaten zu überlagern , die eigentlich den beobachteten Korrelationen in unseren gefälschten Daten ziemlich nahe kommen sollten:

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

Das gibt:

Dichtediagramm

Die überlagerte Kurve ist eindeutig falsch. (Beachten Sie auch, dass die Dichten der y-Achse zwar ungerade sind, aber tatsächlich korrekt sind : Da die x-Werte so klein sind, summiert sich die Fläche auf eins).

Ich erinnere mich (vage), dass die T-Verteilung in diesem Zusammenhang relevant ist, aber ich kann mich nicht darum kümmern, wie man sie richtig parametrisiert. Sind die Freiheitsgrade insbesondere durch die Anzahl der Korrelationen (1000 ^ 2 / 2-500) oder die Anzahl der Beobachtungen, auf denen diese Korrelationen beruhen, gegeben (36)?

In beiden Fällen ist die oben überlagerte Kurve eindeutig falsch.

Ich bin auch verwirrt, weil die Wahrscheinlichkeitsverteilung von Pearsons r begrenzt werden müsste (es gibt keine Werte über (-) 1 hinaus) - aber die t-Verteilung ist nicht begrenzt.

Welche Distribution beschreibt Pearson's in diesem Fall?ρ


Bonus:

Die obigen Daten sind tatsächlich idealisiert: In meiner realen Q-Studie haben Personenvariablen tatsächlich nur sehr wenige Spalten unter einer Normalverteilung, in die ihre Artikelfälle wie folgt sortiert werden können:

q-sort

Tatsächlich handelt es sich bei Personenvariablen tatsächlich um Rangfolge von Artikelfällen , sodass Pearson's nicht anwendbar ist. Als grobe Lösung habe ich mich stattdessen für Spearman's entschieden. Ist die Wahrscheinlichkeitsverteilung für Spearman's ?ρρ


Update : Wenn jemand interessiert ist, ist hier der R-Code, um die fantastische Antwort von @ amoeba unten zu implementieren:

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

Entscheidend sind die pearson.pFunktion und die letzte Ergänzung von ggplot2.

Hier ist das Ergebnis; passt perfekt, wie man erwarten würde:

Geben Sie hier die Bildbeschreibung ein

maxheld
quelle
Auch wie immer wäre es fantastisch, wenn jemand ein "qmethod" -Tag spenden könnte.
Maxheld
3
Eine bestimmte Transformation der empirischen Korrelation hat die Verteilung, nicht die Korrelation selbst. Die Antwort unten von Amöbe nagelt es. t
StasK

Antworten:

11

Im Allgemeinen sind Ihre Fragen in der Regel sehr klar und gut illustriert, gehen jedoch häufig zu weit in die Erläuterung Ihres Themas ("Q-Methodik" oder was auch immer) und verlieren möglicherweise einige Leser auf dem Weg.

In diesem Fall scheinen Sie zu fragen:

Wie ist die Wahrscheinlichkeitsverteilung des Pearson-Korrelationskoeffizienten der Stichprobe ( ) zwischen zwei nicht korrelierten Gaußschen Variablen?n=36

Die Antwort ist leicht zu finden, z. B. in Wikipedia's Artikel über den Pearson-Korrelationskoeffizienten . Die genaue Verteilung kann für jede geschrieben werden und jeden Wert von Bevölkerungskorrelations in Bezug auf die Funktion hypergeometric. Die Formel ist beängstigend und ich möchte sie hier nicht kopieren. In Ihrem Fall von vereinfacht sich dies wie folgt erheblich (siehe denselben Wiki-Artikel):nρρ=0

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

In Ihrem Fall einer zufälligen Matrix ist . Wir können die Formel überprüfen:36×1000n=36

Verteilung der Korrelationskoeffizienten

Hier zeigt die blaue Linie das Histogramm der nicht diagonalen Elemente einer zufällig erzeugten Korrelationsmatrix und die rote Linie zeigt die obige Verteilung. Die Passform ist perfekt.

Beachten Sie, dass die Verteilung möglicherweise Gaußsch erscheint, aber nicht genau Gaußsch sein kann, da sie nur in während die Normalverteilung unendlich unterstützt wird. Ich habe die Normalverteilung mit der gleichen Varianz mit einer schwarzen gestrichelten Linie aufgetragen. Sie können sehen, dass es der roten Linie ziemlich ähnlich ist, aber in der Spitze etwas höher ist.[1,1]


Matlab-Code

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-Korrelationskoeffizient

Mir sind keine theoretischen Ergebnisse zur Verteilung der Spearman-Korrelationen bekannt. In der obigen Simulation ist es jedoch sehr einfach, die Korrelationen von Pearson durch die von Spearman zu ersetzen:

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

und dies scheint die Verteilung überhaupt nicht zu ändern.

Update: @Glen_b wies im Chat darauf hin, dass "die Verteilung nicht dieselbe sein kann, da die Verteilung für den Spearman diskret ist, während die für den Pearson kontinuierlich ist". Dies ist wahr und kann mit meinem Code für kleinere Werte von deutlich gesehen werden . Seltsamerweise überlappt sich das Histogramm perfekt mit dem von Pearson, wenn man einen ausreichend großen Histogrammbehälter verwendet, damit die Diskretion verschwindet. Ich bin mir nicht sicher, wie ich diese Beziehung mathematisch genau formulieren soll.n

Amöbe
quelle
fantastisch, danke für den Zusatz von Spearman. Wird in Zukunft mit der Q-Methodik vorsichtiger sein; Ich kämpfe immer noch darum herauszufinden, wann dieser Q-Move der invertierten PCA principal(t(data.matrix)wichtig ist und wann nicht. Angehängte R-Lösung in obiger Frage; Lassen Sie mich wissen, ob Sie das lieber in Ihrer Antwort haben möchten.
Maxheld
3
Um ganz klar zu sein: Ich habe nichts gegen die Q-Methodik oder damit zusammenhängende Fragen (ich habe auch nichts dafür, ich bin nicht vertraut genug). Aber diese Frage hätte als Zweiliner formuliert werden können; Es ist sehr gut möglich, dass wenn Sie es als solches gepostet hätten, es stärker bewertet und schneller beantwortet worden wäre. Man kann nie vollständig vorhersagen, welche Fragen passieren (oder nicht), um Aufmerksamkeit zu erregen, aber eine lange, komplizierte Formulierung ist meistens nachteilig.
Amöbe
@amoeba ist es möglich, in einem einfachen linearen Modell, in dem die Varianzen von X und Y bekannt sind, in / aus der Korrelation eines Regressionskoeffizienten umzuwandeln?
Sammosummo