- Współczynnik korelacji Pearsona
- Test \(X^2\) Pearsona dla macierzy kontyngencji
- Testy normalności rozkładu
- Test \(X^2\)
- Test Kołmogorowa-Smirnowa
Współczynnik korelacji Pearsona
W celu wyznaczenia istotności zależności pomiędzy dwiema zmiennym często stosuje się współczynnik korelacji Pearsona (Pearson product-moment correlation coefficient), zdefiniowany (populacyjnie) dla dwóch zmiennych losowych \(X\) i \(Y\) jako
\(r_{XY} = \frac{\mathrm{cov}(X,Y)}{\sigma_X \sigma_Y}\),
natomiast w przypadku estymatora próbkowego dla serii \(n\) danych \(x_i\) oraz \(y_i\) jako
\(r = \frac{\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{\sqrt{\sum_{i=1}^n(x_i - \bar{x})^2}\sqrt{\sum_{i=1}^n(y_i - \bar{y})^2}}\).
Współczynnik jest zaimplementowany w R jako funkcja cor().
# PRZYKŁAD 5.1
x <- c(1,2,3,4,5)
y <- c(2,4,1,2,2)
cor(x, y)## [1] -0.2886751W przypadku, gdy mamy do czynienia z korelacją rang, używa się współczynnika Spearmana cor(method="spearman"), który jest de facto współcznnikiem Pearsona dla danych przekształconych w rangi.
# PRZYKŁAD 5.2
cor(x, y, method="spearman")## [1] -0.2236068x; rank(x)## [1] 1 2 3 4 5## [1] 1 2 3 4 5y; rank(y)## [1] 2 4 1 2 2## [1] 3 5 1 3 3cor(rank(x), rank(y))## [1] -0.2236068W przypadku, gdy do funkcji cor() prześlemy ramkę danych, zostaną wyliczone współczynniki korelacji pomiędzy każdą parą zmiennych i przedstawione w posteci macierzy
# PRZYKŁAD 5.3
z <- c(0, 1, 1, 0, 1)
df <- data.frame(x, y, z)
cor(df)## x y z
## x 1.0000000 -0.2886751 0.2886751
## y -0.2886751 1.0000000 0.1666667
## z 0.2886751 0.1666667 1.0000000Warto tu wspomnieć o ciekawym pakiecie corrplot, który udostępnia eponimiczną funkcję do wizualizacji macierzy korelacji. Pakiet należy uprzednio zainstalować za pomocą komendy
install.packages("corrplot")a następnie załadować za pomocą library() i wywołać funkcję corrplot() na utworzonej macierzy korelacji
# PRZYKŁAD 5.4
library(corrplot)## corrplot 0.84 loadedcorrplot(cor(df))
W celu otrzymania pełnych informacji, dotyczących parametrów współczynnika korelacji należy użyć funkcji cor.test()
# PRZYKŁAD 5.5
cor.test(x, y)##
## Pearson's product-moment correlation
##
## data: x and y
## t = -0.52223, df = 3, p-value = 0.6376
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.9332529 0.7964337
## sample estimates:
## cor
## -0.2886751Należy zwrócić tu uwagę na p-wartość (p-value) statystyki. Jest ona bardzo wysoka, zatem nie ma podstaw do odrzucenia hipotezy zerowej, że korelacja pomiędzy zmiennymi jest równa zero.
Test \(X^2\) Pearsona dla macierzy kontyngencji
W przypadku, gdy mamy do czynienia z dwiema zmiennymi jakościowymi i chcemy sprawdzić czy występują one niezależnie, możemy wykonać test \(X^2\) Pearsona, wywoływany funkcją chisq.test()
# PRZYKŁAD 5.6
x <- c(rep("fizyk", 55), rep("humanista", 149))
y <- c(rep("nie pracuje", 7),rep("pracuje", 48),rep("nie pracuje", 45),rep("pracuje", 104))
T <- table(x, y)
T## y
## x nie pracuje pracuje
## fizyk 7 48
## humanista 45 104chisq.test(T)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: T
## X-squared = 5.5711, df = 1, p-value = 0.01826Niska p-wartość wskazuje, że należy odrzucić hipotezę zerową o niezależności zmiennych. Można sprawdzić jak wyglądałaby macierz wartości oczekiwanych (i prawdopodobieństwa) w przypadku zmiennych niezależnych i sprawdzić, czy jest ona równa iloczynom rozkładów brzegowych.
# PRZYKŁAD 5.7
Te <- chisq.test(T)$expected
Te## y
## x nie pracuje pracuje
## fizyk 14.01961 40.98039
## humanista 37.98039 111.01961Te / sum(Te)## y
## x nie pracuje pracuje
## fizyk 0.06872357 0.2008843
## humanista 0.18617839 0.5442138Powyższa tabela mówi o tym, jakie powinny być prawdopodobieństwa wystąpienia danej pary, jeśli zadarzenia (np. fizyk-pracuje etc) byłyby niezależne, czyli występowały losowo. Wartości te można wyznaczyć z naszej tabeli \(T\) obliczając rozkłady brzegowe
| nie pracuje | pracuje | suma | |
|---|---|---|---|
| fizyk | 7 | 48 | 55 |
| humanista | 45 | 104 | 149 |
| suma | 52 | 159 | 204 |
# PRZYKŁAD 5.8
Tx <- apply(T, 1, sum) / sum(T)
Ty <- apply(T, 2, sum) / sum(T)
Tx %o% Ty## nie pracuje pracuje
## fizyk 0.06872357 0.2008843
## humanista 0.18617839 0.5442138Natomiast faktyczne wartości prawdopodobieństw, otrzymane po prostu poprzez podzielenie poszczególnych zliczeń przez całkowitą liczbę przypadków są zupełnie inne, co potwierdza wynik testu.
# PRZYKŁAD 5.9
T / sum(T)## y
## x nie pracuje pracuje
## fizyk 0.03431373 0.23529412
## humanista 0.22058824 0.50980392Testy normalności rozkładu
Powszechnie używanymi testami na to, czy dane pochodzą z rozkładu normalnego są (między innymi) test Shapiro-Wilka oraz test Wilcoxona. W R realizowane są one poprzez funkcje shapiro.test() oraz wilcox.test().
# PRZYKŁAD 5.10
shapiro.test(rnorm(100))##
## Shapiro-Wilk normality test
##
## data: rnorm(100)
## W = 0.99269, p-value = 0.8688shapiro.test(rnorm(100, 0, 5))##
## Shapiro-Wilk normality test
##
## data: rnorm(100, 0, 5)
## W = 0.99038, p-value = 0.695shapiro.test(runif(100, 0, 1))##
## Shapiro-Wilk normality test
##
## data: runif(100, 0, 1)
## W = 0.94216, p-value = 0.0002619wilcox.test(rnorm(100))##
## Wilcoxon signed rank test with continuity correction
##
## data: rnorm(100)
## V = 2555, p-value = 0.9192
## alternative hypothesis: true location is not equal to 0wilcox.test(rnorm(100, 1, 2), mu = 1)##
## Wilcoxon signed rank test with continuity correction
##
## data: rnorm(100, 1, 2)
## V = 2264, p-value = 0.3704
## alternative hypothesis: true location is not equal to 1wilcox.test(rnorm(100, 1, 2), mu = 2)##
## Wilcoxon signed rank test with continuity correction
##
## data: rnorm(100, 1, 2)
## V = 1390, p-value = 9.588e-05
## alternative hypothesis: true location is not equal to 2Test \(X^2\)
Ten test jest chyba najczęściej wykorzystywany do sprawdzania zgodności danych eksperymentalnych z teoretycznymi. W przypadku R realizuje go funkcja chisq.test(x, p), w której jako pierwszy argument podajemy liczbę zliczeń (histogram), a jako \(p\) teoretyczne wartości prawdopodobieństw. Hipoteza zerowa \(H_0\) w tym przypadku jest taka, że liczba zliczeń podzielona przez całkowitą liczbę zliczeń jest równa teoretycznym wartościom prawdopodobieństwa
# PRZYKŁAD 5.11
x <- c(100,50,20,10,5)
p1 <- c(0.5, 0.2, 0.2, 0.05, 0.05)
chisq.test(x, p = p1)##
## Chi-squared test for given probabilities
##
## data: x
## X-squared = 15, df = 4, p-value = 0.004701p2 <- c(0.5, 0.3, 0.1, 0.05, 0.05)
chisq.test(x, p = p2)##
## Chi-squared test for given probabilities
##
## data: x
## X-squared = 3.2883, df = 4, p-value = 0.5108p <- x/sum(x)
chisq.test(x, p = p)##
## Chi-squared test for given probabilities
##
## data: x
## X-squared = 0, df = 4, p-value = 1plot(p)
lines(p1, col="red")
lines(p2, col="blue")
UWAGA 1: jest to test do rozkładów dyskretnych.
UWAGA 2: prawdopodobieństwa muszą się sumować do 1, w przeciwnym wypadku otrzymamy błąd, możliwe jest przeskalowanie poprzez parametr rescale.p=TRUE)
# PRZYKŁAD 5.12
h <- hist(rpois(200, lam = 3), breaks = (0:15)-0.5)
chisq.test(h$counts, p = dpois(h$mids, lam = 3))## Error in chisq.test(h$counts, p = dpois(h$mids, lam = 3)): probabilities must sum to 1.chisq.test(h$counts, p = dpois(h$mids, lam = 3), rescale.p = TRUE)## Warning in chisq.test(h$counts, p = dpois(h$mids, lam = 3), rescale.p =
## TRUE): Chi-squared approximation may be incorrect##
## Chi-squared test for given probabilities
##
## data: h$counts
## X-squared = 8.4756, df = 14, p-value = 0.8631plot(h$mids, h$density)
lines(h$mids, dpois(h$mids, lam = 3))
lines(h$mids, dpois(h$mids, lam = 3), col="red", lwd=2)
Test Kołmogorowa-Smirnowa
W tym przypadku (funkcja ks.test()) możemy:
- stosować test zarówno do rozkładów dyskretnych jak i ciągłych,
- porównywać dwa rozkłady eksperymentalne albo rozkład eksperymentalny z rozkładem teoretycznym.
Test operuje na tzw. statystyce odstępu, badając, czy największa odległość pomiędzy dwiema dystybuantami mieści się w zakresie ww. statystyki. Hipoteza zerowa zakłada, że obie próbki (lub próba i rozkład) pochodzą z tego samego rozkładu.
# PRZYKŁAD 5.13
x <- rpois(100, lam = 5)
y <- rpois(100, lam = 5)
ks.test(x, y)## Warning in ks.test(x, y): p-value will be approximate in the presence of
## ties##
## Two-sample Kolmogorov-Smirnov test
##
## data: x and y
## D = 0.1, p-value = 0.6994
## alternative hypothesis: two-sidedy <- rpois(100, lam = 2)
ks.test(x, y)## Warning in ks.test(x, y): p-value will be approximate in the presence of
## ties##
## Two-sample Kolmogorov-Smirnov test
##
## data: x and y
## D = 0.7, p-value < 2.2e-16
## alternative hypothesis: two-sidedks.test(x, "ppois", lam = 5)## Warning in ks.test(x, "ppois", lam = 5): ties should not be present for the
## Kolmogorov-Smirnov test##
## One-sample Kolmogorov-Smirnov test
##
## data: x
## D = 0.25596, p-value = 4.078e-06
## alternative hypothesis: two-sidedx <- rnorm(100)
y <- rnorm(100)
ks.test(x, y)##
## Two-sample Kolmogorov-Smirnov test
##
## data: x and y
## D = 0.09, p-value = 0.8127
## alternative hypothesis: two-sidedy <- rnorm(100, 0, 2)
ks.test(x, y)##
## Two-sample Kolmogorov-Smirnov test
##
## data: x and y
## D = 0.2, p-value = 0.03663
## alternative hypothesis: two-sidedks.test(x, "pnorm", 0, 1)##
## One-sample Kolmogorov-Smirnov test
##
## data: x
## D = 0.079068, p-value = 0.5594
## alternative hypothesis: two-sided