NMFM301, ZS 2024/25

Cvičenie 11 (praktická časť 5)

(metóda analýzy rozptylu – ANOVA)

Zdrojový soubor pro Knit: Rmd soubor

Na úvod si vyskúšajte v programe RStudio otvoriť podkladový Rmd soubor súbor a skompilovať ho pomocou tlačítka Knit (pre správne fungovanie, resp. zobrazenie kódovania diakritiky, je nutné súbor otvoriť v kódovaní UTF8 – pomocou ponuky “Reopen with encoding” v menu programu). Tento krok síce nie je nutný, ale jeho bezproblemové zvládnutie sa môže ukázať byť prínosné pri spracovaní samostatnej zápočtovej práce budúci týždeň.

Užitočné materiály pre prácu so štatistickým softwarom R

  • Bína, V., Komárek, A. a Komárková, L.: Jak na jazyk R. (PDF súbor)
  • Komárek, A.: Základy práce s R. (PDF súbor)
  • Kulich, M.: Velmi stručný úvod do R. (PDF súbor)
  • De Vries, A. a Meys, J.: R for Dummies. ISBN-13: 978-1119055808 (link na PDF súbor)


Pre účely cvičenia využijeme dva konkrétne datové súbory a tiež niektoré dodatočné funkcie (Rkové príkazy), ktoré do programu R načítame pomocou nasledujúceho kódu (funkcie v rámci skriptu functions.RData budú k dispozícii aj pre samostatnú prácu).

load(url("http://msekce.karlin.mff.cuni.cz/~maciak/NMFM301/cv5_data.RData"))
load(url("http://msekce.karlin.mff.cuni.cz/~maciak/NMFM301/functions.RData"))

K dispozíci máme dva datasety, zam a sex´. Data zam pocházejí ze sociálního průzkumu provedeného v ČR v roce 2004 a datový súbor sex predstavuje maly výber z výsledkov dotazníkového sociálneho šetrenia v Spojených štátoch, ktorého cieľom bolo zistiť mieru sexuálnej aktivity medzi mladými ľudmi.

Zaujímať nás bude niekoľko konkrétnych otázok formulovaných v zmysle štatistického testu. Doležitý je správny výber pravdepodobnostného modelu a exxplicitná formulácia otázky pomocou nulovej a alternatívnej hypotézy. Samotnej voľbe pravdepodobnostného modelu (t.j., konkrétneho štatistického testu) samozrejme predchádza exploratívna analýza uskutočnená pomocou vhodných popisných charakteristík a obrázkov.

Data si prohlédněte poklepáním na příslušný řádek v okénku “Workspace” nebo pomoci příkazů summary(), head(), str(), a pod. Pomocou niekoľkých obrázkov a grafov preskúmajte štruktúru v datach – vždy s ohľadom na kladenú teoretickú otázku.

I. Opakování: Dvouvýběrové testy

Nyní budeme pracovat s datovým souborem sex. Tento soubor obsahuje informace o sexuálním životě 4467 Američanů starších 20 let. Obsahuje veličiny active (byl/a někdy sexuálně aktivní?), veksex (věk sexuální iniciace - pouze u aktivních), pocpart (počet partnerů za celý život), pocpart.rok (počet partnerů za uplynulý rok) a dále demografická data (věk, pohlaví, vzdělání, rodinný status).

Porovnáme nejprve věk při prvním pohlavním styku mužů a žen.

Začneme deskriptivní statistikou (empirické odhady kmonkrétnych teoretických parametrov a ich vhodnou vizualizíciou pomocou grafov a obrázkov). Napr.:

boxplot(veksex~pohl,data=sex, col = "lightblue")

Prípadne pre lepšiu názornosť ten istý boxplot s cenzorovanou osou ‘y’

boxplot(veksex~pohl,data=sex, col = "lightblue", ylim = c(5,30))

Následne spočítáme průměry, výberové mediány a rozptyly pro obě pohlaví (resp. odhadneme neznáme teoretické parametre strednej hodnoty, medianu a rozptylu). Použijeme k tomu funkci tapply (argument na.rm=TRUE říká, že se mají vynechat chybějící hodnoty veličiny veksex).

