NMFM301, ZS 2022/23

Cvičenie 10 (praktické cvičenie 4)

(dvojvyberové testy, kontingenčné tabuľky, ANOVA)

Zdrojový soubor pro Knit: Rmd soubor

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)


Nasledujúce dva príkazy načítajú potrebné data a dodatočné R-kové (príkazy) funkcie, ktoré budú užitočné v ďalšej analýze.

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

Nyní máte k dispozici datové soubory sex a zam. Data sex pocházejí ze zdravotního průzkumu provedeného v USA v roce 2012. Obsahují údaje o respondentech starších 20 let. Data zam pocházejí ze sociálního průzkumu provedeného v ČR v roce 2004.

Samostatne

  • Použijte základné popisné charakteristiky a aplikujte ich na data. Pomocou týchto popisných charakteristik sa pokúste získať základný prehlad o štruktúre dát a jednotlivých nameraných hodnotách.
  • Použijte niektoré grafické nástroje a pomocou nich sa pokúste zmysluplne poskytnúť ucelený obraz o datach, ktoré máte k dispozícii.

  • Data si môžete prezriť aj poklepáním na příslušný řádek v okénku “Workspace” nebo použitím niektorých príkazov v ‘head()’, ‘summary()’, ‘str()’, a alších. Napríklad

    str(sex)
summary(sex)
    summary(sex)
summary(zam)



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 sexuální iniciaci (prvním pohlavním styku) mužů a žen.

Začneme deskriptivní statistikou. Nakreslíme si krabicový diagram věku iniciace pro obě pohlaví

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))

a spočítáme průměry, mediány a rozptyly pro obě pohlaví. 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="Pohlavi",ylab="Prumerny vek iniciace")

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 niekoľko rôzných testovacích postupoch. Každý test ma svoje predpoklady a tiež sílu.

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. 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ý z-test

Lépe by bylo použít dvouvýběrový z-test, který nepožaduje normalitu ani shodnost rozptylů. Jeho 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"])

Welchův test dává prakticky stejné výsledky jako dvouvýběrový t-test se shodnými rozptyly. I zde máme k dispozici 95%-ní interval spolehlivosti pro rozdíl středních hodnot.

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.



II. Odhadování a testování pravděpodobností a šancí

Nyní budeme pracovat s datovým souborem zam. Tento soubor obsahuje data ze sociologického průzkumu provedeného v ČR v roce 2004. Zahrnuje veličiny idno (identifikační kód), tvtot (jak dlouho se obvykle denně díváte na televizi?), gndr (pohlaví), agea (věk), domicil (velikost sídla bydliště), eduyrs (počet let vzdělání), unempl (je nyní nezaměstnaný/á), regioncz (kraj bydliště).

V nasledujúcom sa budeme zabývat nezaměstnaností.

Samostatne

  • Pomocou základných popisných charakteristík sa opäť podívajte na data a pokúste sa popísať základnú štruktúru vzhľadom k nezamestnanosti a alším premenným, ktoré sú v datách k dispozícii.
  • Použijte grafické nástroje a pokuste sa navrhnúť vhodnú vizuálizáciu nezamestnanosti vhľadom k ostatným populačno-sociologickým charakteristikám, ktoré sú v datách obsiahnuté.



Vypišme si, kolik lidí je nezaměstnaných:

table(zam$unempl)
## 
##    0    1 
## 2645  154

Je jich 154 z 2799 respondentů. Odhad pravděpodobnosti, že respondent je nezaměstnaný je \(\widehat{p}_n=154/2799=0.055\). Odhad šance na nezaměstnanost je \(\widehat{p}_n/(1-\widehat{p}_n)=0.0582\). Protože pravděpodobnost nezaměstnanosti je malá, šance je blízko hodnoty pravděpodobnosti.

Budeme testovat, zdali je míra (pravděpodobnost) nezaměstnanosti rovna 0.05. Použijeme funkci propor (není součástí R), která počítá přibližný interval spolehlivosti pro pravděpodobnost a testuje hypotézu \(H_0: p_X=p_0\) dvěma způsoby: jednak přímo pomocí věty 7.1, jednak přes šanci pomocí věty 7.2. Vstupní argumenty jsou tři: počet úspěchů \(X_n\), počet pokusů (velikost populace) \(n\), a hypotetická pravděpodobnost \(p_0\).

