Zastosowanie pakietu R w statystyce medycznej

LABORATORIUM 7

Przed rozpoczęciem zajęć należy zainstalować pakiet car.

install.packages("car")


Dopasowywanie parametrów rozkładów

Za pomocą funkcji fitdistr() (biblioteka MASS) możliwa jest estymacja nieznanych parametrów rozkładów, przy czym postać rozkładu musi być uprzednio znana - może to być albo jeden z wbudowanych rozkładów albo też nasza własna funkcja.

# Przykład 7.1
library(MASS)

df <- data.frame(x = rnorm(200), y = rexp(200, 0.1))
df.fit1 <- fitdistr(df$x, "normal")
df.fit1
##       mean          sd    
##   0.02948909   1.00557069 
##  (0.07110459) (0.05027853)
df.fit2 <- fitdistr(df$y, "exponential")
df.fit2
##       rate    
##   0.096335765 
##  (0.006811967)

Obiekt, który powstaje po dopasowaniu ma następujące istotne pola: estimate, sd oraz vcov. Pierwsze dwa podają odpowiednio oszacowanie parametrów rozkładu oraz odchylenie standardowe, ostatnia to macierz wariancji-kowariancji.

# Przykład 7.2

df.fit1$estimate
##       mean         sd 
## 0.02948909 1.00557069
df.fit1$sd
##       mean         sd 
## 0.07110459 0.05027853
df.fit1$vcov
##             mean          sd
## mean 0.005055862 0.000000000
## sd   0.000000000 0.002527931

Rzecz jasna, nic nie stoi na przeszkodzie, aby sprawdzić hipotezę, czy dane faktycznie mogą pochodzić z rozkładów o wyestymowanych parametrach.

# Przykład 7.3

# Testowanie rozkładu normalnego
mu <- df.fit1$estimate[1]
sigma <- df.fit1$estimate[2]
ks.test(df$x, "pnorm", mu, sigma)
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  df$x
## D = 0.038598, p-value = 0.9269
## alternative hypothesis: two-sided
# Testowanie rozkladu wykladniczego
rate <- df.fit2$estimate
ks.test(df$y, "pexp", rate)
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  df$y
## D = 0.037039, p-value = 0.9467
## alternative hypothesis: two-sided
# Rysowanie
library(ggplot2)
## Registered S3 methods overwritten by 'ggplot2':
##   method         from 
##   [.quosures     rlang
##   c.quosures     rlang
##   print.quosures rlang
ggplot(df) +  geom_histogram(aes(x = x,..density..), fill="blue", colour="black", alpha=0.4) + 
stat_function(fun = dnorm, args = list(mean = mu, sd = sigma), colour="red", size=1.5)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(df) + geom_histogram(aes(x = y,..density..), fill="blue", colour="black", alpha=0.4) + 
stat_function(fun = dexp, args = list(rate = rate), colour="red", size=1.5)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.


Jednoczynnikowa analiza wariancji

Rozważamy pojedynczą zmienną numeryczną lub porządkową i chcielibyśmy się dowiedzieć, czy przeciętna wartość tej zmiennej ulega zmianom w różnych grupach. Jednoczynnikowa analiza wariancji (ANOVA) pozwalająca rozdzielić całkowitą zmienność danych na tę, która może być przypisana różnicom między osobnikami z różnych grup (zmienność międzygrupowa), oraz losowe zmiany między osobnikami wewnątrz każdej grupy (zmienność wewnątrzgrupowa).

\(s^2 = \sum_{i=1}^k \sum_{j=1}^{n_i} (x_{ij}-\overline{x})^2 = \sum_{i=1}^k \sum_{j=1}^{n_i} (x_{ij}-\overline{x_i})^2 + \sum_{i=1}^k n_i (\overline{x_i}-\overline{x})^2 = s_{wew}^2 + s_{pom}^2\)

# Przykład 7.4
plants <- PlantGrowth #30 obserwacji dotyczących zbiorów (w kg) podzielone na trzy grupy (kontrolna i dwie grupy specjalnego traktowania)

# Wizualizacja przy pomocy wykresów skrzynkowych (boxplot)
boxplot(weight~group, plants, ylab="weight", col=c("red","green","blue"))

