NMFM301, ZS 2017/18

Cvičenie 9 (praktické cvičenie 3)

(jednovýberové a dvojvýberové štatistické testy)



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)




Data

V prípade, že používate počítač, ktorý je k dispozícii v posluchárni, tak sa pre lepšie fungovanie programu R najprv presvedčte,
že pracujete so sieťovou verzíou programu R (na disku L:) a váš pracovní direktorář je umístěn na disku H: (nikoli C:).

Budeme potrebovať niekoľko rôznych datových súborov, ktoré spoločne načítate do Rka pomocou príkazu (počítač musí byť pripojený na internet)

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

Následne načítame analogickym príkazom aj dodatočné funkcie, ktoré nám umožnia priamo počítať napr. intervaly spoľahlivosti, alebo robiť niektoré štatistické testy (ktoré v štandardnej inštalácii Rka nie sú k dispozícii)

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

Nyní máte k dispozici datové soubory miry, tlak, sex, deti.

Pomoci příkazů dim(), str() a summary() si prohlédněte rozsah, struktúru a popisné charakteristiky v datovych souborech.

  • Rozsah náhodnych vyběrů a počet sledovaných proměnných:

    str(miry)
    str(tlak)
    str(sex)
    str(deti)
  • Strukturu jednotlivych souborů:

    str(miry)
    str(tlak)
    str(sex)
    str(deti)
  • Základné popisné charakteristiky:

    summary(miry)
    summary(tlak)
    summary(sex)
    summary(deti)



I. Jednovýběrové testy

Jednovýběrový Kolmogorovův-Smirnovův test

Uvažujme náhodný výber \(X_{1}, \dots, X_{n}\) z nejakého spojitého rozdelenia, ktoré je dané distribučnou funkciou \(F_{X}\). 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. Jedná se o neparametrický statistický test, kterého testová statistika je definová

\[ 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 vykreslíme empirickou distribuční funkci a hypotetickou distribuční funkci. Empirická distribuční funkce se počítá funkcí ecdf.

plot(ecdf(miry$vyska),cex.points=0.5, main="")
xpts=seq(135,200,length=500)
lines(xpts,pnorm(xpts,mean=168,sd=10), col = "red")
legend(138, 1, legend = c("Empirická distribuční funkce", "Teoretická distribučni funkce"), col = c("black", "red"), 
       lty = c(1,1), lwd = c(2,2))

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(). Pro help jak zadávat a používat argumenty u této funkce, použite ?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ý).

ks.test(miry$vyska,"pnorm",mean=168,sd=10,exact=FALSE)
## Warning in ks.test(miry$vyska, "pnorm", mean = 168, sd = 10, exact =
## FALSE): ties should not be present for the Kolmogorov-Smirnov test
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  miry$vyska
## D = 0.03703, p-value = 0.4993
## alternative hypothesis: two-sided

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.

Užitočné


Kolmogorovův-Smirnovův test je definovaný pre obecnou nulovú a alternatívnu hypotézu (bez předpokladu normality) a proto je v prípade tejto takto formulovaného testu pomerně slabý v porovnaní so niektorými inými štatistickými testami, ktoré sú špecificky navrhnuté pre testovanie nulovej hypotézy, že náhodný výber pochádza z normálneho rozdelenia.

Samostatne


Pomocou histogramu a teoretickej hustoty normálneho rozdelenia \(N(168, 100)\) porovnajte empirické výšky s rozdelením, ktoré ste v nulovej hypotéze vyššie nezamietli.



Znaménkový test

