Jak ustalić, czy dwie korelacje różnią się znacząco?

9

Chcę ustalić, który z dwóch zestawów danych (B1, B2) lepiej koreluje (gruszki r) z innym zestawem (A). Brakuje danych we wszystkich zestawach danych. Jak mogę ustalić, czy wynikowa korelacja jest znacząco inna, czy nie?

Np. Wartości 8426 są obecne zarówno w A, jak i B1, r = 0,74. 8798 jest obecnych zarówno w A, jak i B2, r = 0,72.

Myślałem, że to pytanie może pomóc, ale nie ma na nie odpowiedzi: Skąd wiedzieć, że jeden system jest znacznie lepszy od drugiego?

— szklanka
źródło

Czy możemy założyć duże

n

$n$ ?

— Firebug

1

@ Firebug n będzie zwykle wynosić od 7000 do 8760.

— Greenglass

6

Czasem może być w stanie osiągnąć to w regresji wielokrotnej, gdzie A jest DV B jest wynik ludzie mają w skali, a C jest kod manekin, który mówi, że jest albo B1 lub B2 lm(A~B+C+B*C). Termin interakcji, B*Cpowie ci, czy korelacje są różne, a proste nachylenia między A i B na obu poziomach C powie ci korelacje.

Nie jest jednak możliwe dopasowanie wszystkich rodzajów porównań między warunkami w tych ramach. cocorPakiet R jest bardzo przydatne i ma bardzo prosty wskaż i kliknij interfejs w internecie. Zauważ, że przy różnych brakujących danych nie masz próbek niezależnych ani zależnych. Użyłbym tutaj listowego usuwania, aby było to proste (a moc nie jest dla ciebie problemem).

— Mark White
źródło

2

Chociaż jest to najkrótsza odpowiedź, link do cocor skierował mnie w stronę potrzebnych informacji. Wielkie dzięki.

— greenglass

15

Och, moc bootstrapu. Spójrzmy na trzy wektory dla ilustracji: $A$ , $B_1$ i $B_2$ gdzie:

do o r (ZA, b_{1}) = 0,92

$Cor(A, B_1) = 0.92$

do o r (ZA, b_{2)}) = 0,86

$Cor(A, B_2) = 0.86$

Celem jest ustalenie, czy korelacja tych dwóch zestawów danych jest znacząco różna. Pobierając próbki bootstrap w taki sposób:

 B <- 10000
 cor1 <- cor2 <- rep(0, B)
 for(i in 1:B){
   samp <- sample(n, n, TRUE)  
   cor1[i] <- cor(A[samp], B1[samp])
   cor2[i] <- cor(A[samp], B2[samp])
 }

Możemy wykreślić rozkłady bootstrap dwóch korelacji:

Możemy również uzyskać 95% przedziały ufności dla $Cor(A, B_i)$ .

95% CI dla $Corr(A, B_1)$ :

(0,897, 0,947)

$(0.897, 0.947)$

95% CI dla $Corr(A, B_2)$ :

(0,810, 0,892)

$(0.810, 0.892)$

Fakt, że przedziały nie nakładają się (ledwo), daje pewne dowody na to, że zaobserwowana różnica korelacji próbek jest rzeczywiście istotna statystycznie.

Jak ameba wskazuje w komentarzach, bardziej „mocny” wynik wynika z uzyskania różnicy dla każdej próbki bootstrap.

95% CI dla różnicy między tymi dwoma wynosi:

(0,019, 0,108)

$(0.019, 0.108)$

Biorąc pod uwagę, że przedział (ledwo) wyklucza 0, mamy podobne dowody jak wcześniej.

Aby poradzić sobie z problemem braku danych, po prostu wybierz próbki bootstrapu z par zawartych w obu zestawach danych.

— knrumsey
źródło

7

Możesz poprawić tę procedurę, obliczając 95% CI dla cor1-cor2 i sprawdzając, czy zawiera zero, czy nie.

— ameba

1

Trafne spostrzeżenie. Dołączę to wkrótce.

— knrumsey

Fajnie, +1.

$\:\:$

— ameba

1

