9 Jednovýběrové neparametrické testy 9.1 Znaménkový jednovýběrový exaktní test Nechť Xi,..., X\n, n > 2 je náhodný výběr ze spojitého rozdělení a nechť xq je konstanta. Na hladině významnosti a testujeme jednu z následujících tří hypotéz oproti příslušné alternativní hypotéze. Hqi : x = xq oproti Hu : x ^ xq (oboustranná alt.) Hq2 : x < xq oproti Hi2 : x > xq (pravostranná alt.) Hqs : x > xq oproti H13 : x, < xq (levostranná alt.) Test nazýváme jednovýběrový znaménkový test o mediánu x. Testovací statistika má tvar n SE = Y/I(Xl-xo>0), (9.1) i=l I(Xi — xq > 0) je indikační funkce, která nabývá hodnoty 1, pokud Xi — xq > 0, a hodnoty 0, pokud Xi — xq < 0. Se udává tedy počet kladných rozdílů Xi — xq.Zů platnosti nulové hypotézy pochází statistika Se z binomického rozdělení s parametry N = m, kde m je počet nenulových rozdílů Xi — x0, a p = i, tj. Se = Y,1^ ~ io > °) ~ Bin m' 9 Kritický obor podle zvolené alternativní hypotézy má tvar Hn:x^x0 W = (-00 ; bm,1/2(a/2) - l) U (&m,i/2(l - a/2); 00) H12 : x > x0 W = 1/2(1 - a); 00) H13 : x < x0 W = (-00; bm^/2{a) - l) kde m je počet nenulových rozdílů Xi — x$, bmi/2(a/2), bmi/2(l — a/2), bmi/2(a) a bmi/2(l — a) jsou kvantily binomického rozdělení, jejichž hodnoty získáme pomocí C:: a implementované funkce qbinom(). Interval spolehlivosti má podle zvolené alternativní hypotézy jeden z následujících tvarů Hn:x^ x0 (d, h) = (x(fc".i/2(«/2)) ; x(6„,1/2(i-«/2)+i)) H12:í> x0 (d, 00) = (x(b".i/2(«)) ; 00) H13:x< x0 (-00, h) = (-00 ; x(b™.i/2(i-«)+i)) kde n je rozsah náhodného výběru, X^ <■■■<. X^ značí vzestupně seřazené hodnoty Xi, i = 1,..., n, a značí fc-tou hodnotu v seřazené posloupnosti X^ < ■ ■ ■ < X^nK p-hodnota má v závislosti na zvolené alternativní hypotéze jeden z následujících tvarů Hn : x, ^ xq p-hodnota = 2min{Pr(5£; < se) , Pr(SE > se)} H\2 : x, > xq p-hodnota = Pr(SE > se) H13 : x < x0 p-hodnota = Pr(5s < sE) kde Se je náhodná veličina pocházející z binomického rozdělení, se je realizace testovací statistiky Se (viz vzorec 9.1), tedy konkrétní číslo, Pr(SE > se) = 1 — Pt(Se < se) = 1 — Pt(Se < se — 1), což vyplývá z faktu, že náhodná veličina Se pochází z binomického (diskrétního) rozdělení (viz kapitola ??), a Pr(SE < se), resp. Pr(SE < se — 1) je distribuční funkce binomického rozdělení, jejíž hodnotu získáme pomocí a implementované funkce pbinom(). Poznámka: Všimněme si, že ve vzorcích intervalu spolehlivosti figuruje rozsah náhodného výběru n, zatímco ve vzorcích hranic kritického oboru a u p-hodnoty pracujeme s počtem nenulových rozdílů m. 1 Příklad 9.1. Znaménkový jednovýběrový exaktní test Mějme datový soubor 21-goldman-measures.csv a proměnnou humer.DR popisující největší délku hlavice ramenní kosti z pravé strany v mm u skeletů z období raného středověku z oblasti Staubing (viz sekce ??). Dále máme k dispozici údaje ze studie (Mali et al.) z roku 1999, v rámci které byla měřena největší délka hlavice ramenní kosti u mužů současné německé populace (mm = 50.0mm, nm = 64). Na hladině významnosti a = 0.05 zjistěte, zda existuje rozdíl mezi největší délkou hlavice ramenní kosti z pravé strany u mužů z raně středověké německé populace a u mužů současné německé populace. Řešení příkladu 9.1 Pomocí příkazu read.delim() načteme datový soubor a pomocí operátoru [] vybereme z datové tabulky údaje o největší délce hlavice ramenní kosti z pravé strany (humer.DR) u mužů (sex == 'm') z oblasti Staubing (pop == 'Staubing'). Z vektoru naměřených hodnot odstraníme chybějící údaje (na.omit()) zjistíme rozsah náhodného výběru (lengthO). 1 data <- read.delim('00-Data//21-goldman-measures.csv', sep = ';') 2 head(data) 3 humer.DRM <- data[data$pop == 'Staubing' k data$sex == 'm', 'humer.DR'] 4 humer.DRM <- na.omit(humer.DRM) 5 n <- length(humer.DRM) # 12 Datový soubor obsahuje údaje o největší délce hlavice ramenní kosti z pravé strany u 12 mužů raně středověké německé populace. Naším úkolem ze zadání je porovnat střední hodnoty dvou německých populací, přičemž u jedné populace (německé raně středověké populace) máme k dispozici naměřené hodnoty. Na základě těchto hodnot můžeme zjistit, zda náhodný výběr pochází z normálního rozdělení, tj. zda náhodná veličina X popisující největší délku hlavice ramenní kosti u mužů raně středověké populace pochází z normálního rozdělení, tj. X ~ N(fi, a2), kde skutečný rozptyl a2 neznáme. Druhá populace (současná německá populace) je reprezentována pouze hodnotou aritmetického průměru (^r = 50.00mm). Řešení příkladu tedy vede na situaci, kdy střední hodnotu jednoho náhodného výběru (jehož skutečnou hodnotu rozptylu neznáme) porovnáváme s konkrétním číslem, tedy na jednovýběrový test o střední hodnotě /i při neznámém rozptylu a2 (viz kapitola ??). Jediným předpokladem k použití tohoto testu je normalita náhodného výběru naměřených délek hlavic ramenních kostí. Před použitím testu tedy tento předpoklad ověříme. Jelikož je rozsah náhodného výběru menší než 30, ověříme předpoklad normality Shapiro-Wilkovým testem (a = 0.05). Grafické ověření provedeme na základě kvantilového diagramu a histogramu superponovaného křivkou normálního rozdělení, jejíž parametry odhadneme pomocí výběrového průměru a výběrového rozptylu (viz obrázek 1). Datový soubor rozdělíme do pěti ekvidistatních intervalů s šířkou 3 mm prostřednictvím stanovených hranic 37, 40,..., 52. Protože p-hodnota = 0.0261 je menší než 0.05, nulovou hypotézu o normalitě dat zamítáme na hladině významnosti a = 0.05. Z histogramu na obrázku 1 vidíme, že naměřené hodnoty jsou vyšikmené doprava s chvostem na levé straně a tvar křivky hustoty normálního rozdělení příliš dobře nekopírují. Na kvantilovém diagramu je zřejmé odchýlení bodů od referenční křivky zejména na levém chvostu. Náhodný výběr největších délek hlavice ramenní kosti mužů z raně středověké německé populace tedy nepochází z normálního rozdělení. Protože náhodný výběr nepochází z normálního rozdělení, nemůžeme hypotézu ze zadání otestovat pomocí parametrického testu o střední hodnotě jj, podle metod uvedených v sekci ??. K testování hypotézy musíme použít neparametrickou alternativu tohoto testu. Zde konkrétně použijeme znaménkový jednovýběrový test. Jelikož rozsah náhodného výběru je menší než 30, použijeme exaktní variantu znaménkového testu. Naším úkolem je zjistit, zda existuje rozdíl mezi největší délkou hlavice ramenní kosti z pravé strany u mužů z raně středověké německé populace a u mužů současné německé populace. Tato věta bude součástí alternativní hypotézy, neboť rozdíl implikuje nerovnost a nerovnost je vždy součástí alternativní hypotézy. Nulovou hypotézu potom stanovíme jako doplněk k alternativní hypotéze. Závěrem poznamenejme, že zatímco u parametrických testů (viz kapitola ??) figuruje v 2 41.25 46.25 51.25 -1.5 -0.5 0.5 1.0 1.5 teoreticky kvantil polomer hlavice ramenni kosti (v mm) Shapiruv test: p-hodnota = 0.0261 Obrázek 1: Histogram a kvantilový diagram poloměru hlavice ramenní kosti na pravé straně u skeletů mužů z raně středověké populace z oblasti Staubing hypotézách pojem střední hodnota, u neparametrických testů tento pojem nahrazujeme neparametrických ekvivalentem, kterým je medián. Testování provedeme v posloupnosti sedmi kroků. 1. Stanovení hypotéz • slovní formulace nulové a alternativní hypotézy Hq : Medián největší délky hlavice ramenní kosti raně středověké německé populace je shodný s mediánem největší délky hlavice ramenní kosti na pravé straně mužů současné německé populace. Hi : Medián největší délky hlavice ramenní kosti raně středověké německé populace není shodný s mediánem největší délky hlavice ramenní kosti na pravé straně mužů současné německé populace. • matematická formulace nulové a alternativní hypotézy Hq : x = xq, kde xq = 50.00 Hi : x, ^ xq, kde xq = 50.00 (oboustranná alternativa) 2. Volba hladiny významnosti • Hladinu významnosti volíme v souladu se zadáním jako a = 0.05. 3. Testování kritickým oborem • Testovací statistika Nejprve vypočítáme vektor rozdílů naměřených hodnot Xi a konstanty xq, tj. Xí—xq, a stanovíme, které rozdíly jsou kladné, (viz tabulka 1). Tabulka 1: Naměřené hodnoty Xi, rozdíly Xi — xq a znaménka těchto rozdílů měření 1 2 3 4 5 6 7 8 9 10 11 12 X, 49.83 47.96 49.49 48.73 49.25 40.90 51.22 51.27 46.16 49.52 50.28 45.60 Xi — x0 -0.17 -2.04 -0.51 -1.27 -0.75 -9.10 1.22 1.27 -3.84 -0.48 0.28 -4.40 +/- - - - - - - + + - - + - Z tabulky 1 vidíme, že celkem tři rozdíly Xi — xq, i = 1,..., 12 jsou kladné. Hodnota testovací statistiky Se, která je definovaná jako počet kladných rozdílů, bude tedy rovná 3. n 13 Se = ^2l(Xi - x0 > 0) = X)/(Xi - 50 > 0) = 3. i=l i=l 3 6 7 8 9 10 11 12 13 14 xO <- 50.00 I <- (numer.DRM > xO) tab <- data.frame(rbind("Xi" = humer.DRM, "Xi-xO" = humer.DRM - xO , "+/-names(tab) <- 1 : 12 # 1 2345678 9 10 11 # Xi 49.83 47.96 49.49 48.73 49.25 40.9 51.22 51.27 46.16 49.52 50.28 # Xi-xO -0.17 -2.O4 -0.51 -1.27 -0.75 -9.1 1.22 1.27 -3.84 -O.48 0.28 # +/- 0.00 0.00 0.00 0.00 0.00 0.0 1.00 1.00 0.00 0.00 1.00 SE <- sum(I) # 3 = D) 12 45. 6 -4.4 0.0 • Kritický obor Z tabulky 1 dále vidíme, že žádný z celkového počtu n = 12 rozdílů není nulový, a tedy počet nenulových rozdílů m = 12. Počet nenulových rozdílů využíváme při stanovení hranic kritického oboru a později také u výpočtu p-hodnoty. W = (-00 ; &m,i/2(a/2) - 1) U <6m,1/2(l - a/2); 00) = (-00 ; 6i2,i/2(0.05/2) - 1) U (612,1/2(1 - 0.05/2); 00) = (-00 ; 6i2,i/2(0.025) - 1) U <6i2,i/2(0.975); 00) = (-00; 3 - 1) U (9; 00) = (-00; 2) U (9; 00) 15 m <- sum(humer.DRM - xO != 0) #12 16 alpha <- 0.05 17 ql <- qbinom(alpha / 2 , m, 1 / 2) - 1 # 2 18 q2 <- qbinomd - alpha / 2, m, 1 / 2) # 9 • Závěr testování Protože realizace testovací statistiky s e = 3 nenáleží do kritického oboru, tj. s e W, Hq nezamítáme na hladině významnosti a = 0.05. 4. Testování intervalem spolehlivosti • Interval spolehlivosti Abychom mohli stanovit hranice 95 % intervalu spolehlivosti, musíme nejprve naměřené hodnoty největší délky hlavice ramenní kosti vzestupně seřadit. To provedeme jednoduše příkazem sort(). Seřazené hodnoty jsou k dispozici v tabulce 2. Tabulka 2: Seřazené hodnoty největší délky hlavice ramenního kloubu pořadí I) 1 2 3 4 5 6 7 8 9 10 11 13 seřazené Xr || 40.90 45.60 46.16 47.96 48.73 49.25 49.49 49.52 49.83 50.28 51.22 51.27 Hranice intervalu spolehlivosti potom tvoří ty hodnoty, které se v seřazeném vektoru hodnot nachází na (6IJ,i/2(a/2))-té pozici a na (&rs,i/2(1 — a/2) + l)-té pozici. Hodnoty &„,i/2(a/2) a 6^,1/2(1 — a/2) + 1 nalezneme opět pomocí funkce qbinomQ. x(612,1/2(0.05/2)) . X(6i2,1/2(1-0.05/2) + !) (46.16; 50.28) 4 19 qbinom(alpha / 2, n, 1/2) #3 20 qbinomd - alpha / 2 , n, 1 / 2) + 1 # 10 21 humer.DRMs [3] # 46.16 22 humer . DRMs [10] # 50.28 • Závěr testování Protože xq = 50.00 náleží do 95% empirického oboustranného intervalu spolehlivosti, tj. xq = 50.00 G IS, Hq nezamítáme na hladině významnosti a = 0.05. 5. Testování p-hodnotou • p-hodnota p-hodnota = 2min{Pr(5E < sE), Pr(se > Se)} = 2mm{Pľ(SE < 3), Pr(SE > 3)} = 2min{Pr(5£; < 3), 1 - Pr(SE < 2)} = 2min{0.07299805, 0.9807129} = 2 x 0.07299805 = 0.1459961 = 0.1460 23 2 * min (pbinom(SE, m, 1/2), 1 - pbinom(SE - 1, m, 1 / 2)) # 0.1459961 • Závěr testování Protože p-hodnota = 0.1460 je větší než a = 0.05, Hq nezamítáme na hladině významnosti a = 0.05. 6. Interpretace výsledků Za základě všech tří typů testování nezamítáme nulovou hypotézu na hladině významnosti a = 0.05. Mezi největší délkou hlavice ramenní kosti na pravé straně u mužů raně středověké a současné populace neexistuje statisticky významný rozdíl. 7. Grafická vizualizace výsledku testování Porovnání náhodného výběru s konstantou xq = 50.00 zobrazíme nejlépe pomocí krabicového diagramu (viz obrázek 2). Poznámka: Znaménkový jednovýběrový exaktní test můžeme provést pomocí funkce SIGN.test() implementované v knihovně BSDA. Vstupními parametry budou vektor reprezentující náhodný výběr (humer.DRM), hodnota parametru xq z nulové hypotézy zadaná argumentem md = 50.00, hodnota hladiny významnosti a zadaná prostřednictvím koeficientu spolehlivosti 1 — a nastavením hodnoty argumentu conf.level = 0.95 a typ zvolené alternativní hypotézy (oboustranná) zadaný pomocí argumentu alternativě = 'two.sideď. Součástí vystupuje hodnota mediánu náhodného 24 BSDA::SIGN.test(humer.DRM, md = 50.0, conf.level = 0.95, 25 alternativě = 'two.sided') # IS interpolovaný výběru medián of x = 49.37, hodnota testovací statistiky s = 3, interpolované hranice 95% Waldova empirického oboustranného intervalu spolehlivosti 46.35145 a 50.23214, které jsou mírně přesnější, než námi stanovené hranice intervalu spolehlivosti (zpřesnění hranic bylo provedeno procesem nazývaným interpolace), a p-hodnota p-value = 0.146. Jediné, co musíme stanovit zvlášť, jsou dolní a horní hranice kritického oboru. 5 55 - rane stredoveká nemecká p. současna nemecká p. 3 50 - 45 - 40 Obrázek 2: Krabicový diagram poloměru hlavice ramenní kosti na pravé straně u skeletů mužů z raně středověké populace z oblasti Staubing One-sample Sign-Test data: humer.DRM s = 3, p-value = 0.146 alternative hypothesis: true median is not equal to 50 95 percent confidence interval: 46.35145 50.23214 sample estimates: median of x 49.37 Achieved and Interpolated Confidence Intervals: Lower Achieved CI Interpolated CI Upper Achieved CI Conf.Level L.E.pt U.E.pt 0.8540 47.9600 49.8300 0.9500 46.3515 50.2321 0.9614 46.1600 50.2800 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 6 Příklad 9.2. Znaménkový jednovýběrový exaktní test Mějme datový soubor 15-anova-means-skull.txt a proměnnou upface.H popisující výšku horní části tváře mužů německé populace (viz sekce ??). Dále máme k dispozici údaje o výšce horní části tváře mužů Cernjachovské kultury na území dnešní Ukrajiny (mck = 70.00 mm, nck = 99). Na hladině významnosti a = 0.10 testujte hypotézu, že výška horní části tváře německé mužské populace je menší nebo rovna výšce horní části tváře mužské populace Cernj achovské kultury. Řešení příkladu 9.2 Úvodu tohoto příkladu jsme se věnovali v příkladu ?? v rámci sekce ??. Zde jsme za základě testování normality pomocí Shapiro-Wilkova testu (p-hodnota = 0.0419) a grafické vizualizace (kombinace histogramu a kvantilového diagramu) dospěli k závěru, že náhodný výběr výšek horní části tváře u 19 mužů německé populace nepochází z normálního rozdělení (pro připomenutí viz obrázek 3). Obrázek 3: Histogram a kvantilový diagram výšky horní části tváře u mužů německé populace Jelikož data nepochází z normálního rozdělení, není možné otestovat hypotézu ze zadání pomocí parametrického testu o střední hodnotě /i, ale bude potřeba použít neparametrickou alternativu tohoto testu. Vzhledem k nízkému rozsahu náhodného výběru použijeme k otestování hypotézy ze zadání exaktní znaménkový jednovýběrový test. Naším úkolem je testovat (nulovou) hypotézu, že výška horní části tváře německé mužské populace je menší nebo rovna výšce horní části tváře mužské populace Cernj achovské kultury. Alternativní hypotézu potom stanovíme jako doplněk k nulové hypotéze. Protože budeme hypotézu ze zadání testovat pomocí neparametrického testu, bude v jejím přesném znění figurovat pojem medián namísto pojmu střední hodnota. 1. Stanovení hypotéz • slovní formulace nulové a alternativní hypotézy H0 : Medián výšky horní části tváře nčmecké mužské populace je menší nebo roven mediánu výšky horní části tváře mužské populace z Cernjachovské kultury. Hi : Medián výšky horní části tváře německé mužské populace je větší než medián výšky horní části tváře mužské populace z Cernjachovské kultury. • matematická formulace nulové a alternativní hypotézy Hq : x < xq, kde xq = 70.00 Hi : x, > xq, kde xq = 70.00 (pravostranná alternativa) 2. Volba hladiny významnosti • Hladinu významnosti volíme v souladu se zadáním jako a = 0.10. 3. Testování kritickým oborem 7 • Testovací statistika Nejprve vypočítáme vektor rozdílů naměřených hodnot Xi a konstanty xq, tj. Xí—xq, a stanovíme, které rozdíly jsou kladné (viz tabulka 3). Tabulka 3: Naměřené hodnoty Xi, rozdíly Xi — xq a znaménka těchto rozdílů měření 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 X, 73 73 67 75 70 62 76 73 71 66 76 73 74 74 70 76 72 69 76 Xi — x0 3 3 -3 5 0 -8 6 3 1 -4 6 3 4 4 0 6 2 -1 6 +/- + + - + 0 - + + + - + + + + 0 + + - + Z tabulky 3 vidíme, že celkem 13 rozdílů Xi — xq, i = 1,..., 19 je kladných. Hodnota testovací statistiky Se, která je definovaná jako počet kladných rozdílů, bude tedy rovná 13. n 19 SE = - x0 > 0) = ^2l(Xi - 70.00 > 0) = 13. i=l i=l 44 xO <- 70.00 45 I <- (upface.HN > xO) 46 tab <- data.frame(rbind("Xi" = upface.HN, "Xi-xO" = upface.HN - xO, "+/-" = I)) 47 names(tab) <- 1 : 19 48 # 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 49 # Xi 73 73 67 75 70 62 76 73 71 66 76 73 74 74 70 76 72 69 76 50 # Xi-xO 3 3 -3 5 0 -8 6 3 1 -4 6 3 4 4 0 6 2 -1 6 51 # +/- 1101001110111101101 52 SE <- sum(I) # 13 • Kritický obor Z tabulky 3 dále vidíme, že z celkového počtu n = 19 rozdílů jsou dva rozdíly nulové, a tedy počet nenulových rozdílů m = 17. Počet nenulových rozdílů využijeme při stanovení hranic kritického oboru a u výpočtu p-hodnoty. W = (bm,i/2(a); 00) = <6i7,i/2(0.10); 00) = (10; 00) 53 m <- sum(upface.HN - xO != 0) # 17 54 alpha <- 0.10 55 q <- qbinomd - alpha, m, 1/2) - 1 # 10 • Závěr testování Protože realizace testovací statistiky s e = 13 náleží do kritického oboru, tj. s e G W, Hq zamítáme na hladině významnosti a = 0.10. 4. Testování intervalem spolehlivosti • Interval spolehlivosti Abychom mohli stanovit hranice 90% intervalu spolehlivosti, musíme nejprve naměřené hodnoty výšky horní části lebky vzestupně seřadit. To provedeme příkazem sort(). Seřazené hodnoty jsou k dispozici v tabulce 4. 8 Tabulka 4: Seřazené hodnoty naměřených výšek horní části tváře pořadí || 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 seřazené X, || 62 66 67 69 70 70 71 72 73 73 73 73 74 74 75 76 76 76 76~ Hranice intervalu spolehlivosti potom tvoří hodnota, která se v seřazeném vektoru hodnot nachází na (&n,i/2(a))-té pozici a nekonečno. (d,QO) = ^(b",l/2(«)) ; = ^(b1B,!/2(0.10)) . ^ = (X^ ; oo) = (71; oo) 56 qbinom(alpha, n, 1/2) #7 57 upface.HNs [7] # 71 • Závěr testování Protože xq = 70.00 nenáleží do 90% empirického levostranného intervalu spolehlivosti, tj. xq = 70.00 ^ IS, Hq zamítáme na hladině významnosti a = 0.10. 5. Testování p-hodnotou • p-hodnota p-hodnota = Pr(SE > sE) = Pľ(SE > 13) = 1 — Pr(5£; < 12) = 1 - 0.9754791 = 0.0245209 = 0.02452 58 1 - pbinom(SE - 1, m, 1/2) # 0. 02452087 • Závěr testování Protože p-hodnota = 0.02452 je menší než a = 0.10, Hq zamítáme na hladině významnosti a = 0.10. 6. Interpretace výsledků Za základě všech tří typů testování zamítáme nulovou hypotézu na hladině významnosti a = 0.10. Výška horní části tváře mužů německé populace je statisticky významně větší než výška horní části tváře mužů Cernjachovské populace z oblasti dnešní Ukrajiny. 7. Grafická vizualizace výsledku testování Významný rozdíl horní části tváře mužů mezi oběma populacemi vizualizujeme krabicovým diagramem (viz obrázek 4). Poznámka: K provedení znaménkového jednovýběrového exaktního testu můžeme použít funkci SIGN.test() z knihovny BSDA. Vstupními parametry budou vektor reprezentující náhodný výběr (upface.HN), hodnota parametru xq z nulové hypotézy (md = 70.00), hodnota hladiny významnosti a zadaná prostřednictvím koeficientu spolehlivosti 1 — a (conf.level = 0.90) a typ zvolené alternativní hypotézy (pravostranná, alternativě = 'greater'). Součástí výstupu je hodnota mediánu náhodného výběru medián of x = 73, hodnota testovací statistiky s = 13, interpolované hranice 90% empirického levostranného intervalu spolehlivosti 71.17133 a Inf a p-hodnota p-value = 0.02452. Jediné, co musíme stanovit zvlášť, je dolní hranice kritického oboru. 9 80 - 75 - nemecká p. cernjachovska k. o 70 65 Obrázek 4: Krabicový diagram výšky horní části tváře u mužů německé populace 59 BSDA::SIGN.test(upface.HN, md = 70.00, conf.level = 0.90, alt = 'greater') # IS interpo Iovany 60 One-sample Sign-Test 61 62 data: upface.HN 63 s = 13, p-value = 0.02452 64 alternative hypothesis: true median is greater than 70 65 90 percent confidence interval: 66 71.17133 Inf 67 sample estimates: 68 median of x 69 73 70 71 Achieved and Interpolated Confidence Intervals: 72 73 Conf.Level L.E.pt U.E.pt 74 Lower Achieved CI 0.8204 72.0000 Inf 75 Interpolated CI 0.9000 71.1713 Inf 76 Upper Achieved CI 0.9165 71.0000 Inf 77 10 Příklad 9.3. Znaménkový jednovýběrový exaktní test Mějme datový soubor 21-goldman-measures.csv a proměnnou tibia.LR popisující největší délku lýtkové kosti z pravé strany v mm u skeletů z období neolitu z oblasti Yoshigo (viz sekce ??). Dále máme k dispozici údaje ze studie (Hasegawa et al.) z roku 2009, v rámci které byla měřena největší délka lýtkové kosti z pravé strany u žen současné japonské populace (uif = 329.40mm, Sf = 17.3mm, rif = 342). Na hladině významnosti a = 0.01 zjistěte, zdaje největší délka lýtkové kosti z pravé strany u žen z neolitické japonské populace významně menší než u žen současné japonské populace. Řešení příkladu 9.3 Načteme datový soubor a pomocí operátoru [] vybereme z datové tabulky údaje o největší délce lýtkové kosti z pravé strany (tibia.LR) u žen (sex == 'f) z oblasti Yoshigo (pop == 'Yoshigo Shell Mounď). Z vektoru naměřených hodnot odstraníme chybějící údaje a zjistíme rozsah náhodného výběru. 78 data <- read.delim('00-Data//21-goldman-measures.csv ' , sep = ';') 79 # head (data) 80 tibia.LRF <- data[data$pop == 'Yoshigo Shell Mounď k data$sex == 'i', 'tibia.LR'] 81 tibia.LRF <- na.omit(tibia.LRF) 82 n <- length(tibia.LRF) # 8 Datový soubor obsahuje údaje o největší délce lýtkové kosti z pravé strany u 8 žen z neolitické japonské populace. V příkladu se zaměříme na porovnání délky lýtkové kosti dvou japonských populací, přičemž u jedné populace (neolitická populace) máme k dispozici naměřené hodnoty. Na základě těchto hodnot můžeme zjistit, zda náhodná veličina X popisující největší délku lýtkové kosti u žen neolitické japonské populace pochází z normálního rozdělení, tj. X ~ N(fi,a2), kde skutečný rozptyl a2 neznáme. Druhá populace (současná japonská populace) je reprezentována pouze hodnotou aritmetického průměru (uif = 329.40mm) a směrodatné odchylky (sf = 17.3mm). Řešení příkladu vede na situaci, kdy střední hodnotu jednoho náhodného výběru porovnáváme s konkrétním číslem, tedy na jednovýběrový test o střední hodnotě jj, při neznámém rozptylu a2. Před použitím parametrického testuje třeba ověřit normalitu náhodného výběru naměřených délek lýtkových kostí. Předpoklad normality ověříme Shapiro-Wilkovým testem (a = 0.05), kvantilovým diagramem a histogramem (viz obrázek 5). Datový soubor rozdělíme do čtyř ekvidistatních intervalů s šířkou 9 mm prostřednictvím stanovených hranic 297, 306,..., 333. 301.5 310.5 319.5 328.5 -1.5 -0.5 0.0 0.5 1.0 1.5 teoreticky kvantil délka lýtkové kosti (v mm) Shapiruv test: p-hodnota = 0.034 Obrázek 5: Histogram a kvantilový diagram délky lýtkové kosti na pravé straně u skeletů žen neolitické populace z oblasti Yoshigo Shell Mound Protože p-hodnota = 0.034 je menší než 0.05, nulovou hypotézu o normalitě dat zamítáme na hladině významnosti 11 a = 0.05. Z histogramu na obrázku 5 vidíme, že naměřené hodnoty charakterem normálního rozdělení příliš nedisponují. Navíc body v kvantilovém diagramu se realizují jve shlucích, kde se střídavě vzdalují a přibližují k referenční přímce. Náhodný výběr největších délek lýtkových kostí žen z neolitické japonské populace nepochází z normálního rozdělení. Protože náhodný výběr nepochází z normálního rozdělení, nemůžeme k testování použít parametrický test o střední hodnotě jj,. Testování tedy provedeme na základě exaktního neparametrického znaménkového jednovýběrového testu. Naším úkolem je zjistit, zda je největší délka lýtkové kosti z pravé strany u žen z neolitické japonské populace významně menší než u žen současné japonské populace. Tato věta bude zněním alternativní hypotézy, neboť v zadání se zaměřujeme na nerovnost, přičemž nikde není zmínka o znění nulové hypotézy. Nulovou hypotézu stanovíme následně jako doplněk k alternativní hypotéze. 1. Stanovení hypotéz • slovní formulace nulové a alternativní hypotézy Hq : Medián největší délky lýtkové kosti z pravé strany u žen z neolitické japonské populace je větší nebo roven mediánu největší délky lýtkové kosti z pravé strany u žen současné japonské populace. Hi : Medián největší délky lýtkové kosti z pravé strany u žen z neolitické japonské populace je menší než medián největší délky lýtkové kosti z pravé strany u žen současné japonské populace. • matematická formulace nulové a alternativní hypotézy H0 : x > x0, kde x0 = 329.40 Hi : x, < xq, kde xq = 329.40 (levostranná alternativa) 2. Volba hladiny významnosti • Hladinu významnosti volíme podle zadání a = 0.01. 3. Testování kritickým oborem • Testovací statistika Nejprve vypočítáme vektor rozdílů naměřených hodnot Xi a konstanty xq, tj. Xi — xq a stanovíme znaménka těchto rozdílů (viz tabulka 6). Tabulka 5: Naměřené hodnoty Xi, rozdíly Xi — xq a znaménka těchto rozdílů měření 1 2 3 4 5 6 7 8 X, 300.0 299.0 323.0 323.0 322.0 301.0 331.5 322.0 X, - x0 -29.4 -30.4 -6.4 -6.4 -7.4 -28.4 2.1 -7.4 +/- - - - - - - + - Z tabulky 6 vidíme, že pouze jeden rozdíl Xi — xq, i = 1,..., 8 je kladný. Hodnota testovací statistiky Se bude tedy rovná 1. n 8 SE = - x0 > 0) = ^2l(Xi - 329.40 > 0) = 1. i=l i=l • Kritický obor Z tabulky 6 dále vidíme, že žádný z celkového počtu n = 8 rozdílů není nulový, a tedy počet nenulových rozdílů m = 8. W = (-oo; &m,i/2(aO - 1) = (-oo;68,1/2(0.01)-l) = (-oo; 1-1) = (-oo; 0) 12 83 xO <- 329.40 84 I <- (tibia.LRF > xO) 85 tab <- data.frame(rbind("Xi" = tibia.LRF, "Xi-xO" = tibia.LRF - xO, "+/-" = I)) 86 names(tab) <- 1 : 8 87 # 12345678 88 # Xi 300.0 299.0 323.0 323.0 322.0 301.0 331.5 322.0 89 # Xi-xO -29.4 -30.4 ~6-4 ~6-4 ~7 ■ 4 -28.4 2-1 ~7■ 4 90 # +/- 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 91 SE <- sum(I) # 1 92 m <- sum(tibia.LRF - xO != 0) #8 93 alpha <- 0.01 94 ql <- qbinom(alpha , m, 1/2) - 1 # 1 • Závěr testování Protože realizace testovací statistiky s e = 1 nenáleží do kritického oboru, tj. s e W, Hq nezamítáme na hladině významnosti a = 0.01. 4. Testování intervalem spolehlivosti • Interval spolehlivosti Abychom mohli stanovit hranice 95% pravostranného intervalu spolehlivosti, musíme nejprve naměřené hodnoty největší délky lýtkové kosti vzestupně seřadit. Seřazené hodnoty jsou k dispozici v tabulce 6. Tabulka 6: Seřazené hodnoty největší délky lýtkové kosti pořadí || 1 2 3 4 5 6 7 8 seřazené Xl || 300.0 299.0 323.0 323.0 322.0 301.0 331.5 322.0 Hranice intervalu spolehlivosti potom tvoří mínus nekonečno a hodnota, která se v seřazeném vektoru hodnot nachází na (&nii/2(l — a) + l)-té pozici. (d,h)= (-oo; X^-1^1"^1)) = (-oo; x(f>s.i/2(i-0.0i+i)^ = (-oo; X 2 340 - 330 - 320 S 310 - 300 14 98 BSDA::SIGN.test(tibia.LRF, md = 329.4, conf.level = 0.99, alt = 'less', exact = F) # interopoIovany IS 99 One-sample Sign-Test 100 101 data: tibia.LRF 102 s = 1, p-value = 0.03516 103 alternative hypothesis: true median is less than 329.4 104 99 percent confidence interval: 105 -Inf 329.8425 106 sample estimates: 107 median of x 108 322 109 110 Achieved and Interpolated Confidence Intervals: 111 112 Conf.Level L.E.pt U.E.pt 113 Lower Achieved CI 0.9648 -Inf 323.0000 114 Interpolated CI 0.9900 -Inf 329.8425 115 Upper Achieved CI 0.9961 -Inf 331.5000 116 15 9.2 Znaménkový jednovýběrový asymptotický test Pro náhodný výběr o rozsahu n > 30 máme možnost použít k otestování nulové hypotézy asymptotickou variantu testu. Tuto variantu nazýváme znaménkovým jedno výběrovým asymptotickým testem. Testovací statistika asymptotického variantu testu má tvar Sa = (9-2) 4 kde Se je testovací statistika definovaná vztahem 9.1 a m je počet nenulových rozdílů Xi — xq. Za platnosti nulové hypotézy pochází statistika S a ze standardizovaného normálního rozdělení, tj. SA = £° N(0,1). V 4 Kritický obor podle zvolené alternativní hypotézy má tvar H n :ž^í0 W = (-oo ; ua/2) U («i_„/2 ; oo) H12 : x, > x0 W = (mi_„ ; oo) íři3 : x < Xq W = (—oo ; ua) kde ua/2, Ui-a/2, ua, «i-a jsou kvantily standardizovaného normálního rozdělení, jejichž hodnoty získáme pomocí Qt a implementované funkce qnorm(). Interval spolehlivosti má podle zvolené alternativní hypotézy jeden z následujících tvarů Hn:x^ x0 (d, h) = (X^-^ ■ x("+1-Cl—/2)) HV2:x> x0 (d, oo) = [X^1-^ ; oo) Hi3:x s a)} Hi2 : x, > xq p-hodnota = Pr(5U > «a) His : x, < xq p-hodnota = Pr(5U < «a) kde Sa je náhodná veličina, sa je realizace testovací statistiky Sa (viz vzorec 9.2), tedy konkrétní číslo, Pr(5^ > s a) = 1 — Pr(5U < s a) = 1 — Pr(5^ < s a), což vyplývá z faktu, že náhodná veličina S a pochází z normálního (spojitého) rozdělení (viz kapitola ??), a Pr(5^ < s a) je distribuční funkce standardizovaného normálního rozdělení, jejíž hodnotu získáme pomocí ^ a implementované funkce pnorm(). Poznámka: Všimněme si, že ve vzorcích intervalu spolehlivosti figuruje rozsah náhodného výběru n, zatímco ve vzorcích testovací statistiky a hranic kritického oboru pracujeme s počtem nenulových rozdílů m. 16 Příklad 9.4. Znaménkový jednovýběrový asymptotický test Mějme datový soubor 19-more-samples-correlations-skull.txt a proměnnou nose.B popisující šířku nosu v mm (viz sekce ??). Dále máme k dispozici údaje o šířce nosu mužů současné malajské populace (mm = 26.90 mm, nm = 45) uveřejněné ve studii (Ibrahim, 2017). Na hladině významnosti a = 0.01 testujte hypotézu o shodě šířky nosu starověké malajské populace a současné malajské populace. Řešení příkladu 9.4 Načteme datový soubor a vybereme z datové tabulky naměřené šířky nosu (nose.B) mužů malajské populace (pop == 'mal'). Nakonec z vektoru naměřených hodnot odstraníme chybějící hodnoty a zjistíme rozsah náhodného výběru. 117 data <- read.delim('00-Data//19-more-samples-correlations-skuli.txt') 118 # head(data) 119 nose.BM <- data[data$pop == 'mal', 'nose.B'] 120 nose.BM <- na.omit(nose.BM) 121 n <- length(nose.BM) # 73 Datový soubor obsahuje údaje o šířce nosu u 73 mužů starověké malajské populace. V příkladu se zaměříme na porovnání šířek nosu dvou malajských populací, přičemž u jedné populace (starověká) máme k dispozici naměřené hodnoty, na základě kterých můžeme zjistit, zda náhodná veličina X popisující šířku nosu u mužů této populace pochází z normálního rozdělení, tj. X ~ 7V(p, tr2), kde skutečný rozptyl a2 neznáme. Druhá populace (současné) je reprezentována pouze hodnotou aritmetického průměru (mm = 26.90mm). Řešení příkladu tedy vede na situaci, kdy střední hodnotu jednoho náhodného výběru porovnáváme s konkrétním číslem, tedy na jednovýběrový test o střední hodnotě p při neznámém rozptylu a2. Nejprve je však potřeba ověřit, zda náhodný výběr šířek nosu u mužů starověké malajské populace pochází z normálního rozdělení. Vzhledem k velkému rozsahu náhodného výběru otestujeme normalitu Lillieforsovým testem (a = 0.05) v kombinaci s kvantilovým diagramem a histogramem. Naměřené hodnoty rozdělíme do sedmi ekvidistantních intervalů se šířkou 1.3 mm prostřednictvím stanovených hranic 20.9, 22.2,..., 30.0 (viz obrázek 7). 21.55 24.15 26.75 29.35 -2 -1 0 1 2 teoreticky kvantil sirka nosu (v mm) Lillieforsuv test: p-hodnota = 0.0094 Obrázek 7: Histogram a kvantilový diagram šířky nosu mužů malajské populace Protože p-hodnota = 0.0094 je menší než 0.05, nulovou hypotézu o normalitě náhodného výběru šířek nosu mužů starověké malajské populace zamítáme na hladině významnosti a = 0.05. Z pohledu na histogram by se mohlo zdát, že naměřené hodnoty kopírují křivku hustoty normálního rozdělení dostatečně. Je však potřeba si uvědomit, že rozsah náhodného výběru je již celkem vysoký a při takovém počtu hodnot by podobnost histogramu s křivkou hustoty měla být mnohem vyšší. Z kvantilového grafu je potom jasně patrné, že vykreslené body se v těsném okolí 17 referenční křivky příliš nepochybují. Náhodný výběr naměřených šířek nosu u mužů starověké malajské populace tedy nepochází z normálního rozdělení. Protože se náhodný výběr neřídí normálním rozdělením, nemůžeme hypotézu ze zadání otestovat pomocí parametrického testu o střední hodnotě jj,. K testování hypotézy použijeme neparametrický znaménkový jednovýběrový test, zde konkrétně jeho asymptotickou variantu, jelikož rozsah náhodného výběru je větší než 30. Připomeňme, že při použití neparametrických tetsů pracujeme s mediány namísto se středními hodnotami. Naším úkolem je otestovat hypotézu o shodě šířky nosu starověké malajské populace a současné malajské populace. Tato věta je zněním nulové hypotézy, jednak proto, že v zadání přímo o nulové hypotéze mluvíme a jednak proto, že shoda implikuje rovnost a rovnost je vždy součástí nulové hypotézy. Zbývá dodefinovat alternativní hypotézu tak, aby byla doplňkem k nulové hypotéze. Testování provedeme v posloupnosti sedmi kroků. 1. Stanovení hypotéz • slovní formulace nulové a alternativní hypotézy Hq : Medián šířky nosu mužů malajské populace je shodný s mediánem šířky nosu mužů současné malajské populace. Hi : Medián šířky nosu mužů malajské populace není shodný s mediánem šířky nosu mužů současné malajské populace. • matematická formulace nulové a alternativní hypotézy Hq : x = xq, kde xq = 26.90 Hi : x, ^ xq, kde xq = 26.90 (levostranná alternativa) 2. Volba hladiny významnosti • Hladinu významnosti podle zadání a = 0.01. 3. Testování kritickým oborem • Testovací statistika Nejprve vypočítáme vektor rozdílů naměřených hodnot Xi a konstanty xq, tj. Xi —xq, následně stanovíme počet kladných rozdílů Se a počet nenulových rozdílů m. n 73 SE = J2I(-X* ~ žo > 0) = Y^I(Xi ~ 26-90 > 0) = 33. i=l i=l 122 123 124 125 xO <- 26.90 I <- (nose.BM > xO) SE <- sum(I) # 33 m <- sum(nose.BM - xO != 0) #73 Všech 73 rozdílů Xi — xq, i = 1,..., 73, je nenulových, tedy m = 73. Počet kladných rozdílů Se = 33. Nyní vypočítáme testovací statistiku asymptotické varianty znaménkového testu podle vzorce 9.2. S a Se 2 33-? 0 V18.25 -3.5 = -0.819288 = -0.8193 33 - 36.5 4.272002 126 SA <- (SE - m / 2) / sqrt (m / 4) # -0.819288 18 • Kritický obor W = ( (-00 ; ua/2) U (u1_a/2 ; 00) (-00 ; «0.01/2) U («1-0.01/2 5 00) (-00 ; M0.005) U («0.99 5 00) (-00; -2.575829) U (2.575829; 00) 127 alpha <- 0.01 128 qnorm(alpha / 2) # -2.575829 129 qnormCl - alpha / 2) # 2.575829 • Závěr testování Protože realizace testovací statistiky sa = —0.819288 nenáleží do kritického oboru, tj. sa <ýz W, Hq nezamítáme na hladině významnosti a = 0.01. 4. Testování intervalem spolehlivosti • Interval spolehlivosti Pomocí příkazu sort() nejprve vzestupně seřadíme naměřené hodnoty Xi, i = 1,..., 73. Hranice 99% intervalu spolehlivosti budou potom (Ci_a/2)-tá a (n + 1 — Ci_a/2)-t& hodnota v seřazeném vektoru naměřených hodnot, kde 130 nose.BMs <- sort(nose.BM) 131 CI <- round(n / 2 - qnorm(l - alpha / 2) * sqrt(n / 4)) # 29 132 C2 <- round (n + 1 - CI) #53 133 nose.BMs [CI] # 25 134 nose.BMs [C2] # 27 • Závěr testování Protože xq = 26.90 náleží do 99% empirického pravostranného intervalu spolehlivosti, tj. xq = 26.90 G IS, Hq nezamítáme na hladině významnosti a = 0.01. = 36.5 — «o.995^18.25 = 36.5 - 2.575829 x 4.272002 = 25.49605 = 25 Interval spolehlivosti má potom následující tvar (25; 27) 19 5. Testování p-hodnotou • p-hodnota p-hodnota = 2min{Pr(5A < s a) , Pr( sa)} = 2min{Pr(5A < -0.819288), Pr(SA > -0.819288)} = 2min{Pr(5A < -0.819288), 1 - Pr(SA < -0.819288) = 2 min{0.206311, 0.7936889} = 2 x 0.2063111 = 0.4126221 = 0.4126 135 p.val <- 2 * min (pnorm(SA), 1 - pnorm(SA)) # 0.4126221 • Závěr testování Protože p-hodnota = 0.4126 je větší než a = 0.01, Hq nezamítáme na hladině významnosti a = 0.01. 6. Interpretace výsledků Za základě všech tří typů testování nezamítáme nulovou hypotézu na hladině významnosti a = 0.01. Mezi šířkou nosu mužů starověké a současné malajské populace není statisticky významný rozdíl. 7. Grafická vizualizace výsledku testování Grafické porovnání šířky nosu obou malajských populací provedeme krabicovým diagramem (viz obrázek 8). 32 - 30 28 - 26 - 24 22 - - staroveká p. O současna p. Obrázek 8: Krabicový diagram šířky nosu mužů malajské populace 20 Příklad 9.5. Znaménkový jednovýběrový asymptotický test Mějme datový soubor 18-more-samples-variances-clavicle.txt a proměnnou cla.L popisující největší délku klíční kosti z pravé strany v mm (viz sekce ??). Dále máme k dispozici údaje o největší délce klíční kosti z pravé strany mužů z populace severní Indie (mj = 148.0 mm, sj = 8.60 mm, rij = 260). Na hladině významnosti a = 0.05 zjistěte, zda je délka klíční kosti z pravé strany u mužů indické populace z Varanasi menší než u mužů ze severní Indie. Řešení příkladu 9.5 Načteme datový soubor a vybereme z datové tabulky naměřené délky klíční kosti (cla.L) mužů populace z Varanasi (pop == 'ind2'). Nakonec z vektoru naměřených hodnot odstraníme chybějící údaje a zjistíme rozsah náhodného výběru. 136 data <- read.delim('00-Data//18-more-samples-variances-clavicle.txt') 137 # head(data) 138 cla.LV <- data[data$pop == 'ind2 ' , 'cla.L'] 139 cla.LV <- na.omit(cla.LV) 140 n <- length(cla.LV) # 81 Datový soubor obsahuje údaje o největší délce klíční kosti z pravé strany u 81 mužů indické populace z Varanasi. V příkladu se zaměříme na porovnání délek klíčních kostí dvou populací, přičemž u jedné populace (indická po-pulcase z Varanasi) máme k dispozici naměřené hodnoty. Na základě těchto hodnot můžeme zjistit, zda náhodná veličina X popisující největší délku klíční kosti u mužů indické populace z Varanasi pochází z normálního rozdělení, tj. X ~ N(p, a2), kde skutečný rozptyl a2 neznáme. Druhá populace (populace ze severní Indie) je reprezentována pouze hodnotou aritmetického průměru (mj = 148.00mm) a směrodatnou odchylkou (sj = 8.60mm). Řešení příkladu tedy vede na situaci, kdy střední hodnotu jednoho náhodného výběru porovnáváme s konkrétním číslem, tedy na jednovýběrový test o střední hodnotě /i při neznámém rozptylu a2 (viz kapitola ??). Nejprve je však potřeba ověřit, zda náhodný výběr délek klíčních kostí u mužů indické populace z Varanasi pochází z normálního rozdělení. Tomuto předpokladu jsme se však už věnovali v rámci příkladu ?? v sekci ??, kde jsme zjistili, že náhodný výběr délek klíčních kostí u mužů indické populace z Varanasi z normálního rozdělení nepochází (pro připomenutí viz obrázek 9). 127 139 151 163 -2 -1 0 1 2 teoreticky kvantil délka klicni kosti z pravé strany (v mm) Lillieforsuv test: p-hodnota =0.0274 Obrázek 9: Histogram a kvantilový diagram největší délky klíční kosti z pravé strany u mužů z populace z Varanasi Protože náhodný výběr nepochází z normálního rozdělení, nemůžeme hypotézu ze zadání otestovat pomocí parametrického testu o střední hodnotě jj,. K testování hypotézy použijeme neparametrický znaménkový jednovýběrový asymptotický test. Naším úkolem je zjistit, zdaje délka klíční kosti z pravé strany u mužů indické populace z Varanasi 21 menší než u mužů ze severní Indie. Tato věta je zněním alternativní hypotézy, neboť v jejím znění není zmínka o (nulové) hypotéze, ani o rovnosti. Zbývá dodefinovat nulovou hypotézu tak, aby byla doplňkem k hypotéze alternativní. 1. Stanovení hypotéz • slovní formulace nulové a alternativní hypotézy H0 : Medián největší délky klíční kosti z pravé strany u mužů indické populace z Varanasi je větší nebo rovný mediánu největší délky klíční kosti z pravé strany u mužů ze severní Indie. Hi : Medián největší délky klíční kosti z pravé strany u mužů indické populace z Varanasi je menší než medián největší délky klíční kosti z pravé strany u mužů ze severní Indie. • matematická formulace nulové a alternativní hypotézy H0 : x > x0, kde x0 = 148.00 Hi : x, < xq, kde xq = 148.00 (oboustranná alternativa) 2. Volba hladiny významnosti • Hladinu významnosti podle zadání a = 0.05. 3. Testování kritickým oborem • Testovací statistika Nejprve vypočítáme vektor rozdílů naměřených hodnot X{ a konstanty ž0, tj. X{ —x0, následně stanovíme počet kladných rozdílů Se a počet nenulových rozdílů m. SE = Y1 J(X» ~ žo > 0) = ]T I(Xi - 148 > 0) = 15. 141 xO <- 148.00 142 I <- (cla.LV > xO) 143 SE <- sum(I) # 15 144 m <- sum(cla.LV - xO != 0) #79 Z celkového počtu 81 rozdílů Xí—xq, i = 1,..., 81, je 79 rozdílů nenulových, tedy m = 79. Počet kladných rozdílů Se = 15. Nyní již můžeme vypočítat testovací statistiku asymptotické varianty znaménkového testu podle vzorce 9.2. SE- f 15-f 0 .- -24.5 Sa = = —ňr = T^7^> = 4^44097 = "5-512931 = "5-5129 145 SA <- (SE - m / 2) / sqrt (m / 4) # -5.512931 • Kritický obor W = (—oo ; ua) = (-oo; «0.05) = (-oo; -1.644854) 22 146 alpha <- 0.05 147 qnorm(alpha) # -1.644854 • Závěr testování Protože realizace testovací statistiky sa = —5.512931 náleží do kritického oboru, tj. sa € W, Hq zamítáme na hladině významnosti a = 0.05. 4. Testování intervalem spolehlivosti • Interval spolehlivosti Pomocí příkazu sort() nejprve vzestupně seřadíme naměřené hodnoty Xi, i = 1,..., 79. Hranice 95% intervalu spolehlivosti budou potom mínus nekonečno a (n + 1 — Ci_a)-tá hodnota v seřazeném vektoru naměřených hodnot, kde 148 cla.LVs <- sort(cla.LV) 149 Cl <- round(n / 2 - qnorm(l - alpha) * sqrt(n / 4) ) #33 150 C2 <- round(n + 1 - Cl) # 49 151 cla.LVs [49] # 142 • Závěr testování Protože xq = 148.00 nenáleží do 95% empirického oboustranného intervalu spolehlivosti, tj. xq = 148.00 ^ IS, Hq zamítáme na hladině významnosti a = 0.05. 5. Testování p-hodnotou = 40.5 - «0.95^20.25 = 40.5 - 1.644854 x 4.5 = 33.09816 = 33 Interval spolehlivosti má potom následující tvar (-00; 27) • p-hodnota p-hodnota = Pr(SA < sA) = Pr(SA < -5.512931) = 1.764536"8 = 1.7645": 8 152 p.val <- pnorm(SA) # 1.764536e-08 23 • Závěr testování Protože p-hodnota = 1.7645-8 je menší než a = 0.05, Hq zamítáme na hladině významnosti a = 0.05. 6. Interpretace výsledků Za základě všech tří typů testování zamítáme nulovou hypotézu na hladině významnosti a = 0.05. Největší délka klíční kosti na pravé straně u mužů indické populace z Varanasi je statisticky významně menší než největší délka klíční kosti na pravé straně u mužů populace ze severní Indie. 7. Grafická vizualizace výsledku testování Rozdíl v délce klíčních kostí zaznamenaný mezi oběma populacemi vizualizujeme pomocí krabicového diagramu (viz obrázek 10). 170 - 160 - - indická p. z Varanasi O indická p. ze severu O o 150 ca 140 H 130 - Obrázek 10: Krabicový diagram délky klíční kosti z pravé strany u mužů indické populace z Varanasi ★ 24 9.3 Znaménkový párový test Nechť (Xi, Yi)T ... (Xn, Yn)T je náhodný výběr z libovolného (ne nutně normálního) dvourozměrného rozdělení Nechť dále Z\,..., Zn, n > 2 je náhodný výběr rozdílů X — Y, tj. Z = (Z\,..., Zn)T, kde Z{ = X{ — Y{, i = 1,..., n, a nechť tento náhodný výběr pochází z libovolného spojitého rozdělení. Konečně, nechť Žq je konstanta. Na hladině významnosti a testujeme jednu z následujících tří hypotéz oproti příslušné alternativní hypotéze. kde ž je medián rozdílů Z\,..., Zn a Žq je konstanta, jejíž hodnotu nejčastěji volíme jako Žq = 0. Tato volba odpovídá hypotéze, že rozdíl mezi mediány náhodných veličin X a Y neexistuje (resp. hypotéze, že medián náhodné veličiny X je menší, resp. větší, než medián náhodné veličiny Y). Vzhledem k tomu, že jde finálně o situaci, kdy medián ž porovnáváme s konstantou testujeme hypotézy o rozdílu mediánů X — Y pomocí exaktní nebo asymptotické varianty znaménkového jednovýběrového testu, analogicky jako je uvedeno v sekcích 9.1 a 9.2. Výše popsaný test, v rámci kterého převádíme problém porovnávání mediánů dvou náhodných veličin X a Y na problém srovnávání mediánu jejich rozdílů Z s konstantou Žq = 0 a následně jej řešíme pomocí exaktní resp. asymptotické varianty znaménkového jednovýběrového testu, nazýváme znaménkový párový test. H01 : ž = ž0 H02 : ž < ž0 H03 : z > ž0 oproti oproti oproti Hu : ž ^ Žq (oboustranná alt.) Hi2 : ž > Žq (pravostranná alt.) His : ž < Žq (levostranná alt.) 25 Příklad 9.6. Znaménkový párový exaktní test Mějme datový soubor 21-goldman-measures.csv obsahující údaje o největší výšce kloubní jamky na pravé straně (acetab.HR) a na levé straně (acetab.HL) u skeletů ze starověké egyptské populace (detaily viz sekce ??). Na hladině významnosti a = 0.05 zjistěte, zda existuje rozdíl mezi výškou kloubní jamky z pravé a levé strany u skeletů mužského pohlaví. Řešení příkladu 9.6 Pomocí příkazu read.delim() načteme datový soubor a pomocí operátoru [] vybereme z datové tabulky údaje o výšce kloubní jamky z pravé strany (acetab.HR), resp. z levé strany (acetabl.HL) u mužů (sex == 'm') ze starověké egyptské populace (pop == 'Dynastie Egyptian, El Hesa'). Z vektoru naměřených hodnot odstraníme pomocí funkce na.omit() chybějící údaje a zjistíme rozsah náhodného výběru (dim()). 153 data <- read.delim('00-Data//21-goldman-measures.csv', sep = ';') 154 # head(data) 155 data.HRL <- data[data$pop == 'Dynastie Egyptian, El Hesa' k data$sex == 'm', 156 c('acetab.HL ' , 'acetab.HR')] 157 data.HRL <- na.omit(data.HRL) 158 dim(data.HRL) # 18x2 Datový soubor obsahuje kompletní údaje o výšce kloubní jamky z pravé i levé strany u 18 mužů ze starověké egyptské populace. Naším úkolem ze zadání je porovnat naměřené hodnoty na pravé a levé straně. Jde tedy o měření stejného znaku (výška kloubní jamky) sledovaného na stejných subjektech (muži), proto použijeme pro tuto situaci párový test. Prvním krokem k použití tohoto testu je vytvoření rozdílů hodnot naměřených na pravé a levé straně. 159 acetab.HR <- data.HRL$acetab.HR 160 acetab.HL <- data.HRL$acetab.HL 161 acetab.HRL <- acetab.HR - acetab.HL Ve druhém kroku je potřeba ověřit normální rozdělení těchto rozdílů. Vzhledem k rozsahu náhodného výběru (n = 18 < 30) otestujeme normalitu rozdílů Shapiro-Wilkovým testem (a = 0.05). Pro potřeby vykreslení histogramu rozdělíme naměřené hodnoty do pěti ekvidistatních intervalů s šířkou 1.15 mm prostřednictvím stanovených hranic —1, 0.15,... ,05.75 (viz obrázek 11). Obrázek 11: Histogram a diagram rozdílů mezi největší výškou kloubní jamky u mužů na pravé a levé straně 26 Protože p-hodnota = 0.0121 je menší než 0.05, nulovou hypotézu o normalitě rozdílů zamítáme na hladině významnosti a = 0.05. Z histogramu na obrázku 11 vidíme, že naměřené hodnoty jsou vyšikmené doleva s odlehlým pozorováním na pravém chvostu. Histogram naměřených hodnot navíc nevykazuje symetrii ani postupné snižování počtu hodnot na levé straně, jak bychom u normálního rozdělení očekávali. Náhodný výběr rozdílů mezi výškou kloubní jamky na pravé a levé straně nepochází z normálního rozdělení. Protože náhodný výběr rozdílů nepochází z normálního rozdělení, nemůžeme hypotézu ze zadání otestovat pomocí parametrického párového testu uvedeného v sekci ??. K testování hypotézy musíme použít neparametrickou alternativu párového testu. S ohledem na nízký rozsah náhodného výběru použijeme exaktní variantu znaménkového párového testu. Naším úkolem je zjistit, existuje rozdíl mezi výškou kloubní jamky z pravé a levé strany u skeletů mužského pohlaví. Tato věta bude součástí alternativní hypotézy, neboť rozdíl implikuje nerovnost a nerovnost je vždy součástí alternativní hypotézy. Nulovou hypotézu stanovíme jako doplněk k tomuto tvrzení. Testování provedeme v posloupnosti sedmi kroků. 1. Stanovení hypotéz • slovní formulace nulové a alternativní hypotézy Hq : Medián rozdílů výšky kloubní jamky na pravé a levé straně u mužů je rovný nule. H\ : Medián rozdílů výšky kloubní jamky na pravé a levé straně u mužů není rovný nule. • matematická formulace nulové a alternativní hypotézy H o : ž = Žq, kde Žq = 0 Hi : ž Žq, kde Žq = 0 (oboustranná alternativa) 2. Volba hladiny významnosti • Hladinu významnosti volíme v souladu se zadáním jako a = 0.05. 3. Testování kritickým oborem • Testovací statistika V úvodu příkladu jsme vytvořili vektor Z i jako vektor rozdílů výšky kloubní jamky na pravé straně Xi a výšky kloubní jamky na levé straně Yi, tj. Zi = Xi — Yi, i = 1,..., 18. Nyní je potřeba stanovit rozdíly Zi — Žq. Protože ale Žq = 0, jsou rozdíly Zi — Žq rovny vektoru rozdílů Zi, tj. Zi — Žq = Zi — 0 = Zi. Původní hodnoty výšky kloubní jamky na pravé straně, resp. na levé straně, rozdíly mezi pravou a levou stranou, rozdíly Zi — Žq a znaménko posledních uvedených rozdílů jsou pro názornost uvedeny v tabulce 7. Tabulka 7: Naměřené hodnoty Xi, Yi, rozdíly Zi = Xi — Yi, rozdíly Z i — žq a znaménka těchto rozdílů měření 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 Xi 50.84 52.36 47.80 51.75 46.06 46.64 50.54 47.03 49.02 47.96 48.25 47.78 49.86 49.56 55.41 52.82 52.47 46.02 Yi 50.97 51.80 48.75 51.40 43.96 44.86 49.64 48.02 49.45 48.52 48.64 48.14 48.13 48.54 50.66 52.10 52.82 45.35 Zi -0.13 0.56 -0.95 0.35 2.10 1.78 0.90 -0.99 -0.43 -0.56 -0.39 -0.36 1.73 1.02 4.75 0.72 -0.35 0.67 Zi - ž0 -0.13 0.56 -0.95 0.35 2.10 1.78 0.90 -0.99 -0.43 -0.56 -0.39 -0.36 1.73 1.02 4.75 0.72 -0.35 0.61 +/- - + - + + + + - - - - - + + + + - + Z tabulky 7 vidíme, že celkem deset rozdílů Zi — žq, i = 1,..., 18 je kladných. Hodnota testovací statistiky Se, která je definovaná jako počet kladných rozdílů, bude tedy rovná 10. n 18 SE = Y,1^ - žo > 0) = ^/(Z2 - 0 > 0) = 10. i=l i=l • Kritický obor Z tabulky 7 dále vidíme, že žádný z celkového počtu n = 18 rozdílů není nulový, a tedy počet nenulových 27 162 zO <- O 163 I <- (acetab.HRL > zO) 164 tab <- data.frame(rbind("Xi" = acetab.HR, "Yi" = acetab.HL, 165 "Zi" = acetab.HRL, "Zi-zO" = acetab.HRL - zO, "+/-" = I)) 166 SE <- sum(I) # 10 rozdílů m = 18. Kritický obor má potom tvar W = (-00 ; 6m,1/2(Q!/2) - 1) U <6m,i/2(l - a/2); oo) = (-00 ; 6i8,i/2(0.05/2) - l) U <618,1/2(1 - 0.05/2); oo) = (-oo ; 6i8,i/2(0.025) - l) U (&iM/2(0.975); oo) = (-oo; 5 - 1) U (13; oo) = (-oo; 4) U (13; oo) 167 m <- sum(acetab.HRL - zO != 0) # 18 168 alpha <- 0.05 169 ql <- qbinom(alpha / 2 , m, 1/2) -1#4 170 q2 <- qbinomd - alpha / 2, m, 1 / 2) # 13 • Závěr testování Protože realizace testovací statistiky s e = 10 nenáleží do kritického oboru, tj. s e W, Hq nezamítáme na hladině významnosti a = 0.05. 4. Testování intervalem spolehlivosti • Interval spolehlivosti Abychom mohli stanovit hranice 95 % intervalu spolehlivosti, musíme nejprve naměřené rozdíly mezi výškami kloubní jamky na pravé a levé straně vzestupně seřadit. To provedeme příkazem sort(). Seřazené hodnoty jsou k dispozici v tabulce 8. Tabulka 8: Seřazené rozdíly mezi výškami kloubní jamky na pravé a levé straně pořadí || 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 seřazené Z; || -0.99 ^095 ^Ôľ56 ^043 1TM ľÔ3fš ^035 ^Ôľl3 Ô35 Ô56 067^ Ô72 Ô9Ô TM L73 L78 210 4.75 Hranice intervalu spolehlivosti potom tvoří ty hodnoty, které se v seřazeném vektoru hodnot nachází na (&IJil/2(a/2))-té pozici a na (&nii/2(l — a/2) + l)-té pozici. (d,h) = (z(iV1 se)} = 2min{Pr(SE < 10), Pľ(SE > 10)} = 2min{Pr(SE < 10), 1 - Pr(SE < 9)} = 2min{0.7596588, 0.4072647} = 2 x 0.4072647 = 0.8145294 = 0.8145 176 2 * min (pbinom(SE, m, 1/2), 1 - pbinom(SE - 1, m, 1 / 2)) # 0.8145294 • Závěr testování Protože p-hodnota 0.8145 je větší než a = 0.05, Hq nezamítáme na hladině významnosti a = 0.05. 6. Interpretace výsledků Za základě všech tří typů testování nezamítáme nulovou hypotézu na hladině významnosti a = 0.05. Mezi největší výškou kloubní jamky na pravé a levé straně u mužů starověké egyptské populace neexistuje statisticky významný rozdíl. 7. Grafická vizualizace výsledku testování Porovnání měření na pravé a levé straně zobrazíme nejlépe pomocí krabicového diagramu. Analogicky jako v sekci ?? můžeme buď sestrojit krabicový diagram rozdílů mezi výškami kloubní jamky na pravé a levé straně a porovnat je s konstantou Žq = 0 (viz obrázek 12 vlevo) nebo sestrojit krabicový diagram zvlášť pro výšky kloubní jamky na pravé straně a zvlášť pro výšky kloubní jamky na levé straně a porovnat tyto diagramy navzájem (viz obrázek 12 vpravo). o 4 - 2 - 0 - -2 - o O 56 - 54 52 - :r 50 - 48 - 46 - 44 pravá strana leva strana Obrázek 12: Krabicový diagram rozdílů mezi největší výškou kloubní jamky u mužů na pravé a levé straně 29 Poznámka: Znaménkový párový test můžeme provést pomocí funkce SIGN.test() implementované v knihovně BSDA. Vstupními parametry budou vektor naměřených hodnot výšek kloubní jamky na pravé straně (acetab.HR), vektor naměřených hodnot výšek kloubní jamky na levé straně (acetab.HL), argument paired = T určující, že oba vektory považujeme za párová pozorování, hodnota hladiny významnosti a zadaná prostřednictvím koeficientu spolehlivosti 1 — a nastavením hodnoty argumentu conf.level = 0.95 a typ zvolené alternativní hypotézy (oboustranná) zadaný pomocí argumentu alternativě = 'two.sideď. Součástí výstupu je hodnota mediánu rozdílů naměřených výšek na 177 BSDA::SIGN.test(acetab.HR, acetab.HL, paired = T, 178 conf.level = 0.95, alternativě = 'two.sided') Dependent - samples Sign-Test data: acetab.HR and acetab.HL S = 10, p-value = 0.8145 alternative hypothesis: true median difference is not equal to 0 95 percent confidence interval: -0.3812269 0.9849076 sample estimates: median of x-y 0.455 Achieved and Interpolated Confidence Intervals: Conf.Level L.E.pt U.E.pt Lower Achieved CI 0.9037 -0.3600 0.9000 Interpolated CI 0.9500 -0.3812 0.9849 Upper Achieved CI 0.9691 -0.3900 1.0200 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 pravé a levé straně medián of x-y = 0.455, hodnota testovací statistiky S = 10, interpolované hranice 95% empirického oboustranného intervalu spolehlivosti -0.3812269 a 0.9849076 pro rozdíl naměřených výšek na pravé a levé straně, které jsou mírně přesnější, než námi stanovené hranice intervalu spolehlivosti (zpřesnění hranic bylo provedeno procesem nazývaným interpolace), a p-hodnota p-value = 0.8145. Jediné, co musíme stanovit zvlášť, jsou dolní a horní hranice kritického oboru. Druhou možností provedení párového testu je opět pomocí funkce SIGN.test(), kde vstupními parametry budou vektor rozdílů naměřených hodnot na pravé a levé straně (acetab.HRL), argument md = 0 určující, že rozdíly porovnáváme s konstantou Žq = 0, hodnota hladiny významnosti a zadaná prostřednictvím koeficientu spolehlivosti 1 — a (conf.level = 0.95) a typ zvolené alternativní hypotézy (alternativě = 'two.sided'). Výstup tohoto příkazu je 197 BSDA::SIGN.test(acetab.HRL, md = 0, conf.level = 0.95, alternativě = 'two.sided') totožný s výše uvedeným výstupem. Záleží tedy na nás, jakou syntaxi k zadání exaktního znaménkového párového testu použijeme. 30 One-sample Sign-Test data: acetab.HRL s = 10, p-value = 0.8145 alternative hypothesis: true median is not equal to 0 95 percent confidence interval: -0.3812269 0.9849076 sample estimates: median of x 0.455 Achieved and Interpolated Confidence Intervals: Lower Achieved CI Interpolated CI Upper Achieved CI Conf.Level 0.9037 0.9500 0.9691 L.E.pt U.E.pt -0.3600 0.9000 -0.3812 0.9849 -0.3900 1.0200 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 31 Příklad 9.7. Znaménkový párový exaktní test Mějme datový soubor 21-goldman-measures.csv obsahující údaje o délkách lýtkových kostí v mm na pravé straně (tibia.LR) a na levé straně (tibia.LL) u skeletů obyvatel útesových obydlí (Cliff Dwellings) v Utahu (detaily viz sekce ??). Na hladině významnosti a = 0.05 testujte nulovou hypotézu, že délka lýtkové kosti z pravé strany u skeletů mužského pohlaví je menší než délka lýtkové kosti z levé strany. Řešení příkladu 9.7 Načteme datový soubor a pomocí operátoru [] vybereme z datové tabulky údaje o délce lýtkové kosti z pravé strany (tibia.LR), resp. z levé strany (tibia.LL) u mužů (sex == 'm') z populace obyvatel útesových obydlí (pop == 'Cliff Dweller'). Z vektoru naměřených hodnot následně odstraníme chybějící hodnoty a zjistíme rozsah náhodného výběru. 216 data <- read.delim('00-Data//21-goldman-measures.csv' , sep = ';') 217 # head(data) 218 data.TRL <- data[data$pop == 'Cliff Dweller' k data$sex == 'm', 219 c('tibia.LR', 'tibia.LL')] 220 data.TRL <- na.omit(data.TRL) 221 dim(data.TRL) # 20 Datový soubor obsahuje kompletní údaje o délce lýtkových kostí u 26 mužů z populace obyvatel útesových obydlí na území Utahu. Naším úkolem ze zadání je porovnat naměřené hodnoty na pravé a levé straně. Jde tedy o měření stejného znaku (délka lýtkové kosti) sledovaného na stejných subjektech (muži). Proto použijeme pro tuto situaci párový test. Prvním krokem k použití párového testu je vytvoření rozdílů hodnot naměřených na pravé a levé straně. 222 tibia.LR <- data.TRL$tibia.LR 223 tibia.LL <- data.TRL$tibia.LL 224 tibia.LRL <- tibia.LR - tibia.LL V druhém kroku je potřeba ověřit normální charakter těchto rozdílů. Tímto předpokladem jsme se zabývali již v rámci příkladu ?? v sekci ??, kde jsme zjistili, že náhodný výběr rozdílů mezi délkami lýtkových kostí na pravé a levé straně nepochází z normálního rozdělení (p-hodnota < 0.0001). Pro připomenutí viz obrázek 13. Obrázek 13: Histogram a diagram rozdílů mezi délkou lýtkové kosti na levé a pravé straně u mužů z populace obyvatel skalních obydlí v Utahu Z histogramu na obrázku 13 vidíme, že rozdíly mezi pravou a levou stranou jsou výrazně vyšikmené doleva s od- 32 lehlým pozorováním na pravé straně. Z kvantilového diagramu potom vidíme, že toho odlehlé pozorování je extrémně vzdálené od referenční přímky. Z příkladu ?? již víme, že toho pozorování narušuje fatálním způsobem normalitu náhodného výběru. V příkladu ?? jsme tuto situaci vyřešili odstraněním odlehlého pozorování a následným provedením parametrického párového testu. Zde budeme pracovat s předpokladem, že naměřené údaje jsou správné a pozorování v souboru ponecháme. Protože však náhodný výběr rozdílů nepochází z normálního rozdělení, nemůžeme hypotézu ze zadání otestovat pomocí parametrického testu, analogicky jako v příkladu ??. K otestování použijeme neparametrický exaktní znaménkový párový test. V zadání příkladu máme explicitně uvedeno, že máme testovat nulovou hypotézu o menší délce lýtkové kosti z pravé strany vzhledem k délce lýtkové kosti z levé strany. To odpovídá tvrzení, že rozdíly získané odečtením naměřených hodnot na levé straně od naměřených hodnot na pravé straně budou menší než nula. Zbývá dodefinovat alternativní hypotézu. 1. Stanovení hypotéz • slovní formulace nulové a alternativní hypotézy H0 : Medián rozdílů délky ramenní kosti na pravé a levé straně u mužů je menší nebo roven nule. Hi : Medián rozdílů délky ramenní kosti na pravé a levé straně u mužů je větší než nula. • matematická formulace nulové a alternativní hypotézy H o : ž < Žq, kde Žq = 0 H\ : ž > Žq, kde Žq = 0 (pravostranná alternativa) 2. Volba hladiny významnosti • Hladinu významnosti volíme v souladu se zadáním jako a = 0.05. 3. Testování kritickým oborem • Testovací statistika V úvodu příkladu jsme vytvořili vektor tibia.LRL (Zi, i = 1, ...,20) jako vektor rozdílů mezi délkami lýtkových kostí na pravé a levé straně. Nyní je potřeba stanovit rozdíly Zi — Žq. Protože ale Žq = 0, jsou rozdíly Z{ — ž0 rovny původnímu vektoru rozdílů Z{. Hodnoty délky lýtkové kosti na pravé straně (Xi), resp. na levé straně (Yí), rozdíly Zi, rozdíly Zi — Žq a znaménko posledních uvedených rozdílů jsou pro názornost uvedeny v tabulce 9. Tabulka 9: Naměřené hodnoty Xi, Yí, rozdíly Zi = Xi — Yí, rozdíly Zi — žq a znaménka těchto rozdílů měření 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Xi 376 370 346 388 380 335 372.5 388 364.0 361 370 396.5 373.0 378 368 370 374 356.0 355 337 Yi 379 370 350 386 380 334 372.0 388 365.5 362 374 396.0 379.5 370 365 372 376 364.5 350 307 Zi -3 0 -4 2 0 1 0.5 0 -1.5 -1 -4 0.5 -6.5 8 3 -2 -2 -8.5 5 30 Zi-ŽQ -3 0 -4 2 0 1 0.5 0 -1.5 -1 -4 0.5 -6.5 8 3 -2 -2 -8.5 5 30 +/- - - - + - + + - - - - + - + + - - - + + Z tabulky 9 vidíme, že celkem osm rozdílů Zi — žq, i = 1,..., 20, je kladných. Hodnota testovací statistiky Se bude tedy rovná 8. n 20 SE = Y,1^ " žo > 0) = ^/(Z2 - 0 > 0) = 8. i=l i=l 225 z0 <- 0 226 I <- (tibia.LRL > z0) 227 tab <- data.frame(rbind("Xi" = tibia.LR, "Yi" = tibia.LL, 228 "Zi" = tibia.LRL, "Zi-zO" = tibia.LRL - zO , "+/-" = I)) 229 names(tab) <- 1 : 20 230 SE <- sum(I) # 8 33 • Kritický obor Z tabulky 9 dále vidíme, že z celkového počtu n = 20 rozdílů jsou tři rozdíly nulové, a tedy počet nenulových rozdílů m = 17. Kritický obor bude potom tvořen hodnotou kvantilu &mii/2(l — a) a nekonečnem. W = (&m,i/2(1 - a); oo) = <6i7,i/2(l-0.05);oo) = <6i7,i/2(0.95); oo) = (12; oo) 231 n <- length(tibia.LRL) 232 m <- sum(tibia.LRL - zO != 0) # 17 233 alpha <- 0.05 234 q <- qbinomd - alpha, m, 1/2) #12 • Závěr testování Protože realizace testovací statistiky se = 8 nenáleží do kritického oboru, tj. se £ W, Hq nezamítáme na hladině významnosti a = 0.05. 4. Testování intervalem spolehlivosti • Interval spolehlivosti Abychom mohli stanovit hranice 95 % intervalu spolehlivosti, musíme nejprve rozdíly mezi naměřenými délkami pravé a levé strany vzestupně seřadit. Seřazené rozdíly jsou k dispozici v tabulce 10. Tabulka 10: Seřazené rozdíly mezi délkami lýtkových kostí na pravé a levé straně pořadí || 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 seřazené Z, || -8.5 ^05 Tči 3~0 1ŠÁ) :2l) ^275 T5 Tč) O0 O0 O0 05 05 LO 275 ŠD 5Á) KO 30ČT Hranice intervalu spolehlivosti potom tvoří hodnota, která se v seřazeném vektoru hodnot nachází na (6n,i/2(a))-té pozici, a nekonečno. (d,h) = (V^.1/^)) ; oo) = (Z(bl7,l/2(0.05)) . ^ = (zW ; oo) = (-2; oo) 235 qbinom(alpha, n, 1/2) #6 236 tibia.LRLs [6] # -2 • Závěr testování Protože Žq = 0 náleží do 95% empirického levostranného intervalu spolehlivosti, tj. Žq = 0 € IS, Hq nezamítáme na hladině významnosti a = 0.05. 5. Testování p-hodnotou • p-hodnota p-hodnota = Pr(5£; > sE) = Pr(SE > 8) = 1 - Pr(5£; < 8) = 1 - Pr^ < 7) = 0.6854706 = 0.6855 34 237 1 - pbinom(SE - 1, m, 1/2) # 0.6854706 • Závěr testování Protože p-hodnota 0.6855 je větší než a = 0.05, Hq nezamítáme na hladině významnosti a = 0.05. Interpretace výsledků Na základě všech tří typů testování na hladině významnosti a = 0.05 nezamítáme nulovou hypotézu, že délka lýtkové kosti z pravé strany u skeletů mužského pohlaví je menší než délka lýtkové kosti z levé strany. Uvědomme si, že zamítnutí nulové hypotézy by vedlo k závěru, že délka lýtkové kosti z pravé strany je statisticky významně větší než délka lýtkové kosti z levé strany. Nezamítnutí nulové hypotézy tedy vede k závěru, že délka lýtkové kosti z pravé strany není statisticky významně větší než délka lýtkové kosti z levé strany. Grafická vizualizace výsledku testování Porovnání měření na pravé a levé straně vizalizujeme pomocí krabicového diagramu, přičemž opět se můžeme rozhodnout, zda dáme přednost diagramu porovnávajícímu rozdíly mezi pravou a levou stranou s konstantou Žq = 0 (obrázek 14 vlevo), nebo diagramu porovnávajícímu naměřené hodnoty na pravé straně s hodnotami naměřenými na levé straně (obrázek 14 vpravo). o ca ti o •a 30 - 20 * 10 - o - -10 - > O 400 380 - 360 - 340 ■° 320 - pravá strana leva strana Obrázek 14: Krabicový diagram rozdílů mezi délkou lýtkové kosti na levé a pravé straně u mužů z populace obyvatel skalních obydlí v Utahu Poznámka: Znaménkový exaktní párový test provedeme pomocí funkce SIGN.test() z knihovny BSDA. Vstupními parametry budou vektor naměřených hodnot délek lýtkové kosti na pravé straně (tibia.LR), vektor naměřených hodnot délek lýtkové kosti na levé straně (tibia.LL), specifikace párového testu (paired = T), hodnota hladiny významnosti a zadaná prostřednictvím koeficientu spolehlivosti 1 — a (conf.level = 0.95) a pravostranný typ alternativní hypotézy (alternativě = 'greater'). Součástí výstupu je hodnota mediánu rozdílů naměřených výšek na pravé a levé straně 238 BSDA::SIGN.test(tibia.LR, tibia.LL, paired = T, conf.level = 0.95, 239 alt = 'greater') # IS interpolovaný medián of x-y = 0, hodnota testovací statistiky s = 8, interpolované hranice 95% empirického levostranného intervalu spolehlivosti -2 a Inf pro rozdíl mezi délkami na pravé a levé straně a p-hodnota p-value = 0.6855. Jediné, co musíme stanovit zvlášť, je dolní hranice kritického oboru. 35 Příklad 9.8. Znaménkový párový exaktní test Mějme datový soubor 21-goldman-measures.csv obsahující údaje o délkách ramenních kostí v mm na pravé straně (humer.LR) a na levé straně (humer.LL) u skeletů z římského pohřebiště v Poundbury (detaily viz sekce ??). Na hladině významnosti a = 0.10 zjistěte, zda je u skeletů ženského pohlaví délka pažní kosti z levé strany menší než délka pažní kosti z pravé strany. Řešení příkladu 9.8 Načteme datový soubor a pomocí operátoru [] vybereme z datové tabulky údaje o délce pažní kosti z levé strany (humer.LL), resp. z pravé strany (humer.LR) u žen (sex == 'f) z římského pohřebiště v Poundbury (pop == 'Poundbury'). Z vektoru naměřených hodnot odstraníme chybějící hodnoty a zjistíme rozsah náhodného výběru. 240 data <- read.delim('00-Data//21-goldman-measures.csv' , sep = ';') 241 # head(data) 242 data.LLR <- data[data$pop == 'Poundbury' k data$sex == 'i', 243 c('numer.LL', 'humer.LR')] 244 data.LLR <- na.omit(data.LLR) 245 dim(data.LLR) # 26x2 Datový soubor obsahuje kompletní údaje o délce pažní kosti z levé i pravé strany u 26 skeletů žen z římského pohřebiště v Poundbury. Naším úkolem ze zadání je porovnat naměřené hodnoty na levé a pravé straně. Jde tedy o měření stejného znaku (délka pažní kosti) sledovaného na stejných subjektech (ženy), proto použijeme pro tuto situaci párový test. Nejprve tedy vytvoříme rozdíly hodnot délek pažních kostí naměřených na levé a pravé straně a následně ověříme normální rozdělení těchto rozdílů. K tomu využijeme Shapirův-Wilkovým testem (a = 0.05), kvantilový diagram a histogram, přičemž naměřené hodnoty rozdělíme do šesti ekvidistatních intervalů s šířkou 5 mm prostřednictvím stanovených hranic —20, —15,..., 10 (viz obrázek ??). 246 numer LL <- data.LLR$humer.LL 247 numer LR <- data.LLR$humer.LR 248 numer LLR < - numer.LL - numer.LR -17.5 -7.5 2.5 7.5 -2 -1 0 1 2 teoreticky kvantil rozdil délky na pravé a leve strane (v mm) Shapiruv test: p-hodnota = 0.0039 Obrázek 15: Histogram a diagram rozdílů mezi délkami těla pažní kosti na pravé a levé straně u skeletů žen z římského pohřebiště v Poundbury Protože p-hodnota = 0.0039 je menší než 0.05, nulovou hypotézu o normalitě rozdílů zamítáme na hladině významnosti 36 a = 0.05. Z histogramu na obrázku ?? vidíme, že naměřené hodnoty jsou výrazně špičatější s mírně prodlouženým pravým chvostem. Z kvantilového diagramu potom vidíme, že hodnoty na chvostech se zjevně odchylují od referenční přímky. Shapiro-Wilkův test vyhodnotil tyto nedostatky jako fatální pro normální charakter náhodného výběru. Náhodný výběr rozdílů mezi délkou pažní kosti na levé a pravé straně nepochází z normálního rozdělení. Protože náhodný výběr rozdílů nepochází z normálního rozdělení, nemůžeme k ověření otázky ze zadání použít parametrický párový test. Použijeme tedy jeho neparametrickou alternativu, a sice exaktní znaménkový párový test. Naším úkolem je zjistit, zda je u skeletů ženského pohlaví délka pažní kosti z levé strany menší než délka pažní kosti z pravé strany. Tato věta je zněním alternativní hypotézy, přičemž nulovou hypotézu dodefinujeme jako doplněk k tomuto tvrzení. 1. Stanovení hypotéz • slovní formulace nulové a alternativní hypotézy Hq : Medián rozdílů délky pažní kosti na levé a pravé straně u žen je větší nebo rovný nule. Hi : Medián rozdílů délky pažní kosti na levé a pravé straně u žen je menší než nula. • matematická formulace nulové a alternativní hypotézy H o : ž > Žq, kde Žq = 0 H\ : ž < Žq, kde Žq = 0 (levostranná alternativa) 2. Volba hladiny významnosti • Hladinu významnosti volíme v souladu se zadáním jako a = 0.10. 3. Testování kritickým oborem • Testovací statistika V úvodu příkladu jsme vytvořili vektor humer.LLR (Zi) jako vektor rozdílů mezi délkou pažní kosti na levé a pravé straně. Nyní je potřeba stanovit rozdíly Zi — Žq a určit počet kladných rozdílů Zi — Žq a počet nenulových rozdílů m. 249 zO <- 0 250 I <- (humer.LLR > zO) 251 tab <- data.frame(rbind("Xi" = humer.LL, "Yi" = humer.LR , 252 "Zi" = humer.LLR, "Zi-zO" = humer.LLR - zO , "+/-" = I)) 253 SE <- sum(I) # 1 254 n <- length(humer.LLR) # 26 255 m <- sum(humer.LLR - zO != 0) #26 Z celkového počtu 26 rozdílů není žádný rozdíl nulový, tj. m tedy 26. Dále pouze jeden rozdíl je kladný, n 26 SE = ^I{Zi - ž0 > 0) = 0 > 0) = 1. i=l i=l • Kritický obor W = ( = ( (-oo; bmA/2(a) - l) (-oo;626,1/2(0.10)-l) (-oo; 10 - 1) (-oo; 9) • Závěr testování Protože realizace testovací statistiky s e = 1 náleží do kritického oboru, tj. s e G W, Hq zamítáme na hladině významnosti a = 0.10. 37 256 alpha <- 0.10 257 q <- qbinom(alpha , m, 1 / 2) - 1 # 6 4. Testování intervalem spolehlivosti • Interval spolehlivosti Abychom mohli stanovit hranice 90 % intervalu spolehlivosti, musíme nejprve rozdíly mezi naměřenými délkami pažní kosti na levé a pravé straně vzestupně seřadit. Hranice intervalu spolehlivosti potom tvoří mínus nekonečno a hodnota, která se v seřazeném vektoru hodnot nachází na (&nii/2(l — a) + l)-té pozici. (d,h) = (-oo; = (-oo; ^(b26,i/2(i-o.io)+i)^ = (-OO; ^(b26,l/2(C).90) + l)^ = (-oo;Z<">) = (-oo; -7) 258 qbinomCl - alpha, n, 1 / 2) + 1 # 17 259 humer . LLRs [17] # -7 • Závěr testování Protože Žq = 0 nenáleží do 90% empirického pravostranného intervalu spolehlivosti, tj. xq = 0 ^ IS, Hq zamítáme na hladině významnosti a = 0.10. 5. Testování p-hodnotou • p-hodnota p-hodnota = Pr(SE < sE) = Pr(SE < 1) = 4.02331"7 260 pbinom(SE, m, 1 / 2) # 4.023314e-07 • Závěr testování Protože p-hodnota = 4.02331-7 je menší než a = 0.10, Hq zamítáme na hladině významnosti a = 0.10. 6. Interpretace výsledků Na základě všech tří typů testování zamítáme nulovou hypotézu na hladině významnosti a = 0.10. Délka pažní kosti na levé straně u žen je statisticky významně větší než na pravé straně. 7. Grafická vizualizace výsledku testování Rozdíl v délce pažní kosti na pravé a levé straně vizualizujeme pomocí krabicového diagramu (viz obrázek 16). Poznámka: Znaménkový párový test můžeme provést pomocí funkce SIGN.test(). Vstupními parametry budou vektor naměřených hodnot délek pažních kostí na levé straně (humer.LL), vektor naměřených hodnot délek pažních kostí na pravé straně (humer.LR), volba párového testu (paired = T), hodnota hladiny významnosti a zadaná prostřednictvím 38 o •a 10 - 0 - -5 - -10 -15 - -20 -25 - B B >^ ř, B 5 o 'B M rú > "u •a 310 - 300 - 290 - 280 - 270 - pravá strana leva strana Obrázek 16: Krabicový diagram rozdílů mezi délkami těla pažní kosti na pravé a levé straně u skeletů žen z římského pohřebiště v Poundbury 261 BSDA::SIGN.test(numer.LL, numer.LR, paired = T, 262 conf.level = 0.90, alternative = less ' ) koeficientu spolehlivosti 1 — a (conf.level = 0.90) a levostranný typ alternativní hypotézy (alternativě = 'less'). Součástí výstupu je hodnota mediánu rozdílů naměřených délek na levé a pravé straně medián of x-y = -7.5, hodnota testovací statistiky S = 1, interpolované hranice 90% empirického pravo intervalu spolehlivosti -Inf a -7 pro rozdíl naměřených délek na levé a pravé straně a p-hodnota p-value = 4.023e-07. Jediné, co musíme stanovit zvlášť, je horní hranice kritického oboru. 39 263 Dependent - samples Sign-Test 264 265 data: humer.LL and humer.LR 266 S = 1, p-value = 4.023e-07 267 alternative hypothesis: true median difference is less than 0 268 90 percent confidence interval: 269 -Inf -7 270 sample estimates: 271 median of x-y 272 -7.5 273 274 Achieved and Interpolated Confidence Intervals: 275 276 Conf.Level L.E.pt U.E.pt 277 Lower Achieved CI 0.8365 -Inf -7 278 Interpolated CI 0.9000 -Inf -7 279 Upper Achieved CI 0.9157 -Inf -7 280 40 Příklad 9.9. Znaménkový párový asymptotický test Máme datový soubor 02-paired-means-clavicle.txt obsahující údaje o hodnotách vertikálního průměru středu délky těla klíční kosti z pravé a levé strany (clavicula) z pohřebiště u Sv. Jakuba v Brně, převážně z období středověku, naměřené jedním výzkumníkem ve dvou opakovaných měřeních (viz sekce ??). Hodnoty naměřené při prvním opakování jsou uloženy v proměnné simd.l, hodnoty naměřené při druhém opakování jsou uloženy v proměnné simd.2. Na hladině významnosti a = 0.05 zjistěte, zda je střední hodnota prvního měření shodná se střední hodnotou druhého měření vertikálního průměru délky těla klíční kosti na pravé straně provedené tímto výzkumníkem. Řešení příkladu ?? Nejprve příkazem read.delim() načteme datový soubor a pomocí operátoru [] z něj vybereme pouze údaje naměřené sledovaným výzkumníkem (sloupce simd.l a simd.2) na pravé straně side == 'R'. Údaje vložíme do proměnné data.l2R. Z datové tabulky následně odstraníme chybějící údaje a zjistíme rozsah náhodného výběru. 281 data <- read.delim('00-Dat 282 #head(data) 283 data.l2R <- data [data$side 284 data.l2R <- na.omit(data . 1 285 dim(data.12R) # 40x2 a//02-paired-means-clavicle.txt') == 'R', cCsimd.ľ, 'simd.2')] 2R) Datový soubor obsahuje 40 opakovaných měření délky těla klíční kosti z pravé strany. Jelikož máme za úkol porovnat opakovaná měření provedená na jednom subjektu, použijeme k tomuto porovnání párový test. Prvním krokem k provedení tohoto testu je vytvoření rozdílů hodnot získaných v prvním a druhém měření. V druhém kroku je potřeba ověřit předpoklad normálního rozdělení těchto rozdílů. 286 simd.lR <- data.12R$simd.1 287 simd.2R <- data.12R$simd.2 288 simd.l2R <- simd.IR - simd.2R Vzhledem k velkému rozsahu náhodného výběru otestujeme normalitu Lillieforsovým testem (a = 0.05) v kombinaci s kvantilovým diagramem a histogramem. Naměřené hodnoty rozdělíme do šesti ekvidistantních intervalů se šířkou 0.17mm prostřednictvím stanovených hranic —0.37, —0.20,..., 0.65 (viz obrázek 17). Obrázek 17: Histogram a diagram rozdílů prvního a druhého měření vertikálního průměru ve středu délky těla klíční kosti na pravé straně získaných jedním výzkumníkem Protože p-hodnota = 0.0001 je menší než 0.05, nulovou hypotézu o normalitě rozdílů zamítáme na hladině významnosti a = 0.05. Histogram naměřených hodnot, stejně jako kvantilový diagram ukazují na vyšikmení naměřených hodnot směrem doleva s odlehlými pozorováními na pravé straně, které jsou pro normální charakter naměřených hodnot 41 fatální. Náhodný výběr rozdílů hodnot získaných v prvním a druhém měření nepochází z normálního rozdělení. Předpoklad normality pro použití parametrického testu není splněn, proto není možné tvrzení ze zadání ověřit pomocí parametrického párovéto testu. Otázku ze zadání tedy ověříme pomocí neparametrické alternativy, a sice asymptotického znaménkového párového testu. Naším úkolem je zjistit, zda je střední hodnota prvního měření shodná se střední hodnotou druhého měření. Toto tvrzení bude zněním nulové hypotézy, neboť shoda implikuje rovnost a rovnost je vždy součástí nulové hypotézy. Dále je potřeba si uvědomit, že možnost testovat hypotézu o středních hodnotách jsme ztratili při rozhodnutí použít k testování neparametrický test. Namísto toho budeme testovat hypotézu o shodě mediánu prvních měření s mediánem druhých měření, což je analogie hypotézy, že medián rozdílů prvního a druhého měření je rovný nule. 1. Stanovení hypotéz • slovní formulace nulové a alternativní hypotézy H0 : Medián rozdílů prvního a druhého měření vertikálního průměru těla klíění kosti na pravé straně měřených jedněm výzkumníkem je rovný nule. H\ : Medián rozdílů prvního a druhého měření vertikálního průměru těla klíění kosti na pravé straně měřených jedním výzkumníkem není rovný nule. • matematická formulace nulové a alternativní hypotézy H o : ž = Žq, kde Žq = 0 Hi : ž < Žq, kde Žq = 0 (levostranná alternativa) 2. Volba hladiny významnosti • Hladinu významnosti podle zadání a = 0.05. 3. Testování kritickým oborem • Testovací statistika V úvodu příkladu jsme vytvořili vektor simd.l2R {Z{) jako vektor rozdílů mezi naměřenými hodnotami v rámci prvního a druhého měření. Nyní je potřeba stanovit rozdíly Z i — Žq a určit počet kladných rozdílů Zi — Žq a počet nenulových rozdílů m. 289 zO <- 0 290 I <- (simd.l2R > zO) 291 SE <- sum(I) # 19 292 n <- length(simd.12R) # 40 293 m <- sum(simd.l2R - zO != 0) #39 Z celkového počtu 40 rozdílů je právě jeden rozdíl nulový, tj. m 39 a 19 rozdílů je kladných, tj. Se = 19. n 40 SE = ^I{Zi - ž0 > 0) = 0 > 0) = 19. i=l i=l Nyní vypočítáme testovací statistiku asymptotické varianty znaménkového testu podle vzorce 9.2. S a 19-f 0 -0.5 = -0.1601282 = -0.16013 39 - 19.5 3.122499 294 SA <- (SE - m / 2) / sqrt (m / 4) # -0.1601282 42 • Kritický obor W = ( = ( (-00 ; ua/2) U (u1_a/2 ; 00) (-00 ; «0.05/2) U («i_o.o5/2 ; 00) (-00 ; M0.025) U («0.975 ; 00) (-00; -1.959964) U (1.959964; 00) 295 alpha <- 0.05 296 qnorm(alpha / 2) # -1.959964 297 qnorm(l - alpha / 2) # 1.959964 • Závěr testování Protože realizace testovací statistiky sa = —0.16013 nenáleží do kritického oboru, tj. sa <ýz W, Hq nezamítáme na hladině významnosti a = 0.05. 4. Testování intervalem spolehlivosti • Interval spolehlivosti Pomocí příkazu sort() nejprve vzestupně seřadíme naměřené hodnoty Zi, i = 1, ...,40. Hranice 95% intervalu spolehlivosti budou potom (Ci_a/2)-tá a (n + 1 — Ci_a/2)-tá hodnota v seřazeném vektoru naměřených hodnot, kde 298 simd.l2Rs <- sort(simd.12R) 299 Cl <- round(n / 2 - qnorm(l - alpha / 2) * sqrt(n / 4)) # 14 300 C2 <- round(n + 1 - Cl) #27 301 simd.12Rs [Cl] # -0.06 302 simd.12Rs[C2] # O.O4 • Závěr testování Protože Žq = 0 náleží do 95% empirického oboustranného intervalu spolehlivosti, tj. Žq = 0 G IS, Hq nezamítáme na hladině významnosti a = 0.05. = 20-«0.975^ = 20 - 1.959964 x 3.162278 = 13.80205 = 14 Interval spolehlivosti má potom následující tvar (-0.06; 0.04) 43 5. Testování p-hodnotou • p-hodnota p-hodnota = 2min{Pr(5'A < s a) , Pt(Sa > sa)} = 2min{Pr(5A < -0.16013), Pr(SA > -0.16013)} = 2min{Pr(5A < -0.16013), 1 - Pr(SA < -0.16013)} = 2min{0.4363901, 0.5636099} = 2 x 0.4363901 = 0.8727801 = 0.8728 303 p.val <- 2 * min(pnorm(SA), 1 - pnorm(SA)) # 0.8727801 • Závěr testování Protože p-hodnota : 0.8727801 je větší než a = 0.05, Hq nezamítáme na hladině významnosti a = 0.05. 6. Interpretace výsledků Za základě všech tří typů testování nezamítáme nulovou hypotézu na hladině významnosti a = 0.05. Mezi prvním a druhém měřením vertikálního průměru délky těla klíční kosti na pravé straně provedenými sledovaným výzkumníkem neexistuje statisticky významný rozdíl. 7. Grafická vizualizace výsledku testování Porovnání naměřených hodnot získaných v rámci prvního a druhého měření zobrazíme pomocí krabicového diagramu (viz obrázek 18). o o 0.4 - 0.2 - 0.0 -0.2 - 3 -0.4 - •a N O o B S 'o •a 14 - 12 - 10 - l.mereni 2.mereni Obrázek 18: Krabicový diagram rozdílů mezi nejvétší výškou kloubní jamky u mužů na pravé a levé straně 44 Příklad 9.10. Znaménkový párový asymptotický test Máme datový soubor 02-paired-means-clavicle.txt obsahující údaje o hodnotách vertikálního průměru středu délky těla klíční kosti z pravé a levé strany (clavicula) z pohřebiště u Sv. Jakuba v Brně, převážně z období středověku, naměřené jedním výzkumníkem ve dvou opakovaných měřeních (viz sekce ??). Hodnoty naměřené při prvním opakování jsou uloženy v proměnné simd.l, hodnoty naměřené při druhém opakování jsou uloženy v proměnné simd.2. Můžeme zjistit, že aritmetický průměr hodnot získaných v rámci prvního měření je větší než aritmetický průměr hodnot získaných v rámmci druhého měření. Na hladině významnosti a = 0.05 zjistěte, zda je střední hodnota prvního měření větší než střední hodnota durhého měření vertikálního průměru délky těla klíční kosti na levé straně provedené tímto výzkumníkem. Řešení příkladu 9.10 Načteme datový soubor a pomocí operátoru [] z něj vybereme pouze údaje naměřené sledovaným výzkumníkem (sloupce simd.l a simd.2) na levé straně side == 'L'. Údaje vložíme do proměnné data.l2L. Z datové tabulky následně odstraníme chybějící údaje a zjistíme rozsah náhodného výběru. 304 data <- read.delim('00-Data//02-paired-means-clavicle.txt') 305 #head(data) 306 data.l2L <- data[data$side == 'L', c('simd.l', 'simd.2')] 307 data.l2L <- na.omit(data.12L) 308 dim(data.l2L) # 40 x 2 Datový soubor obsahuje 40 opakovaných měření délky těla klíční kosti z levé strany. Jelikož chceme porovnat opakovaná měření provedená na jednom subjektu, použijeme k tomuto porovnání párový test. Prvním krokem k provedení tohoto testu je vytvoření rozdílů hodnot získaných v prvním a druhém měření. V druhém kroku je potřeba ověřit předpoklad normálního rozdělení těchto rozdílů. Testováním normality náhodného výběru rozdílů jsme se zabývali v rámci příkladu ?? v sekci ??, přičemž jsme zjistili, že tento náhodný výběr nepochází z normálního rozdělení (p-hodnota = 0.0002), a to ani po odstranění nejvíce odlehlého rozdílu (p-hodnota = 0.0302). Pro připomenutí viz obrázek 19 zobrazující histogram a kvantilový diagram pro původní sadu 40 rozdílů. 309 simd.lL <- data.12L$simd.1 310 simd.2L <- data.12L$simd.2 311 simd.l2L <- simd.lL - simd.2L Obrázek 19: Histogram a diagram rozdílů prvního a druhého měření vertikálního průměru ve středu délky těla klíční kosti na levé straně získaných jedním výzkumníkem Protože předpoklad normality rozdílů prvního a druhého měření není splněn, a k jeho splnění nedošlo ani po 45 odstranění nej odlehlejšího pozorování, vrátíme se zpátky k vektoru obsahujícímu všech 40 rozdílů. Otázku ze zadání ověříme pomocí neparametrické alternativy párového testu, a sice pomocí asymptotického znaménkového párového testu. Naším úkolem je zjistit, zda je střední hodnota prvního měření větší než střední hodnota druhého měření. Toto tvrzení bude zněním alternativní hypotézy, neboť v zadání není o znění nulové hypotézy žádná zmínka. Analogické znění alternativní hypotézy je, že střední hodnota rozdílů, získaných odečtením hodnot naměřených v druhém měření od hodnot naměřených v prvním měření, bude větší než nula. Nulovou hypotézu naformulujeme jako doplněk k této hypotéze. Namísto středních hodnot se opět zaměříme na jejich neparametrické alternativy a sice na mediány, které budou figurovat v konečném znění obou hypotéz. 1. Stanovení hypotéz • slovní formulace nulové a alternativní hypotézy Hq : Medián rozdílů prvního a druhého měření vertikálního průměru těla klíční kosti na levé straně měřených jedním výzkumníkem je menší než nula. Hi : Medián rozdílů prvního a druhého měření vertikálního průměru těla klíční kosti na levé straně měřených jedním výzkumníkem je větší než nula. • matematická formulace nulové a alternativní hypotézy H o : ž < Žq, kde Žq = 0 H\ : ž > Žq, kde Žq = 0 (pravostranná alternativa) 2. Volba hladiny významnosti • Hladinu významnosti volíme podle zadání a = 0.05. 3. Testování kritickým oborem • Testovací statistika V úvodu příkladu jsme vytvořili vektor simd.l2L (Zi) jako vektor rozdílů mezi naměřenými hodnotami v rámci prvního a druhého měření na levé straně. Nyní je potřeba stanovit rozdíly Zi — Žq a určit počet kladných rozdílů Z{ — ž0 a počet nenulových rozdílů m. 312 zO <- 0 313 I <- (simd.l2L > zO) 314 SE <- sum(I) # 23 315 n <- length(simd.12L) # 40 316 m <- sum(simd.l2L - zO != 0) # 38 Z celkového počtu 40 rozdílů jsou dva rozdíly nenulové, tj. m = 38 a 23 rozdílů je kladných, tj. Se = 23. n 40 Se = j2i(-z* " žo > 0) = ^/(Z2 - 0 > 0) = 23. i=l i=l Dále vypočítáme testovací statistiku asymptotické varianty znaménkového testu Sa podle vzorce 9.2. S a 23-f 0 4 = 1.297771 = 1.2978 23 - 19 3.082207 317 SA <- (SE - m / 2) / sqrt(m / 4) # 1.297771 46 • Kritický obor W = ; oo) = («1-0.05; oo) = («0.95 ; oo) = (1.644854; oo) 318 alpha <- 0.05 319 qnormCl - alpha) # 1. 644854 • Závěr testování Protože realizace testovací statistiky sa = 1.297771 nenáleží do kritického oboru, tj. sa W, Hq nezamítáme na hladině významnosti a = 0.05. 4. Testování intervalem spolehlivosti • Interval spolehlivosti Pomocí příkazu sort() nejprve vzestupně seřadíme naměřené hodnoty Zi, i = 1, ...,40. Hranice 95% intervalu spolehlivosti budou potom (Ci_a)-tá hodnota v seřazeném vektoru naměřených hodnot, kde n n °l-« - % _Ml""V 4 40 /4ČJ - y - Mi-o.o5 y y 20 - M0.95VIO 20 - 1.644854 x 3.162278 14.79851 = 15 a plus nekonečno. Interval spolehlivosti má potom následující tvar (d,h) = (a-(ci —); (X<15) ; 00) (-0.01; 00) 320 simd.l2Ls <- sort(simd.12L) 321 Cl <- round(n / 2 - qnorm(l - alpha) * sqrt(n / 4)) #14 322 simd.12Ls[Cl] # -0.01 • Závěr testování Protože Žq = 0 náleží do 95% empirického oboustranného intervalu spolehlivosti, tj. Žq = 0 G IS, Hq nezamítáme na hladině významnosti a = 0.05. 5. Testování p-hodnotou 47 • p-hodnota p-hodnota = Pr(SA > sA) = Pr(SA > 1.297771) = 1 - Pr(SA < 1.297771) = 0.09718296 = 0.09718 323 p.val <- 1 - pnorm(SA) # 0.09718296 • Závěr testování Protože p-hodnota = 0.09718 je větší než a = 0.05, H0 nezamítáme na hladině významnosti a = 0.05. 6. Interpretace výsledků Za základě všech tří typů testování nezamítáme nulovou hypotézu na hladině významnosti a = 0.05. Hodnoty získané v prvním měření nejsou statisticky významně větší než hodnoty získané v druhém měření. 7. Grafická vizualizace výsledku testování Porovnání naměřených hodnot získaných v rámci prvního a druhého měření zobrazíme pomocí krabicového diagramu (viz obrázek 20). o 'o 1.5 - 1.0 0.5 0.0 - -0.5 - o B 14 12 10 l.mereni 2.mereni Obrázek 20: Krabicový diagram rozdílů prvního a druhého měření vertikálního průměru ve středu délky těla klíční kosti na levé straně získaných jedním výzkumníkem 48 Příklad 9.11. Znaménkový párový asymptotický test Máme datový soubor 03-paired-means-clavicle2.txt obsahující údaje o délkách klíční kosti (clavicula) z pravé strany (length.R) a levé strany (length.L) z anglického souboru dokumentovaných skeletů (Parsons, 1916, viz soubor D-03-paired-means-clavicle2). Na hladině významnosti a = 0.05 testujte hypotézu, že střední hodnota délky klíční kosti u žen z pravé strany je větší než střední hodnota délky klíční kosti u žen z levé strany. Řešení příkladu 9.11 Nejprve načteme datový soubor a pomocí operátoru [] z něj vybereme pouze údaje o délce klíční kosti z pravé strany (sloupce length.R resp. z levé strany (length.L) u žen sex == 'f. Údaje vložíme do proměnné length.RLF. Z datové tabulky následně odstraníme chybějící hodnoty a zjistíme rozsah náhodného výběru. 324 data <- read.delim('00-Data//03-paired-means-clavicle2.txt') 325 #head(data) 326 data.RLF <- data[data$sex == 'i', c('length.R ' , 'length.L')] 327 data.RLF <- na.omit(data.RLF) 328 dim(data.RLF) # 50 Datový soubor obsahuje údaje o délce klíční kosti z pravé a levé strany u 50 mužů. Úkolem ze zadání je porovnat naměřené hodnoty na pravé a levé straně. Jde tedy o měření stejného znaku (délka klíční kosti) sledovaného na stejných subjektech (ženy), proto použijeme na tuto situaci párový test. Prvním krokem k tohoto testu je vytvoření rozdílů hodnot naměřených na pravé a levé straně. V druhém kroku je potřeba ověřit předpoklad normality těchto rozdílů. Předpoklad normality ověříme Lillieforsovým testem (a = 0.05) v kombinaci s kvantilovým diagramem a histogramem (viz obrázek 21). Datový soubor rozdělíme do sedmi ekvidistatních intervalů s šířkou 2mm prostřednictvím stanovených hranic —8, —6,..., 6. Obrázek 21: Histogram a kvantilový diagram rozdílů délky těla klíčních kostí na levé a pravé strané žen Protože p-hodnota = 0.0481 je menší než 0.05, nulovou hypotézu o normalitě dat, zamítáme na hladině významnosti a = 0.05. Z histogramu vidíme, že rozdíly pravé a levé strany vykazují plošší trend než očekávaná křivka hustoty normálního rozdělení. Taktéž zde pozorujeme mírné vyšikmení hodnot směrem doleva s odlehlými pozorováními na pravé straně. Lillieforsův test tyto nedostatky vyhodnotil vzhledem k poměrně vysokému rozsahu náhodného výběru jako fatální pro normální charakter náhodného výběru. Náhodný výběr rozdílů délek klíčních kostí z levé a pravé strany tedy nepochází z normálního rozdělení. Předpoklad normality pro použití parametrického testu není splněn, proto není možné tvrzení ze zadání ověřit pomocí parametrického párového testu. Otázku ze zadání tedy ověříme pomocí neparametrického asymptotického 49 znaménkového párového testu. Naším úkolem je tetsovat hypotézu, že střední hodnota délky klíční kosti u žen z pravé strany je větší než střední hodnota délky klíční kosti u žen z levé strany. Toto tvrzení je zněním nulové hypotézy, protože pojmem hypotéza zmíněným v zadání je vždy myšlena nulová hypotéza. Dále je potřeba si uvědomit, že možnost testovat hypotézu o středních hodnotách jsme ztratili při rozhodnutí použít k testování neparametrický test. Pojem střední hodnoty nahradíme ve slovní formulaci hypotéz pojmem medián. 1. Stanovení hypotéz • slovní formulace nulové a alternativní hypotézy Hq : Medián rozdílů délky klíční kosti z pravé a z levé strany u žen je větší nebo roven nule. H\ : Medián rozdílů délky klíční kosti z pravé a z levé strany u žen je menší než nula. • matematická formulace nulové a alternativní hypotézy H o : ž > Žq, kde Žq = 0 H\ : ž < Žq, kde Žq = 0 (levostranná alternativa) 2. Volba hladiny významnosti • Hladinu významnosti podle zadání zvolíme jako a = 0.05. 3. Testování kritickým oborem • Testovací statistika V úvodu příkladu jsme vytvořili vektor length.RLF (Zi) jako vektor rozdílů mezi naměřenými hodnotami na pravé a levé straně. Nyní je potřeba stanovit rozdíly Zi — Žq a určit počet kladných rozdílů Zi — žq a počet nenulových rozdílů m. 329 zO <- 0 330 I <- (length.RLF > zO) 331 SE <- sum(I) # 19 332 n <- length(length.RLF) # 50 333 m <- sum(length.RLF - zO != 0) # 43 Z celkového počtu 50 rozdílů je sedm rozdílů rovných nule, tj. počet nenulových rozdílů m = 43 a 19 rozdílů je kladných, tj. S e = 19. 50 SE = Y,1^ " žo > 0) = - 0 > 0) = 19. i=l i=l Nyní vypočítáme testovací statistiku asymptotické varianty znaménkového testu podle vzorce 9.2. SE-f 19-f 0 .- -2.5 Sa = r-2 =-t==~ = t7,—^^10.75 = ——— = -0.7624929 = -0.7625 ^ /43 19 - 21.5 3.278719 334 SA <- (SE - m / 2) / sqrt(m / 4) # -0.7624929 • Kritický obor W = (—oo ; ua) = (-oo; «0.05) = (-oo; -1.644854) 50 335 alpha <- 0.05 336 qnorm(alpha) # -1.644854 • Závěr testování Protože realizace testovací statistiky sa = —0.16013 nenáleží do kritického oboru, tj. sa <ýz W, Hq nezamítáme na hladině významnosti a = 0.05. 4. Testování intervalem spolehlivosti • Interval spolehlivosti Pomocí příkazu sort() nejprve vzestupně seřadíme naměřené hodnoty Zi, i = 1, ...,50. Hranice 95% pravostranného intervalu spolehlivosti budou potom mínus nekonečno a (n + 1 — Ci_a)-tá hodnota v seřazeném vektoru naměřených hodnot, kde = 25 — M0.95 v 12.5 = 25 - 1.644854 x 3.535534 = 19.18456 = 19 Interval spolehlivosti má potom následující tvar (d,h) = (-00; a- 5 - 0 - -5 - -10 - l.mereni 2.mereni Obrázek 22: Krabicový diagram rozdílů délky těla klíčních kostí na levé a pravé strané žen 52 9.4 Wilcoxonův jednovýběrový exaktní test Nechť Xi,..., Xin, n > 2 je náhodný výběr ze spojitého rozdělení s hustotou f(x), která je symetrická okolo mediánu x, a nechť xq je konstanta. Na hladině významnosti a testujeme jednu z následujících tří hypotéz oproti příslušné alternativní hypotéze. Hqi : x = xq oproti Hu : x ^ xq (oboustranná alt.) Hq2 : x < xq oproti H12 : x > xq (pravostranná alt.) Hqs : x > xq oproti H13 : x < xq (levostranná alt.) Test nazýváme Wilcoxonův jednovýběrový exaktní test o mediánu x. Testovací statistika má tvar SE = ^i?í/(Xí -x0 > 0) (9.3) i=l kde m je počet nenulových rozdílů Xi — xq, Rí je pořadí absolutní hodnoty rozdílu i-té náhodné veličiny Xi a konstanty xq, tj. Ri je pořadí seřazených rozdílů \Xi — xq\, i = 1,..., m, I(Xi — xq > 0) je indikační funkce, která nabývá hodnoty 1, pokud Xi — xq > 0, a hodnoty 0 jinak. Statistika Se je potom součet pořadí absolutních hodnot rozdílů \Xi — xq\, pro něž je rozdíl Xi — xq kladný. Kritický obor podle zvolené alternativní hypotézy má tvar Hn-.x^ío W = (-00; sm(a/2) - 1) U (sm(l - a/2); 00) H12 : x > x0 W = {sm(l - a); 00) íři3 : x < Xq W = (—00 ; sm(a) — 1) kde sm(a/2), sm(l — a/2), sm{a) a sm(l — a) jsou tabelované kvantily pro jednovýběrový Wilcoxonův test, jejichž hodnoty získáme pomocí softwaru CŠT a implementované funkce qsignrank(). Interval spolehlivosti má podle zvolené alternativní hypotézy jeden z následujících tvarů Hu:xr x0 (d, h) = (j/(s™(«/2)) . y(sm{í-a/2) + í)^ H12:x> x0 (d, 00) = (l/(s-(«)) ; 00) H13:í< x0 (-00, h) = (-00 ; V^1-"^) kde < < ■ ■ ■ < (2+ ') značí poslupnost vzestupně seřazených "1<-Tľ+1-) Walshových průměrů ííiiíil r»(m+l)^ ^^oln^n^ot irmotnr.no oor-c 7amrr.h - noioiiuvjui jJL umel u -2~ i = 1,..., m, j = 1,..., m, j < i a značí fc-tý seřazený Walshův průměr. Posloupnost Walschových průměrů získáme příkazem owa z knihovny NSM3. p-hodnota má v závislosti na zvolené alternativní hypotéze jeden z následujících tvarů Hu : x, ^ xq p-hodnota = 2min{Pr(5£; < se) , Pr(5£; > se)} H12 : x > Íq p-hodnota = Vt{Se > se) H13 : x, < xq p-hodnota = Pr(5s < se) kde Se je náhodná veličina, se je realizace testovací statistiky Se (viz vzorec 9.3), tedy konkrétní číslo, a Pr(5E < se) je distribuční funkce tabelovaného rozdělení pro jednovýběrový Wilcoxonův test, jejíž hodnotu získáme pomocí <® a implementované funkce psignrankQ. 53 Příklad 9.12. Wilcoxonův jednovýběrový exaktní test Mějme datový soubor 21-goldman-measures.csv a proměnnou humer.DR popisující největší délku hlavice ramenní kosti z pravé strany v mm u skeletů z období raného středověku z oblasti Staubing (viz sekce ??). Dále máme k dispozici údaje ze studie (Mali et al.) z roku 1999, v rámci které byla měřena největší délka hlavice ramenní kosti u mužů současné německé populace (mm = 50.0mm, nm = 64). Na hladině významnosti a = 0.05 zjistěte, zda existuje rozdíl mezi největší délkou hlavice ramenní kosti z pravé strany u mužů z raně středověké německé populace a u mužů současné německé populace. Řešení příkladu 9.12 Zadání příkladu je shodné se zadáním příkladu 9.1. Datový soubor obsahující údaje o největší délce hlavice ramenní kosti z pravé strany u 12 mužů nepochází z normálního rozdělení (viz kvantilový diagram, obrázek 23 vlevo), proto na prozkoumání otázky ze zadání použijeme neparametrický jednovýběrový test. Tentokrát bychom však chtěli nulovou hypotézu otestovat pomocí Wilcoxonova jednovýběrového testu. K použití tohoto testu však musí být splněn předpoklad symetrického rozdělení naměřených hodnot okolo mediánu. Ten ověříme pomocí testu symetrie (příkaz symmetry.test() z knihovny lawstat) uvedeného v sekci ?? (hladinu významnosti zvolíme a = 0.05). Graficky vizualizujeme míru symetrie naměřených hodnot okolo mediánu pomocí histogramu superponovaného jádrovým odhadem křivky hustoty (viz kapitola ??) se zvýrazněnou hodnotou mediánu pomocí vertikální přímky (viz obrázek 23 vpravo). -1.5 -0.5 0.5 1.0 1.5 41.25 46.25 51.25 teoreticky kvantil polomer hlavice ramenniho klubu (v mm) Shapiruv test: p-hodnota = 0.0261 Test symetrie: p-hodnota = 0.0518 Obrázek 23: Kvantilový diagram (vlevo) a histogram naměřených hodnot poloměru hlavice ramenniho kloubu na pravé straně u skeletů mužů ze starověké populace z oblasti Staubing superponovaný křivkou jadrového odhadu hustoty (vpravo) Protože p-hodnota = 0.0518 je větší než 0.05, nezamítáme hypotézu o symetrii okolo mediánu na hladině významnosti a = 0.05. Z histogramu na obrázku 23 vidíme, že data nejsou ukázkově symetrická okolo mediánu, jsou vyšikmená doprava s prodlouženým levým koncem. S ohledem na nízký rozsah náhodného výběru jsou však tyto prohřešky vzhledem k symetrii ještě přijatelné. Naměřené hodnoty největší délky hlavice ramenní kosti z pravé strany tedy považujeme za symetricky rozdělené okolo mediánu. Protože hypotéza o symetrii data okolo mediánu nebyla zamítnuta, můžeme k otestování otázky ze zadání použít Wilcoxonův jednovýběrový parametrický test. Analogicky jako v příkadu 9.1 provedeme testování v posloupnosti sedmi kroků. 1. Stanovení hypotéz • slovní formulace nulové a alternativní hypotézy Hq : Medián největší délky hlavice ramenní kosti raně středověké německé populace je shodný s mediánem největší délky hlavice ramenní kosti na pravé straně mužů současné německé populace. Hi : Medián největší délky hlavice ramenní kosti raně středověké německé populace není shodný s mediánem největší délky hlavice ramenní kosti na pravé straně mužů současné německé populace. 54 • matematická formulace nulové a alternativní hypotézy Hq : x = xq, kde xq = 50.00 Hi : x, ^ xq, kde xq = 50.00 (oboustranná alternativa) 2. Volba hladiny významnosti • Hladinu významnosti volíme v souladu se zadáním jako a = 0.05. 3. Testování kritickým oborem • Testovací statistika Nejprve vypočítáme vektor rozdílů naměřených hodnot Xi a konstanty xq, tj. Xi — xq, stanovíme, které rozdíly jsou kladné, a zjistíme, kolik rozdílů je nenulových (m). Dále vypočítáme absolutní hodnoty nenulových rozdílů, tj. \Xi — xq\, i = 1,..., m a stanovíme pořadí těchto rozdílů v absolutních hodnotách (viz tabulka 11). Tabulka 11: Naměřené hodnoty Xi, rozdíly Xi — xq, znaménka těchto rozdílů, absotulní hodnoty nenulových rozdílů \Xi — xq\ a pořadí nenulových rozdílů v absolutních hodnotách Ri měření 1 2 3 4 5 6 7 8 9 10 11 12 49.83 47.96 49.49 48.73 49.25 40.9 51.22 51.27 46.16 49.52 50.28 45.6 Xi — Xq t -0.17 -2.04 -0.51 -1.27 -0.75 -9.1 1.22 1.27 -3.84 -0.48 0.28 -4.4 1 \Xl - x0\ 0.17 2.04 0.51 1.27 0.75 9.1 + 1.22 + 1.27 3.84 0.48 + 0.28 4.4 1.00 9.00 4.00 7.50 5.00 12.0 6.00 7.50 10.00 3.00 2.00 11.0 Z tabulky 1 vidíme, že všechny rozdíly jsou neulové, tedy m = 12. Dále vidíme, že celkem tři rozdíly Xi — xq, i = 1, • • •, 12 jsou kladné, přičemž pořadí těchto rozdílů v absolutních hodnotách jsou 2, 6 a 7.5. Hodnota testovací statistiky Se, která je definovaná jako součet pořadí absolutních hodnot kladných rozdílů, bude tedy rovná 15.5. 12 SE = RíI(Xí - ío > 0) = 2 + 6 + 7.5 = 15.5 i=l 342 xO <- 50 343 I <- (numer.DRM - xO > 0) 344 Ri <- rank(abs(numer.DRM-x0)) 345 tab <- rbind("Xi" = numer.DRM, "Xi-xO" = numer.DRM - xO, 346 "I" = I, "|Xi-x0|" = abs(numer.DRM - xO), 347 "Ri" = Ri) 348 tab <- data.frame(tab) 349 names(tab) <- 1 : 12 350 # 1 2 3 4 5 6 7 8 9 10 11 12 351 # Xi 49 83 47. 96 49 49 48. 73 49 25 40 9 51 22 51 . 27 46 16 49 52 50 28 45. 6 352 # Xi-xO -0 17 -2. 04 -0 51 -1. 27 -0 75 -9 1 1 22 1 . 27 -3 84 -0 48 0 28 -4.4 353 # I 0 00 0. 00 0 00 0. 00 0 00 0 0 1 00 1 . 00 0 00 0 00 1 00 0. 0 354 # IXi-xO 1 0 17 2. 04 0 51 1 . 27 0 75 9 1 1 22 1 . 27 3 84 0 48 0 28 4.4 355 # Ri 1 00 9. 00 4 00 7. 50 5 00 12 0 6 00 7. 50 10 00 3 00 2 00 11.0 356 SE <- sum(Ri * I) # 15.5 • Kritický obor 55 W = (-00; sm(a/2) - 1) U (sm(l - a/2); 00) = (-00; si2(0.05/2) - 1) U (s12(l - 0.05/2); 00) = (-00 ; si2(0.025) - 1) U (si2(0.975); 00) = (-00; 14 - 1) U (64; 00) = (-00; 13) U (64; 00) 357 m <- sum(humer.DRM - xO != 0) 358 alpha <- 0.05 359 qsignrank(alpha / 2, m) - 1 # 13 360 qsignrank (1 - alpha / 2, m) #64 • Závěr testování Protože realizace testovací statistiky sE = 15.5 nenáleží do kritického oboru, tj. sE £ W, H0 nezamítáme na hladině významnosti a = 0.05. 4. Testování intervalem spolehlivosti • Interval spolehlivosti Abychom mohli stanovit hranice 95 % intervalu spolehlivosti, musíme nejprve vypočítat T"(™+1) Wal-shových průměrů x'+x?; i = l,...,m a ty následně vzestupně seřadit. Walshovy průměry získáme pomocí funkce owa z knihovny NSM3. Prvním argumentem této funkce bude vektor nul o délce m = 12, druhým argumentem bude vektor naměřených délek ramenních kostí z pravé strany humer.DRM. Získané Walschovy průměry seřadíme příkazem sort(). 361 walsh <- NSM3::owa(rep(0, 12), humer.DRM)$owa 362 V <- sort(walsh) 363 # [ U 40 900 43 250 43 530 44 430 44 815 45 075 45 195 45 210 45 365 45 590 45 600 364 # [12] 45 880 46 060 46 085 46 160 46 780 47 060 47 165 47.425 47.445 47. 545 47. 560 365 # [23] 47 705 47. 715 47. 825 47. 840 47. 94O 47 960 47 995 48 220 48 345 48 4IO 48 435 366 # [34] 48 605 48 690 48 715 48 725 48 730 48 74O 48 895 48 990 49 110 49 120 49 125 367 # [45] 49 250 49 280 49 370 49 385 49 490 49 505 49 505 49 520 49 54O 49 590 49 615 368 # [56] 49 660 49 675 49 765 49 830 49 885 49 900 49 975 50 000 50 055 50 235 50 260 369 # [67] 50 280 50 355 50 370 50 380 50 395 50 525 50 550 50 750 50 775 51 220 51 245 370 # [78] 51 270 Hranice intervalu spolehlivosti potom tvoří ty hodnoty, které se v seřazeném vektoru Walshových průměrů Vs nachází na (sm(oí/2))-té pozici a na (sm(l — a/2) + l)-té pozici. Hodnoty sm(a/2) asra(l- a/2) + 1 nalezneme opět pomocí funkce qsignrankQ. = (y(sm(a/2)) . ^(Sm(l-a/2) + l)^ = ^(si2(0.05/2)). ^(s12(l-0.05/2) + l)^ V(14) ; V{e4)^j = (46.085; 50.235) • Závěr testování Protože xq = 50.00 náleží do 95% empirického oboustranného intervalu spolehlivosti, tj. xq = 50.00 G IS, Hq nezamítáme na hladině významnosti a = 0.05. 56 371 qsignrank(alpha / 2, n) #14 372 qsignrank (1 - alpha / 2 , n) + 1 # 64 373 V [14] # 46. 085 374 V [65] # 50.235 5. Testování p-hodnotou Príslušnou p-hodnotu vypočítáme pomocí vzorce 2 min{Pr(5£; < s e) , Pt(Se > se)}- Zde si uvědomme, že realizace testovací statistiky se = 15.5. Zároveň Se je diskrétni náhodná veličina. Z vlastností pravděpodobnostní funkce diskrétních náhodných veličin víme, že ~Pt(Se < 15.5) = Pr(5e < 15) a Pr(SE > 15.5) = Pr(SE > 16) = 1 - Pr(SE < 16) = Pr(SE < 15). Viz kapitola ??. • p-hodnota p-hodnota = 2min{Pr(5£; < se) , Pt(Se > se)} = 2mm{Pľ(SE < 15), 1 - Pr(SE < 15)} = 2min{0.03198242, 0.9680176} = 2 x 0.03198242 = 0.06396484 = 0.06396 375 2 * min (psignrank(15, m), 1 - psignrank(15, m)) # 0.06396484 • Závěr testování Protože p-hodnota = 0.06396 je větší než a = 0.05, Hq nezamítáme na hladině významnosti a = 0.05. 6. Interpretace výsledků Za základě všech tří typů testování nezamítáme nulovou hypotézu na hladině významnosti a = 0.05. Mezi největší délkou hlavice ramenní kosti na pravé straně u mužů raně středověké a současné populace neexistuje statisticky významný rozdíl. Ke stejnému závěru jsme dospěli také v příkladu 9.1. 7. Grafická vizualizace výsledku testování Grafické porovnání náhodného výběru s konstantou xq = 50.00 bychom provedli analogicky jako v příkladu 9.1 pomocí krabicového diagramu (viz obrázek 2). Poznámka: Wilcoxonův jednovýběrový exaktní test můžeme provést pomocí funkce wilcox.test(). Vstupními parametry budou vektor reprezentující náhodný výběr (humer.DRM), hodnota parametru xq z nulové hypotézy zadaná argumentem mu = 50.00, požadavek na výpočet hranic intervalu spolehlivosti zadaný nastavením argumentu conf.int = T, hodnota hladiny významnosti a zadaná prostřednictvím koeficientu spolehlivosti 1 — a nastavením hodnoty argumentu conf.level = 0.95, typ zvolené alternativní hypotézy (oboustranná) zadaný pomocí argumentu alternativě = 'two.sideď a argument correct = F, který zakazuje aplikování spojité korekce na výsledné intervaly spolehlivosti a p-hodnotu. 376 wilcox.test(humer.DRM, mu = 50.0, conf.int = T, conf.level = 0.95, 377 alternativě = 'two.sided', correct = F) # interpolovaný IS i p-hodnota Součástí výstupu je hodnota mediánu náhodného výběru medián of x = 48.79584, hodnota testovací statistiky V = 15.5, interpolované hranice 95% Waldova empirického oboustranného intervalu spolehlivosti 46.085 a 50.235 a interpolovaná p-hodnota p-value = 0.07108 adaptovaná na situaci, kdy jsou některé rozdíly \Xi —xq \ shodné (viz tabulka 11, kde ix4 —xq\ = \Xr —Íq| = 1.27). Jediné, co musíme stanovit zvlášť, jsou dolní a horní hranice kritického oboru. 57 Wilcoxon signed rank test data: humer.DRM V = 15.5, p-value = 0.06515 alternative hypothesis: true location is not equal to 50 95 percent confidence interval: 46.15993 50.05504 sample estimates: (pseudo)median 48 .79584 378 379 380 381 382 383 384 385 386 387 388 Poznámka: Výstup funkce wilcox.test() v tomto případě provází varovná hláška Warning messages: In wilcox.test.default (humer.DRM, mu = 50, conf.int = T, conf.level = 0.95,: cannot compute exact p-value with ties. Tato hláška nás upozorňuje právě na výskyt duplicitních rozdílů X4 — xq\ = \X? — xq\ = 1.27 a na použití modifikovaného postupu na výpočet p-hodnoty. 58 Příklad 9.13. Wilcoxonův jednovýběrový exaktní test Mějme datový soubor 15-anova-means-skull.txt a proměnnou upface.H popisující výšku horní části tváře mužů německé populace (viz sekce ??). Dále máme k dispozici údaje o výšce horní části tváře mužů Cernjachovské kultury na území dnešní Ukrajiny (mck = 70.00 mm, nck = 99). Na hladině významnosti a = 0.10 testujte hypotézu, že výška horní části tváře německé mužské populace je menší nebo rovna výšce horní části tváře mužské populace Cernj achovské kultury. Řešení příkladu 9.13 Zadání příkladu je shodné se zadáním příkladu 9.2. Datový soubor obsahující údaje o výšce horní části tváře 19 mužů německé populace nepochází z normálního rozdělení (p-hodnota = 0.0419, viz obrázek 9.13), proto na prozkoumání otázky ze zadání použijeme neparametrický jednovýběrový test. Opět bychom chtěli nulovou hypotézu otestovat pomocí Wilcoxonova jednovýběrového testu, k čemuž je potřeba ověřit předpoklad symetrického rozdělení naměřených hodnot okolo mediánu. Ten ověříme pomocí testu symetrie (a = 0.05). Graficky vizualizujeme míru symetrie naměřených hodnot okolo mediánu pomocí histogramu superponovaného jádrovým odhadem křivky hustoty (viz kapitola ??) se zvýrazněnou hodnotou mediánu pomocí vertikální přímky (viz obrázek 24). Obrázek 24: Kvantilový diagram (vlevo) a histogram naměřených hodnot výšky horní části tváře u mužů německé populace superponovaný křivkou jadrového odhadu hustoty (vpravo) Protože p-hodnota = 0.0682 je větší než 0.05, hypotézu o symetrickém rozdělení náhodného výběru výšek horní části tváře mužů okolo mediánu nezamítáme na hladině významnosti a = 0.05. Z grafu 24 není symetrie okolo mediánu příliš zjevná, hodnoty jsou vyšikmené směrem doprava. U náhodného výběru většího rozsahu by to byl problém, vhledem k nízkému rozsahu náhodného výběru je však test symetrie schovívavější a hypotéza o symetrii nebyla ještě zamítnuta. Náhodný výběr výšek horní části tváře u mužů považujeme tedy za symetricky rozdělený okolo mediánu. Protože hypotéza o symetrii data okolo mediánu nebyla zamítnuta, můžeme k otestování otázky ze zadání použít Wilcoxonův jednovýběrový parametrický test. Analogicky jako v příkadu 9.2 testujeme (nulovou) hypotézu v posloupnosti sedmi kroků. 1. Stanovení hypotéz • slovní formulace nulové a alternativní hypotézy Hq : Medián výšky horní části tváře německé mužské populace je menší nebo roven mediánu výšky horní části tváře mužské populace z Cernjachovské kultury. Hi : Medián výšky horní části tváře německé mužské populace je větší než medián výšky horní části tváře mužské populace z Cernjachovské kultury. • matematická formulace nulové a alternativní hypotézy Hq : x < xq, kde xq = 70.00 Hi : x, > xq, kde xq = 70.00 (pravostranná alternativa) 59 2. Volba hladiny významnosti • Hladinu významnosti volíme v souladu se zadáním jako a = 0.10. 3. Testování kritickým oborem • Testovací statistika Nejprve vypočítáme vektor rozdílů naměřených hodnot Xi a konstanty xq, tj. Xi — xq, a stanovíme, které rozdíly jsou kladné. Dále vypočítáme absolutní hodnoty těchto rozdílů, tj. \Xi — xq\ a stanovíme pořadí rozdílů v absolutních hodnotách (viz tabulka 12). Tabulka 12: Naměřené hodnoty Xi, rozdíly Xi — xq, znaménka těchto rozdílů, absotulní hodnoty rozdílů \Xi — xq\ a pořadí rozdílů v absolutních hodnotách Ri měření l 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 73 73 67 75 70.0 62 76.0 73 71.0 66 76.0 73 74 74 70.0 76.0 72 69.0 76.0 Xi — Xq t 3 3 -3 5 0.0 -8 6.0 3 1.0 -4 6.0 3 4 4 0.0 6.0 2 -1.0 6.0 1 \Xl - x0\ + 3 + 3 3 + 5 _ 8 + 6.0 + 3 + 1.0 4 + 6.0 + 3 + 4 + 4 _ + 6.0 + 2 1.0 + 6.0 Ri 6 6 6 12 - 17 14.5 6 1.5 10 14.5 6 10 10 14.5 - 3 1.5 14.5 Z tabulky 12vidíme, že z celkového počtu n = 19 rozdílů jsou dva rozdíly nulové, a tedy počet nenulových rozdílů m = 17. Dále vidíme, že celkem 13 rozdílů Xi —xq, i = 1,..., 17 je kladných, přičemž pořadí těchto rozdílů v absolutních hodnotách jsou 6, 6,12,14.5, 6,1.5,14.5, 6,10,10,14.5, 3,14.5. Hodnota testovací statistiky Se bude tedy rovná součtu těchto pořadí. 17 SE = RíI(Xí - ío > 0) = 6 + 6 + 12 + 14.5 + 6 + 1.5 + 14.5 + 6 + 10 + 10 + 14.5 + 3 + 14.5 = 118.5 i=l 389 # upface.HNI <- upface.HN 390 xO <- 70.00 391 upface.HN <- upface.HN[-which(upf 392 m <- length(upface.HN) 393 394 I <- (upface.HN - xO > 0) 395 Ri <- rank(abs(upface.HN-xO)) 396 tab <- rbind("Xi" = upface.HN, 397 "I" = I, "|Xi-xO|" = 398 "Ri" = Ri) 399 tab <- data.frame (tab) 400 nanes (tab) <- 1 : 17 401 # 1 2 3 4 5 6 7 402 # Xi 73 73 67 75 62 76 0 73 403 # Xi- xO 3 3 -3 5 -8 6 0 3 404 # I 1 1 0 1 0 1 0 1 405 # IXi -xO / 3 3 3 5 8 6 0 3 406 # Ri 6 6 6 12 17 H 5 6 407 SE <- sum(Ri * I) # 118.5 ace.HN == xO)] Xi-xO" = upface.HN - xO, abs(upface.HN - xO), 8 9 10 11 12 13 H 15 16 17 71 0 66 76. 0 73 74 74 76. 0 72 69. 0 76. 0 1 0 -4 6. 0 3 4 4 6. 0 2 -1. 0 6. 0 1 0 0 1. 0 1 1 1 1 . 0 1 0. 0 1 . 0 1 0 4 6. 0 3 4 4 6. 0 2 1 . 0 6. 0 1 5 10 14.5 6 10 10 14. 5 3 1 . 5 14.5 • Kritický obor 60 w oo; sm(a/2) - 1) U (sm(l - a/2); oo) ■oo; si7(0.10/2) - 1) U (s17(l - 0.10/2); oo) oo; si7(0.05) - 1) U (si7(0.95) ; oo) oo; 42 - 1) U (111; oo) oo; 41) U (111; oo) 408 m <- sum (upf ace . HN - xO != 0) 409 alpha <- 0.10 410 qsignrank(alpha / 2, m) - 1 # 41 411 qsignrank(1 - alpha / 2, m) # 111 • Závěr testování Protože realizace testovací statistiky sE = 118.5 náleží do kritického oboru, tj. sE £ W, Hq zamítáme na hladině významnosti a = 0.10. 4. Testování intervalem spolehlivosti • Interval spolehlivosti Abychom mohli stanovit hranice 90 % levostranného intervalu spolehlivosti, musíme nejprve vypočítat m(™+1) Walshových průměrů x'+xi; i = l,...,m. Walshovy průměry získáme pomocí funkce owa z knihovny NSM3. Prvním argumentem této funkce bude vektor nul o délce m = 17, druhým argumentem bude vektor naměřených výšek horní části tváře upface.HN. Získané Walschovy průměry seřadíme příkazem sort(). 412 walsh <- NSM3::owa(rep(0, m) , upface.HN)$owa 413 V <- sort(walsh) Hranice intervalu spolehlivosti potom tvoří hodnota, která se v seřazeném vektoru Walshových průměrů V nachází na (sm(a))-té pozici a nekonečno. Hodnotu sm(a) nalezneme pomocí funkce qsignrankQ. (71; 00) 414 alpha <- 0.10 415 qsignrank(alpha , m) # 49 416 V[49] • Závěr testování Protože xq = 70.00 nenáleží do 90% empirického levostranného intervalu spolehlivosti, tj. xq = 70.00 ^ IS, Hq zamítáme na hladině významnosti a = 0.10. 61 Testování p-hodnotou Příslušnou p-hodnotu vypočítáme pomocí vzorce ~Pt(Se > se)}- Zde si uvědomme, že realizace testovací statistiky se = 118.5. Zároveň Se je diskrétní náhodná veličina. Z vlastností pravděpodobnostní funkce diskrétních náhodných veličin víme, že Pr(SE > 118.5) = Pr(SE > 119) = 1 - Pr(SE < 119) = 1 - Pr(SE < 118). Viz kapitola ??. • p-hodnota p-hodnota = Pr(SE > sE) = Pľ(SE > 118.5) = 1 - Pr(SE < 118) = 1 - 0.9776154 = 0.02238464 = 0.02238 417 1 - psignrank(118, m) # 0.02238464 • Závěr testování Protože p-hodnota = 0.02238 je menší než a = 0.10, Hq zamítáme na hladině významnosti a = 0.10. 6. Interpretace výsledků Za základě všech tří typů testování zamítáme nulovou hypotézu na hladině významnosti a = 0.10. Výška horní části tváře mužů německé populace je statisticky významně větší než výška horní části tváře mužů Cernjachovské populace z oblasti dnešní Ukrajiny. Ke stejnému závěru jsme dospěli také v příkladu 9.2. 7. Grafická vizualizace výsledku testování Významný rozdíl horní části tváře mužů mezi oběma populacemi vizualizujeme analogicky jako v příkladu 9.2 krabicovým diagramem (viz obrázek 4). Poznámka: K provedení Wilcoxonova jednovýběrového exaktního testu použijeme funkci wilcox.test(). Vstupními parametry budou původní vektor délky n = 19 reprezentující náhodný výběr (upface.HN), hodnota parametru xq z nulové hypotézy (mu = 70.00), argument conf.int = T zadávající požadavek na výpočet intervalu spolehlivosti, hodnota hladiny významnosti a zadaná prostřednictvím koeficientu spolehlivosti l — a (conf.level = 0.90), typ zvolené alternativní hypotézy (pravostranná, alternativě = 'greater') a argument correct = F, který zakazuje aplikování spojité korekce na výsledné intervaly spolehlivosti a p-hodnotu. 418 data <- read.delim('00-Data//15-anova-means-skuli.txt ' ) 419 upface.HN <- data[data$pop == 'nem', 'upface.H'] 420 upface.HN <- na.omit(upface.HN) 421 wilcox.test(upface.HN, mu = 70.00, conf.int = T, conf.level = 0.90, alt = 'g', 422 correct = F) # IS i p-hdonota interpolované Wilcoxon signed rank test data: upface.HN V = 118.5, p-value = 0.02286 alternative hypothesis: true location is greater than 70 90 percent confidence interval: 71.00003 Inf sample estimates: (pseudo)median 72.99994 423 424 425 426 427 428 429 430 431 432 433 Součástí výstupu je hodnota mediánu náhodného výběru medián of x = 72.99994, hodnota testovací statistiky s = 118.5, interpolované hranice 90% empirického levostranného intervalu spolehlivosti 71.00003 a Inf a interpolovaná p-hodnota p-value = 0.02286. Jediné, co musíme stanovit zvlášť, je dolní hranice kritického oboru. ★ 62 Příklad 9.14. Wilcoxonův jednovýběrový exaktní test Mějme datový soubor 21-goldman-measures.csv a proměnnou tibia.LR popisující největší délku lýtkové kosti z pravé strany v mm u skeletů z období neolitu z oblasti Yoshigo (viz sekce ??). Dále máme k dispozici údaje ze studie (Hasegawa et al.) z roku 2009, v rámci které byla měřena největší délka lýtkové kosti z pravé strany u žen současné japonské populace (uif = 329.40mm, Sf = 17.3mm, rif = 342). Na hladině významnosti a = 0.01 zjistěte, zdaje největší délka lýtkové kosti z pravé strany u žen z neolitické japonské populace významně menší než u žen současné japonské populace. Řešení příkladu 9.14 Zadání příkladu je shodné se zadáním příkladu 9.3. Datový soubor obsahující údaje o největší délce lýtkové kosti z pravé strany nepochází z normálního rozdělení (p-hodnota = 0.034, viz obrázek 9.3), proto na prozkoumání otázky ze zadání použijeme neparametrický jednovýběrový test. Opět bychom chtěli nulovou hypotézu otestovat pomocí Wilcoxonova jednovýběrového testu, k čemuž je potřeba ověřit předpoklad symetrického rozdělení naměřených hodnot okolo mediánu. Ten ověříme pomocí testu symetrie (a = 0.05) a histogramem superponovaným jádrovým odhadem křivky hustoty se zvýrazněnou hodnotou mediánu pomocí vertikální přímky (viz obrázek 25). teoreticky kvantu délka lýtkové kosti (v mm) Shapiruv test: p-hodnota = 0.034 Test symetrie: p-hodnota = 0.0357 Obrázek 25: Kvantilový diagram (vlevo) a histogram naměřených hodnot délky lýtkové kosti na pravé straně u skeletů žen populace z oblasti Youshigo Shell Mound superponovaný křivkou jadrového odhadu hustoty (vpravo) Protože p-hodnota = 0.0357 je menší než 0.05, zamítáme hypotézu o symetrii naměřených hodnot okolo mediánu na hladině významnosti a = 0.05. Z grafu 25 vidíme, že data jsou nesymetrická okolo mediánu, a to natolik že i přes nízký rozsah náhodného výběru dochází k zamítnutí hypotézy o symetrii. Naměřené hodnoty největší délky lýtkové kosti z pravé strany u žen nejsou symetricky rozdělené okolo mediánu. Předpoklad pro použití Wilcoxonova testu není splněn a proto tento test nemůžeme v tomto případě použít. Musíme se tedy spokojit s výsledkem znaménkového testu použitého v rámci příkladu 9.3. 63 9.5 Wilcoxonův jednovýběrový asymptotický test Pro náhodné výběry o rozsazích n > 30 máme možnost použít k otestování nulové hypotézy asymptotickou variantu testu. Tuto variantu nazýváme Wilcoxonův jednovýběrový asymptotický test o mediánu x. Testovací statistika asymptotického testu má tvar o _ m{m+l) S a = i~ (9.4) 24 kde Se je statistika definovaná vztahem 9.3 a m je počet nenulových rozdílů Xi — xq. Za platnosti nulové hypotézy pochází statistika S a ze standardizovaného normálního rozdělení, tj. o _ m(m+l) SA = 4 = ~° ÍV(0,1). ' m(m+í)(2m+í) 24 Kritický obor podle zvolené alternativní hypotézy má tvar H n :ž^í0 W = (-oo ; ua/2) U («i_„/2 5 oo) H12 : x, > x0 W = (mi_„ ; oo) íři3 : x < Xq W = (—oo ; ua) kde ua/2, Wi-a/2i ua, wi-a jsou kvantily standardizovaného normálního rozdělení, jejichž hodnoty získáme pomocí Qt a implementované funkce qnorm(). Interval spolehlivosti má podle zvolené alternativní hypotézy jeden z následujících tvarů Hn:x^ x0 (d, h) = (v(Cí-°; V^"/^) HV2:x> x0 (d, oo) = (V^1-^ ; oo) íři3 : í < x0 (-oo, h) = (-oo ; V^Ca)) i i _ m(m+l) m(m+l)(2m+l) ^ _ m(m+l) m(m+l)(2m+l) ^ _ m(m+l) K(le L^i_a/2 — -j--ul-a/2\]-2i-' °«/2 — -4--Ua/2\J-24-' ~ -4-- ui_a^m{m+1}2^±^, Ca = ^±±i - ua^m{m+1,^žm±^, <■■■< V^21^) značí posloupnost vzestupně seřazených T"(™+1) Walshových průměrů (x<-+x^ ; i = 1,...,™, j = l,...,m, j < i a V^k\ značí fc-tý seřazený Walshův průměr. Posloupnost Walshových průměrů získáme příkazem owa z knihovny NSM3. p-hodnota má v závislosti na zvolené alternativní hypotéze jeden z následujících tvarů Hn : x, ^ xq p-hodnota = 2min{Pr(5A < s a) , Pí(5a > s a)} H12 : x > xq p-hodnota = Pr(5U > «a) His : x, < xq p-hodnota = Pr(5U < «a) kde S a je náhodná veličina, s a je realizace testovací statistiky S a (viz vzorec 9.4), tedy konkrétní číslo, a Pr(5^ < sa) je distribuční funkce standardizovaného normálního rozdělení, jejíž hodnotu získáme pomocí <® a implementované funkce pnormQ. 64 Příklad 9.15. Wilcoxonův jednovýběrový asymptotický test Mějme datový soubor 18-more-samples-variances-clavicle.txt a proměnnou cla.L popisující největší délku klíční kosti z pravé strany v mm (viz sekce ??). Dále máme k dispozici údaje o největší délce klíční kosti z pravé strany mužů z populace severní Indie (mj = 148.0mm, sj = 8.60mm, rij = 260). Na hladině významnosti a = 0.01 testujte hypotézu o shodě délky klíční kosti z pravé strany u mužů indické populace z Varanasi a u mužů ze severní Indie. Řešení příkladu 9.15 Zadání příkladu je shodné se zadáním příkladu ??. Datový soubor obsahující údaje o největší délce klíční kosti z pravé strany mužů indické populace z Varanasi nepochází z normálního rozdělení (p-hodnota = 0.0274, viz obrázek ??), proto na prozkoumání otázky ze zadání použijeme neparametrický jednovýběrový test. Nulovou hypotézu bychom chtěli otestovat pomocí Wilcoxonova jednovýběrového testu, k čemuž je potřeba ověřit předpoklad symetrického rozdělení naměřených hodnot okolo mediánu. Ten ověříme pomocí testu symetrie (a = 0.05). Graficky vizualizujeme míru symetrie naměřených hodnot okolo mediánu pomocí histogramu superponovaného jádrovým odhadem křivky hustoty (viz kapitola ??) se zvýrazněnou hodnotou mediánu pomocí vertikální přímky (viz obrázek ??). 434 data <- read.delim('00-Data//18-more-samples-variances-clavicle.txt') 435 cla.LV <- data[data$pop == 'ind2' , 'cla.L'] 436 cla.LV <- na.omit(cla.LV) 437 n <- length(cla.LV) Obrázek 26: Míra symetrie naměřených hodnot délky klíční kosti z pravé strany u mužů indické populace z Varanasi okolo mediánu 438 wilcox.test(cla.LV, mu = 148, alt = 'ť, conf.level = 0.95, conf.int = T, correct = T) ★ 65 Wilcoxon signed rank test with continuity correction data: cla.LV V = 428, p-value = 1.805e-08 alternative hypothesis: true location is not equal to 148 95 percent confidence interval: 139.0000 142.5001 sample estimates: (pseudo)median 140.75 439 440 441 442 443 444 445 446 447 448 449 66 Příklad 9.16. Wilcoxonův jednovýběrový asymptotický test Řešení příkladu 9.16 450 data <- read.delim('00-Data//19-more-samples-correlations - skull.txt') 451 nose.BM <- data[data$pop == 'mal', 'nose.B'] 452 nose.BM <- na.omit(nose.BM) 453 n <- length(nose.BM) 21.55 24.15 26.75 29.35 sirka nosu (v mm) Test symetrie: p-hodnota = 0.7182 Obrázek 27: Míra symetrie naměřených hodnot šířky nosu mužů malajské populace okolo mediánu 454 wilcox . test (nose . BM , mu = 148, alt = 'less', exact = F) ★ 67 455 Wilcoxon signed rank test with continuity correction 456 457 data: nose.BH 458 V = 0, p-value = 5.003e -14 459 alternative hypothesis: true location is less than 148 460 68 9.6 Wilcoxonův párový test Nechť (Xi, Yi)T ... (Xn, Yn)T je náhodný výběr z libovolného (ne nutně normálního) dvourozměrného rozdělení Nechť dále Z\,..., Zn, n > 2 je náhodný výběr rozdílů X — Y, tj. Z = (Z\,..., Zn)T, kde Z{ = X{ — Y{, i = 1,..., n, a nechť tento náhodný výběr pochází ze spojitého rozdělení s hustotou f(x), která je symetrická okolo mediánu ž. Konečně, nechť ž0 je konstanta. Na hladině významnosti a testujeme jednu z následujících tří hypotéz oproti příslušné alternativní hypotéze. kde ž je medián rozdílů Z\,..., Zn a ž0 je konstanta, jejíž hodnotu nejčastěji volíme jako ž0 = 0. Tato volba odpovídá hypotéze, že rozdíl mezi mediány náhodných veličin X &Y neexistuje (resp. hypotéze, že medián náhodné veličiny X je menší, resp. větší, než medián náhodné veličiny Y). Vzhledem k tomu, že jde finálně o situaci, kdy medián ž porovnáváme s konstantou žoj testujeme hypotézy o rozdílu mediánů X — Y pomocí exaktní nebo asymptotické varianty Wilcoxonova jednovýběrového testu, analogicky jako je uvedeno v sekcích 9.4 a 9.5. Výše popsaný test, v rámci kterého převádíme problém porovnávání mediánů dvou náhodných veličin X a Y na problém srovnávání mediánu jejich rozdílů Z s konstantou ž0 = 0 a následně jej řešíme pomocí exaktní resp. asymptotické varianty Wilcoxonova jednovýběrového testu, nazýváme párový Wilcoxonův test. H01 : ž = ž0 H02 : ž < ž0 H03 : z > ž0 oproti oproti oproti Hu : ž ^ Žq (oboustranná alt.) H12 : ž > Žq (pravostranná alt.) His : ž < Žq (levostranná alt.) 69 Příklad 9.17. Wilcoxonův párový exaktní test Řešení příkladu 9.17 461 data <- read.delim('00-Data//21-goldman-measures.csv' , sep = ';') 462 # head(data) 463 data.HLF <- data [data$pop == 1Poundbury (English Roman)' & data$sex == 'f' , c('LHML1 RHML ' ) ] Error in '[.data.frame'(data, data$pop selected 1 Poundburyu(EnglishuRoman)" & : undefined columns 464 465 data.HLF <- na.omit(data.HLF) 70 Error in na.omit(data.HLF): object 5 data.HLF5 not found length.HRLF <- data.HLF$RHML - data.HLF$LHML Error in evalCexpr, envir, enclos): object 5 data.HLF5 not found n <- length(length.HRLF) # 26 Error in evalCexpr, envir, enclos): object 5 length.HRLF5 not found Error in hist(length.HRLF, prob = T, breaks = seq(-10, 20, by = 5), main = "", : object 5 length. HRLF5 not found Error in boxCbty = "o"): plot.new has not been called yet Error in axis(l, seq(-7.5, 17.5, by = 5)): plot.new has not been called yet Error in axis(2, las = 1): plot.new has not been called yet Error in density(length.HRLF): object 3 length.HRLF3 not found Error in mean (length . HRLF ) : object 3 length . HRLF 3 not found Error in mtext (" rozdiludelkyuiiaupraveuauleveuStraiieu Cvumm) " , side = 1, : plot.new has not been called yet Error in na.omit(x): object 31ength.HRLF3 not found Error in mtext(bquote(paste("Testusymetrie:up_hodnotau=u". •(round(as.p, : plot.new has not been called yet Error in strwidth(legend, units = "user", cex = cex, font = text.font): plot.new has not been called yet uilcox.test(length.HRLF , mu = 0, alt = 'two.sided5, exact = F) Error in wilcox . test (length . HRLF , mu = 0, alt = "two.sided", exact = F) : object 3 length . HRLF 3 not f ound 71 Příklad 9.18. Wilcoxonův párový exaktní test Řešení příkladu 9.18 483 data <- read.delim('00-Data//21-goldman-measures.csv' , sep = ';') 484 # head(data) 485 data.ARLM <- data[data$pop == 'Dynastic Egyptian, El Hesa' & data$sex == 'm', c('LAcH: 'RAcH ' )] Error in '[.data.frame'(data, data$pop == selected 1 DynasticuEgyptian,uEluHesa" undefined columns 486 487 data.ARLM <- na.omit(data.ARLM) 72 Error in na.omit(data.ARLM): object 3 data.ARLM3 not found length.ARLM <- data.ARLM$RAcH - data.ARLM$LAcH Error in eval(expr, envir, enclos): object 3 data.ARLM3 not found n <- length(length.ARLM) # 18 Error in evalCexpr, envir, enclos): object 3 length.ARLM3 not found Error in hist(length.ARLM, prob = T, breaks = seq(-l, 5.75, by = 1.15), : object 3 length.ARLM3 not f ound Error in boxCbty = "o"): plot.new has not been called yet Error in axis(l, seq(-0.425, 5.175, by = 1.15)): plot.new has not been called yet Error in axis(2, las = 1): plot.new has not been called yet Error in density(length.ARLM): object 3 length.ARLM3 not found Error in mean(length.ARLM): object 3 length.ARLM3 not found Error in mt ext ( " rozdil ua • u vy sky unauprave ua-ule ve u str ane u ( vum) " , side = 1, : plot.new has not been called yet Error in na.omit(x): object 31ength.ARLM3 not found Error in mtext(bquote(paste("Testusymetrie:up_hodnotau=u". •(round(as.p, : plot.new has not been called yet Error in strwidth(legend, units = "user", cex = cex, font = text.font): plot.new has not been called yet uilcox.test(length.ARLM, mu = 0, alt = 'greater5, exact = F) Error in wilcox . test (length. ARLM , mu = 0, alt = "greater", exact = F): object 3 length. ARLM3 not f ound 73 Příklad 9.19. Wilcoxonův párový exaktní test Řešení příkladu 9.19 505 data <- read.delim('00-Data//21-goldman-measures.csv' , sep = ';') 506 # head(data) 507 data.TRLM <- data[data$pop == 'Cliff Dweller' & data$sex == 'm', c('LTML' , 'RTML') ] Error in '[.data.frame'(data, data$pop == "CliffuDweller" & data$sex == undefined columns selected 509 data.TRLM <- na.omit(data.TRLM) * 74 Error in na.omit(data.TRLM): object 3 data.TRLM3 not found 510 511 length.TRLM <- data.TRLM$RTML - data.TRLM$LTML Error in evalCexpr, envir, enclos): object 3 data.TRLM3 not found 512 513 n <- length(length.TRLM) # 20 Error in evalCexpr, envir, enclos): object 3 length.TRLM3 not found 514 Error in hist(length.TRLM, prob = T, breaks = seq(-8.5, 30, by = 7.7) f ound object 3 length.TRLM3 not 515 Error in boxCbty = "o"): plot.new has not been called yet 516 Error in axis(l, seq(-4.65, 26.15, by = 7.7)): plot.new has not been called yet 517 Error in axis(2, las = 1): plot.new has not been called yet 518 Error in density(length.TRLM): object 3 length.TRLM3 not found 519 Error in mean(length.TRLM): object 3 length.TRLM3 not found 520 Error in mtext ( " rozdiludelkyuiiaupraveuauleveuStraiieu Cvumm) " , side = 1. called yet plot.new has not been 521 Error in na.omit(x): object 3 1ength.TRLM3 not found 522 Error in mtext(bquote(paste("Testusymetrie:up-hodnotau = u" called yet (round(as.p, : plot.new has not been 523 Error in strwidth(legend , units = called yet 'user", cex = cex, font = text.font): plot.new has not been 524 525 uilcox.test(length.TRLM, mu = 0, alt = 'less5, exact = F) Error in wilcox . test (length. TRLM , mu = 0, alt = "less", exact = F) : object 3 length. TRLM 3 not found 526 75 Příklad 9.20. Wilcoxonův párový asymptotický test Řešení příkladu 9.20 527 data <- read.delim('00-Data//02-paired-means-clavicle.txt') 528 # head(data) 529 data.l2R <- data[data$side == 'R', c('simd.l', 'simd.2')] 530 data.l2R <- na.omit(data.12R) 531 simd.l2R <- data.12R$simd.1 - data.12R$simd.2 532 n <- length(simd . 12R) # 40 T T T T median T T" -0.285 0.055 0.395 rozdil prvního a druhého mereni (v mm) Test symetrie: p-hodnota = 0.4316 Obrázek 28: Míra symetrie rozdílu prvního a druhého měření vertikálního průměru ve středu délky těla klíční kosti na pravé straně získaných jedním výzkumníkem 533 wilcox.test(simd.12R, mu = 0, alt = 'two.sided', exact = F) ★ 76 4489999999945 Wilcoxon signed rank test with continuity correction data: simd.12R V = 366.5, p-value = 0.7482 alternative hypothesis: true location is not equal to 0 534 535 536 537 538 539 77 Příklad 9.21. Wilcoxonův párový asymptotický test Řešení příkladu 9.21 540 data <- read.delim('00-Data//02-paired-means-clavicle.txt') 541 # head(data) 542 data.l2L <- data[data$side == 'L', c('simd.l', 'simd.2')] 543 data.l2L <- na.omit(data.12L) 544 simd.l2L <- data.12L$simd . 1 - data.12L$simd . 2 545 simd.l2L <- simd.12L [-36] 546 n <- length(simd.12L) # 39 -0.445 -0.025 0.395 rozdil prvního a druhého mereni (v mm) Test symetrie: p-hodnota = 0.5608 Obrázek 29: Míra symetrie rozdílu prvního a druhého měření vertikálního průměru ve středu délky těla klíční kosti na levé straně získaných jedním výzkumníkem po odstranvení odlehlé hodnoty 547 wilcox.test(simd.12L, mu = 0, alt = 'greater', exact = F) ★ 78 Wilcoxon signed rank test with continuity correction data: simd.12L V = 442, p-value = 0.08722 alternative hypothesis: true location is greater than 0 548 549 550 551 552 553 79 Příklad 9.22. Wilcoxonův párový asymptotický test Řešení příkladu 9.22 554 data <- read.delim('00-Data//03-paired-means-clavicle2.txt') 555 # head(data) 556 data.RL <- data[data$sex == 'f', c('length.R ' , 'length.L')] 557 data.RL <- na.omit(data.RL) 558 length.RL <- data.RL$length.R - data.RL$length.L 559 n <- length(length.RL) # 50 -7 -5-3-11 3 5 rozdil delek na pravé a leve strane (v mm) Test symetrie: p-hodnota = 0.1745 Obrázek 30: Míra symetrie rozdílu délky těla klíčních kostí na levé a pravé strané žen 560 wilcox.test(length.RL, mu = 0, alt = 'less', exact = F) ★ 80 Wilcoxon signed rank test with continuity correction data: length.RL V = 125.5, p-value = 0.039 alternative hypothesis: true location is less than 0 561 562 563 564 565 566 81 9.7 Pořadový exaktní test o nezávislosti Nechť (Xi, Yi)T ... (Xn, Yn)T je náhodný výběr z dvourozměrného rozdělení. Na hladině významnosti a testujeme jednu z následujících tří hypotéz oproti příslušné alternativní hypotéze. Hqi : p = 0 oproti Hu : p ^ 0 (oboustranná alt.) í^02 : p < 0 oproti Hi2 : p > 0 (pravostranná alt.) Hq3 : p > 0 oproti H13 : p < 0 (levostranná alt.) Test nazýváme pořadovým testem o nezávislosti. Testovací statistika má tvar SE = Rs = l- V(J?2 - Qx)2, (9.5) n(nz — 1) -f—^ kde n je rozsah náhodného výběru, Ri je pořadí náhodné veličiny Xi a Qi je pořadí nádhodné veličiny Y{. Všimněme si, že statistika Se je současně Spearmanovým korelačním koeficientem R$, resp. Spearmanův korelační koeficient je používán jako testovací statistika pořadového testu o nezávislosti. Kritický obor podle zvolené alternativní hypotézy má tvar Hn : p = 0 W=(-l; r„(a/2)> U (r„(l - a/2); 1) H12:p>Q W=(-l;rn(a)) H13:p<0 W= (r„(l-a); 1) kde rIj(a/2), r„(l — a/2), r„(l — a) a r„(l — a/2) jsou kvantily pro pořadový test o nezávislosti, jejichž hodnoty získáme pomocí softwaru CĚt a funkce qSpearman() z knihovny SuppDists. Tvary intervalů spolehlivosti si pro pořadový test o nezávislosti neuvádíme, p-hodnota má v závislosti na zvolené alternativní hypotéze jeden z následujících tvarů Hu : x ^ xq p-hodnota = 2min{Pr(5£; < se) , Pr(SE > se)} H12 : x, > xq p-hodnota = Pr(SE > se) H13 : x, < xq p-hodnota = Pr(SE < se) kde Se je náhodná veličina, sE je realizace testovací statistiky Se (viz vzorec 9.5), tedy konkrétní číslo, Vt(Se > se) = 1 — Pt(Se < se) = 1 — Pt(Se < se _ l)j coz vyplývá z faktu, že náhodná veličina Se se řídí tabelo-vaným (diskrétním) rozdělením pro pořadový test o nezávislosti, a Pr(SE < se) je distribuční funkce tabelovaného rozdělení pro pořadový test o nezávislosti, jejíž hodnotu získáme pomocí 9 a implementované funkce pSpearman() z knihovny SuppDists. 82 Příklad 9.23. Pořadový exaktní test o nezávislosti Řešení příkladu 9.23 567 data <- read.delim('00-Data//21-goldman-measures.csv' , sep = ';') 568 head(data) 569 data <- data [data$sex == ' m' & data$pop == 'Hawikuh', cCRAcH", ' RIBL ') ] Error in '[.data.frame'Cdata, data$sex == "m" & data$pop == "Hawikuh", : undefined columns 570 selected 571 data <- na.omit(data) 572 573 MVN::mvn(data, mvnTest = 'mardia')$multi S.'i Error in colHeans(x, na.rm = TRUE): must be numeric 574 575 # sikmost: O.O4II66 # spicatost: 0.9761901 576 MVN::mvn(data, mvnTest = 'hz')$multi # O.O436649 Error in cov(data): is.numeric(x) || is.logicalCx) is not TRUE 577 578 MVN::mvn(data, mvnTest = 'royston')$multi # 0.2507807 Warning in mean.defaultCx): argument is not numeric or logical: returning NA 579 Error in x - mean(x): non-numeric argument to binary operator 580 581 acetab.HRM <- data$RAcH 582 iblade.LRM <- data$RIBL 583 n <- length(iblade.LRM) # 14 584 alpha <- 0.05 585 rS <- cor(acetab.HRM , iblade.LRM. method = 'spearman') # 0.4856559 Error in cor(acetab.HRM, like 'x' iblade.LRM, method = "spearman"): supply both and or a matrix- 586 587 library(SuppDists) 588 qSpearman(alpha / 2, n) # -0.5252747 589 qSpearman(1 - alpha / 2, n) # 0.5384615 590 591 2 * min(pSpearman(rS, n), 1 - pSpearman(rS, n)) # 0.08218591 Error in pSpearman(rS, n): object 'rS' not found 592 84 Příklad 9.24. Pořadový exaktní test o nezávislosti Řešení příkladu 9.24 593 data <- read.delim('00-Data//21-goldman-measures.csv' , sep = ';') 594 head(data) 595 data <- data[data$sex == 'm' & data$pop == 'Cliff Dweller', c('LTML', 'LRML')] Error in '[.data.frame'(data, data$sex == selected data$pop == "CI iffuDweller" , : undefined columns 596 597 data <- na.omit(data) 598 MVN::mvn(data, mvnTest 'mardia')$mult i 85 Error in colHeans(x, na.rm = TRUE): must be numeric 599 600 # sikmost: 0.000229396 # spicatost : 0.0056161 601 MVN::mvn(data, mvnTest = 'hz')$multi # 0.009711723 Error in cov(data): is.numeric(x) || is.logicalCx) is not TRUE 602 603 MVN::mvn(data, mvnTest = 'royston')$multi # 0.07266684 Warning in mean.defaultCx): argument is not numeric or logical: returning NA 604 Error in x - mean(x): non-numeric argument to binary operator 605 606 607 608 609 610 611 612 613 614 radius.LLM <- data$LRML tibia.LLM <- data$LTML n <- length(tibia.LLM) alpha <- 0.05 rS <- cor(radius.LLM # 17 tibia.LLM, method = 'spearman') # 0.8725396 library(SuppDists) qSpearmand - alpha, 1 - pSpearman (rS , n) n) # 0.4166667 # 2.981788e-06 86 9.8 Pořadový asymptotický test o nezávislosti Pro náhodný výběr o rozsahu n > 30 máme možnost použít k otestování nulové hypotézy asymptotickou variantu testu. Tuto variantu nazýváme pořadový asymptotický test o nezávislosti. Stejně jako u exaktního testu si pro asymptotický test uvedeme pouze postup testování pomocí kritického oboru a p-hodnoty. Testovací statistika asymptotického testu má tvar SA = rsVr^l (9.6) kde r s je testovací statistika definovaná vztahem 9.5 a n je rozsah náhodného výběru. Za platnosti nulové hypotézy pochází statistika S a ze standardizovaného normálního rozdělení, tj. SA = rsVrT^l ~° JV(0,1). Kritický obor podle zvolené alternativní hypotézy má tvar Hn:xjíx0 W = (-oo; ua/2) U (u1_a/2 ; oo) Hi2 : x, > x0 W = ; oo) ifi3 : x < Xq W = (—oo ; ua) kde ua/2, Ui-a/2, ua, «i-a jsou kvantily standardizovaného normálního rozdělení, jejichž hodnoty získáme pomocí Ol a implementované funkce qnorm(). p-hodnota má v závislosti na zvolené alternativní hypotéze jeden z následujících tvarů Hu : x ^ xq p-hodnota = 2min{Pr(5A < s a) , P?(Sa > sa)} H\2 : x, > xq p-hodnota = Pr(5^ > s a) His : x, < xq p-hodnota = Pr(5^ < s a) kde Sa je náhodná veličina, sa je realizace testovací statistiky Sa (viz vzorec 9.6), tedy konkrétní číslo, Pr(5^ > s a) = 1 — Pr(5^ < s a), což vyplývá z faktu, že náhodná veličina S a pochází z normálního (spojitého) rozdělení (viz kapitola ??), a Pr(5^ < s a) je distribuční funkce standardizovaného normálního rozdělení, jejíž hodnotu získáme pomocí ^ a implementované funkce pnormQ. 87 Příklad 9.25. Pořadový asymptotický test o nezávislosti Řešení příkladu 9.25 ★ 615 data <- read.delim('00-Data//16-anova-head.txt ' ) 616 head(data) 617 data <- data [data$sex == 'f', cChead.L', 'bizyg.W')] 618 data <- na.omit(data) 619 MVN::mvn(data, mvnTest = 1mardia1)$multi 620 # sikmost: 0.0427523 # spicatost: O.545548O 621 MVN::mvn(data, mvnTest = 'hz')$multi # 0.09943109 622 MVN::mvn(data, mvnTest = 'royston')$multi # 0.04040737 623 624 head.LF <- data$head.L 625 bizyg.WF <- data$bizyg.W 626 n <- length(head.LF) # 100 627 alpha <- 0.05 628 629 630 631 632 633 rS <- cor(head.LF, bizyg.WF, method = 'spearman sA <- rS * sqrt(n - 1) # 1.703704 qnorm(alpha / 2) # -1.959964 qnorm(l - alpha / 2) # 1.959964 2 * min(pnorm(sA), 1 - pnorm(sA)) # 0.08843632 ') # 0.1712287 88 Příklad 9.26. Pořadový asymptotický test o nezávislosti Řešení příkladu 9.26 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 data <- read.delim('00-Data//16-anova-head.txt1) head(data) data <- data[data$sex == 'm: data <- na.omit(data) cObigo.W , ' bizyg.W')] MVN::mvn(data, mvnTest # skimost: 0.002218132 MVN::mvn(data, mvnTest MVN::mvn(data , mvnTest = 'mardia')$multi # spicatost: 0.05247948 = 'hz')$multi # 0.006343478 = 'royston')$multi # 0.00644503 bigo.WM <- data$bigo.W bizyg.WM <- data$bizyg.W n <- length(bigo.WM) # alpha <- 0.05 rS <- cor(bigo.WM, bizyg sA <- rS * sqrt (n - 1) # qnorm(l - alpha) # 1.644854 1 - pnorm(sA) # 3.185348e-05 75 .WM, method = 3.998642 'spearman') # 0.4648327 * 89