12 Dvouvýběrové neparametrické testy 12.1 Wilcoxonův dvouvýběrový exaktní test Nechť Xn,..., Xini, ri\ > 2 je náhodný výběr ze spojitého rozdělení s distribuční funkcí F\{x) a X21,..., X2„2, U'2 > 2 je na něm nezávislý náhodný výběr ze spojitého rozdělení s distribuční funkcí F'2(x). O distribučních funkcích Fi(x) a F2(x) předpokládáme, že se navzájem liší pouze posunutím. Nechť dále í\ je medián náhodného výběru In,..., X\ni a X2 je medián náhodného výběru X21, ■ ■ ■, X2n2 a 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,\ — ÔC2 = Íq oproti Hu : í\ — X2 7^ ío (oboustranná alt.) H02 : x,\ — ÔC2 < xq oproti H12 : í\ — X2 > xq (pravostranná alt.) Hqs : x,\ — ÔC2 > xq oproti H13 : í\ — X2 < xq (levostranná alt.) Test nazýváme Wilcoxonův dvouvýběrový exaktní test o rozdílu mediánů x,\ — xi. Testovací statistika má tvar n 2 kde ri2 je rozsah druhého náhodného výběru a Sj, j = 1,..., «2, je pořadí náhodných veličin X21, ■ ■ ■, X2n2 druhého náhodného výběru v seřazeném vektoru všech ri\ +ri2 hodnot z obou náhodných výběrů. Statistika W je tedy součet pořadí hodnot X21, ■ ■ ■, ^2n2i v seřazeném vektoru všech ri\ + ri2 hodnot. Kritický obor podle zvolené alternativní hypotézy má tvar H11 : íi - x2 ^ x0 W = (-00; wnun2(a/2)) U (wnuJl2(l - a/2); 00) H12 : Xx - x2 > x0 W = {wnun,2(l - a) ; 00) H13 : xi - Í2 < x0 W = (-00; wnun,2(a)) kde wnii„2(a/2), wni^n2(l — a/2), wni^n2(a), wni^n2(l — a) jsou tabelované kvantily pro dvouvýběrový Wilcoxonův test, jejichž hodnoty získáme příkazem qwilcox(). Interval spolehlivosti má podle zvolené alternativní hypotézy jeden z následujících tvarů Hu : xx - Í2 ^ x0 (d, h) = (ř7(^i,"2(«/2)) ; í/(^1,n2(i-a/2)+i)) H12 :x1-x2>x0 (d, 00) = (f/(^i."2(«)) ; 00) H13 :x1-x2 se)} H\2 ■ x\ — Í2 > Íq p-hodnota = Yt{Se > se) = 1 — Pr(5.E < se) H13 : í\ — x,2 < xq p-hodnota = Pr(SE < se) kde Se je náhodná veličina, se je realizace testovací statistiky Se (viz vzorec 12.1), tedy konkrétní číslo, a Pt(Se < se) je distribuční funkce tabelovaného rozdělení pro Wilcoxonův dvouvýběrový exaktní test, jejíž hodnotu získáme pomocí ^1 a implementované funkce pwilcoxQ. 1 1 # ZKONTROLOVAT BARVY 2 Příklad 12.1. Wilcoxonův dvouvýběrový exaktní test Řešení příkladu 12.1 2 data <- read.delim('00-Data//18-more-samples-variances-clavicle.txt') 3 head(data) 4 cla.LIl <- data [data$pop == 'indl' , 'cla.L'] 5 cla.LIl <- na.omit(cla.LI1) 6 cla.LI2 <- data [data$pop == 'ind2' , 'cla.L'] 7 cla.LI2 <- na.omit(cla.LI2) 8 nl <- length(cla.LIl) 9 n2 <- length(cla.LI2) 10 par(mar = c(5, 4, 1, 2), family = 'Times') 11 qqnorm(cla.LI1, main = '', xlab = '', pen = 21, col = 'red', bg = 'bisque', 12 yla-b = 'vyberovy kvantil', las = 1) 13 mtext('teoreticky kvantil', side = 1, line = 2) 14 mtext(bquote(paste('Lillieforsuv test: p-hodnota = ', .(round(nortest::1 illie.test(cla. LI1)$p.val , 4) ) ) ) , 15 side = 1, line = 3.2) 16 qqline(cla.LI1 , lwd = 2, col = 'darkred ' ) 17 18 qqnorm(cla.LI2, main = '', xlab = '', pen = 21, col = 'dodgerblue3', bg = 'mintcream', 19 yla-b = 'vyberovy kvantil', las = 1) 20 mtext('teoreticky kvantil', side = 1, line = 2) 21 mtext(bquote(paste('Lillieforsuv test: p-hodnota = ', .(round(nortest::1 illie.test(cla. LI2)$p.val , 4) ) ) ) , 22 side = 1, line = 3.2) 23 qqline(cla.LI2, lwd = 2, col = 'darkblue') ★ 3 teoreticky kvantil teoreticky kvantil Lillieforsuv test: p-hodnota = 0.0956 Lillieforsuv test: p-hodnota = 0.0274 Obrázek 1: Kvantilové diagramy délky klíční kosti z levé strany u mužů indické populace z Amritsaru (vlevo) a indické populace z Varanasi (vpravo) 125 137 149 161 120 130 140 150 160 170 délka leve klicni kosti (v mm) délka leve klicni kosti (v mm) KS-test: p-hodnota = 0.919 Obrázek 2: Porovnání histogramů a graf distribučních funkcí délky klíční kosti z levé strany u mužů indické populace z Amritsaru a indické populace z Varanasi 4 24 par(mar = c(5, 4, 1, 2), family = 'Times') 25 hist(cla.LIl, prob = T, breaks = seq(122, 170, by = 6), main = '', axes = F, 26 col = rgb(0.9, 0, 0, 0.2), density = 80, xlim = c(120, 170), 27 xlab = 'delka leve klicni kosti (v mm)', ylim = c(0, 0.085), 28 yla-b = 'relativni cetnost') 29 hist(cla.LI2, prob = T, breaks = seq(124, 166, by = 6), main = '', axes = F, 30 col = rgb(0, 0, 0.9, 0.2), density = 80, add = T) 31 box(bty = 'o') 32 axis(l, seq(125, 167, by = 6)) 33 axis (2, las = 1) 34 lines(density(cla.LI2), lwd = 2, col = 'blue') 35 lines(density(cla.LI1), lwd = 2, col = 'red') 36 legend('topleft', fill = c(rgb(0.9, 0, 0, 0.2), rgb(0, 0, 0.9, 0.2)), 37 legend = c('p. z Amritsaru', 'p. z Varanasi'), bty = 'n') 38 39 plot(ecdf(cla.LI1), col = 'red', main = '', xlim = c(120, 170), ylab = 'F(x)', 40 xlab = '', las = 1, cex = 0.6) 41 plot(ecdf(cla.LI2), add = T, col = 'blue', cex = 0.6) 42 ks.val <- round(ks.test(cla.LI2-mean(cla.LI2) , cla.LI1 -(mean(cla.LI1)))$p.val, 4) 43 mtext('delka leve klicni kosti (v mm)', side = 1, line = 2.1) 44 mtext(bquote(paste('KS-test: p-hodnota = ', .(ks.val))), side = 1, line = 3.3) 45 legend('bottomright' , col = c('red' , 'blue'), lty = 1, pch = 20, 46 legend = c('p. z Amritsaru', 'p. z Varanasi'), bty = 'n') 5 Příklad 12.2. Wilcoxonův dvouvýběrový exaktní test Řešení příkladu 12.2 47 data <- read.delim('00-Data//15-anova-means - skull.txt ' ) 48 head(data) 49 upface.HN <- data [data$pop == 'nem' , 'upface.H'] 50 upface.HN <- na.omit(upface.HN) 51 upface.HC <- data[data$pop == 'cin', 'upface.H'] 52 upface.HC <- na.omit(upface.HC) 53 nl <- length(upface.HN) 54 n2 <- length(upface.HC) 55 par(mar = c(5, 4, 1, 2), family = 'Times') 56 qqnorm(upface.HN, main = '', xlab = '', pch = 21, col = 'red', bg = 'bisque', 57 yla-b = 'vyberovy kvantil', las = 1) 58 mtext('teoreticky kvantil', side = 1, line = 2) 59 mtext(bquote(paste('Shapiro-Wilkuv test: p-hodnota = ', .(round(Shapiro.test(upface.HN)$ p.val, 4)))), 60 side = 1, line = 3.2) 61 qqline(upface.HN , lwd = 2, col = 'darkred') 62 63 qqnorm(upface.HC, main = '', xlab = '', pch = 21, col = 'dodgerblue3', bg = 'mintcream', 64 yla-b = 'vyberovy kvantil', las = 1) 65 mtext('teoreticky kvantil', side = 1, line = 2) 66 mtext(bquote(paste('Shapiro-Wilkuv test: p-hodnota = ', .(round(Shapiro.test(upface.HC)$ p.val, 4)))), 67 side = 1, line = 3.2) 68 qqline(upface.HC , lwd = 2, col = 'darkblue') ★ 6 -2-1012 -2-1012 teoreticky kvantil teoreticky kvantil Shapiro-Wilkuv test: p-hodnota = 0.0419 Shapiro-Wilkuv test: p-hodnota = 0.0513 Obrázek 3: Kvantilové diagramy výšky horní části tváře mužů německé populace (vlevo) a čínské populace (vpravo) Obrázek 4: Porovnání histogramů a graf distribučních funkcí výšky horní části tváře mužů německé a čínské populace 7 69 par(mar = c(5, 4, 1, 2), family = 'Times') 70 hist(upface.HN, prob = T, breaks = seq(61, 76, by = 3), main = '', axes = F, 71 col = rgb(0.9, 0, 0, 0.2), density = 80, ylim = c(0, 0.13), xlim = c(57, 81), 72 xlab = 'vyska horni casti tvare (v mm)', 73 yla-b = 'relativni cetnost') 74 75 hist(upface.HC, prob = T, breaks = seq(59, 79, by = 4), main = '', axes = F, 76 col = rgb(0, 0, 0.9, 0.2), density = 80, add = T) 77 box(bty = 'o') 78 axis(l, seq(53.5, 80.5, by = 3)) 79 axis (2, las = 1) 80 lines(density(upface.HC), lwd = 2, col = 'blue') 81 lines(density(upface.HN), lwd = 2, col = 'red') 82 legend('topleft', fill = c(rgb(0.9, 0, 0, 0.2), rgb(0, 0, 0.9, 0.2)), 83 legend = c('nemecka p.', 'cinska p.'), bty = 'n') 84 85 plot (ecdf (upf ace . HN) , col = 'red', main = '', xlim = c (57 , 81), ylab = 'F(x)', 86 xlab = '', las = 1, cex = 0.6) 87 plot (ecdf (upf ace . HC) , add = T, col = 'blue', cex = 0.6) 88 ks.val <- round(ks.test(upface.HC - mean(upface.HC), upface.HN - (mean(upface.HN)))$p. val, 4) 89 mtext('vyska horni casti tvare', side = 1, line = 2.1) 90 mtext(bquote(paste('KS-test: p-hodnota = ', .(ks.val))), side = 1, line = 3.3) 91 legend('topleft', col = c('red', 'blue'), lty = 1, pch = 20, 92 legend = c('nemecka p.', 'cinska p.'), bty = 'n') 8 Příklad 12.3. Wilcoxonův dvouvýběrový exaktní test Řešení příkladu 12.3 93 94 95 96 97 98 99 100 data <- read head(data) intorb.BP <-intorb.BP <-intorb.BB <-intorb.BB <-nl < n2 < .delim('00-Data//19-more-samples-correlations-skull.txt') data [data$pop == 'per' na.omit(intorb.BP) data[data$pop == 'ban' na.omit(intorb.BB) length(intorb.BP) # 46 length(intorb.BB) # 14 intorb.B'] intorb.B'] 101 par(mar = c(5, 4, 1, 2), family = 'Times') 102 qqnorm(intorb.BP, main = '', xlab = '', pch = 21, col = 'red', bg = 'bisque', 103 yla-b = 'vyberovy kvantil', las = 1) 104 mtext('teoreticky kvantil', side = 1, line = 2) 105 mtext(bquote(paste('Lillieforsuv test: p-hodnota = ', .(round(nortest::1 illie.test( intorb.BP)$p.val, 4)))), 106 side = 1, line = 3.2) 107 qqline(intorb.BP , lwd = 2, col = 'darkred') 108 109 qqnorm(intorb.BB, main = '', xlab = '', pch = 21, col = 'dodgerblue3', bg = 'mintcream', 110 ylab = 'vyberovy kvantil', las = 1) 111 mtext('teoreticky kvantil', side = 1, line = 2) 112 mtext(bquote(paste('Shapiro-Wilkuv test: p-hodnota = ', .(round(Shapiro.test(intorb.BB)$ p.val, 4)))), 113 side = 1, line = 3.2) 114 qqline(intorb.BB , lwd = 2, col = 'darkblue') ★ 9 -2-10 1 2 teoreticky kvantil Lillieforsuv test: p-hodnota = 0.0023 -1 0 1 teoreticky kvantil Shapiro-Wilkuv test: p-hodnota = 0.4537 Obrázek 5: Porovnání histogramů a graf distribučních funkcí interorbitální šířky mužů peruánské a bantuské popu lace 0.30 0.25 0.20 0.15 0.10 0.05 0.00 □ peruánská p ■ bantuska p. 0h 1.0 H 0.8 0.6 -0.4 -0.2 0.0 peruánská p. bantuska p. interorbitalni sirka (v mm) "i-1-1-r 15 20 25 30 interorbitalni sirka (v mm) KS-test: p-hodnota = 0.2958 Obrázek 6: Porovnání histogramů a graf distribučních funkcí interorbitální šířky mužů peruánské a bantuské popu lace 10 115 par(mar = c(5, 4, 1, 2), family = 'Times') 116 hist(intorb.BP, prob = T, breaks = seq(17, 29, by = 2), main = '', axes = F, 117 col = rgb(0.9, 0, 0, 0.2), density = 80, xlim = c(15, 32), ylim = c(0, 0.33), 118 xlab = 'interorbitalni sirka (v mm)', 119 yla-b = 'relativni cetnost') 120 hist(intorb.BB, prob = T, breaks = seq(20, 30, by = 2), main = '', axes = F, 121 col = rgb(0, 0, 0.9, 0.2), density = 80, add = T) 122 box(bty = 'o') 123 axis(l, seq(12, 36, by = 2)) 124 axis(2, las = 1) 125 lines(density(intorb.BP), lwd = 2, col = 'red') 126 lines(density(intorb.BB), lwd = 2, col = 'blue') 127 legend('topleft', fill = c(rgb(0.9, 0, 0, 0.2), rgb(0, 0, 0.9, 0.2)), 128 legend = c('peruanska p.', 'bantuska p.'), bty = 'n') 129 130 plot (ecdf (intorb . BP) , col = 'red', main = '', xlim = c(15, 32), ylab = 'F(x)', 131 xlab = '', las = 1, cex = 0.6) 132 plot (ecdf (intorb . BB) , add = T, col = 'blue', cex = 0.6) 133 ks.val <- round(ks.test(intorb.BP - mean(intorb.BP), intorb.BB - (mean(intorb.BB)))$p. val, 4) 134 mtext('interorbitalni sirka (v mm)', side = 1, line = 2.1) 135 mtext(bquote(paste('KS-test: p-hodnota = ', .(ks.val))), side = 1, line = 3.3) 136 137 legend('topleft', col = c('red', 'blue'), lty = 1, pch = 20, 138 legend = c('peruanska p.', 'bantuska p.'), bty = 'n') 11 12.2 Wilcoxonův dvouvýběrový asymptotický test Pro náhodné výběry o rozsazích ri\ > 30 a ni > 30 máme možnost použít k otestování nulové hypotézy asymptotickou variantu testu. Tuto variantu nazýváme Wilcoxonův dvouvýběrový asymptotický test o rozdílu mediánů x\ — obi. Testovací statistika asymptotického testu má tvar a _ "2("i+"2 + l) SA = °E (12.2) 'n1Ti2("i+Ti2 + l) 12 kde Se je testovací statistika definovaná vztahem 12.1, n\ je rozsah prvního náhodného výběru, ni je rozsah druhého náhodného výběru. Za platnosti nulové hypotézy pochází statistika U a ze standardizovaného normálního rozdělení, tj- a _ "2("i+"2+l) S a = e, 2 =- ~° N(0,1). / nin2(ni+n2 + l) 12 Kritický obor podle zvolené alternativní hypotézy má tvar Hn :x1-x2^x0 W = (-oo ; ua/2) U («i_a/2 ; oo) Hn : xi - xi > x0 W = (mi_„ ; oo) íři3 : X\ — Xi < Xq W = (—oo ; ua) kde ua/i, Ui_a/i, ua, «i_a jsou kvantily standardizovaného normálního rozdělení, jejichž hodnoty získáme pomocí Oř a implementované funkce qnorm(). Interval spolehlivosti má podle zvolené alternativní hypotézy jeden z následujících tvarů Hn :xi-x2^x0 (d, h) = {U<~c"^ ; [/("i^+i-C^)) H12 ■ íi — X'i > Xq (d, oo) = ([/(c=) ; oo) Hia :x1-x2 značí vzestupně seřazené rozdíly Yj — Xi, i = 1,..., n\, j = 1,..., rig, a značí a;-tý seřazený rozdíl. p-hodnota má v závislosti na zvolené alternativní hypotéze jeden z následujících tvarů Hn : í\ — x,2 7^ xq p-hodnota = 2min{Pr(5A < s a) , ~P*(Sa > sa)} H\2 : x\ — Í2 > Íq p-hodnota = Yt{Sa > s a) = 1 — Pr(5U < s a) His : í\ — x,2 < xq p-hodnota = Pr(5^ < s a) kde S a je náhodná veličina, s a je realizace testovací statistiky S a (viz vzorec 12.2), tedy konkrétní číslo, a Pr(5U < sa) je distribuční funkce standardizovaného normálního rozdělení, jejíž hodnotu získáme pomocí ď a implementované funkce pnormQ. 12 Příklad 12.4. Wilcoxonův dvouvýběrový asymptotický test Řešení příkladu 12.4 139 data <- read.delim('00-Data//16- anova-head.txt') 140 head(data) 141 head.LMSa <- data[data$sex == 'm ' & data$sexor == ' sa ' , 'head.L'] 142 head.LMSa <- na.omit(head.LMSa) 143 head.LMOp <- data [data$sex == 'm ' & data$sexor == ' op ' , 'head.L'] 144 head.LMOp <- na.omit(head.LMOp) 145 nl <- length(head.LMSa) 146 n2 <- length(head.LMOp) 147 par(mar = c(5, 4, 1, 2), family = 'Times') 148 qqnorm(head.LMSa, main = '', xlab = '', pch = 21, col = 'red', bg = 'bisque', 149 yla-b = 'vyberovy kvantil', las = 1) 150 mtext('teoreticky kvantil', side = 1, line = 2) 151 mtext(bquote(paste('Shapiro-Wilkuv test: p-hodnota = ', .(round(Shapiro.test(head.LMSa)$ p.val, 4)))), 152 side = 1, line = 3.2) 153 qqline(head.LMSa , lwd = 2, col = 'darkred') 154 155 qqnorm(head.LMOp, main = '', xlab = '', pch = 21, col = 'dodgerblue3', bg = 'mintcream', 156 yla-b = 'vyberovy kvantil', las = 1) 157 mtext('teoreticky kvantil', side = 1, line = 2) 158 mtext(bquote(paste('Lillieforsuv test: p-hodnota = ', .(round(nortest::1 illie.test(head. LM0p)$p.val, 4)))) , 159 side = 1, line = 3.2) 160 qqline(head.LMOp , lwd = 2, col = 'darkblue') ★ 13 -1.5 -0.5 0.5 1.0 1.5 -2-10 1 2 teoreticky kvantil teoreticky kvantil Shapiro-Wilkuv test: p-hodnota = 0.4092 Lillieforsuv test: p-hodnota = 0.043 Obrázek 7: Kvantilové diagramy délky hlavy mužů orientovaných výhradně heterosexuálně (vlevo) a mužů orientovaných jinak než výhradně heterosexuálně (vpravo) 177 185 193 201 209 217 180 190 200 210 220 délka hlavy (v mm) délka hlavy (v mm) KS-test: p-hodnota = 0.9826 Obrázek 8: Porovnání histogramů a graf distribučních funkcí délky hlavy mužů orientovaných výhradně heterosexuálně a mužů orientovaných jinak než výhradně heterosexuálně 14 161 par(mar = c(5, 4, 1, 2), family = 'Times') 162 hist(head.LMSa, prob = T, breaks = seq(187, 203, by = 4), main = '', axes = F, 163 col = rgb(0, 0, 0.9, 0.2), density = 80, xlim = c(178, 217), 164 xlab = 'delka hlavy (v mm)', ylim = c(0, 0.12), 165 yla-b = 'relativni cetnost') 166 167 hist(head.LMOp, prob = T, breaks = seq(180, 215, by = 5), main = '', axes = F, 168 col = rgb(0.9, 0, 0, 0.2), density = 80, add = T) 169 box(bty = 'o') 170 axis(l, seq(169, 221, by = 4)) 171 axis(2, las = 1) 172 lines(density(head.LMSa), lwd = 2, col = 'blue') 173 lines(density(head.LMOp), lwd = 2, col = 'red') 174 legend('topleft', fill = c(rgb(0.9, 0, 0, 0.2), rgb(0, 0, 0.9, 0.2)), 175 legend = c('heterosexualni o.', 'jina o.'), bty = 'n') 176 177 plot(ecdf(head.LMSa), col = 'red', main = '', xlim = c(178, 220), 178 ylab = 'F(x)', xlab = '', las = 1, cex = 0.6) 179 plot(ecdf(head.LMOp), add = T, col = 'blue', cex = 0.6) 180 ks.val <- round(ks.test(head.LMOp - mean(head.LMOp), head.LMSa - (mean(head.LMSa)))$p. val, 4) 181 mtext('delka hlavy (v mm)', side = 1, line = 2.1) 182 mtext(bquote(paste('KS-test: p-hodnota = ', .(ks.val))), side = 1, line = 3.3) 183 legend('bottomright', col = c('red', 'blue'), lty = 1, pch = 20, 184 legend = c('heterosexualni o.', 'jina o.'), bty = 'n') 15 Příklad 12.5. Wilcoxonův dvouvýběrový asymptotický test Řešení příkladu 12.5 185 data <- read.delim('00-Data//05-one - sample-correlation - skull-mf.txt') 186 head(data) 187 skull.pHM <- data[data$sex == 'm', 'skull.pH ' ] 188 skull.pHM <- na.omit(skull.pHM) 189 skull.pHF <- data [data$sex == 'f, 'skull.pH ' ] 190 skull.pHF <- na.omit(skull.pHF) 191 nl <- length(skull.pHM) # 216 192 n2 <- length(skull.pHF) # 107 193 par(mar = c(5, 4, 1, 2), family = 'Times') 194 qqnorm(skull.pHM, main = '', xlab = '', pch = 21, col = 'dodgerblue3', bg = 'mintcream', 195 yla-b = 'vyberovy kvantil', las = 1) 196 mtext('teoreticky kvantil', side = 1, line = 2) 197 mtext(bquote(paste('Lillieforsuv test: p-hodnota = ', .(round(nortest::1 illie.test(skull .pHM)$p.val, 4)))) , 198 side = 1, line = 3.2) 199 qqline(skull.pHM , lwd = 2, col = 'darkblue') 200 201 qqnorm(skull.pHF, main = '', xlab = '', pch = 21, col = 'red', bg = 'bisque', 202 ylab = 'vyberovy kvantil', las = 1) 203 mtext('teoreticky kvantil', side = 1, line = 2) 204 mtext(bquote(paste('Lillieforsuv test: p-hodnota = ', .(round(nortest::1 illie.test(skull .pHF)$p.val, 4)))) , 205 side = 1, line = 3.2) 206 qqline(skull.pHF, lwd = 2, col = 'darkred') ★ 16 -3-2-10123 -2-1012 teoreticky kvantil teoreticky kvantil Lillieforsuv test: p-hodnota = 0.0329 Lillieforsuv test: p-hodnota = 0.2933 Obrázek 9: Kvantilové diagramy největší výšky mozkovny mužů (vlevo) a žen (vpravo) 118.5 127.5 136.5 145.5 120 125 130 135 140 145 150 nejvetsi vyska mozkovny (v mm) nejvetsi vyska mozkovny (v mm) KS-test: p-hodnota = 0.806 Obrázek 10: Porovnání histogramů a graf distribučních funkcí největší výšky mozkovny mužů a žen 17 207 par(mar = c(5, 4, 1, 2), family = 'Times') 208 hist(skull.pHM, prob = T, breaks = seq(123, 150, by = 3), main = '', axes = F, 209 col = rgb(0, 0, 0.9, 0.2), xlim = c(117, 152), density = 80, ylim = c(0, 0.10), 210 xlab = 'nejvetsi vyska mozkovny (v mm) ' , 211 yla-b = 'relativní četnost') 212 213 hist(skull.pHF, prob = T, breaks = seq(119, 143, by = 3), main = '', axes = F, 214 col = rgb(0.9, 0, 0, 0.2), density = 80, 215 xlab = 'vyska horni časti tvare (v mm) ' , 216 ylab = 'relativní četnost', add = T) 217 box(bty = 'o') 218 axis(l, seq(115.5, 157.5, by = 3)) 219 axis(2, las = 1) 220 lines(density(skull.pHM), lwd = 2, col = 'blue') 221 lines(density(skull.pHF), lwd = 2, col = 'red') 222 legend('topleft', fill = c(rgb(0, 0, 0.9, 0.2), rgb(0.9, 0, 0, 0.2)), 223 legend = c('muzi', 'zeny'), bty = 'n') 224 225 plot(ecdf(skull.pHM), col = 'blue', main = '', xlim = c(119, 150), ylab = 'F(x)', xlab = J J 226 cex = 0.6, las = 1) 227 plot (ecdf (skull . pHF) , add = T, col = 'red', cex = 0.6) 228 ks.val <- round(ks.test(skull.pHM - mean(skull.pHM), skull.pHF - (mean(skull.pHF)))$p. val, 4) 229 mtext('nejvet si vyska mozkovny (v mm) ' , side = 1, line = 2.1) 230 mtext(bquote(paste('KS-test: p-hodnota = ', .(ks.val))), side = 1, line = 3.3) 231 legend('bottomright', col = c('blue', 'red'), lty = 1, pch = 20, 232 legend = c('muzi', 'zeny'), bty = 'n') 18 Příklad 12.6. Wilcoxonův dvouvýběrový asymptotický test Řešení příkladu 12.6 233 data <- read.delim('00-Data//10-two - samples-means-birth.txt') 234 head(data) 235 birth.WO <- data [data$o.sib.N == '0', 'birth.W'] 236 birth.WO <- na.omit(birth.WO) 237 birth.Wl <- data [data$o.sib.N == '1', 'birth.W'] 238 birth.Wl <- na.omit(birth.Wl) 239 nl <- length(birth.WO) 240 n2 <- length(birth.Wl) 241 par(mar = c(5, 4, 1, 2), family = 'Times') 242 qqnorm(birth.W0, main = '', xlab = '', pch = 21, col = 'red', bg = 'bisque', 243 yla-b = 'vyberovy kvantil', las = 1) 244 mtext('teoreticky kvantil', side = 1, line = 2) 245 mtext(bquote(paste('Lillieforsuv test: p-hodnota = ', .(round(nortest::1 illie.test(birth .W0)$p.val , 4) ) ) ) , 246 side = 1, line = 3.2) 247 qqline(birth.W0, lwd = 2, col = 'darkred') 248 legend('topleft', fill = c(rgb(0, 0, 0.9, 0.2), rgb(0.9, 0, 0, 0.2)), 249 legend = c('muzi', 'zeny'), bty = 'n') 250 qqnorm(birth.Wl, main = '', xlab = '', pch = 21, col = 'dodgerblue3', bg = 'mintcream', 251 yla-b = 'vyberovy kvantil', las = 1) 252 mtext('teoreticky kvantil', side = 1, line = 2) 253 mtext(bquote(paste('Lillieforsuv test: p-hodnota = ', .(round(nortest::1 illie.test(birth .Wl)$p.val , 4) ) ) ) , 254 side = 1, line = 3.2) 255 qqline(birth.Wl, lwd = 2, col = 'darkblue') ★ 19 Obrázek 11: Kvantilové diagramy porodní hmotnosti novorozenců s žádným starším sourozencem (vlevo) a s jedním starším sourozencem (vpravo) 1112.5 2471.5 3830.5 1000 2000 3000 4000 5000 porodní hmotnost (v g) porodní hmotnost (v g) KS-test: p-hodnota = 0.8326 Obrázek 12: Porovnání histogramů a graf distribučních funkcí porodní hmotnosti novorozenců s žádným starším sourozencem a s jedním starším sourozencem 20 256 par(mar = c(5, 5, 1, 2), family = 'Times') 257 hist(birth.WO, prob = T, breaks = seq(886, 4963, by = 453), main = '', axes = F, 258 col = rgb(0.9, 0, 0, 0.2), density = 80, xlim = c(800, 5000), 259 xlab = 'porodni hmotnost (v g)', ylim = c(0, 0.001), 260 ylab = '') 261 hist(birth.Wl, prob = T, breaks = seq(986, 4973, by = 443), main = '', axes = F, 262 col = rgb(0, 0, 0.9, 0.2), density = 80, add = T) 263 box(bty = 'o') 264 axis(l, seq(1112.5, 4736.5, by = 453)) 265 axis(2, las = 1, at = seq(0, 0.0010, by = 0.0002), labels = c('0.0000', '0.0002', ' 0.0004', '0.0006', '0.0008', '0.0010')) 266 lines(density(birth.W0), lwd = 2, col = 'red') 267 lines(density(birth.Wl), lwd = 2, col = 'blue') 268 mtext('relativni cetnost', side = 2, line = 4) 269 legend('topleft', fill = c(rgb(0, 0, 0.9, 0.2), rgb(0.9, 0, 0, 0.2)), 270 legend = c('jeden starsi s.', 'zadny starsi s.'), bty = 'n') 271 272 plot(ecdf(birth.W0), col = 'red', main = '', xlim = c(800, 5000), 273 ylab = 'F(x)', xlab = '', cex = 0.6, las = 1) 274 plot(ecdf(birth.Wl), add = T, col = 'blue', cex = 0.6) 275 ks.val <- round(ks.test(birth.Wl - mean(birth.Wl) , birth.W0 - (mean(birth.W0)))$p.val , 4) 276 mtext('porodni hmotnost (v g)', side = 1, line = 2.1) 277 mtext(bquote(paste('KS-test: p-hodnota = ', .(ks.val))), side = 1, line = 3.3) 278 legend('bottomright', col = c('blue', 'red'), lty = 1, pch = 20, 279 legend = c('jeden starsi s.', 'zadny starsi s.'), bty = 'n') 21 12.3 Mann-Whitneyův dvouvýběrový test Nechť Xn,..., Xini, ri\ > 2 je náhodný výběr ze spojitého rozdělení s distribuční funkcí F\{x) a X21, ■ ■ ■, X2n2, U'2 > 2 je na něm nezávislý náhodný výběr ze spojitého rozdělení s distribuční funkcí F'2(x). O distribučních funkcích Fi(x) a F'2(x) předpokládáme, že se navzájem liší pouze posunutím. Nechť dále í\ je medián náhodného výběru In,..., Xíni a í2 Je medián náhodného výběru X2i, ■ ■ ■, X2„2 a x0 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,\ — ÔC2 = xq oproti Hu : í\ — X2 7^ ío (oboustranná alt.) H02 : x,\ — ÔC2 < xq oproti H12 : íi — X2 > xq (pravostranná alt.) Hqs : x,\ — ÔC2 > xq oproti H13 : í\ — X2 < xq (levostranná alt.) Test nazýváme Mann-Whitneyův dvouvýběrový test o rozdílu mediánů x\ — ž2. Testovací statistika má tvar U = SE + ^^1, (12.3) kde Se je testovací statistika definovaná vztahem 12.1, n2 je rozsah druhého náhodného výběru. Za platnosti nulové hypotézy pochází statistika U a ze standardizovaného normálního rozdělení, tj. U = SE + Mn2 + 1) Homiy Hn : ži - Í2 ^ x0 W= (-00; wnun2(a/2) - 22Ím±ll^ y ^ni>„2(l - a/2) - 22ÍI^±Íl; ^ H12 : ži - x2 > x0 W = (wnun2(l - a) - "'("2'+1) ; H13 : xi - x2 < x0 W = (-00; wnun2(a) - "2("22+1)\> kde wni^n2(a/2), wnitTl2(l — a/2), wni^n2(a), wnitTl2(l — a) jsou tabelované kvantily pro dvouvýběrový Wilcoxonův test, jejichž hodnoty získáme příkazem qwilcox(). Interval spolehlivosti má podle zvolené alternativní hypotézy jeden z následujících tvarů ifn : íi - Í2 ^ io (d, h) = (u^-^)-23^) ; u^^d-^+i-^^l^ H12 :x1-x2>x0 (d, 00) = (V(^i."2(")-"2("22+1)) ; 00) H13 :x1-x2 «)} H12 : x,\ — ÔC2 > xq p-hodnota = Pr(č7 > u) = 1 — Pr(č7 < u) H13 : x,\ — ÔC2 < xq p-hodnota = Pr(č7 < u) kde U je náhodná veličina, u je realizace testovací statistiky U (viz vzorec 12.3), tedy konkrétní číslo, a Pr(č7 < u) je distribuční funkce tabelovaného rozdělení pro Wilcoxonův dvouvýběrový exaktní test, jejíž hodnotu získáme pomocí ď a implementované funkce pwilcoxQ. 22 Příklad 12.7. Mann-Whitneyův dvouvýběrový test Řešení příkladu 12.7 ★ 280 data <- read.delim('00-Data//18-more-samples-variances-clavicle.txt') 281 head(data) 282 cla.Lil <- data[data$pop == 'indl' , 'cla.L'] 283 cla.Lil <- na.omit(cla.LI1) 284 cla.LI2 <- data[data$pop == 'ind2' , 'cla.L'] 285 cla.LI2 <- na.omit(cla.LI2) 286 287 wilcox.test(cla.Lil, cla.LI2, alt = 't')$p.val # Z 23 Příklad 12.8. Mann-Whitneyův dvouvýběrový test Řešení příkladu 12.8 ★ 288 data <- read.delim('00-Data//15-anova-means - skull.txt ' ) 289 head(data) 290 upface.HN <- data[data$pop 291 upface.HN <- na.omit(upface 292 upface.HC <- data [data$pop 293 upface.HC <- na.omit(upface 294 295 wilcox.test(upface.HN, upface.HC, alt = 'g')$p.val # N == 'nem', 'upface.H'] . HN) == 'cin', 'upface.H'] . HC) 21 Příklad 12.9. Mann-Whitneyův dvouvýběrový test Řešení příkladu 12.9 ★ 296 data <- read.delim('00-Data//19-more-samples-correlations - skull.txt') 297 head(data) 298 intorb.BP <- data[data$pop == 'per', 'intorb.B'] 299 intorb.BP <- na.omit(intorb.BP) 300 intorb.BB <- data [data$pop == 'ban', 'intorb.B'] 301 intorb.BB <- na.omit(intorb.BB) 302 303 wilcox.test(intorb.BP, intorb.BB, alt = '1')$p.val # Z 25