# NMFM301 # Praktické cvičení 3 -- Instruktáž ************** ## Příprava k práci Přesvědčte se, že pracujete se síťovou verzí R (na disku L:) a váš pracovní direktorář je umístěn na disku H: (nikoli C:). Načtěte data pro dnešní cvičení a potřebné funkce: ```{r} load(url("http://www.karlin.mff.cuni.cz/~pesta/NMFM301/cv3.RData")) load(url("http://www.karlin.mff.cuni.cz/~pesta/NMFM301/functions.RData")) ``` Nyní máte k dispozici datové soubory `miry`, `tlak`, `sex`, `deti`. Data si prohlédněte poklepáním na příslušný řádek v okénku "Workspace". Vypište si strukturu dat ```{r} str(miry) str(tlak) str(sex) str(deti) ``` a základní popisné statistiky ```{r} summary(miry) summary(tlak) summary(sex) summary(deti) ``` ## Jednovýběrový Kolmogorovův-Smirnovův test Kolmogorovův-Smirnovův test testuje hypotézu $H_0: F_X(x)=F_0(x) \quad\forall x\in\mathbb{R}$ proti alternativě $H_1: \exists x\in\mathbb{R}: F_X(x)\neq F_0(x)$, kde $F_0$ je nějaká pevně specifikovaná spojitá distribuční funkce (bez neznámých parametrů). Jeho testová statistika je \[ K_n=\sup_{x\in\mathbb{R}}\left|\widehat{F}_{n}(x)-F_0(x)\right|, \] kde $\widehat{F}_{n}$ je empirická distribuční funkce náhodného výběru $X_1,\ldots,X_n$. Hypotézu zamítáme pro velké hodnoty $K_n$ nebo $\sqrt{n}K_n$. Používáme asymptotické kritické hodnoty vypočtené numericky. Otestujeme hypotézu, že výška obyvatel USA zaznamenaná v datovém souboru `miry` se řídí rozdělením N(168,100). Nejprve spočítáme a namalujeme empirickou distribuční funkci a hypotetickou distribuční funkci. Empirická distribuční funkce se počítá funkcí `ecdf`. ```{r} plot(ecdf(miry$vyska),cex.points=0.5, main="Empiricka distribucni funkce vysky") xpts=seq(135,200,length=500) lines(xpts,pnorm(xpts,mean=168,sd=10)) ``` Funkce `lines` domaluje na předchozí obrázek propojené hodnoty hypotetické distribuční funkce spočítané v bodech `xpts`. Provedeme Kolmogorovův-Smirnovův test pomocí funkce `ks.test`. Jako první argument má testovaný náhodný výběr, další tři argumenty specifikují hypotetickou distribuční funkci, poslední argument říká, že máme počítat asymptotický test (nikoli přesný). ```{r} ks.test(miry$vyska,"pnorm",mean=168,sd=10,exact=FALSE) ``` R nás varuje, že v datech jsou shodné hodnoty (ties), čímž porušujeme předpoklady. Výsledek testu tedy nemusí být spolehlivý. Nyní budeme tento problém ignorovat. Hodnota testové statistiky $K_n$ byla 0.037, p-hodnota vyšla 0.4993. Na hladině $\alpha=0.05$, nemůžeme zamítnout hypotézu, že výška má rozdělení N(168,100). Nyní opakujte celý postup pro výšku mužů (`miry$vyska[miry$pohl=="Male"]`) a otestujte, zdali má výška mužů normální rozdělení se střední hodnotou 174 a směrodatnou odchylkou 7.5. ## 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: ```{r} sign.test(miry$vaha,78) ``` 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*. ```{r} binom.test(sum(miry$vaha>78),length(miry$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č asi? ## Párové testy Nyní budeme pracovat s datovým souborem `tlak`, který obsahuje 2309 pozorování. Nalezneme v něm dvě měření systolického tlaku v mm Hg (`sys.tl.1` a `sys.tl.2`) vykonaná na téže osobě v určitém časovém intervalu (několik minut) a dvě měření diastolického tlaku (`dia.tl.1` a `dia.tl.2`). Dále máme informace o věku (veličina `vek`, od 50 let výš) a pohlaví (veličina `pohl`) pacientů. Nejprve spočítejme průměrný systolický tlak při prvním a druhém měření ```{r} mean(tlak$sys.tl.1) mean(tlak$sys.tl.2) ``` a namalujme krabicový diagram obou měření. Takto se malují krabicové diagramy dvou veličin (nebo dvou skupin) vedle sebe: ```{r} boxplot(c(tlak$sys.tl.1,tlak$sys.tl.2)~rep(1:2,nrow(tlak)),main="Systolicky tlak", names=c("1. mereni","2. mereni")) ``` Prohlédněme si ještě histogram rozdílů mezi prvním a druhým měřením (s 24 intervaly o délce 2 mm Hg): ```{r} hist(tlak$sys.tl.1-tlak$sys.tl.2,breaks=24) ``` Řekli byste podle obrázků, že mezi prvním a druhým měřením tlaku je nějaký systematický rozdíl? Otestujme si, zdali se systematicky liší střední hodnota prvního a druhého měření systolického tlaku. Problém je třeba řešit párovým testem, protože se jedná o náhodný výběr dvojic pozorování $(X_i,Y_i)$ učiněných na tom samém pacientovi $i\in\{1,\ldots,n\}$. Párové testy jsou vlastně jednovýběrové testy aplikované na rozdíly $Z_i=X_i-Y_i$. Ukážeme si na tomto příkladě několik různých párových testů. ### Párový t-test Párový t-test spočítáme příkazem ```{r} t.test(tlak$sys.tl.1,tlak$sys.tl.2,paired=TRUE) ``` T-statistika nabyla hodnoty více než 14, proto s velikou rezervou zamítáme hypotézu o rovnosti středních hodnot. Můžeme tedy tvrdit, že střední hodnota systolického tlaku mezi prvním a druhým měřením klesla. ### Párový znaménkový test Párový znaménkový test testuje rovnost středních hodnot, jestliže rozdělení rozdílů mezi prvním a druhým měřením tlaku je symetrické. Pomocí histogramu rozdílů (viz výše) můžeme neformálně posoudit, zdali je tento předpoklad splněn. Myslíte, že ano? Párový znaménkový test provedeme jako jednovýběrový znaménkový test aplikovaný na rozdíly mezi prvním a druhým měřením tlaku. ```{r} sign.test(tlak$sys.tl.1-tlak$sys.tl.2,m0=0) ``` Vidíme, že i znaménkový test zamítá hypotézu o rovnosti středních hodnot. Počet lidí, jimž se mezi oběma měřeními tlak snížil, byl příliš velký. Poznáte z výstupu funkce `sign.test`, kolik jich bylo? ### Párový Wilcoxonův test I párový Wilcoxonův test testuje rovnost středních hodnot za předpokladu symetrie rozdělení rozdílů mezi prvním a druhým měřením tlaku. Počítá se funkcí `wilcox.test`, argumenty vypadají podobně jako u párového t-testu. ```{r} wilcox.test(tlak$sys.tl.1,tlak$sys.tl.2,paired=TRUE,correct=FALSE) ``` Wilcoxonův test také zamítá hypotézu o rovnosti středních hodnot. ## Dvouvýběrové testy Nyní budeme pracovat s datovým souborem `sex`. Tento soubor obsahuje informace o sexuálním životě 4467 Američanů. 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í ```{r} boxplot(veksex~pohl,data=sex) ``` 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`). ```{r} tapply(sex$veksex,sex$pohl,mean,na.rm=TRUE) tapply(sex$veksex,sex$pohl,median,na.rm=TRUE) tapply(sex$veksex,sex$pohl,var,na.rm=TRUE) ``` Ještě můžeme namalovat průměry věku iniciace a intervaly spolehlivosti pro střední věk iniciace pro obě pohlaví. ```{r} plotmeans(veksex~pohl,data=sex,xlab="Pohlavi",ylab="Prumerny vek iniciace") ``` Nakonec porovnáme empirické distribuční funkce věku iniciace obou pohlaví. ```{r} 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? Nyní přikročíme k testování několika různými metodami. ### 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`. ```{r} 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. ```{r} 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. ### 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í. ```{r} 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.