1 Test parametryczny - Analiza wariancji

  • Celem analizy wariancji (ANOVA) jest zazwyczaj testowanie istotności różnic pomiędzy średnimi w wielu grupach.

  • służą do oceny czy średnie wartości cechy Y (cecha ilościowa) różnią się istotnie pomiędzy grupami wyznaczonymi przez zmienną X.

  • W analizie jednoczynnikowej badamy zależność pomiędzy cechą Y a jedną zmienną jakościową.

  • Hipoteza zerowa - średnie w grupach są jednakowe: H0: m1=m2=…mk

  • Hipoteza alternatywna: H1: co najmniej dwie średnie różnią się między sobą.

  • analiza wariancji nie pozwala stwierdzić pomiędzy którymi grupami występują różnice.

  • aby określić między którymi grupami występują różnice trzeba wykonać porównania wielokrotne (“post-hoc”)

1.1 Podstawowe założenia analizy wariancji:

  1. Zmienna w każdej z analizowanych grup ma rozkład zbliżony do rozkładu normalnego.
  2. Wariancje w grupach są do siebie podobne.
  3. Próby zostały losowo wybrane z populacji.

Przykład: Czy średnie ceny mieszkań różnią się w poszczególnych dzielnicach?

W przykładzie chcemy sprawdzić czy istnieją (istotne statystycznie) różnice pomiędzy średnimi cenami mieszkań w 3 dzielnicach Wrocławia. Przykład wykorzystuje dane z podręcznika „Przewodnik po pakiecie R” (Przemysław Biecek).

1.1.1 Dane

Dane dotyczą cen 200 różnych mieszkań. Dane zawierają następujące kolumny:

  • cena – cena mieszkania
  • pokoi – liczba pokoi w mieszkaniu,
  • powierzchnia – powierzchnia mieszkania
  • dzielnica – zmienna tekstowa przyjmująca wartości: Biskupin, Śródmieście, Krzyki.
  • typ.budynku – wieżowiec, niski blok (4-piętrowy blok), kamienica.
library(Przewodnik)
data(mieszkania)
daneMieszkania <- mieszkania
names(daneMieszkania)
## [1] "cena"         "pokoi"        "powierzchnia" "dzielnica"    "typ.budynku"
head(daneMieszkania)
##       cena pokoi powierzchnia   dzielnica typ.budynku
## 1 226298.2     4         76.4      Krzyki   wiezowiec
## 2 128121.4     1         19.4    Biskupin   wiezowiec
## 3 118333.4     1         21.5 Srodmiescie   wiezowiec
## 4 202064.2     4         65.2 Srodmiescie   wiezowiec
## 5 130303.0     1         21.4      Krzyki   wiezowiec
## 6 156241.9     2         35.8    Biskupin   wiezowiec

1.1.2 Statystyki w grupach

Statystyki opisowe dla całego zbioru danych

summary(daneMieszkania)
##       cena            pokoi       powierzchnia         dzielnica 
##  Min.   : 83280   Min.   :1.00   Min.   :17.00   Biskupin   :65  
##  1st Qu.:143304   1st Qu.:2.00   1st Qu.:31.15   Krzyki     :79  
##  Median :174935   Median :3.00   Median :43.70   Srodmiescie:56  
##  Mean   :175934   Mean   :2.55   Mean   :46.20                   
##  3rd Qu.:208741   3rd Qu.:3.00   3rd Qu.:61.40                   
##  Max.   :295762   Max.   :4.00   Max.   :87.70                   
##      typ.budynku
##  kamienica :61  
##  niski blok:63  
##  wiezowiec :76  
##                 
##                 
## 
library(dplyr)
by_dzielnica <- group_by(daneMieszkania, dzielnica)

stat_cena_dzielnice <- summarise(by_dzielnica, 
                            ilosc = n(),
                            srednia = mean(cena), 
                            wariancja = var(cena), 
                            min = min(cena), 
                            max = max(cena))
stat_cena_dzielnice
## # A tibble: 3 x 6
##   dzielnica   ilosc srednia   wariancja     min     max
##   <fct>       <int>   <dbl>       <dbl>   <dbl>   <dbl>
## 1 Biskupin       65 189494. 1652233491. 120290. 295762.
## 2 Krzyki         79 168173. 2037842235.  83280. 252418.
## 3 Srodmiescie    56 171143. 1574695352. 103749. 245846.