# Obliczmy średnie i wariancje w każdej grupie
grupy <- levels(plants$group)
srednie <- sapply(grupy, function(g) mean(plants$weight[plants$group==g]));srednie
##  ctrl  trt1  trt2 
## 5.032 4.661 5.526
plants$mu <- rep(srednie, each=10)
plants$s <- (plants$weight-plants$mu)^2
s.wew <- sum(plants$s); s.wew
## [1] 10.49209
# Obliczmy całkowitą średnią i wariancję międzygrupową
k <- 3
n <- rep(10,k)
N <- sum(n)
mu.tot <- mean(plants$weight)
s.pom <- sum(n*(srednie-mu.tot)^2); s.pom
## [1] 3.76634
# Obliczmy stosunek wariancji międzygrupowej do wariancji wewnątrzgrupowej
s.wew <- s.wew/(N-k)
s.pom <- s.pom/(k-1)
s.pom/s.wew
## [1] 4.846088
# Wykonujemy jednoczynnikową analizę wariancji przy pomocy funkcji aov()
plants.aov <- aov(weight~group, plants)
summary(plants.aov)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## group        2  3.766  1.8832   4.846 0.0159 *
## Residuals   27 10.492  0.3886                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

P-wartość wynosi 0.0159, co oznacza, że na poziomie istotności \(\alpha = 0.05\) można odrzucić hipotezę zerową o tym, że średnie w grupach są takie same. Aby się dowiedzieć, w której grupie średnia jest inna można wykonać test Tukey HSD (Tukey Honest Significant Differences), jako argument podając wynik testu ANOVA.

TukeyHSD(plants.aov)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = weight ~ group, data = plants)
## 
## $group
##             diff        lwr       upr     p adj
## trt1-ctrl -0.371 -1.0622161 0.3202161 0.3908711
## trt2-ctrl  0.494 -0.1972161 1.1852161 0.1979960
## trt2-trt1  0.865  0.1737839 1.5562161 0.0120064

Założenia stojące za analizą wariancji

  1. W badanej populacji zmienna w każdej grupie ma rozkład normalny.

  2. Wariancje we wszystkich grupach są takie same.

Równość wariancji wielu grup można sprawdzić testem Levene'a lub Bartletta. Test Levene'a jest mniej wrażliwy na odchylenia od rozkładu normalnego. Wykorzystamy funkcję leveneTest() z pakietu car (należy go uprzednio zainstalować). Hipoteza zerowa testu Levene'a mówi, że wariancje poszczególnych grup nie różnią się istotnie od siebie.

# Przykład 7.5
library(car)
leveneTest(weight ~ group, plants)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  1.1192 0.3412
##       27


Test Welcha

W sytuacji, gdy nie jest spełniony warunek równości wariancji między różnymi grupami, możemy posłużyć się testem Welcha, który w pakiecie R jest realizowany przez funkcję oneway.test(). Alternatywnie, można wykonać test T-Studenta parami przy użyciu funkcji pairwise.t.test()

# Przykład 7.6
oneway.test(weight ~ group, plants)
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  weight and group
## F = 5.181, num df = 2.000, denom df = 17.128, p-value = 0.01739
pairwise.t.test(plants$weight, plants$group)
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  plants$weight and plants$group 
## 
##      ctrl  trt1 
## trt1 0.194 -    
## trt2 0.175 0.013
## 
## P value adjustment method: holm


Test Kruskala-Wallisa

Nieparametryczną alternatywą jednoczynnikowej analizy wariancji jest test Kruskalla-Wallisa, który jest rozszerzeniem testu sumy rang Wilcoxona. Test ten nie zakłada równości wariancji ani normalności rozkładów. Hipotezą zerową \(H_0\) jest równość dystrybuant rozkładów w porównywanych populacjach.

# Przykład 7.7
kruskal.test(weight ~ group, plants)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  weight by group
## Kruskal-Wallis chi-squared = 7.9882, df = 2, p-value = 0.01842


Zadanie punktowane

Wczytaj zbiór danych iris z biblioteki datasets i zapoznaj się z nim. Następnie sprawdź (wykonując testy normalności i równości wariancji międzygrupowej), która cecha kwiatów spełnia najlepiej założenia jednoczynnikowej analizy wariancji. Dla tej cechy wykonaj ANOVA, skomentuj wynik analizy i wykonaj wykres pudełkowy.


##              Df Sum Sq Mean Sq F value Pr(>F)    
## Species       2  11.35   5.672   49.16 <2e-16 ***
## Residuals   147  16.96   0.115                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1