tapply(sex$veksex,sex$pohl,mean,na.rm=TRUE)
##     Male   Female 
## 17.35301 17.90112
tapply(sex$veksex,sex$pohl,median,na.rm=TRUE)
##   Male Female 
##     17     17
tapply(sex$veksex,sex$pohl,var,na.rm=TRUE)
##     Male   Female 
## 18.23625 15.54390

Ještě můžeme namalovat průměry věku iniciace a intervaly spolehlivosti pro střední věk iniciace pro obě pohlaví.

plotmeans(veksex~pohl,data=sex,xlab="Pohlaví",ylab="Průmerný věk sexuální iniciace", ylim = c(17, 18.2))

Nakonec porovnáme empirické distribuční funkce věku iniciace obou pohlaví.

plot(ecdf(sex$veksex[sex$pohl=="Male"]),col="blue",cex.points=0.5,verticals=TRUE,
     main="Sexualni iniciace")
lines(ecdf(sex$veksex[sex$pohl=="Female"]),col="red",cex.points=0.5,verticals=TRUE) 
legend(40,0.4,lty=1,col=c("blue","red"),legend=c("Muzi","Zeny"))

Nasvědčují podle vás popisné statistiky a grafy tomu, že se věk iniciace mužů a žen může lišit?

Formálne rozhodnutie o tom, či sa vek sexuálnej iniciácie u mužov a žien líši, je ale potrebné urobiť na základe korektného štatistického testu - dvojvýberového testu. K dispozícii je ale samozrejme množstvo rôzných modelov – štatistických testov, resp. testovacích postupoch. Pritom každý model, resp. štatistický test, má svoje predpoklady, určitú silu voči konkrétnym alternatívam a rôzne výhody, či nevýhody.

Výber konkrétneho testu (pravdepodobnostného modelu) by mal byť preto dostatočne dobre zdôvodnený, vysvetlený.

Diskutujte, z akých predpokladov nasledujúce testy vychádzajú a akú nulovú a alternatívnu hypotézu sú ideálne. Ktorý z uvedených testov považujete ako najlepší na daný problém a prečo?

Dvouvýběrový Kolmogorovův-Smirnovův test

Dvouvýběrový Kolmogorovův-Smirnovův test testuje rovnost distribučních funkcí věku iniciace mužů a žen. Jeho testová statistika je \[ K_{n,m}=\sup_{x\in\mathbb{R}}\left|\widehat{F}_{X}(x)-\widehat{F}_{Y}(x)\right|. \] Výsledky nám spočítá funkce ks.test.

ks.test(sex$veksex[sex$pohl=="Male"],sex$veksex[sex$pohl=="Female"])

Opět dostáváme varování o výskytu shodných hodnot. Věk iniciace zaokrouhlený na kalendářní rok fakticky není spojitá veličina, proto bychom Kolmogorovův-Smirnovův test raději neměli používat. Nicméně, p-hodnota je velice malá, kdybychom tomuto testu věřili navzdory porušeným předpokladům, zamítli bychom hypotézu o rovnosti distribučních funkcí s velkou rezervou.

Dvouvýběrový t-test

Proveďme nyní dvouvýběrový t-test na rovnost středních hodnot. Tento test ale predpokladá normálne rozdelenie a navyše aj zhodnosť rozptylov. Testová statistika pro hypotézu \(H_0: \mu_X=\mu_Y\) je \[ T_{n,m}=\sqrt{\frac{mn}{n+m}}\,\frac{\overline{X}_{n}-\overline{Y}_{m}}{S_{n,m}}, \] kde \[ S^2_{n,m}=\frac{n-1}{n+m-2}S^2_{X}+\frac{m-1}{n+m-2}S^2_{Y} \] je nestranný odhad společného rozptylu spočítaný z obou výběrů. Tento test předpokládá normalitu a shodné rozptyly. Ani jeden předpoklad není nejspíš splněn, i když porušení normality by při takhle velkém počtu pozorování vadit nemělo.

t.test(sex$veksex[sex$pohl=="Male"],sex$veksex[sex$pohl=="Female"],var.equal=TRUE)

I t-test zamítá rovnost středních hodnot u mužů a u žen s velkou rezervou.

Všimněte si, že výstup funkce t.test zahrnuje 95%-ní interval spolehlivosti pro rozdíl středního věku iniciace mužů a žen.

