NMFM301, ZS 2025/26

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ť užitočné pri spracovaní samostatnej zápočtovej práce budúci týždeň (hlavne v prípade študentov, ktorý sa rozhodnú zadanie spracovať v zdrojovom súbore samotného zadania).

Užitočné materiály pre prácu so štatistickým programom 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 tohto cvičenia využijeme dva konkrétne datové súbory a niektoré dodatočné funkcie (resp. 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 – ich účelom je hlavne zjednodušiť a prípadne zrýchliť 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/functions2.RData"))

K dispozíci sú dva datové súbory, 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 menší pod-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, ktoré sa dajú formulovať v zmysle použitia určitého štatistického testu. Doležitý je správny výber pravdepodobnostného modelu a explicitná formulácia otázky pomocou nulovej a alternatívnej hypotézy, testovej štatistiky a celkového rozdnutia (t.j., zamítam/nezamítam nulovou hypotézu).

Samotnej voľbe pravdepodobnostného modelu (t.j., konkrétneho štatistického testu) samozrejme predchádza exploratívna analýza pomocou vhodných popisných charakteristík (t.j., odhadov konkrétnych, ale neznámych parametrov) a obrázkov.

Datové súbory zam a sex si prohlédněte poklepáním na příslušný řádek v okénku “Workspace” (v RStudio) nebo pomoci příkazů summary(), head(), str(), a pod. (v konzole programu R). Pomocou niekoľkých obrázkov a grafov preskúmajte charakter a štruktúru dat – 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ě \(n = 4467\) Američanů starších 20 let. Obsahuje veličiny/proměnné active (byl někdy respondent sexuálně aktivní?), veksex (věk sexuální iniciace – pouze u aktivních), pocpart (počet partnerů v průběhu života), pocpart.rok (počet partnerů za uplynulý rok) a dále demografická data (věk, pohlaví, vzdělání, rodinný status).

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

Začneme deskriptivní statistikou (empirické odhady konkrétnych teoretických parametrov0 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))

Boxplot zobrazuje robustné charakteristiky spočítané z náhodného výberu – odhad neznámeho (teoretického) mediánu v celej populaci, odhad prvního a třetího kvartilu (t.j., okraje boxplot krabice) nebo informaci o prípadných odlehlých pozorovaniach. Takýto boxplot ale môžeme doplniť aj o odhady menej robustných charakteristík – napr. odhad teoretickej strednej hodnoty:

boxplot(veksex~pohl,data=sex, col = "lightblue", ylim = c(5,30), xlab = "Pohlavie respondenta", ylab = "Vek [v rokoch]")
lines(c(0.5, 1.5), rep(mean(sex$veksex[sex$pohl == "Male"], na.rm = T), 2), col = "red", lwd = 3)
lines(c(1.5, 2.5), rep(mean(sex$veksex[sex$pohl == "Female"], na.rm = T), 2), col = "red", lwd = 3)



Průměry v závislosti na pohlaví, prípadne výberové mediány, rozptyly, či směrodatné chyby pro obě pohlaví spočítame (resp. odhadneme neznáme teoretické parametre strednej hodnoty, medianu, rozptylu a směrodatnej chyby) aj pomocou funkcie 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
tapply(sex$veksex,sex$pohl,sd,na.rm=TRUE)
##     Male   Female 
## 4.270392 3.942575

Ještě můžeme vykřeslit průměry věku iniciace a příslušné intervaly spolehlivosti (pro \(\alpha = 0.05\), což je defaultní nastavení funkce) pro střední věk iniciace pro obě pohlaví.

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

Na pôvodnej škále osi \(y\) (ako to bolo v prípade boxplotov) je tento rozdíl ale v podstatě nepatrný – napriek tomu sa ale asi bude jednat o statisticky významný rozdíl, kedže intervaly spolehlivosti se vzájemne nepřekrývají.

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

Môžeme ješte porovnat empirické distribuční funkce pro věk iniciace v obou skupinách:

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 (myslené statisticky průkazně líšit)?