Dziękuję za tę odpowiedź. To było bardzo interesujące i pouczające, mimo że nie tego użyłem.

— klepsydra

7

Załóżmy, że transformacja Fishera: $r_1'=\tanh^{-1}(r_1)$ i $r_2'=\tanh^{-1} \left(r_2\right)$ . Lub w równoważny i być może jaśniejszy sposób (dzięki @dbwilson !), $r_1'={1\over2}\ln\left({1+r_1\over1-r_1}\right)$ i $r_2'={1\over2}\ln\left({1+r_2\over1-r_2}\right)$ .

Wynika z tego, że ze względu na fakt, że zmienne transformowane przez Fishera są teraz rozkładem normalnym, a suma rozkładów losowych rozkładów normalnych jest nadal rozkładem normalnym:

z = \frac{r_{1}^{'} - r_{2)}^{'}}{S.} \sim N. (0, 1)

$z={r_1'-r_2'\over S}\sim N(0,1)$ Z

S. = \sqrt{{S.}_{1}^{2)} + {S.}_{2)}^{2)}} = \sqrt{\frac{1}{n_{1} - 3)} + \frac{1}{n_{2)} - 3)}}

$S=\sqrt{S_1^2+S_2^2}=\sqrt{{1\over n_1-3}+{1\over n_2-3}}$

Więc testujesz hipotezę zerową $H_0:z=0$ poprzez uzyskanie $P(z\neq0)=2\cdot P(Z>|z|)$ .

W porównaniu do zwykłego $t$ -testuj, zauważ, że nie mogliśmy użyć $t$ -statystyki tak łatwo, zobacz jaki jest rozkład różnicy rozkładów dwóch t , więc należy wziąć pod uwagę stopnie swobody dostępne w obliczeniach, tzn. zakładamy $n$ wystarczająco duży, aby normalne przybliżenie mogło być rozsądnie odpowiednie $t$ Statystyka.

-

Po komentarzu @Josha możemy nieco uwzględnić możliwość współzależności między próbkami (pamiętaj, że obie korelacje zależą od rozkładu A). Bez zakładania niezależnych próbek i stosowania nierówności Cauchy'ego-Schwarza możemy uzyskać następującą górną granicę (patrz: Jak znaleźć standardowe odchylenie różnicy między dwoma średnimi? ):

S. \leq {S.}_{1} + {S.}_{2)}

$S\leq S_1+S_2$

S. \leq \sqrt{\frac{1}{n_{1} - 3)}} + \sqrt{\frac{1}{n_{2)} - 3)}}

$S\leq {\sqrt{1\over n_1-3}+\sqrt{1\over n_2-3}}$

— Firebug
źródło

2

To byłoby moje zalecenie, ale alternatywną formułą transformacji F Fishera jest z = 0,5 * ln ((1 + r) / (1-r)). Zrób to dla każdego ri postępuj jak wyżej.

— dbwilson

@dbwilson O tak (+1), są równoważne, dodam twoją sugestię, aby była bardziej zrozumiała dla szerszego grona odbiorców.

— Firebug

Czy ta formuła nie zakłada niezależności między

r_{1}

$r_1$ i

r_{2}

$r_2$ ? Myślę, że nie są ...

— Josh

6

Edytowane po pomocnych opiniach Marka White'a (dziękuję!)

Jedną z opcji jest obliczenie obu zależności (B1 z A i B2 z A) w jednym modelu, który również szacuje różnicę między nimi. Jest to łatwe do osiągnięcia dzięki regresji wielokrotnej . Uruchomiłbyś model z A jako zmienną zależną, a następnie jedną zmienną ciągłą ze wszystkimi wynikami dla B1 i B2, zmienną kategoryczną wskazującą, która to zmienna (B1 lub B2) i interakcję między nimi. R:

> set.seed(24601)
> 
> library(tidyverse)
> library(mvtnorm)
> cov <- matrix(c(1, .4, .16,.4, 1, .4, .16, .4, 1), ncol=3, byrow=TRUE)
> mydata <- rmvnorm(n=100, sigma = cov)
> colnames(mydata) = c("A", "B1", "B2")
> head(mydata)
              A         B1         B2