Dvouvýběrový t-test pre nezhodné rozptyly a z-test

Lépe by bylo použít dvouvýběrový z-test, který nepožaduje rovnosť rozptylov. Príslušná testová statistika je \[ Z_{n,m}=\frac{\overline{X}_{n}-\overline{Y}_{m}}{\sqrt{S^2_X/n+S^2_Y/m}}. \] Tento test dostaneme funkcí t.test, v níž vynecháme argument var.equal=TRUE (požadující shodnost rozptylů). Kritické hodnoty a p-hodnoty se počítají Welchovou aproximací.

t.test(sex$veksex[sex$pohl=="Male"],sex$veksex[sex$pohl=="Female"])



V prípade, že není opodstatnený ani predpoklad normality, ale lze předpokládat konečnost druhých momentov (t.j. , možnpsť aplikácie centrálnej limitnej vety), tak je možné využiť tzv. z-test, ktorý je založený na rovnakej testovej štatistike, ale príslušný kritický obor a p-hodnota sú počítane pouze približne, resp. asymptoticky, s využitím normálneho rozdelenia \(N(0,1)\).

Samostatne

Pomocou programu R aplikujte dvojvýberový t-test pre obe varianty – za predpokladu zhodných rozptylov aj bez tohto predpoladu. Výsledky následne porovnajte aj s asymptotickým z-testom.


Dvouvýběrový Wilcoxonův test

Dvouvýběrový Wilcoxonův test testuje shodnost rozdělení (za předpokladu, že platí model posunutí v poloze), obecně jde o test hypotézy \(H_0: P[X_i\lt Y_j]=1/2\). Tento test počítá funkce wilcox.test. Testová statistika je \[ W^*_{n,m}=\sum_{i=1}^n\sum_{j=1}^m \mathbb{I}_{\{X_i\lt Y_j\}}, \] tj. odpovídá Mann-Whitneyho formulaci Wilcoxonova testu. Kritické hodnoty a p-hodnoty se počítají normální aproximací.

wilcox.test(sex$veksex[sex$pohl=="Male"],sex$veksex[sex$pohl=="Female"],correct=FALSE)

Argument correct=FALSE zabraňuje provedení tzv. opravy na spojitost, kterou jsme do testové statistiky nezařazovali, ale R ji normálně provádí. Také Wilcoxonův test zamítá nulovou hypotézu. Výstup Wilcoxonova testu však neudává žádný odhad (bodový ani intervalový) rozdílu mezi porovnávanými populacemi.

Samostatne

Na konkrétnom datovom súbore si pripomeňte základné jednovýberové testy, parametrické aj neparametrické. Pripomeňte si potrebné teoretické základy a tiež príncíp, ako ich aplikovať na konkrétny test nulovej hypotézy.


Dvouvýběrový Kuiperův test

