Zastosowanie pakietu R w statystyce medycznej

LABORATORIUM 6


Test \(\chi^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 6.1
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.004701
p2 <- 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.5108
p <- x/sum(x)
chisq.test(x, p = p)
## 
##  Chi-squared test for given probabilities
## 
## data:  x
## X-squared = 0, df = 4, p-value = 1
plot(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 6.2
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 = 6.5595, df = 14, p-value = 0.9504
plot(h$mids, h$density)
lines(h$mids, dpois(h$mids, lam = 3), col="red", lwd=2)


Test Kołmogorowa-Smirnowa

W tym przypadku (funkcja ks.test()) możemy:

  1. stosować test zarówno do rozkładów dyskretnych jak i ciągłych,
  2. porównywać dwa rozkłady eksperymentalne albo rozkład eksperymentalny z rozkładem teoretycznym (ale tylko ciągłym).

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 6.3
x <- 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-sided
y <- rnorm(100, 0, 2)
ks.test(x, y)
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.21, p-value = 0.02431
## alternative hypothesis: two-sided
ks.test(x, "pnorm", 0, 1)
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  x
## D = 0.083154, p-value = 0.4938
## alternative hypothesis: two-sided
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.08, p-value = 0.9062
## alternative hypothesis: two-sided
y <- 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.61, p-value < 2.2e-16
## alternative hypothesis: two-sided
# Wynik poniższego testu będzie niepoprawny, ponieważ rozkład Poissona jest rozkładem dyskretnym
ks.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.20049, p-value = 0.0006449
## alternative hypothesis: two-sided


Testy 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 6.4
shapiro.test(rnorm(100))
## 
##  Shapiro-Wilk normality test
## 
## data:  rnorm(100)
## W = 0.99156, p-value = 0.7886
shapiro.test(rnorm(100, 0, 5))
## 
##  Shapiro-Wilk normality test
## 
## data:  rnorm(100, 0, 5)
## W = 0.98079, p-value = 0.153
shapiro.test(runif(100, 0, 1))
## 
##  Shapiro-Wilk normality test
## 
## data:  runif(100, 0, 1)
## W = 0.94573, p-value = 0.000439
wilcox.test(rnorm(100))
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  rnorm(100)
## V = 2700, p-value = 0.5485
## alternative hypothesis: true location is not equal to 0
wilcox.test(rnorm(100, 1, 2), mu = 1)
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  rnorm(100, 1, 2)
## V = 2548, p-value = 0.9383
## alternative hypothesis: true location is not equal to 1
wilcox.test(rnorm(100, 1, 2), mu = 2)
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  rnorm(100, 1, 2)
## V = 1116, p-value = 1.28e-06
## alternative hypothesis: true location is not equal to 2


Współczynnik korelacji Pearsona

Do kwantifykacji 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}}\).

Funkcją obliczającą współczynnik korelacji w R jest cor().

# Przykład 6.5
x <- c(1,2,3,4,5)
y <- c(2,4,1,2,2)
cor(x, y)
## [1] -0.2886751

W przypadku, gdy mamy do czynienia z korelacją rang, używa się współczynnika Spearmana cor(method="spearman"), który jest de facto współczynnikiem Pearsona dla danych przekształconych w rangi.

# Przykład 6.6
cor(x, y, method="spearman")
## [1] -0.2236068
x; rank(x)
## [1] 1 2 3 4 5
## [1] 1 2 3 4 5
y; rank(y)
## [1] 2 4 1 2 2
## [1] 3 5 1 3 3
cor(rank(x), rank(y))
## [1] -0.2236068

Jeśli do funkcji cor() przekażemy ramkę danych, zostaną wyliczone współczynniki korelacji pomiędzy każdą parą zmiennych i przedstawione w posteci macierzy

# Przykład 6.7
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.0000000

Warto 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 6.8
library(corrplot)
## corrplot 0.84 loaded
corrplot(cor(df))

W celu otrzymania pełnych informacji, dotyczących parametrów współczynnika korelacji należy użyć funkcji cor.test()

# Przykład 6.9
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.2886751

Należ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 \(\chi^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 \(\chi^2\) Pearsona, wywoływany funkcją chisq.test()

# Przykład 6.10
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     104
chisq.test(T)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  T
## X-squared = 5.5711, df = 1, p-value = 0.01826

Niska 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 6.11
Te <- chisq.test(T)$expected
Te
##            y
## x           nie pracuje   pracuje
##   fizyk        14.01961  40.98039
##   humanista    37.98039 111.01961
Te / sum(Te)
##            y
## x           nie pracuje   pracuje
##   fizyk      0.06872357 0.2008843
##   humanista  0.18617839 0.5442138

Powyż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 6.12
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.5442138

Natomiast 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 6.13
T / sum(T)
##            y
## x           nie pracuje    pracuje
##   fizyk      0.03431373 0.23529412
##   humanista  0.22058824 0.50980392


Zadanie punktowane

Wczytaj dane z pliku http://www.if.pw.edu.pl/~paluch/MSR/data/zad6.txt zawierającego 200 liczb losowych z pewnego rozkładu. Zweryfikuj następujące hipotezy:

  1. dane pochodzą z rozkładu normalnego (użyj testu normalności),

  2. dane pochodzą z rozkładu Poissona o parametrze \(\lambda=4\) (użyj testu \(\chi^2\)),

  3. dane pochodzą z rozkładu wykładniczego o parametrze \(\lambda=0.25\) (użyj test Kołmogorowa-Smirnowa).

Następnie wykonaj histogram w oparciu o te dane, umieść empiryczną gęstość prawdopodobieństwa na wykresie w postaci punktów i teoretyczną gęstość prawdopodobieństwa w postaci linii ciągłej.