Formálne rozhodnutie o tom, či sa vek sexuálnej iniciácie u mužov a žien líši, je potrebné urobiť na základe korektného štatistického testu – v tomto prípade vhodného dvojvýberového testu. K dispozícii je ale samozrejme množstvo rôzných modelov (štatistických testov), resp. testovacích postupoch. Každý model, resp. štatistický test, má ale svoje (teoretické) predpoklady, určitú silu voči konkrétnym alternatívam a prípadne rôzne ďalšie výhody a nevýhody. Toto všetko by malo byť pri voľbe konkrétneho testu nejak zohľadnené. Vo väčšine prípadov sa ale nedá jednoznačne povedať, ktorý model, resp. štatistický test, je tym najlepším, resp. jediným správnym modelom.

Výber konkrétneho testu by mal byť preto dostatočne dobre zdôvodnený.

Za akých predpokladov sú formulované nasledujúce testy? Ako je definovaná nulová a alternatívni hypotéza? Ktorý z uvedených testov považujete pre tento konkrétny problém jako najvhodnější? Vysvětlete proč?

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ů. Testová štatistika má za platnosti nulovej hypotézy študentovo \(t\) rozdelenie s \(n + m - 2\) stup+nami voľnosti. Štandardne ale tento test předpokládá normalitu jednotlivých náhodných výberov 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 (test by fungoval alespoň asymptoticky). Test spočítame pomocou funkcie \(t.test()\).

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

Všimněte si, že výstup funkce t.test zahrnuje i 95%-ní interval spolehlivosti pro rozdíl středního věku iniciace mužů a žen. To je užitočné hlavne pre interpretáciu výsledku testu v prípade, že nulové hypotéza je zamietnutá,

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

Mono by bylo lépe využí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}}. \] Takáto testová statistika má za platnosti nulovej hypotézy opět studentovo \(t\) rozdělení, ale počet stupňu volnosti není \(n + m - 2\), ale stupn+e volnosti jsou spočítaný pomoci tzv. Welchovej aproximace. Tento test získame pomoci funkce t.test(), v níž vynecháme argument var.equal=TRUE (požadující shodnost rozptylů).

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žnosť využitia centrálnej limitnej vety), tak je možné aplikovať tzv. z-test, ktorý je založený na rovnakej testovej štatistike, ale príslušný kritický obor a p-hodnota sú počítané 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. Výsledky vhodne interpretujte.


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 prípadne pokúste nájsť ďalšie alternatívy štatistických testov pre daný dvojvýberový problém. Pomocou programu R tento test skúste implementovať a opět porovnat výsledky.



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 – teda napr. \(K \in \mathbb{N}\), pro \(K \geq 2\)). Slouží teda k rozdělení jedinců do \(K\) disjunktných skupin, které pak pomoci ANOVA vzájemně porovnávame. V podstate sa jedná o tzv. K-výběrový problém. Metóda ANOVA teda 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 na 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?

Príklad 1: Je úroveň vzdelania rovnaká v mestách a na venkově?

V tejto úlohe budeme používať datový súbor zam, ktorý predstavuje podsúbor jedincov pocházejích 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ť napr. 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")

Metóda ANOVA ale vo svojej základnej verzii predstavuje zobecnenie dvojvýberového \(t\)-testu (t.j., za predpokladu normality náhodných výberov a taktiež za predpokladu zhodnosti rozptylov). V podstate sa jedna o test nulovej hypotézy \[ H_0: \mu_1 = \mu_2 = \dots = \mu_K \]

oproti alternatíve, že alespoň jedna z uvedených rovností neplatí. Test je založený na ``datech’’, ktoré predstavujú \(K\) nezávislých náhodných výberov \(\{X_{k, 1}, \dots, X_{k, n_k}\}_{k = 1}^K\), kde \(n_k \in \mathbb{N}\) je velikost (rozsah) \(k\)-teho náhodného výberu a plati, že \(X_{k, i} \sim N(\mu_k, \sigma^2)\), pro \(i = 1, \dots, n_k\) a \(k = 1, \dots, K\). Parametre \(\mu_{k} \in \mathbb{R}\) predstavují neznáme střední hodnoty.