1.1.3 Wizualizacja danych w grupach

ggplot(daneMieszkania, aes(x = dzielnica, y = cena)) + geom_boxplot() + theme_bw()

1.1.4 Sprawdzenie założeń analizy wariancji

Testowanie podobieństwa do rozkładu normalnego (Test Shapiro-Wilka)

Jednym z założeń analizy wariancji jest zgodność danych z rozkładem normalnym. Do testowania podobieństwa do rozkładu normalnego można wykorzystać testu Shapiro-Wilka. Alternatywą test Shapiro-Wilka jest test Kołmogorowa-Smirnowa.

Test Shapiro-Wilka (w R: shapiro.test()):

  • Służy do testowania podobieństwa rozkładu danej zmiennej do rozkładu normalnego (sprawdzamy czy interesujące nas zmienne mają rozkłady zbliżone do rozkładu normalnego)

  • Hipoteza zerowa: rozkład naszej zmiennej jest zbliżony do rozkładu normalnego

  • UWAGA: Istotny wynik testu Shapiro-Wilka świadczy o tym, że rozkład zmiennej obserwowanej NIE JEST podobny do rozkładu normalnego.

Sprawdzenie zgodności rozkładu normalnego dla całego zbioru danych

Wykresy ilustrujące rozkład cen mieszkań we Wrocławiu (odpowiednio histogram, wykres gęstości, wykres kwantylowy)

p1 <- ggplot(daneMieszkania, aes(x=cena)) + geom_histogram() + theme_bw()
p2 <- ggplot(daneMieszkania, aes(x=cena)) + geom_density() + theme_bw()
p3 <- ggplot(daneMieszkania, aes(sample=cena)) + geom_qq() + theme_bw()
gridExtra::grid.arrange(p1, p2, p3, ncol=3)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Wyniki testu Shapiro-Wilka dla cen mieszkań we Wrocławiu 

shapiro.test(daneMieszkania$cena)
## 
##  Shapiro-Wilk normality test
## 
## data:  daneMieszkania$cena
## W = 0.98791, p-value = 0.08729

Interpretacja: Zerową hipotezą dla tego testu jest, że dane pochodzą z rozkładu normalnego. Zatem jeśli p-value jest poniżej istotnego progu (zazwyczaj 0,05), wtedy hipoteza zerowa jest odrzucana i przyjmowana jest hipoteza alternatywa. Hipotezą alternatywną jest, że dane nie pochodzą z rozkładu normalnego.

W przykładzie p-value 0,08 > p=0,05 (założony próg) zatem nie ma podstaw do odrzucenia hipotezy zerowej (dane pochodzą z rozkładu normalnego).

Sprawdzenie zgodności rozkładu normalnego dla poszczególnych dzielnic

Przy analizie wariancji założenie o zgodności danych z rozkładem normalnym powinno być spełnione w każdej analizowanej grupie (w tym przykładzie w każdej dzielnicy).

Wykresy

library(ggplot2)
p <- ggplot(daneMieszkania, aes(x=cena)) + geom_density() + theme_bw()
p + facet_grid(.~dzielnica)

library(ggplot2)
p <- ggplot(daneMieszkania, aes(sample=cena)) + geom_qq() + theme_bw()
p + facet_grid(.~dzielnica)

Test Shapiro-Wilka dla każdej dzielnicy

dzielnica <- unique(daneMieszkania$dzielnica)
test.out <- c()#wartości p value dla testu 
for(d in dzielnica) {
    df <- daneMieszkania[daneMieszkania$dzielnica==d,]
    st <- shapiro.test(df$cena) #przypisanie wyników testu do zmiennej. Wyniki testu są przechowywane w liście
    test.out <- c(test.out, st$p.value) #st$p.value pozwala na odczytanie tylko poziomu istotności 
}
names(test.out) <- dzielnica
test.out
##      Krzyki    Biskupin Srodmiescie 
##  0.08893684  0.09416835  0.09193965

Interpretacja wyników: Dla każdej z dzielnic Wrocławia otrzymano p-wartość większą od 0.05. Nie ma zatem podstaw do odrzucenia hipotezy zerowej o zgodności z rozkładem normalnym.

Założenie o zgodności z rozkładem normalnym jest spełnione dla każdej dzielnicy oraz dla całego zbioru danych.