Existuje samozrejme množstvo dalších variant dvojvýberových testov, ktoré možu byt opoužité pre testovanie zmienenej nulovej a alternatívnej hypotézy. Jednou z možnosti je napr. tzv. Kuiperův test, ktorý vychádza z podobného princípu porovnania empirických distribučných funkcií, ako je to u Kolmogorov-Smirnovho testu (více napr. stručne tu https://en.wikipedia.org/wiki/Kuiper%27s_test).

V programe R je tento test implementovaný v knižnici kuiper.2samp’ (pre inštaláciu knižnice nutné použíť príkaz install.packages("kuiper.2samp’")).

library("kuiper.2samp’")
kuiper.2samp(sex$veksex[sex$pohl=="Male"],sex$veksex[sex$pohl=="Female"])

Samostatne

Pomocou internetu sa pokústa nájsť aspoň jednu ďalšiu alternatívu pre dvojvýberový problem a pomocou programu R tento test implementujte a porovnajte výsledky s ostatnými testami.



II. Metóda analýzy rozptylu (ANOVA)

Metóda analýzy rozptylu (z anglického The Analysis of Variance - ANOVA) nám umožňuje ověřit, zda na nějaku (spojitou) hodnotu náhodné veličiny má statisticky významný vliv hodnota jinej náhodnej veličiny (diskrétnej). Táto diskrétní veličina (resp. charakteristika, nebo znak) musí nabývat pouze konečného počtu možných hodnot (nejméně dvou). Slouží teda k rozdělení jedinců do skupin, které pak pomoci ANOVA vzájemně porovnávame.

Z určitého pohľadu teda metóda ANOVA predstavuje zobecnenie dvojvýberových testovacích problémov pre viacvýberové problémy.

Samostatne

Pripomeňte asi teoretické základy metódy analýzy rozptylu - ANOVA. Na základe akých predpokladov metóda funguje a pro testy akých hypotéz je navrhnutá? Ako by ste metódu ANOVA aplikovali a niektorý z uvažovaných datových súborov?

  • Opäť sa na problém podívajte prostredníctvom základných popisných charakteristík, následne pomocou grafických a vizuálnych nástrojov a nakoniec formálne udelajte štatistický test s využitim metody analýzy rozptylu - ANOVA.
  • V kontexte daného problému interpretujte výsledok testu a náležite vysvetlite.
  • Ako by vyzeral analogický problem ale s využitím neparametrického postupu v prípade testovania podobnej nulovej hypotézy? Aké sú predpoklady takého modelu? Líšia sa závery paramtrického a neprametrického testu?

Je úroveň vzdelania rovnaká v mestách a na venkově?

V tejto úlohe budeme používať datový súbor zam, které pocházejí ze sociálního průzkumu provedeného v ČR v roce 2004. Úroveň vzdelania v tejto úlohe budeme sledovať pomocou premennej eduyrs, ktorá zaznamenáva počet rokov strávených štúdiom u každého subjektu zvlášť. Na jednotlivé populácie definované miestom trvalého bydliska sa možeme pozrieť pomocou boxplotu.

levels(zam$domicil) <- c("Big City", "Suburbs", "Town/City", "Village", "Countryside")
cols <- rainbow(5, start = 0, end = 0.1)
boxplot(zam$eduyrs ~ zam$domicil, col = cols, main = "Dlzka vzdelania v rokoch")

Prvé štyri poplulácie sa zdajú byť rovnaké čo sa týka školskej dochádzky (v rokoch), ale vidiecké lokality vykazujú v priemere o niečo nižiu mieru vzdelanosti medzi svojími obyvateľmi (t.j. v primere menší počet rokov trávených na škole).

Pomocou metódy analýzy rozptylu sa pokúsime zistiť, či je tento rozdiel štatisticky významný, alebo nie. Využijeme k tomu najprv príkaz aov() (pre podrobnú implementáciu príkazu použijte help, t.j. ?aov). štandardne dostupných v programe R.

aov(zam$eduyrs ~ zam$domicil)
## Call:
##    aov(formula = zam$eduyrs ~ zam$domicil)
## 
## Terms:
##                 zam$domicil Residuals
## Sum of Squares       319.82  15219.46
## Deg. of Freedom           4      2794
## 
## Residual standard error: 2.333922
## Estimated effects may be unbalanced

Vyššie uvedený príkaz porovnajte ešte s výstupom nižšie, ktorý získame dvoj-kombináciou iných dvoch príkazov – anova() a lm():

anova(lm(zam$eduyrs ~ zam$domicil))
## Analysis of Variance Table
## 
## Response: zam$eduyrs
##               Df  Sum Sq Mean Sq F value    Pr(>F)    
## zam$domicil    4   319.8  79.955  14.678 7.183e-12 ***
## Residuals   2794 15219.5   5.447                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dokážete v jednotlivých výstupoch identifikovať jednotlivé súčty štvorcov a uvedené čísla interpretovať?

Ako je formulovaná nulová hypotéza, ktorú testujeme? Aký je výsledok testu a aká je interpretácia vzhľadom k pôvodnej otázke? Je možné z výsledkov testu určiť, kde a ako sa jednotlivé (sub-)populácie od seba líšia?

Je rozdíl ve vzdělání mužov a žien?

Odpoveď budeme opäť zisťovať pomocou metódy analýzy rozptylu. Najprv sa ale pozrieme na obe populácie (mužov a ženy) pomocou grafických nástrojov. Využijeme boxplot

boxplot(zam$eduyrs ~ zam$gndr, col = c("lightblue","pink"), main = "Dlzka vzdelania v rokoch")

prípadne histogram

par(mfrow = c(1,2))
hist(zam$eduyrs[zam$gndr == "Male"], col = "lightblue", xlab = "Roky stravene na skole", main = "Males", breaks = 10)
hist(zam$eduyrs[zam$gndr == "Female"], col = "pink", xlab = "Roky stravene na skole", main = "Females", breaks = 10)

Možeme skúsiť porovnať aj empirický odhad hustoty (neparametrický) s odhadom hustoty založenom na predpoklade, že rozdelenia sú oba normálne, a líšía sa iba vrámci parametru strednej hodnoty (homoskedasticita), prípadne aj rozptylu (heteroskedasticita). Analogické porovanie je možné urobiť aj vzhľadom k distribučným funciám.

Histogramy a hustoty:

par(mfrow = c(1,2))
hist(zam$eduyrs[zam$gndr == "Male"], col = "lightblue", xlab = "Roky stravené vzdelávaním", main = "Males", freq = F, ylim = c(0, 0.35), xlim = c(5,25), breaks = 20)
lines(density(zam$eduyrs[zam$gndr == "Male"], adjust = 3), col = "red", lwd = 2)
lines(dnorm(seq(0,30, length = 1000), mean(zam$eduyrs[zam$gndr == "Male"]), sd(zam$eduyrs)) ~ seq(0,30, length = 1000), col = "blue", lwd = 1, lty = 2)
lines(dnorm(seq(0,30, length = 1000), mean(zam$eduyrs[zam$gndr == "Male"]), sd(zam$eduyrs[zam$gndr == "Male"])) ~ seq(0,30, length = 1000), col = "green", lwd = 1, lty = 2)
legend(15, 0.35, c("Empirical", "Homoscedastic", "Heteroscedastic"), col = c("red", "blue", "green"), lty = c(1,2,2), cex = 0.8)

hist(zam$eduyrs[zam$gndr == "Female"], col = "pink", xlab = "Roky stravené vzdelávaním", main = "Females", freq = F, ylim = c(0,0.35), xlim = c(5,25), breaks = 20)
lines(density(zam$eduyrs[zam$gndr == "Female"], adjust = 3), col = "red", lwd = 2)
lines(dnorm(seq(0,30, length = 1000), mean(zam$eduyrs[zam$gndr == "Female"]), sd(zam$eduyrs)) ~ seq(0,30, length = 1000), col = "blue", lwd = 1, lty = 2)
lines(dnorm(seq(0,30, length = 1000), mean(zam$eduyrs[zam$gndr == "Female"]), sd(zam$eduyrs[zam$gndr == "Female"])) ~ seq(0,30, length = 1000), col = "green", lwd = 1, lty = 2)
legend(15, 0.35, c("Empirical", "Homoscedastic", "Heteroscedastic"), col = c("red", "blue", "green"), lty = c(1,2,2), cex = 0.8)

Empirické a teoretické distribučné funkcie:

par(mfrow = c(1,2))
plot(ecdf(zam$eduyrs[zam$gndr == "Male"]), col = "lightblue", xlab = "Roky stravené vzdelávaním", main = "Males", xlim = c(5, 25))
lines(pnorm(seq(0,30, length = 1000), mean(zam$eduyrs[zam$gndr == "Male"]), sd(zam$eduyrs)) ~ seq(0,30, length = 1000), col = "blue", lwd = 2, lty = 1)
lines(pnorm(seq(0,30, length = 1000), mean(zam$eduyrs[zam$gndr == "Male"]), sd(zam$eduyrs[zam$gndr == "Male"])) ~ seq(0,30, length = 1000), col = "green", lwd = 2, lty = 2)
legend(17, 0.25, c("Empirical", "Homoscedastic", "Heteroscedastic"), col = c("red", "blue", "green"), lty = c(1,2,2), cex = 0.8)

plot(ecdf(zam$eduyrs[zam$gndr == "Female"]), col = "pink", xlab = "Roky stravené vzdelávaním", main = "Females", xlim = c(5,25))
lines(pnorm(seq(0,30, length = 1000), mean(zam$eduyrs[zam$gndr == "Female"]), sd(zam$eduyrs)) ~ seq(0,30, length = 1000), col = "blue", lwd = 2, lty = 1)
lines(pnorm(seq(0,30, length = 1000), mean(zam$eduyrs[zam$gndr == "Female"]), sd(zam$eduyrs[zam$gndr == "Female"])) ~ seq(0,30, length = 1000), col = "green", lwd = 2, lty = 2)
legend(17, 0.25, c("Empirical", "Homoscedastic", "Heteroscedastic"), col = c("red", "blue", "green"), lty = c(1,2,2), cex = 0.8)



Pomocou vhodného štatistického testu sa pozrieme, či je rozdiel, ktorý prostredníctvom histogramov sledujeme, statistický významný, alebo je zanedbateľný. Test založíme na vyhodnotení rozptylu v jednotlivých skupinách - tzv. metoda analýzy rozptylu - ANOVA.

anova(lm(zam$eduyrs ~ zam$gndr))
## Analysis of Variance Table
## 
## Response: zam$eduyrs
##             Df Sum Sq Mean Sq F value    Pr(>F)    
## zam$gndr     1    120 120.049  21.777 3.207e-06 ***
## Residuals 2797  15419   5.513                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Keďže porovnávame iba dve populácie, môžeme použiť aj klasický dvojvýberový postup.

t.test(zam$eduyrs ~ zam$gndr, paired = F, var.equal = T)
## 
##  Two Sample t-test
## 
## data:  zam$eduyrs by zam$gndr
## t = 4.6665, df = 2797, p-value = 3.207e-06
## alternative hypothesis: true difference in means between group Male and group Female is not equal to 0
## 95 percent confidence interval:
##  0.2408303 0.5898852
## sample estimates:
##   mean in group Male mean in group Female 
##             12.41004             11.99468

Ktoré výsledky sú vo výstupe analýzy rozptylu a t.testu rovnaké a prečo? Aká je interpretácia výsledkov vzhľadom k pôvodnej otázke?

Samostatne

  • Metóda analýzy rozptylu pracuje s konkrétny štatistickým modelom. Ako tento model vyzerá a na akých predpokladoch je založeny?
  • Aké následky má nesplnenie jednotlivých predpokladov modelu ANOVA na výsledky a interpretáciu výsledkov získaných v teste pomocou metody analýzy rozptylu?



Neparametrické verze pro analýzu rozptylu

Neparametrická verzia metódy analýzy rozptylu – porovnávanie niekoľkých skupín – je implementována v programe R pomocou funkcie kruskal.test(). Pripomeňte si, s akým modelom (t.j., predpokladmi) Kruskal-Wallisov test pracuje. Vpodstate sa jedná o test nulovej hypotézy

\[ H_{0}:~F_1(x) = F_2(x) = \dots = F_K(x), \quad \textrm{pre všetky hodnoty $x \in \mathbb{R}$} \] na základe (vzájomne nezávislých) náhodných výberov \(X_{k 1}, \dots, X_{k n_k} \sim F_k\), kde \(k = 1, \dots, K \in \mathbb{N}\) a \(n_k \in \mathbb{N}\) (teda pre aspoň \(K \geq 2\) rôznych skupín).



Alternatvívna hypotéza je formulovaná tak, že aspoň jeden náhodný výber sa od ostatných líší – mysleno v zmysle toho, že distribučná funkcia (aka rozdelenie) niektorého náhodného výberu sa líši od rozdelenia ostatných náhodných vyberov. Uvedomte si, že za platnosti nulovej hypotézy tvoria všetky pozorovania t.j. \(X_{1,1}, \dots, X_{K,n_K}\) dohromady jeden náhodný výber (t.j., súbor nezávislých a rovnako rozdelených náhodných veiičín).



Implementácia v programe R:

kruskal.test(zam$eduyrs ~ zam$domicil)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  zam$eduyrs by zam$domicil
## Kruskal-Wallis chi-squared = 59.495, df = 4, p-value = 3.703e-12



Samostatne

  • Aplikujte neparametrický postup - tzv. Kruskal-Wallisov test na oba príklady uvedené vyššie a podrobne interpretujte výsledky.
  • Líšia sa výsledky testov založených na metode ANOVA a na neparmetrickom postupe pomocou Kruskal-Wallisovho testu?
    Aká je interpretácia výsledkov vzhľadom k celkovému kontextu experimentu?