Jak badać zależności pomiędzy parą zmiennych?

Testów weryfikujących strukturę zależności pomiędzy parą zmiennych jest co niemiara. W języku polskim całkiem bogata kolekcja testów jest przedstawiona w książce Ryszarda Magiery ,,Modele i metody statystyki Matematycznej''.

Zamiast jednak opisywać wszystkie testy statystyczne (co jest niemożliwe) poniżej skupimy się na przedstawieniu tych najbardziej popularnych, tworzących swego rodzaju szkielet, który można rozbudowywać w różnych kierunkach.

Jak weryfikować niezależność dwóch zmiennych jakościowych?

Problem: Mamy dwie zmienne jakościowe (np. kolor oczu i kolor włosów), ciekawi jesteśmy czy jest pomiędzy nimi zależność.

Model: Przyjmijmy, że obserwujemy dwuwymiarową zmienną losową (X,Y)(X,Y). XX przyjmuje wartości ze zbioru (x1,...,xr)(x_1, ..., x_r) a YY przyjmuje wartości ze zbioru (y1,...,ys)(y_1, ..., y_s).

Brzegowy rozkład zmiennej X określmy jako P(X=xi)=pi.P(X = x_i) = p_{i.} a zmiennej YY jako P(Y=Yj)=p.jP(Y = Y_j) = p_{.j}. Łączny rozkład zmiennych (X,Y)(X,Y) określimy jako P(X=xi,Y=yj)=pijP(X=x_i, Y=y_j) = p_{ij}.

Hipoteza: Określmy hipotezę zerową (niezależności), dla każdego i,ji,j

H0:pij=pi.p.j H_0: p_{ij} = p_{i.}p_{.j}

Hipoteza alternatywna jest taka, że dla dowolnego i,ji,j ta równość nie zachodzi.

Statystyka testowa:

Statystyka testowa jest oparta o tablicę kontyngencji (tablicę zliczeń).

Y \ X x1x_1 x2x_2 xsx_s
y1y_1 n11n_{11} n12n_{12} ... n1sn_{1s}
... ... ...
yry_r nr1n_{r1} nr2n_{r2} ... nrsn_{rs}

Tablica kontyngencji opisuje obserwowane w danych liczebności wystąpienia poszczególnych kombinacji wartości. Statystyka testowa porównuje te wartości z wartościami oczekiwanymi, przy założeniu prawdziwej hipotezy zerowej.

Ogólna postać statystyki testowej w teście χ2\chi^2 to

T=i(OiEi)2Ei T = \sum_i \frac{(O_i - E_i)^2}{E_i}

gdzie OiO_i oznacza obserwowane liczebności, EiE_i oczekiwane liczebności dla i-tej komórki.

Ta generyczna postać statystyki testowej, w tym konkretnym przypadku przyjmuje postać (indeksy z kropkami to sumy po wierszach/kolumnach)

T=i=1rj=1s(nijni.n.j/n)2ni.n.j/n.. T = \sum_{i=1}^r \sum_{j=1}^s \frac{(n_{ij} - n_{i.}n_{.j}/n)^2}{n_{i.}n_{.j}/n_{..}}

Ta statystyka testowa ma asymptotyczny rozkład χ(r1)(s1)2\chi^2_{(r-1)(s-1)}.

Za obszar odrzucenia przyjmuje się przedział postaci [c,][c, \infty]. Przykładowa gęstość rozkładu χ42\chi^2_4 z zaznaczonym kwartylem rzędu 0.950.95.

library(ggplot2)
x <- seq(0,15,0.01)
(q <- qchisq(0.95, 4))
## [1] 9.487729
df <- data.frame(x, d = dchisq(x,4), rejection = x>q )
ggplot(df, aes(x, y=d, fill=rejection)) + geom_area() + theme(legend.position="none") + scale_fill_manual(values=c("grey","red3"))

plot of chunk unnamed-chunk-1

Przykład

Rozważmy następującą tabelę liczebności. Opisuje ona występowanie różnych kolorów oczu / włosów.

(tab <- HairEyeColor[,,1])
##        Eye
## Hair    Brown Blue Hazel Green
##   Black    32   11    10     3
##   Brown    53   50    25    15
##   Red      10   10     7     7
##   Blond     3   30     5     8
# archivist::aread("pbiecek/Przewodnik/arepo/a3666b4084daa5db9251dc36e3286298")

Aby przeprowadzić test χ2\chi^2 można wykorzystać funkcję chisq.test(). Wyznacza ona zarówno macierz oczekiwanych częstości, statystykę testową jak i wartość p.

wynik <- chisq.test(tab)

wynik$p.value
## [1] 4.447279e-06
wynik$statistic
## X-squared 
##  41.28029
wynik$expected
##        Eye
## Hair       Brown     Blue     Hazel     Green
##   Black 19.67025 20.27240  9.433692  6.623656
##   Brown 50.22939 51.76703 24.089606 16.913978
##   Red   11.94265 12.30824  5.727599  4.021505
##   Blond 16.15771 16.65233  7.749104  5.440860
wynik$observed
##        Eye
## Hair    Brown Blue Hazel Green
##   Black    32   11    10     3
##   Brown    53   50    25    15
##   Red      10   10     7     7
##   Blond     3   30     5     8