Testowanie założenia o jednorodności wariancji.

  • Testy te sprawdzają czy dwie różne próbki mają równe wariancje (homogeniczność wariancji, jednorodność wariancji).

  • Hipotezą zerową jest, że wariancje są równe w poszczególnych grupach.

  • Jeśli poziom istotności p-value spada poniżej istotnego progu (zazwyczaj 0,05) wtedy możemy odrzucić hipotezę zerową i przyjąć hipotezę altenatywną (wariancje nie są równe).

W przypadku testowania założenia o homogeniczności (jednorodności) wariancji dla wielu grup można użyć jednego z dwóch testów:

  • Test Barletta (Barlett’s test) - służy do testowania czy próby pochodzą z populacji o takiej samej wariancji; jest czuły na odchylenia od normalności.

  • Test Levene’a - jest alternatywą dla testu Barlett’a; jest mniej czuły na odchylenia od normalności. Jeśli mamy silne dowody, że nasze dane mają rozkład normalny, lub zbliżony do normalnego, to test Barletta da lepsze wyniki.

bartlett.test(daneMieszkania$cena, g=daneMieszkania$dzielnica)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  daneMieszkania$cena and daneMieszkania$dzielnica
## Bartlett's K-squared = 1.3075, df = 2, p-value = 0.5201
library(car)
leveneTest(daneMieszkania$cena, g=daneMieszkania$dzielnica)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   2  0.6934 0.5011
##       197

W obu przypadkach wartość p-value jest większa od 0.05 (0.5201 dla testu Barletta oraz 0.5011 dla testu Levene). Zatem nie ma podstaw do odrzucenia hipotezy zerowej o równości wariancji w analizowanych grupach. Założenie o równości wariancji jest spełnione.

Równoliczność grup.

Przy stosowaniu analizy wariancji grupy powinny mieć także taką samą (lub bardzo zbliżoną liczność). W przypadku danych dotyczących cen mieszkań we Wrocławiu założenie to nie jest spełnione. Mimo tego, analiza wariancji może być zastosowana.

table(daneMieszkania$dzielnica)
## 
##    Biskupin      Krzyki Srodmiescie 
##          65          79          56

1.1.5 Jednoczynnikowa analiza wariancji – ANOVA

Do oceny różnic między średnimi cenami wykorzystamy analizę wariancji.

CenaDzielnica <- aov(cena~dzielnica, data = daneMieszkania)
summary(CenaDzielnica)
##              Df    Sum Sq   Mean Sq F value  Pr(>F)   
## dzielnica     2 1.800e+10 8.998e+09   5.046 0.00729 **
## Residuals   197 3.513e+11 1.783e+09                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dla każdej zmiennej i dla reszt wyznacza się następujące wartości:

  • Df – liczba stopni swobody. Dla zmiennej jakościowej jest ona równa k-1, gdzie k to liczba poziomów danej zmiennej.
  • Sum Sq – suma kwadratów wartości wyjaśnionych przez daną zmienną. Wartość wykorzystana do wykonania testu F.
  • Mean Sq – średnia suma kwadratów.
  • F value – wartość statystyki testowej dla testu F., weryfikującego hipotezę o równości średniej wartości cechy w podgrupach danej zmiennej jakościowej.
  • Pr(>F) – wartość p wyznaczona dla testu F

Istotność dla poszczególnych zmiennych jest oznaczana na dwa sposoby. Podawana jest p-wartość dla testu istotności F, ale również zaznaczany jest kod istotności (po prawej stronie p-wartości). Stosowane jest następujące kodowanie:

*** (wysoka istotność) – p wartość jest mniejsza od 0,001, ** p-wartość jest większa od 0,001 ale mniejsza od 0,01, * p-wartość jest większa od 0,01 ale mniejsza od 0,05 (ta wartość jest najczęściej wybierana za graniczą wartość istotności statystycznej), . (kropka) p-wartość jest mniejsza od 0,1.

Brak kodu oznacza brak podstaw do odrzucenia hipotezy zerowej na poziomie istotności 0,1.

W omawianym przykładzie p-wartość równa się 0,007294 a zatem na poziomie istotności 0,01 skłonni jesteśmy uznać, że średnie ceny mieszkań różnią się dla przynajmniej jednej pary. Aby sprawdzić, które pary różnią się między sobą wykonuje się testy post-hoc (tzw. Testy wielokrotnych porównań).

