NMFM301, ZS 2024/25
Cvičenie 10 (praktická časť 4)
(testy proporcií, kontingenčné tabuľky a \(\chi^2\) testy)
Zdrojový soubor pro Knit: Rmd soubor
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)
Úvod/Intro
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_data.RData"))
load(url("http://msekce.karlin.mff.cuni.cz/~maciak/NMFM301/functions.RData"))Po načítaní súboru cv4_data.RData máme v programe R k
dispozici tri datové soubory deti, kraje a
zam. Data deti udávají počty dětí narozených v
Československu v jednotlivých měsících roku 1970 (zdroj: Anděl,
Matematická statistika, SNTL 1977). Datový súbor kraje
pochádza z veřejných databází Českého statistického úřadu. Nakoniec
datový súbor zam pochádza ze sociálního průzkumu
provedeného v ČR v roce 2004.
Samostatne
- Použijte základné popisné charakteristiky a pomocou týchto popisných charakteristik sa pokúste získať základný prehlad o štruktúre jednotlivých datových súborov a nameraných hodnotách (t.j., jednoduchá exploratívna analýza).
- Použijte niektoré grafické nástroje a pomocou nich sa pokúste zmysluplne poskytnúť ucelený obraz o datach, ktoré máte k dispozícii.
-
Využijte príkazy
dim(), ‘head()’, ‘summary()’, ‘str()’, a ďalšíe, Zamyslite sa nad vhodnou grafickou reprezentáciou diskretných premenných (t.j., absolútných a relatívných počtov v rôznych kategóriach).
I. Odhadování a testování pravděpodobností a šancí
Ako prvý využijeme datový súbor zam, ktorý 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 v jednotlivých regiónoch ČR..
- 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é.
Začneme s celkovým pohľadom na počet zamestnaných a nezaměstnaných ľudi obsiahnutých v prieskume:
table(zam$unempl)##
## 0 1
## 2645 154resp. relatívne počty, ktoré mávajú často viac výpovednú hodnotu (ale maly by byť doplnené informáciou o celkovom počte pozorovaní, t.j., \(n = 2799\)):
format(100 * table(zam$unempl)/nrow(zam), digits = 2, trim = T)##
## 0 1
## "94.5" "5.5"Nezamestnaných je teda 154 z 2799 respondentov a teda príslušný odhad
pravděpodobnosti, že jedinec v populácii (uvedomte si rozdiel,
medzi jedincom v populácii a respondentom v náhodnom výbere) 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.
Uvedomte si rozdiel medzi pojmami pravdepodobnosť, šanca a pomer šancí (z angl. probability, odds a odds ratio).
Budeme testovat, zdali je míra (pravděpodobnost) nezaměstnanosti
rovna 0.05. Použijeme funkci propor (není součástí
štandardnej inštalácie programu R, ale je obsiahnuta v
súbore functions.RData, ktorý bol načítaný v úvode), 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 (asymptotická normalita odhadu \(\widehat{p}_n\)), jednak přes šanci pomocí
Věty 7.2 (delta veta pre transformáciu parametru \(p\) na logaritmus šance, t.j., \(\theta = \log \frac{p}{1 - p}\)). 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.223Hypoté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.05501965Vý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.05501965Znížovanie rozsahu \(n \in \mathbb{N}\) a vyšetrovanie sily jednotlivých testov
Nyní budeme zkoumat nezaměstnanost mezi některými subpopulacemi (což vlastne povede ke zmenšování rozsahu uvažovaného náhodného výběru a následne menší sile jednotlivých testov).
Napr., zameriame sa pouze na zaměstnanost/nezam+estnanost ve Zlínském kraji (pozor na správny rozsah zkoumané populace – t.j., zmena hodnoty \(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.479binom.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.03539823Nezamě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.895Tady už je rozsah výběru tak malý (pouze 4 ,,úspěchy’’), že bychom rozhodně neměli používat asymptotické metody. Porovnanie s presným testom (dole) ale ukazuje, že i asymptotické metódy stále dávajú použitelné výsledky (celkový počet pozorovaní je totiž \(4 + 71 = 75\), takže celkom postačujúci počet pre aproximáciu pomocou centrálnej limitnej vety – i keď zastsúpenie v jednotlivých kategóriach je dosť nevyvážené). Aké to má teoretické následky/praktické ná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.05333333Nezamě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.63binom.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.03125V 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.
II. Kontingenčne tabuľky a \(\boldsymbol{\chi^2}\) testy nezávislosti a dobré zhody
V tejto časti budeme používať datovy súbor deti, který
udáva počty dětí narozených v Československu v jednotlivých měsících
roku 1970. Pomocou grafických nástrojov bežne dostupných v programe R sa
na data podívajte a pokuste se data zobrazit vždy tak, aby výsledný graf
poskytoval určité možnosti jak nahlédnout do struktúry, které se týka
daný statistický test.
Rodí se děti rovnoměrně během roku?
Datový súbor obsahuje len 12 pozorovaní, ktoré predstavujú súhrnné počty narodených detí pre každý mesiac v roku 1970.
barplot(deti$poc.deti,names.arg = c("Leden", " Únor", "Březen", "Duben", "Květen", "Červen", "Červenec", "Srpen", "Žárí", "Říjen", "Listopad", "Prosinec"), col = "lightblue", cex.names = 0.8, ylim = c(0,25000))
abline(mean(deti$poc.deti), 0, col = "red", lwd = 2, lty =2)
text(13.5, 21800, "Celkový priemer", col = "red")
Ověříme, jestli se děti rodí rovnoměrně během roku. Pokud by tomu tak bylo, pravděpodonosti jednotlivých měsíců by byly přímo úměrné počtu dní v měsíci. Zapíšeme tedy počty dní v měsících do vektoru a ověříme, že součet dá 365 (rok 1970 nebyl přestupný). Vieme teda rozdelenie pravdepodobnosti, ktoré zodpovedá nulovej hypotéze. Z matematického/štatistického hľadiska preto využijeme tzv. \(\chi^2\) test dobrej zhody.
dni <- c(31,28,31,30,31,30,31,31,30,31,30,31)
sum(dni)## [1] 365Nyní napočítáme pravděpodobnosti měsíců za nulové hypotézy.
psti <- dni/sum(dni)
psti## [1] 0.08493151 0.07671233 0.08493151 0.08219178 0.08493151 0.08219178
## [7] 0.08493151 0.08493151 0.08219178 0.08493151 0.08219178 0.08493151Popripade tie ist pravdepodobnosti s presnosťou na štyri desatinné miesta:
round(psti, 4)## [1] 0.0849 0.0767 0.0849 0.0822 0.0849 0.0822 0.0849 0.0849 0.0822 0.0849
## [11] 0.0822 0.0849Odhadnuté pravdepodobnosti z datového súboru (opäť na štyri desatinne miesta):
round(deti$poc.deti / sum(deti$poc.deti), 4)## [1] 0.0838 0.0790 0.0902 0.0902 0.0915 0.0865 0.0845 0.0805 0.0829 0.0793
## [11] 0.0741 0.0775
Budeme testovať nulovô hypotézu, že vektor pravděpodobností
multinomického rozdělení, které vygenerovalo počty narozených dětí, je
roven vektoru pravdepodobnosti, ktorý máme uložený v R objekte
psti. Formálne (matematicky) zapísané, zaujíma nas test
nulovej hypotézy \[
H_0:~\boldsymbol{p} = \boldsymbol{p}_0,
\] kde \(\boldsymbol{p} = (p_1, \dots,
p_{12})^\top \in [0,1]^{12}\) je vektor pravdepodobnosti, že sa
náhodne vybrané dieťa narodí v niektorom z 12tich mesiacov. Samozrejme
platí, že \(\sum_{i = 1}^{12} p_i =
1\), inak by sa totiž nejednalo o pravdepodobnostné rozdelenie.
Pripomeňte si teoretické základy \(\chi^2\) testu dobrej zhody.
V programe R oužijeme funkci chisq.test. Její první
argument je vektor pozorovaných četností a druhý, argument
p, obsahuje hypotetické pravděpodobnosti - vektor
pravdepodobnosti v multinomickom rozdelení za platnosti nulovej
hypotézy.
chisq.test(deti$poc.deti,p=psti)##
## Chi-squared test for given probabilities
##
## data: deti$poc.deti
## X-squared = 1004.1, df = 11, p-value < 2.2e-16Testová statistika je velmi velká (referenční rozdělení \(\chi^2_{11}\) má střední hodnotu 11), p-hodnota je takřka nulová. Hypotézu tedy s velkou rezervou zamítáme.
Podívejme se blíže na způsob výpočtu testové statistiky. Testová štatistika \(\chi^2\) testu dobrej zhody, je definovaná predpisom \[ \chi^2=\sum_{k=1}^K\frac{(X_k-np^0_k)^2}{np^0_k}. \]
Pre lepšie pochopenie pomůže nám funkce kuchej.chisq
(není součástí R).
kuchej.chisq(deti$poc.deti,psti)## skupina cetnost pst ocek rozdil statistika
## 1 1 21182 0.08493151 21465.59 -2.835890e+02 3.7465892
## 2 2 19960 0.07671233 19388.27 5.717260e+02 16.8591929
## 3 3 22787 0.08493151 21465.59 1.321411e+03 81.3453998
## 4 4 22805 0.08219178 20773.15 2.031849e+03 198.7378661
## 5 5 23120 0.08493151 21465.59 1.654411e+03 127.5099237
## 6 6 21859 0.08219178 20773.15 1.085849e+03 56.7592636
## 7 7 21367 0.08493151 21465.59 -9.858904e+01 0.4528084
## 8 8 20357 0.08493151 21465.59 -1.108589e+03 57.2530136
## 9 9 20946 0.08219178 20773.15 1.728493e+02 1.4382453
## 10 10 20037 0.08493151 21465.59 -1.428589e+03 95.0762005
## 11 11 18728 0.08219178 20773.15 -2.045151e+03 201.3484323
## 12 12 19592 0.08493151 21465.59 -1.873589e+03 163.5331734
## 13 Celkem 252740 1.00000000 252740.00 1.818989e-11 1004.0601088V prvním sloupečku je číslo \(k\)
nebo název skupiny (kategorie). Ve druhém je pozorovaná četnost \(X_k\). Dále následuje hypotetická
pravděpodobnost \(p^0_k\) a očekávaná
četnost (kdyby hypotéza platila) \(np^0_k\). V posledních dvou sloupečcích je
rozdíl \(X_k-np^0_k\) mezi pozorovanou
a očekávanou četností a příspěvek \((X_k-np^0_k)^2/(np^0_k)\) do testové
statistiky. Vidíme, které kategorie (měsíce) nejvíce přispěly k hodnotě
1004.06, která vedla k zamítnutí hypotézy. Byly to duben a květen, kdy
se narodilo mnohem více dětí, než kolik by mělo, kdyby se děti rodily
rovnoměrně, a dále listopad a prosinec, kdy se naopak narodilo dětí
méně.
Je poměr nezamestnaných stejný u mužů a u žen?
Funkciu chisq.test je možné použíť aj k testovaniu
nulovej hypotézy o nezávislosti dvoch premenných – tzv. \(\chi^2\) test nezávislosti. Ak budeme pre
jednoduchosť uvažovať pouze dve binárne premenné – napr. náhodnú
veličinu \(X_1\) (ktorá nabýva hodnot 1
a 2 s pravdepodobnostou \(p_11\) a
\(p_12\)) a druhú náhodnú veličinu
\(X_2\) (ktorá nabýva hodnot 1 a 2 s
pravdepodobnostou \(p_21\) a \(p_22\)), tak nulovú hypotézu nezávislosti
lze formálne (matematicky) zapísať takto:
\[ H_0: p_{i j} = p_{i \bullet} p_{\bullet j} \qquad \textrm{pre všetky $i, j \in \{1,2\}$.} \] Samozrejme platí, že \(p_{11} + p_{12} = 1\) a \(p_{21} + p_{22} = 1\) (v opačnom prípade by sa nejednalo o rozdelenie pravdepodobnosti) a tiež
\[ p_{i \bullet} = p_{i 1} + p_{i 2} \qquad \textrm{a} p_{\bullet j} = p_{1 j} + p_{2 j} \qquad \textrm{opäť pre všetky $i,j \in \{1,2\}$.} \]
Použijeme opäť datový súpbor zam a uvažovať budeme mieru
zamestnanosti/nezamestnanosti medzi mužmi a ženami. Pozorovaná
(empirická) kontingenčná tabuľka (pre nezamestnanosť mužov a žien)
je
tbl <- table(zam$unempl, zam$gndr)
tbl##
## Male Female
## 0 1238 1407
## 1 57 97Pomocou funkcie barplot() možeme na data nahliadnúť
pomocou obrázku
barplot(tbl, cex.names = 0.8, horiz = T, col = c("darkblue","red"), xlim = c(0,1600))
Keďže nás ale zaujíma pomer, lepšie bude pozrieť sa na data prostredníctvom relatívnych čísel:
(tbl_rep<- 100 * cbind(tbl[,1]/sum(tbl[,1]), tbl[,2]/sum(tbl[,2])))## [,1] [,2]
## 0 95.598456 93.550532
## 1 4.401544 6.449468barplot(tbl_rep, cex.names = 0.8, horiz = T, col = c("darkblue","red"), xlim = c(0,100))
Spočítajte relatívne četnosti a pomocou funkcie
chisq.test aplikujte test nezávislosti medzi pohlavím a
mierou nezamestnanosti.
chisq.test(tbl)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tbl
## X-squared = 5.2261, df = 1, p-value = 0.02225Aká je výsledná interpretácia testu?
Je úmrtnost ve věku 15-19 let stejná ve všech krajích?
Nyní použijeme data kraje a zaujímať nás bude miera
úmrtnosti v jednotlivých krajoch ČR.
kraje## kraj obyv naroz posl zemr.15.19 obyv.15.19
## 1 Hlavn\xed m\xecsto Praha 1246780 14233 24 9 47373
## 2 St\xf8edo\xe8esk\xfd kraj 1291816 14483 25 27 61552
## 3 Jiho\xe8esk\xfd kraj 636611 6672 12 22 31887
## 4 Plze\xf2sk\xfd kraj 572687 5785 11 8 26926
## 5 Karlovarsk\xfd kraj 301726 2831 5 13 15177
## 6 \xdasteck\xfd kraj 826764 8246 14 15 42749
## 7 Libereck\xfd kraj 438594 4609 8 7 22035
## 8 Kr\xe1lov\xe9hradeck\xfd kraj 552946 5489 11 12 27779
## 9 Pardubick\xfd kraj 516440 5405 10 8 26552
## 10 Kraj Vyso\xe8ina 511207 5166 11 12 27476
## 11 Jihomoravsk\xfd kraj 1168650 12385 23 17 55909
## 12 Olomouck\xfd kraj 637609 6319 12 9 31777
## 13 Zl\xednsk\xfd kraj 587693 5512 12 10 29521
## 14 Moravskoslezsk\xfd kraj 1226602 11820 22 19 63552Budeme testovat, zdali všechny kraje měly stejné riziko úmrtí mladých
lidí ve věku 15-19 let. Ve sloupečku zemr.15.19 vidíme
počty zemřelých v této věkové kategorii v jednotlivých krajích (v roce
2012). Nesmíme však zapomenout, že kraje jsou různě velké a úmrtnost
musíme posuzovat úměrně k počtu obyvatel v dané věkové kategorii. Tyto
počty jsou uvedeny ve sloupečku obyv.15.19.
Napočítáme nejprve pravděpodobnosti, které by odpovídaly stejné úmrtnosti ve všech krajích.
psti <- kraje$obyv.15.19/sum(kraje$obyv.15.19)
psti## [1] 0.09283999 0.12062752 0.06249106 0.05276866 0.02974337 0.08377804
## [7] 0.04318344 0.05444034 0.05203571 0.05384653 0.10956856 0.06227548
## [13] 0.05785425 0.12454705Nyní provedeme test hypotézy, že rozdělení počtu zemřelých do krajů se řídí právě těmito pravděpodobnostmi a vykucháme výpočet testové statistiky.
chisq.test(kraje$zemr.15.19,p=psti)##
## Chi-squared test for given probabilities
##
## data: kraje$zemr.15.19
## X-squared = 27.376, df = 13, p-value = 0.01105kuchej.chisq(kraje$zemr.15.19,psti,skup=kraje$kraj)## skupina cetnost pst ocek rozdil
## 1 Hlavn\xed m\xecsto Praha 9 0.09283999 17.453919 -8.453919e+00
## 2 St\xf8edo\xe8esk\xfd kraj 27 0.12062752 22.677973 4.322027e+00
## 3 Jiho\xe8esk\xfd kraj 22 0.06249106 11.748319 1.025168e+01
## 4 Plze\xf2sk\xfd kraj 8 0.05276866 9.920508 -1.920508e+00
## 5 Karlovarsk\xfd kraj 13 0.02974337 5.591753 7.408247e+00
## 6 \xdasteck\xfd kraj 15 0.08377804 15.750271 -7.502709e-01
## 7 Libereck\xfd kraj 7 0.04318344 8.118487 -1.118487e+00
## 8 Kr\xe1lov\xe9hradeck\xfd kraj 12 0.05444034 10.234784 1.765216e+00
## 9 Pardubick\xfd kraj 8 0.05203571 9.782713 -1.782713e+00
## 10 Kraj Vyso\xe8ina 12 0.05384653 10.123148 1.876852e+00
## 11 Jihomoravsk\xfd kraj 17 0.10956856 20.598889 -3.598889e+00
## 12 Olomouck\xfd kraj 9 0.06227548 11.707791 -2.707791e+00
## 13 Zl\xednsk\xfd kraj 10 0.05785425 10.876599 -8.765994e-01
## 14 Moravskoslezsk\xfd kraj 19 0.12454705 23.414845 -4.414845e+00
## 15 Celkem 188 1.00000000 188.000000 8.881784e-16
## statistika
## 1 4.09471059
## 2 0.82370304
## 3 8.94570219
## 4 0.37179053
## 5 9.81483197
## 6 0.03573948
## 7 0.15409449
## 8 0.30445078
## 9 0.32486544
## 10 0.34797223
## 11 0.62877181
## 12 0.62626095
## 13 0.07064952
## 14 0.83241457
## 15 27.37595759Testová statistika vyšla 27.38. Porovnává se s kvantilem rozdělení \(\chi^2_{13}\) (14 krajů - 1), které má střední hodnotu 13. Její hodnota stačí na zamítnutí nulové hypotézy (p-hodnota 0.011). Tedy zamítáme hypotézu, že všechny kraje měly stejné riziko úmrtí mladých lidí ve věku 15-19 let.
Které kraje se nejvíce liší? Prakticky celou pozorovanou hodnotu
testové statistiky vytvořily pouze tři kraje: Praha, kde byla úmrtnost
nižší než jinde, a Jihočeský a Karlovarský kraj, kde byla úmrtnost
naopak velmi vysoká. Máte-li v těchto dvou krajích mladší kamarády,
radši jim poraďte, aby se raději přestěhovali buď do Prahy anebo alespoň
do kraje Ústeckého nebo Moravskoslezského.
Samostatne
Otázka, zda je úmrtnosť lidí ve věku 15-19 let stejná ve všech krajích je vpodstatě ekvivalentní s otázkou, zda lze považovat úmrtnost a kraj za nezávisle veličiny. Formálne to vypadá na použití statistického testu pro testování nezávislosti - \(\chi^2\) testu nezávislosti. Je možné tento test aplikovat na data z datového souboru
zam? Pokud ano, udělejte to, pokud ne, vysvětlete
proč.
III. Inferencia pre lineárne kombinácie pravdepodobnosti
Budeme pracovat s daty zam a zkoumat velikosti městské
vs. venkovské populace v ČR v roce 2004. Velikost bydliště respondenta
je uvedena ve veličině zam$domicil. Funkce
table nám spočítá, kolik respondentů se nachází v
existujících kategoriích této veličiny.
tbl <- table(zam$domicil)
tbl##
## A big city Suburbs or outskirts of big city
## 519 138
## Town or small city Country village
## 1360 776
## Farm or home in countryside
## 6Celkový počet respondentů je
sum(tbl)## [1] 2799Vizuálne môžeme data reprezentovať napr. pomocou koláčového grafu (tzv. piechart):
library("colorspace")
pie(table(zam$domicil), col = diverging_hcl(5, c = 100, l = c(50, 90)))
Vhodným modelem pro tuto tabulku je multinomické rozdělení s \(K=5\) kategoriemi, celkovým počtem pozorování \(n=2799\) a neznámým vektorem pravděpodobností \(p\). Ten můžeme odhadnout vektorem relativních četností \(\widehat{p}_n\):
tbl/sum(tbl)##
## A big city Suburbs or outskirts of big city
## 0.185423365 0.049303323
## Town or small city Country village
## 0.485887817 0.277241872
## Farm or home in countryside
## 0.002143623Mohli bychom využít i funkce prop.table, která počítá
totéž:
prop.table(tbl)##
## A big city Suburbs or outskirts of big city
## 0.185423365 0.049303323
## Town or small city Country village
## 0.485887817 0.277241872
## Farm or home in countryside
## 0.002143623Otestujme nejprve, zdali je pravda, že procento lidí žijících na venkově je o 10 větší než procento lidí žijících ve velkoměstech. Venkovská populace odpovídá posledním dvěma kategoriím tabulky. Velkoměstská populace je v první kategorii tabulky. Jde tedy o to, zdali součet posledních dvou pravděpodobností, tj. \(p_4+p_5\), je o 0.1 větší než pravděpodobnost \(p_1\).
Stanovíme vektor \(c=(-1,0,0,1,1)^T\) a budeme testovat
hypotézu \(H_0: c^T p=0.1\). Použijeme
funkci multitest, která není součástí R, ale
implementuje vzorec (7.4) z větníku (str. 68-69). Prvním argumentem této
funkce je realizace vektoru s multinomickým rozdělením, druhý argument
specifikuje vektor \(c\) pro lineární
kombinaci parametrů. Třetí argument uvádí hodnotu lineární kombinace
\(c^T p\) za nulové hypotézy.
multitest(tbl,comb=c(-1,0,0,1,1),c0=0.1)## Warning in c(-1, 1) * d: Recycling array of length 1 in vector-array arithmetic is deprecated.
## Use c() or as.vector() instead.##
## Testovani linearnich kombinací pravdepodobnosti
##
## Multinomicky vektor: 519 138 1360 776 6
##
## Odhadnute pravdepodobnosti: 0.18500 0.04930 0.48600 0.27700 0.00214
## H0: - p1 + p4 + p5 = 0.1
##
## Odhad - p1 + p4 + p5 = 0.09396
##
## Testova statistika = -0.4731, p-hodnota = 0.636
##
## Priblizny 95%-ni interval spolehlivosti:
## ( 0.0689 , 0.119 )Odhad požadované lineární kombinace je 0.09396, což je blízko 0.1. Testová statistika (7.4) je -0.4731. Test má p-hodnotu 0.636, takže nulovou hypotézu nelze zamítnout. Interval spolehlivosti říká, že rozdíl mezi procentem lidí žijících na venkově a procentem lidí žijících ve velkoměstech je s pravděpodobností 0.95 někde mezi 6.9 a 11.9.
Podobně můžeme otestovat hypotézu, že počet lidí žijících v menších městech je dvakrát větší než počet lidí žijících ve velkoměstech a jejich bezprostředním okolí.
multitest(tbl,comb=c(1,1,-1/2,0,0),c0=0)## Warning in c(-1, 1) * d: Recycling array of length 1 in vector-array arithmetic is deprecated.
## Use c() or as.vector() instead.##
## Testovani linearnich kombinací pravdepodobnosti
##
## Multinomicky vektor: 519 138 1360 776 6
##
## Odhadnute pravdepodobnosti: 0.18500 0.04930 0.48600 0.27700 0.00214
## H0: p1 + p2 -0.5 p3 = 0
##
## Odhad p1 + p2 -0.5 p3 = -0.008217
##
## Testova statistika = -0.7285, p-hodnota = 0.466
##
## Priblizny 95%-ni interval spolehlivosti:
## ( -0.0303 , 0.0139 )Sami si promyslete, proč je vektor \(c\) nastaven takto a co znamenají výsledky
testu.
Samostatne
- Pomocou základných popisných charakteristík a pomocou rôznych grafických nástrojov sa podívajte na štruktúru dat a navrhnite vlastnú lineárnu kombináciu pre vektor neznámich pravdepodobnosti v multinomickom rozdelení a test pre nulovosť vami zvolenej lineárnej kombinácie pomocou programu R spočítajte.
- Interpretujte výsledky testu vrámci kontextu celého problému a dát, ktoré máte k dispozícii.