NMFM301, ZS 2023/24

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

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

Zdrojový soubor pro Knit: Rmd soubor

Užitočné materiály pre prácu so štatistickým programom 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)


Cieľom tohto cvičenia je aplikácia (prostredníctvom programu R) niektorých jednovýberových a dvojvýberových štatistických testov na reálne data. Pre tieto účely využijeme postupne niekoľko rôznych datových súborov.

Data

Datové súbory spoločne načítate do programu R pomocou príkazu (počítač musí byť pripojený na internet)

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


Analogickým príkazom môžeme načítať aj niekoľko pomocných funkcii, ktoré umožňujú priamo spočítať niektoré užitočné charakteristiky – napr. intervaly spoľahlivosti, prípadne niektoré základné štatistické testy. Konkrétne sa jedná o funkcie ci.asym(), ci.t(), invalid(), kuchej.chisq(), multitest(), plotCI(), plotmeans(), print.multitest(), print.propor(), propor() a sign.test(). 

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


K dispozici sú štyri rôzne datové soubory, pomenované ako miry, tlak, sex a deti.

Pomoci příkazů dim(), str(), head() a summary() si zistíme rozsah, struktúru a formát dat a tiež niektoré základné popisné charakteristiky (prvého rádu).

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

    dim(miry)
dim(tlak)
dim(sex)
dim(deti)

  • Štruktúra jednotlivych premenných v daných súboroch:

    str(miry)
str(tlak)
str(sex)
str(deti)

  • Niekoľko prvých pozorovaní z každého datového súboru:

    head(miry)
head(tlak)
head(sex)
head(deti)

  • Základné popisné charakteristiky (bez charakteristík druhého rádu – t.j., napr. rozptyl, smerodajná chyba, atď.):

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

  • Odhadnutá variabilita (resp. variančná-kovariančná matica) dat (popisné charakteristiky druhého rádu):

    var(miry[,c("vaha", "vyska")])
var(tlak[,-c(1,3)])
var(sex[,c("vek", "veksex", "pocpart", "pocpart.rok")], na.rm = T)
var(deti[, c("vek", "poc.deti")])



V nasledujúcich častiach sa postupne zameriame na niektoré základne jednovýberové, dvojvýberové párové a dvojvýberové nepárové štatistické testy. Jednotlivé testy budú ilustrované na príkladoch naviazaných a uvedené štyri datové súbory.

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() doplní 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.

  • Pre vizuálne overenie normality nejakého náhodného výberu je samozrejme vhodnejšie použiť histogram a nie empiricku distribučnú funkciu. Pre porovnanie, nagenerujte nahodný výber z ľubovolného rozdelenia (ideálne s nulovou strednou hodnotou) a porovnajte empirickú distribučnú funkciu spočítanú z tohto náhodného výberu s distribučnou funkciou štandardného normálneho rozdelenia. Analogické porovnanie urobte aj pre histogram a teoreticku hustotu.

  • Silnejší (v tomto prílade aj vhodnejší) test na testovanie normality náhodného výberu je napr. Shapiro-Wilkův test - v programe R implementovaný pod funkciou shapiro.test(). Vyskúšajte nasledujúci príkaz:

    shapiro.test(miry$vyska)
    ## 
##  Shapiro-Wilk normality test
## 
## data:  miry$vyska
## W = 0.99282, p-value = 0.01698

  • Iný test, ktorý testuje normalitu náhodného výberu na základe porovnania empirickej šikmosti a špičatosti s teoretickými hodnotami v normálnom rozdelení je Jarque-Bera test, ktorá je v Rku implementovaný v knižnici ‘tseries’ (inštalácia knižnice pomocou príkazu install.packages("tseries")).

    library(tseries)
jarque.bera.test(miry$vyska)
    ## 
##  Jarque Bera Test
## 
## data:  miry$vyska
## X-squared = 5.5187, df = 2, p-value = 0.06333

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 Kolmogorov-Smirnovovho nezamietli.
  • Pokúste sa navrhnúť malú simulačnú študiu, v ktorej porovnáte silu jednotlivých testov uvedených vyššie.

  • Jak je to v prípade porovnania požadovanej hladiny u jednotlivých testov?

    N <- c(10, 100, 1000, 5000)

