NMFM301, ZS 2015/16

Cvičenie 8 (praktické cvičení 2 – Instruktáž)


Kvantilový diagram

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í?

Úkol:

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

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

  2. Nakreslete normální kvantilový diagram pro tato data.

  3. 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í?

Práce s daty

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í?

Intervaly spolehlivosti pro střední hodnotu

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

Úkol:

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í?

Aplikace na datech

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

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?