NMFM301, ZS 2019/20

Cvičenie 8 (praktické cvičenie 2)

(popisné charakteristiky a štatistické testy)

Zdrojový soubor pro Knit: Rmd soubor

Cieľom druhého praktického cvičenia je práca s reálnymi datovými súbormi, jednoduchá analýza dat v programe R a aplikácia základných štatistických metód (napr. počítanie jednoduchých výberových charakteristík, konštrukcia intervalov spoľahlivosti pre neznámy parameter v modeli, alebo test hypotézy o neznámom parametri) pomocou štatistického softwaru R.

V programe R je k dispozícii niekolko datových súborov, ktorých zoznam získame pomocou príkazu data(). Ďalšie datové súbory je možné získať nainštalovaním dodatočných R knižníc, tzv. libraries. Okrem týchto preddefinovaných (vzorových) datových súborov je možné do programu R načítať vlastné súbory (napr. xls súbor, alebo csv a txt tabulky, alebo aj načítať data priamo z databázy SQL, alebo z internetového zdroja).

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)


I. Kvantilový diagram

Pre ľubovoľný náhodný výběr \(X_1,\ldots,X_n\) z nejakého rozdělení \(F_X\) používame kvantilový diagram ako neformálnu grafickú metódu slúžiacu k porovnaniu neznámeho rozdělení \(F_X\) s hypotetickým rozdělením \(F_0\) (napr. v prípade potreby overiť predpoklad normality je \(F_0\) štandardne Gaussovo rozdelenie \(N(0,1)\)).

Na os \(y\) sa vykreslia pořádkové statistiky \(X_{(i)}\), pro \(i = 1, \dots, n\). Na os \(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 hodnota argumentu/parametru \(p\).
  • Vygenerujte 50 pozorování z rozdělení \(\Gamma(a,p)\), kde a=0.4 a p=0.5. Použijte help k funkci rgamma() pro správnou parametrizaci Gama rozdělení.
  • Nakreslete normální kvantilový diagram pro takto vygenerované 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í?



II. Práce s daty

V nasledujúcom budeme pracovať so skutočnými datami, ktoré pochádzajú z reálneho experimentu. Datovy súbor je možne stiahnuť na adrese http://msekce.karlin.mff.cuni.cz/~maciak/NMFM301/cv2.RData, prípadne je možné použiť link na data a data priamo načítať do programu R pomocou príkazu.

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

Momentálne máme v programe R k dispozíci datový soubor s názvom cv2. Stručny náhľad na štruktúru datového súboru získame pomocou príkazu

head(cv2)
##         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
## 3591 65919 Female Elementary unfinished             Widowed  49.4 145.8

Data mají celkovo 500 pozorování (řádků) a 6 veličin (sloupců). Túto informáciu zjistíme pomocou funkcie 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")= 'omit' Named int  13 32 48 75 83 93 103 105 119 145 ...
##   ..- attr(*, "names")= chr  "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



Užitočné

Několik užitočných funkcii, ktoré uľahčia prácu v programe R (napr. funkcie plotCI,plotmeans,ci.asym,ci.t,sign.test).

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



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

Použijte iné grafické nástroje v programe R a pokúste za podívat na data alternatívnym spôsobom: napr. pomocou histogramu (príkaz hist()), alebo pomocou odhadu hustoty (príkaz density()), alebo empirickej distribučnej funkcie (príkaz ecdf()). Jednotlivé empirické charakteristiky vykreslite na vhodnom obrázku a porovnajte s príslušným protejškem založeným na predpoklade normality.

III. 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.604417 2.007028 2.409640

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.93
(prum.sirka <- mean(vs.ci[3,]-vs.ci[1,]))
## [1] 0.9073477

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í?
  • Vygenerujte 100 náhodných výběrů o rozsahu 10 z exponenciálního rozdělení Exp(0.25). Zjistěte průměrnou šířku intervalu. Kolik jich pokrylo skutečnou střední hodnotu? Jak se to liší od onormálního rozdělení? [Pokuste se o vysvětlení] Nakreslete obrázek zachycující všech 100 intervalů a skutečnou střední hodnotu.



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

Úkol:

  • Sestrojte intervaly spolehlivosti pro střední výšku ženatých/vdaných lidí a lidí žijících odděleně. Který z nich je delší?
    Cím je způsobena rozdílnost jejich délky?



IV. Štatistické testy

V nasledujúcej časti sa budeme podrobnejšie venovať niekoľkým štatistických testom, ktorými budeme testovať skutočnu hodnotu parametru (resp. parametrov).

Znaménkový test

Znaménkový test testuje hypotézu \(H_0: m_X=m_0\) proti alternativě \(H_1: m_X\neq m_0\), kde \(m_X\) je výberový médian z nejakého spojitého rozdelenia. Test je založen na testovej 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))

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

\(t\)-test

V prípade, že náhodný výber pochádza z normálneho rozdelenia (rozdelenie je symetrické), teoretická hodnota mediánu je totožná s teoretickou strednou hodnotou. V prípade, že nas zaujíma test ohľadom skutočnej hodnoty mediánu, môžeme za týchto predpokladov použíť aj jednovýberový \(t\)-test (hodnota parametru rozptylu je neznáma) - resp. testovať nulovú hypotézu o skutočnej hodnote parametru strednej hodnoty. Implementácia testu je pomocou príkazu t.test().

t.test(cv2$vaha, mu = 78)
## 
##  One Sample t-test
## 
## data:  cv2$vaha
## t = 3.8617, df = 499, p-value = 0.0001274
## alternative hypothesis: true mean is not equal to 78
## 95 percent confidence interval:
##  79.87366 83.75474
## sample estimates:
## mean of x 
##   81.8142

Ako vysvetliť rozpor medzi výsledkom znamienkového testu a klasického \(t\)-testu?

hist(cv2$vaha, col = "lightblue", xlab = "Hmotnosť/váha", breaks = 20, main = "")



Úkol:

  • Otestujte asymptotickým znaménkovým testem hypotézu, že medián výšky je 167 cm (celkově, mezi muži a mezi ženami). Interpretujte výsledek testů.
  • Otestujte nulovú hypotézu pre skutočnú hodnotu mediánu a tiež pre skutočnú hodnotu strednej hodnoty za použitia znamiekového testu a \(t\)-testu, pričom bude uvažovať samostatne populáciu mužov a žien. Pre nulovú hypotézu zvoľte vhodnú hodnotu, ktorú budete testovať.
  • Niekolko dalších príkladov k samostatnému precvičovaniu: zadanie a zdrojový Rmd súbor.



Užitočné

V programe R je defaultne implementoványch niekoľko funkcii so základnými štatistickými testami, napr. test parametru binomického rozdelenia binom.test() (znamienkový test), alebo jedno/dvoj výberové testy wilcox.test(), t.test(), ks.test(), fisher.test(), var.test(), prípadne viacvýberové testy kruskal.test(), anova() a mnoho ďalších.