2 Porównanie parami średnich wielu prób – test i wykres Tukey’a

  • Wynikiem testu Tukeya są różnice między średnimi oraz 95% przedział ufności dla tych różnic dla każdej pary porównań.
  • Wynik jest pokazywany w postaci numerycznej i graficznej.
  • Szukamy takich par, dla których przedział ufności średnich nie przecina 0, wskazując na statystycznie istotne różnice między grupami.

Graficzna reprezentacja wyniku – przedziały ufności dla różnic w średnich cenach mieszkań w poszczególnych dzielnicach.

posthoc <- TukeyHSD(x=CenaDzielnica, 'dzielnica', conf.level=0.95)
posthoc
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = cena ~ dzielnica, data = daneMieszkania)
## 
## $dzielnica
##                            diff       lwr        upr     p adj
## Krzyki-Biskupin      -21321.019 -38021.10 -4620.9333 0.0081457
## Srodmiescie-Biskupin -18350.541 -36532.88  -168.2053 0.0473579
## Srodmiescie-Krzyki     2970.478 -14450.28 20391.2340 0.9145465
par(mar=c(5.1,10.1,4.1,4.1)) # zmiana marginesów
plot(posthoc, las=1)

Wnioski z analizy danych:

  • Wnioskiem z przeprowadzenia testu post-hoc jest stwierdzenie, że ceny mieszkań w dzielnicy Krzyki i Śródmieście są porównywalne.

  • Ceny mieszkań w dzielnicy Krzyki oraz Biskupin różnią się istotnie.

Jak zmienia się cena mieszkań dla poszczególnych typów budynku (oblicz odpowiednie statystyki i wykonaj odpowiednie wizualizacje). Sprawdź założenia analizy wariancji (rozkład normalny, równość wariancji, równość grup). Wykonaj odpowiednie testy i wizaulizacje. Czy cena mieszkań w różnych typach budynków zmienia się istotnie? (Wykonaj analize wariancji oraz testy post-hoc). Wskaż dla jakiej pary różnice są istotne statystycznie.

3 Test nieparametryczny dla wielu prób

  • test Kruskala-Wallisa
  • test Friedmana

3.1 Test Kruskala-Wallisa

Nieparametryczną alternatywą dla jednoczynnikowej analizy wariancji stanowi test Kruskala-Wallisa. Za pomocą tego testu sprawdzamy, czy n niezależnych próbek pochodzi z tej samej populacji, bądź z populacji z taką samą medianą. Poszczególne próbki nie muszą mieć takiej samej liczebności. Test ten jest rozszerzeniem testu Two sample Wilcoxon dla 3 i więcej grup.

Przykład: Czy średnie wartości mpg różnią się pomiędzy samochodami z Europy, Japonii i USA?

3.1.1 Dane

Zbiór danych w pliku cardata1.csv zawiera informacje nt. 154 samochodów wyprodukowanych w Europie, Japoni oraz USA.

Zmienne:

  • mpg – ilość mil jaką można przejechać na jednym galonie paliwa.
  • cylinders – ilość cylindrów w silniku
  • displace –
  • horsepower – moc silnika w koniach mechanicznych.
  • Accel - przyspieszenie
  • Year – rok produkcji
  • Weight – waga samochodu
  • Origin – pochodzenie
  • Price - cena
  • Make - marka
cardata1 <- read.csv('dane/cardata1.csv')
summary(cardata1)
##       mpg          cylinders        displace       horsepower    
##  Min.   :15.50   Min.   :3.000   Min.   : 70.0   Min.   : 48.00  
##  1st Qu.:22.55   1st Qu.:4.000   1st Qu.: 98.0   1st Qu.: 70.00  
##  Median :28.90   Median :4.000   Median :134.0   Median : 85.00  
##  Mean   :28.79   Mean   :4.844   Mean   :153.7   Mean   : 88.86  
##  3rd Qu.:34.27   3rd Qu.:6.000   3rd Qu.:179.0   3rd Qu.:102.25  
##  Max.   :46.60   Max.   :8.000   Max.   :360.0   Max.   :165.00  
##                                                  NA's   :4       
##      accel            year           weight        origin         
##  Min.   :11.20   Min.   :78.00   Min.   :1755   Length:154        
##  1st Qu.:14.70   1st Qu.:79.00   1st Qu.:2144   Class :character  
##  Median :15.85   Median :80.00   Median :2618   Mode  :character  
##  Mean   :16.28   Mean   :79.94   Mean   :2672                     
##  3rd Qu.:17.75   3rd Qu.:81.00   3rd Qu.:3068                     
##  Max.   :24.80   Max.   :82.00   Max.   :4360                     
##                                                                   
##      make              model               price      
##  Length:154         Length:154         Min.   : 1900  
##  Class :character   Class :character   1st Qu.: 3281  
##  Mode  :character   Mode  :character   Median : 4225  
##                                        Mean   : 4610  
##                                        3rd Qu.: 5469  
##                                        Max.   :15475  
##