Uvažujme opäť náhodný výber \(X_{1}, \dots, X_{n}\) z nejakého spojitého rozdelenia, ktoré je dané distribučnou funkciou \(F_{X}\). 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 teoretický médian rozdelenia \(F_{X}\) a \(m_0\) je nejaká predom zvolená konstanta. 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(miry$vaha,78)
##                n              Y_n Test. statistika    Krit. hodnota 
##      500.0000000      267.0000000        1.5205262        1.9599640 
##        P-hodnota 
##        0.1283788

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(miry$vaha>78),length(miry$vaha))
## 
##  Exact binomial test
## 
## data:  sum(miry$vaha > 78) and length(miry$vaha)
## number of successes = 267, number of trials = 500, p-value =
## 0.1399
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.4891844 0.5784114
## sample estimates:
## probability of success 
##                  0.534

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?

Otázky a úlohy


  • Pomocou jednovýberového \(t\)-testu otestujte, že stredná hodnota váhy je 78kg a stredn hodnota výšky je 168cm.
  • Porovnajte výsledky predchádzajúcich dvoch testov s výsledkami \(t\)-testov a vysvetlite.
  • Aké sú predpoklady jednovýberového \(t\)-testu? Sú tieto predpoklady splnené?



II. Dvouvýběrové 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ů.

Párový t-test

Nejprve spočítejme průměrný systolický tlak při prvním a druhém měření

mean(tlak$sys.tl.1)
## [1] 132.1369
mean(tlak$sys.tl.2)
## [1] 130.0017

a vykreslime krabicový diagram obou měření. Takto se malují krabicové diagramy dvou veličin (nebo dvou skupin) vedle sebe:

boxplot(c(tlak$sys.tl.1,tlak$sys.tl.2)~rep(1:2,nrow(tlak)),main="Systolický tlak",
        names=c("1. měření","2. měření"), col = "lightblue")

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

hist(tlak$sys.tl.1-tlak$sys.tl.2,breaks=24, col = "lightblue", main = "Rozdíl systolických tlaků")

Ř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 spočítáme v programu R pomoci příkazu (podívejte se na help ?t.test pro pochopení správne implementace).

t.test(tlak$sys.tl.1,tlak$sys.tl.2,paired=TRUE)
## 
##  Paired t-test
## 
## data:  tlak$sys.tl.1 and tlak$sys.tl.2
## t = 14.269, df = 2308, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  1.841697 2.428550
## sample estimates:
## mean of the differences 
##                2.135123

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.

sign.test(tlak$sys.tl.1-tlak$sys.tl.2,m0=0)
##                n              Y_n Test. statistika    Krit. hodnota 
##     2.309000e+03     1.272000e+03     4.890530e+00     1.959964e+00 
##        P-hodnota 
##     1.005650e-06

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.

wilcox.test(tlak$sys.tl.1,tlak$sys.tl.2,paired=TRUE,correct=FALSE)
## 
##  Wilcoxon signed rank test
## 
## data:  tlak$sys.tl.1 and tlak$sys.tl.2
## V = 1371100, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0

Wilcoxonův test také zamítá hypotézu o rovnosti středních hodnot.

Otázky a úlohy


  • Uveďte předpoklady k párovému \(t\)-testu, párovému znaménkovému testu a k párovému Wilcoxonovu testu.
  • Který z uvedených testů je podle vás nejlepší pro zrovnání dvou systolických tlaků?



III. Dvouvýběrové testy \(|\) nepárové

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í

boxplot(veksex~pohl,data=sex, pch = 21, bg = "pink", col = "lightblue")

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

tapply(sex$veksex,sex$pohl,mean,na.rm=TRUE)
##     Male   Female 
## 17.35301 17.90112
tapply(sex$veksex,sex$pohl,median,na.rm=TRUE)
##   Male Female 
##     17     17
tapply(sex$veksex,sex$pohl,var,na.rm=TRUE)
##     Male   Female 
## 18.23625 15.54390

Ještě můžeme namalovat průměry věku iniciace a intervaly spolehlivosti pro střední věk iniciace pro obě pohlaví.

plotmeans(veksex~pohl,data=sex,xlab="Pohlavi",ylab="Prumerny vek iniciace")

Nakonec porovnáme empirické distribuční funkce věku iniciace obou pohlaví.

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