Jak weryfikować niezależność dwóch zmiennych binarnych?

Specyficzną wersją testu na niezależność dwóch zmiennych jakościowych jest test dla dwóch zmiennych binarnych. Zamiast wykorzystywać w tym przypadku asymptotyczny rozkład statystyki testowej można badać dokładny rozkład statystyki testowej. Stąd też nazwa testu - dokładny test Fishera.

Statystyka testowa jest oparta o tablicę kontyngencji

Y \ X x1x_1 x2x_2
y1y_1 n11n_{11} n12n_{12}
y2y_2 n21n_{21} n22n_{22}

o rozkładzie hipergeometrycznym.

Jeżeli badamy zależnośc pomiędzy parą zmiennych binarnych to zalecane jest użycie tego testu. Umożliwia on również weryfikowanie hipotez kierunkowych (a więc częstsze/rzadsze niż przypadkowe współwystępowanie x2y2x_2y_2).

Przykład

Ograniczmy nadanie kolorów oczu do niebieskie/brązowe a włosów do czarne / blond.

tab22 <- tab[c(1,4),c(1,2)]

fisher.test(tab22)
## 
##     Fisher's Exact Test for Count Data
## 
## data:  tab22
## p-value = 7.391e-09
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##    6.673302 168.168624
## sample estimates:
## odds ratio 
##   27.39659

Jak weryfikować niezależność dwóch zmiennych ilościowych?

Dla zmiennych ilościowych zależność może przybierać bardzo różną postać. Możemy obserwować zależność w kwadratach (czyli korelacje wartości skranych bez względu na ich znak), uwikłanie zmiennych, wiele możliwych odstępst od niezależności.

Z powodu łatwości interpretacji, najczęściej w pierwszym kroku interesują nas dwa rodzaje zależności: liniowa oraz monotoniczna. Do ich badania najczęściej wykorzystuje się testy na współczynnik korelacji Pearsona i Spearmana.

Dwuwymiarowy rozkład normalny

Model: Przyjmijmy, że obserwujemy dwuwymiarową zmienną losową z dwuwymiarowego rozkładu normalnego (X,Y)N(μ,Σ)(X,Y) \sim \mathcal N(\mu, \Sigma). Gdzie Σ\Sigma to macierz kowariancji, element poza przekątną oznaczny przez σ12\sigma_{12}.

Hipoteza:

H0:σ12=0. H_0: \sigma_{12} = 0. HA:σ12<>0. H_A: \sigma_{12} <> 0.

Statystyka testowa:

Statystyka testowa oparta jest o współczynnik korelacji Pearsona

ρ^=i(xix¯)(yiy¯)i(xix¯)2i(yiy¯)2 \hat \rho = \frac{\sum_i (x_{i} - \bar x)(y_i - \bar y)}{\sqrt{\sum_i (x_{i} - \bar x)^2\sum_i (y_{i} - \bar y)^2}}

Rozkład statystyki ρ^\hat \rho jest określony na odcinku [-1,1]. Aby ułatwić jej analizę stosowane jest następujące przekształcenie

T=n2ρ^1ρ^2. T = \sqrt{n-2}\frac{\hat \rho}{\sqrt{1 - \hat \rho^2}}.

Po takim przekształceniu statystyka TT ma rozkład tn2t_{n-2} i w oparciu o niego konstruowany jest obszar krytyczny. Dla dwustronnej hipotezy alternatywnej jest to obszar (,c][c,)(-\infty, -c] \cup [c, \infty).

x <- seq(-5,5, .01)
(q <- qt(0.975, 5))
## [1] 2.570582
df <- data.frame(x, d = dt(x,5), rejection = abs(x)>q )
ggplot(df, aes(x, y=d, fill=rejection)) + geom_bar(stat="identity") + theme(legend.position="none") + scale_fill_manual(values=c("grey","red3"))

plot of chunk unnamed-chunk-5

Czasem weryfikowana jest też inna hipoteza zerowa. Jeżeli chceny zbadać, czy korelacja jest istotnie różna (np. istotnie większa) od określonej wartości ρ0\rho_0 to interesuje nas raczej hipoteza.

Hipoteza:

H0:σ12=ρ0. H_0: \sigma_{12} = \rho_0. HA:σ12<>ρ0. H_A: \sigma_{12} <> \rho_0.

Statystyka testowa:

W tym przypadku stosuje się inną transformację, tzw. transformację Fishera.

U=1/2ln(1+ρ^1ρ^). U = 1/2 \ln(\frac{1+\hat \rho}{1-\hat \rho}).

Przy prawdziwej hipotezie zerowej ta statystyka ma asymptotycznie rozkład normalny N(1/2ln(1+ρ01ρ0)+ρ0/(2(n1)),1/(n3))\mathcal N (1/2\ln(\frac{1+\rho_0}{1-\rho_0}) + \rho_0/(2(n-1)), 1/(n-3)).

