Mějme náhodný výběr \(X_1,\ldots,X_n\) z rozdělení \(F_X\). Kvantilový diagram je grafická metoda sloužící k neformálnímu porovnání rozdělení \(F_X\) s hypotetickým rozdělením \(F_0\).
Na ose \(y\) se kreslí pořádkové statistiky \(X_{(i)}\). Na ose \(x\) se vynášejí aproximace střední hodnoty pořádkových statistik rozdělení \(F_0\), tj. \(F_0^{-1}\bigl(\frac{i}{n+1}\bigr)\) [proč zrovna tohle?].
Jestliže data pocházejí z rozdělení \(F_0\), vykreslené body leží přibližně na přímce \(y=x\). Jestliže data pocházejí z lineární transformace rozdělení \(F_0\), vykreslené body stále leží přibližně na přímce (ale jiné než \(y=x\)). Jestliže data pocházejí z rozdělení odlišného od \(F_0\), vykreslené body tvoří rostoucí křivku. Z tvaru této křivky dokážeme poznat, v čem se rozdělení dat liší od \(F_0\).
V R se dá porovnávat rozdělení dat s rozdělením N(0,1) pomocí funkce qqnorm
. Argumentem funkce je vektor obsahující náhodný výběr.
Vygenerujme 50 pozorování z rozdělení N(10,4) a namalujme normální kvantilový diagram:
x <- rnorm(50,10,2)
qqnorm(x, pch = 21, bg = "red")
qqline(x)
Funkce qqline
namaluje přímku procházející dolním a horním výběrovým kvartilem dat, která usnadnuje čtení grafu.
Nyní vygenerujme 50 pozorování z rozdělení Exp(0.4) a namalujme opět normální kvantilový diagram:
x <- rexp(50,0.4)
qqnorm(x, pch = 21, bg = "red")
qqline(x)
Poznali byste, že tato data nepocházejí z normálního rozdělení?
Vytvořte normální kvantilový diagram pro data z rozdělení \(\Gamma(a,p)\) a sledujte, co se děje, když se zvětšuje argument \(p\).
Vygenerujte 50 pozorování z rozdělení \(\Gamma(a,p)\), kde a=0.4 a p=0.5.
Musíte dát pozor na správnou parametrizaci Gama rozdělení; použijte volání rgamma(50,rate=0.4,shape=0.5)
(argument rate
udává hodnotu parametru a, argument shape
udává hodnotu parametru p, oba musíte pojmenovat).
Nakreslete normální kvantilový diagram pro tato data.
Zvětšujte hodnotu parametru p na 1, 2, 20, 50 a sledujte, co se děje s kvantilovým diagramem. Jaké to má vysvětlení?
Nyní se naučíme pracovat s datovými soubory. Načtěte data z pracovního direktoráře nebo přímo z webu příkazem
load(url("http://msekce.karlin.mff.cuni.cz/~maciak/NMFM301/cv2.RData"))
?load
Nyní máte k dispozici datový soubor nazvaný cv2
.
Načteme ješte několik funkcí, s kterými dnes budeme pracovat (plotCI
,plotmeans
,ci.asym
,ci.t
,sign.test
).
load(url("http://msekce.karlin.mff.cuni.cz/~maciak/NMFM301/functions.RData"))
Data mají 500 pozorování (řádků) a 6 veličin (sloupců). Jejich velikost zjistíme funkcí dim
:
dim(cv2)
## [1] 500 6
Pozorování tvoří náhodný výběr 500 obyvatel USA. V jednotlivých veličinách o nich máme následující údaje:
Název | Význam |
---|---|
id |
ID jedince |
pohl |
Pohlaví |
vzdelani |
Vzdělání |
stav |
Rodinný stav |
vaha |
Hmotnost [kg] |
vyska |
Výška [cm] |
Data si prohlédněte poklepáním na příslušný řádek v okénku “Workspace”. Strukturu dat vypíšeme pomocí funkce str
:
str(cv2)
## 'data.frame': 500 obs. of 6 variables:
## $ id : num 68573 71157 65182 69993 64705 ...
## $ pohl : Factor w/ 2 levels "Male","Female": 1 2 1 2 2 2 1 1 1 1 ...
## $ vzdelani: Factor w/ 5 levels "Elementary unfinished",..: 2 4 2 2 4 1 3 4 2 4 ...
## $ stav : Factor w/ 6 levels "Married","Widowed",..: 1 1 6 5 1 2 3 6 1 1 ...
## $ vaha : num 57.1 94.5 71.6 118.4 76.7 ...
## $ vyska : num 165 163 170 157 149 ...
## - attr(*, "na.action")=Class 'omit' Named int [1:385] 13 32 48 75 83 93 103 105 119 145 ...
## .. ..- attr(*, "names")= chr [1:385] "30" "62" "87" "130" ...
Základní popisné statistiky všech veličin získáme pomocí funkce summary
:
summary(cv2)
## id pohl vzdelani
## Min. :62172 Male :259 Elementary unfinished: 40
## 1st Qu.:64728 Female:241 Elementary : 82
## Median :67260 High school :100
## Mean :67133 College unfinished :155
## 3rd Qu.:69597 College :123
## Max. :71897
## stav vaha vyska
## Married :268 Min. : 35.30 Min. :139.3
## Widowed : 34 1st Qu.: 65.67 1st Qu.:160.4
## Divorced : 45 Median : 78.60 Median :167.8
## Separated : 18 Mean : 81.81 Mean :168.1
## Single : 97 3rd Qu.: 94.08 3rd Qu.:175.3
## Living with partner: 38 Max. :188.50 Max. :204.5
S datovým souborem se zachází podobně jako s maticí. Můžeme vypsat jeden nebo více řádků pomocí hranaté závorky:
cv2[1:5,]
## id pohl vzdelani stav vaha vyska
## 6126 68573 Male Elementary Married 57.1 165.0
## 8610 71157 Female College unfinished Married 94.5 162.7
## 2886 65182 Male Elementary Living with partner 71.6 170.3
## 7500 69993 Female Elementary Single 118.4 156.6
## 2430 64705 Female College unfinished Married 76.7 148.9
K jednotlivým veličinám máme přístup buď pomocí čísla sloupce (cv2[,2]
je vektor hodnot pohlaví) nebo lépe pomocí názvu veličiny napsané za názvem dat a oddělené znakem \(, např `cv2\)pohl`. Výšky prvních pěti obyvatel vypíšeme takto:
cv2$vyska[1:5]
## [1] 165.0 162.7 170.3 156.6 148.9
Hmotnosti lidí se základním vzděláním získáme pomocí podmínky v hranaté závorce (rovnítko musí být zdvojeno):
cv2$vaha[cv2$vzdelani=="Elementary"]
## [1] 57.1 71.6 118.4 77.6 82.8 62.1 35.3 73.2 130.0 67.0 61.0
## [12] 83.0 78.7 109.0 104.1 66.6 87.1 84.5 137.4 91.2 62.1 63.0
## [23] 124.7 87.2 84.4 86.6 82.4 70.0 76.7 83.7 113.8 90.4 62.1
## [34] 90.6 84.0 64.6 81.0 167.8 71.3 129.2 73.4 83.1 69.8 51.9
## [45] 77.2 88.3 109.3 83.0 78.2 87.9 75.0 94.5 78.1 140.1 101.4
## [56] 51.8 87.4 81.1 70.0 109.7 69.1 73.3 64.5 57.7 110.2 76.1
## [67] 60.2 143.5 62.2 90.2 59.6 88.6 68.2 70.3 104.5 54.5 88.1
## [78] 86.4 65.6 110.1 75.9 64.5
Zkusme se nyní podívat na normální kvantilový diagram výšky.
qqnorm(cv2$vyska, pch = 21, bg = "red")
qqline(cv2$vyska)
Obecně se věří, že výška v populaci má normální rozdělení. Nasvědčuje tomu tento obrázek?
Podívejme se ještě na výšky mužů a žen zvlášť.
qqnorm(cv2$vyska[cv2$pohl=="Male"], pch = 21, bg = "red")
qqline(cv2$vyska[cv2$pohl=="Male"])
qqnorm(cv2$vyska[cv2$pohl=="Female"], pch = 21, bg = "red")
qqline(cv2$vyska[cv2$pohl=="Female"])
Jsou výšky pro každé pohlaví normální?
Funkce ci.t
počítá 95% interval spolehlivosti pro střední hodnotu založený na t-rozdělení: \[
\left(\overline{X}_{n}-\frac{S_n}{\sqrt{n}}t_{n-1}\Bigl(1-\frac{\alpha}{2}\Bigr),
\ \overline{X}_{n}+\frac{S_n}{\sqrt{n}}t_{n-1}\Bigl(1-\frac{\alpha}{2}\Bigr)\right).
\] (pro \(\alpha=0.05\)). Je to přesný interval pro normální rozdělení a přibližný pro jakékoli rozdělení s konečným druhým momentem.
smp <- rnorm(20,2,1)
prum <- mean(smp)
ci.t(smp)
## CI.low Mean CI.up
## 1.372024 1.942404 2.512784
Výsledkem funkce ci.t
je vektor o třech složkách: dolní konec intervalu, průměr, horní konec intervalu.
Vygenerujme nyní 100 náhodných výběrů o rozsahu 20 se střední hodnotou 2. Nejprve nastavíme požadované vstupy:
nobs <- 20
nvyb <- 100
str.h <- 2
Nyní vygenerujeme 2000 veličin z rozdělení N(2,1) a uspořádáme je do matice 20 x 100. Ve sloupcích bude 100 výběrů o rozsahu 20 z rozdělení N(2,1).
vybery <- rnorm(nobs*nvyb,mean=str.h,sd=1)
data.mat <- matrix(vybery,nrow=nobs,ncol=nvyb)
Na každý sloupec (výběr) spočítáme interval spolehlivosti pro střední hodnotu funkcí ci.t
:
vs.ci <- apply(data.mat,2,ci.t)
Zjistíme, kolik intervalů pokrylo skutečnou střední hodnotu (2) a jaká byla průměrná šířka intervalu. Kolik intervalů by ji teoreticky mělo pokrýt?
(pokryti <- sum(vs.ci[1,]<str.h & str.h<vs.ci[3,])/nvyb)
## [1] 0.94
(prum.sirka <- mean(vs.ci[3,]-vs.ci[1,]))
## [1] 0.9016771
Nakonec nakreslíme všech 100 spočítaných intervalů do obrázku. Červená čára značí skutečnou střední hodnotu. Vidíte, které intervaly ji nepokryly?
plotCI(vs.ci[2,],uiw=vs.ci[3,]-vs.ci[2,],liw=vs.ci[2,]-vs.ci[1,],
gap=0.15,sfrac=0.002,ylab="Int. spol. pro stredni hodnotu",xlab="Vyber")
abline(h=str.h,col="red")
Nyní zvětšujte rozsah výběru. Generujte 100 výběrů z N(2,1) o rozsahu 100 a 1000. Jak se mění šířka intervalů? Jak se mění pokrytí?
Sestrojíme interval spolehlivosti pro střední hmotnost mužů a žen:
ci.t(cv2$vaha[cv2$pohl=="Male"])
## CI.low Mean CI.up
## 84.11979 86.67181 89.22384
ci.t(cv2$vaha[cv2$pohl=="Female"])
## CI.low Mean CI.up
## 73.77435 76.59378 79.41320
Co tyto intervaly znamenají? Umíte je interpretovat?
Chceme-li získat obrázky těchto intervalů, použijeme funkci plotmeans
:
plotmeans(vaha~pohl,data=cv2,connect=FALSE, ylim = c(70,90))
Dá se říci, že muži a ženy váží stejně?
Zkuste namalovat intervaly spolehlivosti pro hmotnost pro různé stupně vzdělání.
Znaménkový test testuje hypotézu \(H_0: m_X=m_0\) proti alternativě \(H_1: m_X\neq m_0\). Je založen na statistice \[ Y_n=\sum_{i=1}^n \mathbb{I}_{(0,\infty)}(X_i-m_0). \] Můžeme jej provést přesně (\(Y_n\) má za hypotézy binomické rozdělení) nebo asymptoticky pomocí statistiky \[ \frac{2}{\sqrt{n}}Y_n-\sqrt{n}, \] která má za \(H_0\) přibližně normované normální rozdělení.
Asymptotický test je implementován funkcí sign.test
. Jejím prvním argumentem je náhodný výběr, druhým argumentem je hypotetický medián \(m_0\). Např. test hypotézy, že medián hmotnosti je 78 kg:
sign.test(cv2$vaha,78)
## n Y_n Test. statistika Krit. hodnota
## 500.0000000 254.0000000 0.3577709 1.9599640
## P-hodnota
## 0.7205148
Funkce sign.test
počítá celkový počet pozorování n, hodnotu \(Y_n\), hodnotu asymptotické testové statistiky, kritickou hodnotu pro testování a p-hodnotu.
Otestujte nyní asymptotickým znaménkovým testem hypotézy, že medián hmotnosti mužů (žen) je 78 kg.
Přesný asymptotický test lze získat pomocí funkce binom.test
(jde vlastně o test na parametr binomického rozdělení). Do této funkce musíme zadat hodnotu \(Y_n\) a hodnotu n.
binom.test(sum(cv2$vaha>78),length(cv2$vaha))
##
## Exact binomial test
##
## data: sum(cv2$vaha > 78) and length(cv2$vaha)
## number of successes = 254, number of trials = 500, p-value =
## 0.7543
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
## 0.4632436 0.5526613
## sample estimates:
## probability of success
## 0.508
Otestujte nyní přesným znaménkovým testem hypotézy, že medián hmotnosti mužů (žen) je 78 kg. Liší se nějak záSadně výsledky přesného testu od testu asymptotického? Proč asi?