3.1.2 Wizualizacja danych

Dane zawierają parametry samochodów z Europy, USA, Japonii.

ggplot(cardata1, aes(x = origin, y = mpg)) + geom_boxplot() + theme_bw()

p <- ggplot(cardata1, aes(x=mpg)) + geom_histogram(binwidth = 5) + theme_bw()
p + facet_grid(.~origin)

3.1.3 Założenia: rozkład normalny i równość wariancji

Testowanie podobieństwa do rozkładu normalnego dla całego zbioru danych:

shapiro.test(cardata1$mpg)
## 
##  Shapiro-Wilk normality test
## 
## data:  cardata1$mpg
## W = 0.9731, p-value = 0.004127

Testowanie podobieństwa do rozkładu normalnego dla poszczególnych grup:

grupa <- unique(cardata1$origin)
test.out <- c()#wartości p value dla testu 
for(i in grupa) {
    df <- cardata1[cardata1$origin==i,]
    st <- shapiro.test(df$mpg) #przypisanie wyników testu do zmiennej. Wyniki testu są przechowywane w liście
    test.out <- c(test.out, st$p.value) #st$p.value pozwala na odczytanie tylko poziomu istotności 
}
names(test.out) <- grupa
test.out
##      Europa         USA     Japonia 
## 0.296997525 0.001946911 0.373110870

Założenie o zgodności z rozkładem normalnym nie jest spełnione dla USA oraz dla całego zbioru danych.

Testowanie założenia o jednorodności wariancji

library(car)
leveneTest(mpg~origin, cardata1)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value  Pr(>F)  
## group   2  3.6899 0.02725 *
##       151                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Hipotezę zerową o równości wariancji można odrzucić przy założeniu progu (p-wartość = 0.027).

Liczność grup

table(cardata1$origin)
## 
##  Europa Japonia     USA 
##      25      44      85

Podsumowanie: Założenia o zgodności z rozkładem normalnym, jednorodności wariancji oraz równoliczności grup nie są spełnione. Ponadto jedna z grup (Europa) ma tylko 25 obserwacji. Zostanie zatem zastosowana nieparametryczna alternatywa ANOVA – test Kruskal-Wallis

3.1.4 Test Kruskala-Wallisa

kruskal.test(mpg ~ origin, data=cardata1)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  mpg by origin
## Kruskal-Wallis chi-squared = 44.126, df = 2, p-value = 2.619e-10

Interpretacja: p-value poniżej założonego progu – mediany nie są równe.

3.1.5 Porównania wielokrotne

library(FSA)
dunnTest(cardata1$mpg, cardata1$origin, method = "bonferroni")
##         Comparison          Z      P.unadj        P.adj
## 1 Europa - Japonia -0.7664177 4.434278e-01 1.000000e+00
## 2     Europa - USA  4.1614688 3.162073e-05 9.486218e-05
## 3    Japonia - USA  6.1316118 8.699314e-10 2.609794e-09

Jeśli wartość P.adj jest mniejsza od założonego progu (zazwyczaj 0,05) to hipotezę zerową o braku różnic między grupami można odrzucić na rzecz hipotezy alternatywnej (grupy różnią się).

Pytanie: Pomiędzy, którymi grupami istnieją statystycznie istotne różnice?

Czy oczekiwana dalsza długość życia w 2007 roku różni się pomiędzy kontynentami? Oblicz statystyki w grupach (oczekiwaną długość trwania życia dla poszczególnych kontynentów) oraz wykonaj odpowiednią wizualizację. Używając testu Kruskala - Wallisa sprawdź, czy różnice są istotne statystycznie? Pomiędzy, którymi kontynentami odnotowano różnice.