x <- sum(zam$unempl)
n <- nrow(zam)
propor(x,n,p0=0.05)
## 
## Testovani pravdepodobnosti
## 
## 154 uspechu z 2799 pokusu
## 
## Odhadnuta pravdepodobnost: 0.055 
## Odhadnuta sance          : 0.0582 
## Odhadnuta log-sance      : -2.84 
## 
## 95%-ni as. interval spolehlivosti pro pravdapodobnost - metoda 1:
## ( 0.0466 , 0.0635 )
## 
## As. test hypotezy p=0.05 proti p!=0.05 - typ 1:
## Testova statistika 1.165, p-hodnota 0.244
## 
## 95%-ni as. interval spolehlivosti pro pravdepodobnost - metoda 2:
## ( 0.0472 , 0.0641 )
## 
## As. test hypotezy p=0.05 proti p!=0.05 - typ 2:
## Testova statistika 1.218, p-hodnota 0.223

Hypotézu, že pravděpodobnost nezaměstnanosti je 0.05, nemůžeme zamítnout. Obě metody dávají velmi podobné výsledky jak co se týče výsledku testu, tak co se týče intervalu spolehlivosti pro parametr \(p_X\). Statistický program R nám dává i možnost otestovat hypotézu přesným testem založeným na binomickém rozdělení a získat přesný interval spolehlivosti. Toto provádí funkce binom.test. Její vstupní argumenty jsou stejné, jako u funkce propor (jen poslední argument se jmenuje p místo p0).

binom.test(x,n,p=0.05)
## 
##  Exact binomial test
## 
## data:  x and n
## number of successes = 154, number of trials = 2799, p-value = 0.2245
## alternative hypothesis: true probability of success is not equal to 0.05
## 95 percent confidence interval:
##  0.04686276 0.06412182
## sample estimates:
## probability of success 
##             0.05501965

Její výsledky jsou velmi podobné oběma metodám z funkce propor, neboť pracujeme s dosti velkým počtem pozorování.

Součástí základní distribuce R je funkce prop.test, která provádí asymptotický test metodou 2, ale uvádí čtverec testové statistiky a používá asymptotické rozdělení \(\chi^2_1\):

prop.test(x,n,p=0.05,correct=FALSE)
## 
##  1-sample proportions test without continuity correction
## 
## data:  x out of n, null probability 0.05
## X-squared = 1.4848, df = 1, p-value = 0.223
## alternative hypothesis: true p is not equal to 0.05
## 95 percent confidence interval:
##  0.04716603 0.06409301
## sample estimates:
##          p 
## 0.05501965

Nyní budeme zkoumat nezaměstnanost mezi některými subpopulacemi (což povede ke zmenšování rozsahu výběru).

Nezaměstnanost ve Zlínském kraji: [Nesmíme zapomenout správně měnit rozsah celé zkoumané populace \(n\).]

x <- sum(zam$unempl[zam$regioncz=="Zlin Reg."])
n <- sum(zam$regioncz=="Zlin Reg.")
propor(x,n,p0=0.05)
## 
## Testovani pravdepodobnosti
## 
## 4 uspechu z 113 pokusu
## 
## Odhadnuta pravdepodobnost: 0.0354 
## Odhadnuta sance          : 0.0367 
## Odhadnuta log-sance      : -3.31 
## 
## 95%-ni as. interval spolehlivosti pro pravdapodobnost - metoda 1:
## ( 0.00133 , 0.0695 )
## 
## As. test hypotezy p=0.05 proti p!=0.05 - typ 1:
## Testova statistika -0.84, p-hodnota 0.401
## 
## 95%-ni as. interval spolehlivosti pro pravdepodobnost - metoda 2:
## ( 0.0133 , 0.0905 )
## 
## As. test hypotezy p=0.05 proti p!=0.05 - typ 2:
## Testova statistika -0.7083, p-hodnota 0.479
binom.test(x,n,0.05)
## 
##  Exact binomial test
## 
## data:  x and n
## number of successes = 4, number of trials = 113, p-value = 0.6645
## alternative hypothesis: true probability of success is not equal to 0.05
## 95 percent confidence interval:
##  0.009727554 0.088156443
## sample estimates:
## probability of success 
##             0.03539823