[1,] -0.1046382  0.6031253  0.5641158
[2,] -1.9303293 -0.7663828 -0.7921836
[3,]  0.1244192 -0.4413581 -1.2376256
[4,] -3.2822601 -1.2512055 -0.5586773
[5,] -0.9543368 -0.1743740  1.1884185
[6,] -0.4843183 -0.2612668 -0.7161938

Oto korelacje z danych, które wygenerowałem:

> cor(mydata)
           A        B1        B2
A  1.0000000 0.4726093 0.3043496
B1 0.4726093 1.0000000 0.3779376
B2 0.3043496 0.3779376 1.0000000
>

Zmiana formatu danych w celu zaspokojenia potrzeb modelu (przeformatowanie na „długie”):

> mydata <- as.data.frame(mydata) %>% 
+   gather("var", "value", B1, B2)
>

Oto model:

podsumowanie (lm (wartość A ~ * var, data = moje dane))

Call:
lm(formula = A ~ value * var, data = mydata)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.89310 -0.52638  0.02998  0.64424  2.85747 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.09699    0.09014  -1.076    0.283    
value        0.47445    0.09305   5.099 8.03e-07 ***
varB2       -0.10117    0.12711  -0.796    0.427    
value:varB2 -0.13256    0.13965  -0.949    0.344    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.891 on 196 degrees of freedom
Multiple R-squared:  0.158, Adjusted R-squared:  0.1451 
F-statistic: 12.26 on 3 and 196 DF,  p-value: 2.194e-07

Wyniki tutaj (z moich opracowanych danych) sugerują, że istnieje znaczący związek między B1 i A (test współczynnika „wartości”, ponieważ B1 jest grupą odniesienia dla współczynnika „var”), ale że różnica związek między B1 z A i B2 z A nie jest istotny (test współczynnika „wartość: varB2”).

Jeśli wolisz myśleć w kategoriach korelacji niż współczynników regresji, po prostu ujednolicaj wszystkie swoje zmienne (A, B1 i B2) przed uruchomieniem modelu, a otrzymane współczynniki regresji zostaną ustandaryzowane (nie do końca to samo korelacja zerowego rzędu, ale znacznie bliższa interpretacji).

Zauważ również, że ograniczy to analizę tylko do przypadków, które mają zarówno B1, jak i B2 ( usunięcie listowe ). Tak długo, jak pozostawia Ci to wystarczającą ilość danych, abyś nie był niedostatecznie zasilony, i tak długo, jak brakujące dane są losowo tracone (lub wystarczająco mała część wszystkich danych, aby nie mieć większego znaczenia, nawet jeśli brakuje ich losowo), to w porządku.

Ograniczenie analizy do tego samego zestawu danych w celu oszacowania efektów zarówno dla B1, jak i B2 (zamiast używania nieco innych zestawów danych w oparciu o różne wzorce braków) ma tę zaletę, że nieco interpretuje różnicę między korelacjami bardziej proste. Jeśli obliczysz korelacje osobno dla każdej z nich, a następnie przetestujesz różnicę między nimi , napotkasz problem polegający na tym, że dane leżące u ich podstaw są nieco inne w każdym przypadku --- każda widoczna różnica może wynikać z różnic w próbkach tak samo jak różnic w rzeczywistych relacjach między zmiennymi.

— Rose Hartman
źródło

2

Nie jest to przypadek, że lm(A~B1*B2)będzie sprawdzić, czy korelacja pomiędzy B1i A zależy od czyjegoś B2wynik ? Ten termin interakcji nie sprawdza, czy korelacje są różne; sprawdza, czy dwa predyktory oddziałują na siebie. Możesz stworzyć sztuczny kod, Cktóry koduje, czy skala dla Bjest, B1czy nie B2. To by powiedziało ci, że korelacja pomiędzy Bi Azależy od tego, czy jest, B1czyli od B2tego, czy korelacje są różne.

— Mark White

1

@MarkWhite O rany, masz całkowitą rację! Dzięki za złapanie tego. Yikes! Będę edytować, aby to naprawić.

— Rose Hartman