power <- NULL
set.seed(1234)
for (n in N){
  p1 <- p2 <- p3 <- NULL
for (i in 1:100){
  x <- rchisq(n, df = 50)
  p1 <- c(p1, as.numeric(ks.test(x,"pnorm",mean=mean(x),sd=sd(x),exact=FALSE)$p.value < 0.05))
  p2 <- c(p2, as.numeric(shapiro.test(x)$p.value < 0.05))
  p3 <- c(p3, as.numeric(jarque.bera.test(x)$p.value < 0.05))
}
power <- rbind(power, c(n, mean(p1), mean(p2), mean(p3)))
}

plot(power[,4] ~ log10(power[,1]), col = "red", type = "l", pch = 21, bg = "red", xlab = "Logarimus rozsahu náhodného výberu", ylab = "Empirická síla testu")
points(power[,4] ~ log10(power[,1]), pch = 21, bg = "red")
points(power[,3] ~ log10(power[,1]), pch = 21, bg = "blue")
points(power[,2] ~ log10(power[,1]), pch = 21, bg = "green")

lines(power[,3] ~ log10(power[,1]), col = "blue")
lines(power[,2] ~ log10(power[,1]), col = "green")
abline(1,0, col = "black", lty = 3)

legend(1, 1, legend = c("Shapiro-Wilk test", "Jarque-Beta test", "Kolmogorov-Smirnov test"), col = c("blue", "red", "green"), lty= rep(1,3))



Jednovýberový t-test (z-test)

Tento test je určený pre testovanie nulovej hypotézy o strednej hodnote náhodného výberu. Matematicky formálne môže test zapísať následovne:


Pre náhodný výber \(X_1, \dots, X_n\) z normálneho rozdelenia \(N(\mu, \sigma^2)\), kde \(\mu \in \mathbb{R}\) je neznámy parameter strednej hodnoty a \(\sigma^2 > 0\) je známy rozptyl, testujeme nulovú hypotézu \[ H_{0}: \mu = \mu_0, \] oproti obecnej alternatíve \[ H_{1}: \mu \neq \mu_0, \] kde \(\mu_0 \in \mathbb{R}\) je predom zvolená pevna hodnota (stredná hodnota \(X_i\) za platnosti nulovej hypotézy). Za predpokladu známeho rozptylu \(\sigma^2 > 0\) je jednovýberový t-test založený na testovej štatistike \[ T = \sqrt{n} \frac{\overline{X}_n - \mu_0}{\sigma}, \] ktorá ma za platnosti nulovej hypotézy \(H_0\) štandardné normálne rozdelenie \(N(0,1)\). Nulovú hypotézu zamietame pre veľké hodnoty testovej štatistiky, konkrétne pre \(|T| > u_{1 - \alpha/2}\). Symbol \(\overline{X}_n\) označuje výberový priemer, \(u_{1 - \alpha/2}\) je príslušný \(1 - \alpha/2\) kvantil štandardného normálneho rozdelenia a \(\alpha \in (0,1)\). Jedná sa o tzv. z-test (presný).

V praxi ale častejšie nastáva situácia, že teoretický rozptyl náhodných veličín \(X_1, \dots, X_n\) je neznámy a preto je nutné použíť výberový rozptyl \(S_n^2 = \frac{1}{n - 1}\sum_{i = 1}^n (X_i - \overline{X}_n)^2\), pričom platí, že pri vhodnom pre3k8lovan9 sa jedná o náhodnú veličinu s \(\chi^2\) rozdelením s \(n - 1\) stupňami voľnosti, teda \[ (n - 1) \frac{S_n^2}{\sigma^2} \sim \chi_{n - 1}^2. \] Z definície studentovho t rozdelenia dostaneme, že prislušná testová štatistika \[ T = \sqrt{n} \frac{\overline{X}_n - \mu_0}{S_n} \] má studentovo t rozdelenie s (n - 1) stupňami voľnosti. Nulovú hypotézu zamietame pre príliš veľké, resp, príliš malé hodnoty \(T\), presnejšie ak \(|T| > q_{n - 1}(1 - \alpha/2)\), kde \(q_{n - 1}(1 - \alpha/2)\) je príslušný \(1 - \alpha/2\) kvantil studentovho t rozdelenia s \(n - 1\) stupňami voľnosti. Jedná sa o tzv. jednovýberový t-test (presný).



Pomocou programu R, s využitím predpokladu normality náhodného výberu (výška ľudí v populácii) môžeme otestovať nulovú hypotézu, že stredná hodnota výšky je 168 cm:

t.test(miry$vyska, mu = 168) 
## 
##  One Sample t-test
## 
## data:  miry$vyska
## t = -0.22664, df = 499, p-value = 0.8208
## alternative hypothesis: true mean is not equal to 168
## 95 percent confidence interval:
##  167.0099 168.7853
## sample estimates:
## mean of x 
##  167.8976

a porovnajte výsledok s

print(mean(miry$vyska))
## [1] 167.8976
print(var(miry$vyska))
## [1] 102.0674
print(sd(miry$vyska))
## [1] 10.10284
print(sqrt(var(miry$vyska)))
## [1] 10.10284
print(T <- sqrt(nrow(miry)) * (mean(miry$vyska) - 168)/sd(miry$vyska))
## [1] -0.2266425
print(2 * pt(- abs(T), df = nrow(miry) - 1))
## [1] 0.8207946



Užitočné

Niekedy nie je reálne/rozumné predpokládať normálne rozdelenie náhodného výberu. Avšak ak je možné predpokladať aspoň konečný druhý moment, možeme využíť asymptotický z-test, ktorý má testovú štatistiku \[ T = \sqrt{n} \frac{\overline{X}_n - \mu_0}{S_n}. \]

  • Testová štatistika má za platnosti nulovej hypotézy \(H_0\) asymtoticky normálne rozdelenie \(N(0,1)\). Nulovú hypotézu preto zamietame ak \(|T| > u_{1 - \alpha /2}\).

  • Použitie príkazu t.test() nie je uplne správne. Prákaz totíž spočíta správnu hodnotu testovej štatistiky, ale nesprávnu \(p\)-hodnotu. Porovnajte:

    print(2 * pnorm(- abs(T)))
print(2 * pt(- abs(T), df = nrow(miry) - 1))
  • Pomocou helpu k funkcii t.test() zistite, ako testovať a v programe R implementovať štatistický test pre jednostrannú alternatívu.



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

  • Porovnajte výsledky nerobustného t testu a robusnejšej varianty , teda znamienkového testu. Vyskúštajte rúzne nulové a alternatívne hypotézy, prípadne simulujte náhodný výber z rúznych rozdelení.
  • Pomocou analogickej simulácie ako v prípade testovania normality sa podívajte na sílu znaménkového testu a klasického t-testu (resp. z-testu).



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ů", xlab = "Rozdiel systolického a diastolického tlaku")

Ř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 = 1371104, 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ů?
  • Pomocou jednoduchej simulace porovnajte sílu vyššie uvedených troch testov.



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

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.

ks.test(sex$veksex[sex$pohl=="Male"],sex$veksex[sex$pohl=="Female"])
## Warning in ks.test(sex$veksex[sex$pohl == "Male"], sex$veksex[sex$pohl == : p-
## value will be approximate in the presence of ties
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  sex$veksex[sex$pohl == "Male"] and sex$veksex[sex$pohl == "Female"]
## D = 0.084467, p-value = 5.431e-06
## alternative hypothesis: two-sided

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.

t.test(sex$veksex[sex$pohl=="Male"],sex$veksex[sex$pohl=="Female"],var.equal=TRUE)
## 
##  Two Sample t-test
## 
## data:  sex$veksex[sex$pohl == "Male"] and sex$veksex[sex$pohl == "Female"]
## t = -3.9956, df = 3591, p-value = 6.583e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.8170783 -0.2791568
## sample estimates:
## mean of x mean of y 
##  17.35301  17.90112

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

t.test(sex$veksex[sex$pohl=="Male"],sex$veksex[sex$pohl=="Female"])
## 
##  Welch Two Sample t-test
## 
## data:  sex$veksex[sex$pohl == "Male"] and sex$veksex[sex$pohl == "Female"]
## t = -3.9985, df = 3577.5, p-value = 6.503e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.8168818 -0.2793533
## sample estimates:
## mean of x mean of y 
##  17.35301  17.90112

Welchův test dává prakticky stejné výsledky jako dvouvýběrový t-test se shodnými rozptyly.

Otázky a úlohy

  • V čem spočíva základný rozdíl mezi dvouvýběrovým párovým a dvouvýběrovým nepárovým testem?
  • Jaké jsu předpoklady jednotlivých dvouvýběrových nepárových testů, které jsme využili?
  • Který z uvedených testů je podle vás nejlepší k zrovnání věku prvního sexu mužů a žen?
  • Pomoci jednoduchých simulaci se podívejte na sílu a hladinu uveděných testů.