Nezaměstnanost na vesnicích ve Zlínském kraji:

x <- sum(zam$unempl[zam$regioncz=="Zlin Reg." & zam$domicil=="Country village"])
n <- sum(zam$regioncz=="Zlin Reg." & zam$domicil=="Country village")
propor(x,n,p0=0.05)
## 
## Testovani pravdepodobnosti
## 
## 4 uspechu z 75 pokusu
## 
## Odhadnuta pravdepodobnost: 0.0533 
## Odhadnuta sance          : 0.0563 
## Odhadnuta log-sance      : -2.88 
## 
## 95%-ni as. interval spolehlivosti pro pravdapodobnost - metoda 1:
## ( 0.00248 , 0.104 )
## 
## As. test hypotezy p=0.05 proti p!=0.05 - typ 1:
## Testova statistika 0.1285, p-hodnota 0.898
## 
## 95%-ni as. interval spolehlivosti pro pravdepodobnost - metoda 2:
## ( 0.0202 , 0.134 )
## 
## As. test hypotezy p=0.05 proti p!=0.05 - typ 2:
## Testova statistika 0.1324, p-hodnota 0.895

Tady už je rozsah výběru tak malý (pouze 4 ,,úspěchy’’), že bychom rozhodně neměli používat asymptotické metody. Ale, jak ukazuje přesný test (dole), metoda 2 stále dává použitelné výsledky.

binom.test(x,n,0.05)
## 
##  Exact binomial test
## 
## data:  x and n
## number of successes = 4, number of trials = 75, p-value = 0.7899
## alternative hypothesis: true probability of success is not equal to 0.05
## 95 percent confidence interval:
##  0.01472076 0.13096081
## sample estimates:
## probability of success 
##             0.05333333

Nezaměstnanost mezi ženami na vesnicích ve Zlínském kraji:

x <- sum(zam$unempl[zam$regioncz=="Zlin Reg." & zam$domicil=="Country village" & zam$gndr=="Female"])
n <- sum(zam$regioncz=="Zlin Reg." & zam$domicil=="Country village"  & zam$gndr=="Female")
propor(x,n,p0=0.05)
## 
## Testovani pravdepodobnosti
## 
## 1 uspechu z 32 pokusu
## 
## Odhadnuta pravdepodobnost: 0.0312 
## Odhadnuta sance          : 0.0323 
## Odhadnuta log-sance      : -3.43 
## 
## 95%-ni as. interval spolehlivosti pro pravdapodobnost - metoda 1:
## ( -0.029 , 0.0915 )
## 
## As. test hypotezy p=0.05 proti p!=0.05 - typ 1:
## Testova statistika -0.6096, p-hodnota 0.542
## 
## 95%-ni as. interval spolehlivosti pro pravdepodobnost - metoda 2:
## ( 0.00438 , 0.191 )
## 
## As. test hypotezy p=0.05 proti p!=0.05 - typ 2:
## Testova statistika -0.4818, p-hodnota 0.63
binom.test(x,n,0.05)
## 
##  Exact binomial test
## 
## data:  x and n
## number of successes = 1, number of trials = 32, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.05
## 95 percent confidence interval:
##  0.0007908686 0.1621709942
## sample estimates:
## probability of success 
##                0.03125

V tomto extrémním případě už metoda 1 zcela selhává - interval spolehlivosti zahrnuje i záporná čísla a jeho horní okraj je špatný ve srovnání s přesným intervalem.

Máte-li čas, otestujte, zdali pravděpodobnost, že občan má aspoň 12 let školní docházky je rovna jedné polovině. Proveďte to v celé populaci, mezi ženami, mezi ženami staršími 50 let, a mezi ženami staršími 50 let žijícími na vesnici.

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?