NMFM301, ZS 2022/23

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)


Data

Budeme používať niekoľko 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"))

Analogickým príkazom načítame aj pomocné funkcie, ktoré umožňunjú spočítať priamo niektoré užitočné charakteristiky – napr. intervaly spoľahlivosti, prípadne niektoré základné štatistické testy.

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:

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



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é?
  • 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.



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