Znając rozkład możemy już w prosty sposób zbudować test statystyczny. Tę transformację można wykorzystać również do konstrukcji testu dla równości dwóch współczynników korelacji.

Korelacja rang

Założenie o dwuwymiarowym rozkładzie normalnym jest silnie ograniczające. Co prawda test Pearsona stosuje się nawet jeżeli zmienna nie ma rozkładu normalnego, ale jedynie do niego zbliżony. Wciąż jednak tego typu test wykrywa jedynie liniowe zależnosci.

Dlatego często stosowanym testem dla zbioru hipotez jest test korelacji Spearmana. Jest ona w stanie identyfikować zależności monotoniczne. Ideę testu Spearmana można streścić w określeniu: badanie korelacji rang.

Model: Przyjmijmy, że obserwujemy dwuwymiarową zmienną losową (X,Y)(X,Y) o rozkładzie ciągłym.

Oznaczmy dodatkowo ri=rank(Xi)r_i = rank(X_i), si=rank(Yi)s_i = rank(Y_i).

Hipoteza:

H0:ρ(ri,si)=0. H_0: \rho(r_i, s_i) = 0. HA:ρ(ri,si)<>0. H_A: \rho(r_i, s_i) <> 0.

Statystyka testowa:

Statystyką testową jest korelacja Pearsona ale liczona dla rang, a nie oryginalnych obserwacji. Ponieważ średnia ranga to (n+1)/2(n+1)/2 więc otrzymujemy

ρ=i(ri(n+1)/2)(si(n+1)/2)i(ri(n+1)/2)2i(si(n+1)/2)2. \rho = \frac{\sum_i (r_{i} - (n+1)/2)(s_i - (n+1)/2)}{\sqrt{\sum_i (r_{i} - (n+1)/2)^2\sum_i (s_{i} - (n+1)/2)^2}}.

Po prostych przekształceniach otrzymujemy

ρ=16(risi)2n(n21). \rho = 1 - \frac{6 \sum(r_i - s_i)^2}{n(n^2-1)}.

Rozkład tej statystyki można tablicować dla małych nn. Asymptotycznie ma ona rozkład normalny z wariancją 1/(n1)1/(n-1). Ale w implementacji najczęściej stosuje się podobną transformację co w przypadku testu Pearsona, czyli

T=n2ρ^1ρ^2. T = \sqrt{n-2}\frac{\hat \rho}{\sqrt{1-\hat \rho^2}}.

Asymptotycznie ta statystyka ma rozkład tn2t_{n-2}.

A co jeżeli interesuje mnie inna zależność?

Opisane powyżej testy badają dwa rodzaje zależności - liniowy i monotoniczny.

Jeżeli interesuje nas inna klasa zależności to możliwy wybór jest albo przez badanie funkcji łączących (tzw. kopule) albo przez badanie zmiennych jakościowych. Każdą zmienną ilosciową możemy zdyskretyzować, dzieląc ją na pewną liczbę podprzedziałów i analizując zależność pomiedzy przedziałami (testem χ2\chi^2 lub analizą korespondencji).

Przykład

Korzystając ze zbioru danych koty_ptaki sprawdzimy czy jest zależność pomiędzy długością kota a jego wagą.

library(PogromcyDanych)
head(koty_ptaki[,1:3])
##   gatunek waga dlugosc
## 1  Tygrys  300     2.5
## 2     Lew  200     2.0
## 3  Jaguar  100     1.7
## 4    Puma   80     1.7
## 5 Leopard   70     1.4
## 6  Gepard   60     1.4
ggplot(koty_ptaki, aes(waga, dlugosc)) + geom_point(size=3) + geom_smooth(method="lm", color="grey", se=FALSE)

plot of chunk z1

Korelację Pearsona (liniową) badalibyśmy testem cor.test()

cor.test(koty_ptaki$waga, koty_ptaki$dlugosc)
## 
##     Pearson's product-moment correlation
## 
## data:  koty_ptaki$waga and koty_ptaki$dlugosc
## t = 6.4499, df = 11, p-value = 4.743e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6633850 0.9666441
## sample estimates:
##       cor 
## 0.8893128

Korelację Spearmana (monotoniczną) badalibyśmy testem cor.test(). Współczynnik korelacji rang jest prawie równy 1.

cor.test(koty_ptaki$waga, koty_ptaki$dlugosc, method="spearman")
## Warning in cor.test.default(koty_ptaki$waga, koty_ptaki$dlugosc, method =
## "spearman"): Cannot compute exact p-value with ties
## 
##     Spearman's rank correlation rho
## 
## data:  koty_ptaki$waga and koty_ptaki$dlugosc
## S = 4.5155, p-value = 3.401e-10
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.9875947

Po dyskretyzacji badalibyśmy np. następującą tabelę (widzimy silną zależność nawet bez testu).

table(cut(koty_ptaki$waga, c(0,20,110,400)), 
      cut(koty_ptaki$dlugosc, c(0, 1, 1.9, 3)))
##            
##             (0,1] (1,1.9] (1.9,3]
##   (0,20]        5       0       0
##   (20,110]      0       5       0
##   (110,400]     0       0       3