V boxplote vyššie príslušné odhady neznámych stredných hodnôt \(\mu_1, \dots, \mu_K\) ale nevidíme. Môžeme ich ale jednoducho získať (spolu s príslušnými odhadmi směrodatných chýb) napr. pomocou príkazov:

tapply(zam$eduyrs, zam$domicil,mean)
##    Big City     Suburbs   Town/City     Village Countryside 
##    12.66859    12.72464    12.20735    11.73840    11.50000
tapply(zam$eduyrs, zam$domicil,sd)
##    Big City     Suburbs   Town/City     Village Countryside 
##    2.484313    2.599180    2.328049    2.182007    2.880972



Prvé štyri poplulácie sa zdajú byť pomerne rovnaké čo sa týka školskej dochádzky (v rokoch je to rádovo 12 rokov), ale vidiecké lokality vykazujú v priemere o niečo nižiu mieru vzdelanosti (približne o pol roka) a zároveň sledujeme najväčšiu variabilitu (a aj vychýlenie smerom k nižším hodnotám). Pripomeňme, že klasická metóda ANOVA predpokladá rovnosť teoretických rozptylov.

Tieto zistenia budeme nateraz ignorovať (metóda ANOVA totíž funguje aj na asymptotickom princípe a s modifikáciou pre rôzne rozptyly) a pomocou analýzy rozptylu sa pokúsime zistiť, či je rozdiel, ktorý medzi jednotlivými skupinami sledujeme, š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ý 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



Samostatne

Pokúste sa zreprodukovať jednotlivé hodnoty, ktoré sú obsiahnuté vo výstupoch a pripomeňte si základne teoretické princípy metódy ANOVA.

Ako je definovaná testová štatistika? Aké rozdelenie za platnosti nulovej hypotézy testová štatistika má? 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?

manuálna rekoknštrukcia testu:

X1 <- zam$eduyrs[zam$domicil == "A big city"] ### 1. nahodny vyber 
X2 <- zam$eduyrs[zam$domicil == "Suburbs or outskirts of big city"] ### 2. nahodny vyber 
X3 <- zam$eduyrs[zam$domicil == "Town or small city"] ### 3. nahodny vyber 
X4 <- zam$eduyrs[zam$domicil == "Country village"] ### 4. nahodny vyber 
X5 <- zam$eduyrs[zam$domicil == "Farm or home in countryside"] ### 5. nahodny vyber

Príslušné odhady stretných hodnôt pre jednotlivé skupiny aj celkový súbor:

mean1 <- mean(X1)
mean2 <- mean(X2)
mean3 <- mean(X3)
mean4 <- mean(X4)
mean5 <- mean(X5)

meanTotal <- mean(zam$eduyrs)

Počty pozorovaní v jednotlivých náhodných výberoch:

n1 <- length(X1)
n2 <- length(X2)
n3 <- length(X3)
n4 <- length(X4)
n5 <- length(X5)

N <- n1 + n2 + n3 + n4 + n5

Pomocou vyššie spočítaných hodnôt už môžeme zrekonštruovať jednotlivé hodnoty vo výstupoch príkazov aov(...) a anova(lm(...)):

### Sum of squares: 
(SS <- sum(c(n1, n2, n3, n4, n5) * (c(mean1, mean2, mean3, mean4, mean5) - meanTotal)^2))
## [1] NaN
### Residuals: 
(SE <- sum((X1 - mean1)^2) + sum((X2 - mean2)^2) + sum((X3 - mean3)^2) + sum((X4 - mean4)^2) + sum((X5 - mean5)^2))
## [1] 0
### Mean Squares: 
SS / 4
## [1] NaN
SE / 2794
## [1] 0

Testová štatistika \(F = 79.955\) je definovaná ako podiel stredných štvorcov, teda hodnoty \(79.955\) a hodnoty \(5.447\) a má Fisherovo \(F\) rozdelenie so \(4\) a \(2794\) stupňami voľnosti. Uvedená \(p\)-hodnota je teda spočítaná ako:

1 - pf(14.678, df1 = 4, df2 = 2794)
## [1] 7.185696e-12



Príklad 2: 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, 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?