Testy středních hodnot

Často srovnáváme skupiny, které vykazují jistou variabilitu. Rádi bychom ovšem zjistili, zda dochází „v průměru“ k nějakému rozdílu mezi skupinami. K tomu slouží testy na porovnávání středních hodnot.

Dobře udělaný exploratorní graf nám můžeme na otázku v podstatě rovnou odpovědět.

library(ggplot2)
ggplot(ToothGrowth,
    aes(y = len, x = supp, color = as.factor(dose), group = dose, shape = supp)) +
    geom_point(size=3) +
    stat_summary(fun.y = mean, geom = "line") +
    stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = 0.1)

Tolik zajímavých informací v jednom zobrazení! Vliv potravinových doplňků na růst zubů u morčat. Supp: typ doplňku (OJ = pomerančový džus, VC = vitamin C), dose: dávka (0.5mg, 1mg a 2mg), len: velikost zubů. Přestože je veliká variabilita, tak čáry spojující středy stejných dávek obou doplňků ukazují, jak by na data nahlížel t-test. Dopředu je vidět, jak by asi dopadl. Znázorněny jsou také 95% intervaly spolehlivosti odhadu středních hodnot. V malých dávkách je džus lepší, v největší dávce jsou oba doplňky srovnatelné.

Nicméně teprve statistický test dokáže odpověď opravdu kvantifikovat, bere v úvahu míru nejistoty.

Diagram pro výběr správného testu

Jak zjistit závislost výběrů? Např. disfluence v řeči střízliví vs. po požití alkoholu.

  1. Ti samí za prvé před a za druhé po => závislé skupiny (páry). Zdá se lepší pro srovnání, ale těžší zajistit, aby v 2. testu již nevěděli, co je čeká (systematická chyba).

  2. Jedna skupina střízlivá, jiná skupina alkohol => nezávislé výběry (jiní lidé).

Párové testy jsou silnější. Zěměřují se totiž přímo na rozdíly hodnot vždy v daném páru, takže jim nevadí vysoká míra variability mezi jednotlivými subjekty. Oproti tomu testy nezávislých skupin nazírají na skupinu jako na jeden celek a zaměřují se pouze na její celkový střed. Odlehlé body jedné skupiny tak mohou výsledky nepříjemně ovlivnit.

ggplot(sleep, aes(x = group, y = extra)) +
    geom_point(size=5, pch=21, fill="salmon", alpha=.5)

ggplot(sleep, aes(x = group, y = extra, group = ID, fill=ID)) +
    geom_point(size=5, pch=21, alpha=.5) +
    geom_line(size=1, aes(color=ID))

Nalevo nahlížíme na data jako na nezávislé skupiny, nezbyde nám tedy než porovnat jen celkové průměry skupin, které vzhledem k variabilitě mohou mít široké konfidenční intervaly.

Napravo říkáme, že skupiny jsou spárované, a je jasně vidět stoupající trend většiny bodů i přes velkou variabilitu subjektů.

T-testy

T-testy slouží pro metrická data s normálním rozdělením, porovnávají průměry. Zjišťujeme, zda rozdíl průměru od dané hodnoty je jen dílem náhody. Fungují na principu intervalových odhadů. Řadí se do kategorie parametrických testů, protože pracují s parametry normálního rozdělení.

  • Jednovýběrový. Porovnáváme průměr jedné skupiny s konkrétním číslem \(\mu_0\). H0: \(\mu = \mu_0\). Oboustranný H1: \(\mu \neq \mu_0\). Jednostranný buď H1: \(\mu > \mu_0\) nebo H1: \(\mu < \mu_0\).

    Jednostranné testy se používají méně často a jen tam, kde si jsme velmi jisti, že jen jeden směr dává smysl.

  • Dvouvýběrový párový. Zjišťujeme, zda průměr rozdílů spárovaných výběrů je jen dílem náhody. H0: \(\Delta \mu _{xy} = 0\). Oboustranný H1: \(\Delta \mu _{xy} \neq 0\). Jednostranný buď H1: \(\Delta \mu _{xy} > 0\) nebo H1: \(\Delta \mu _{xy} < 0\).

  • Dvouvýběrový test nezávislých výběrů. Porovnáváme průměry dvou nezávislých skupin. H0: \(\mu_x = \mu_y\). Oboustranný H1: \(\mu_x \neq \mu_y\). Jednostranný buď H1: \(\mu_x > \mu_y\) nebo H1: \(\mu_x < \mu_y\).

    Nejdříve je nutno otestovat pomocí F-testu, zda mají výběry stejný rozptyl. Pak teprve zvolíme vhodnou variantu t-testu (dvouvýběrový t-test s rovností rozptylů / bez rovnosti rozptylů). F-test pro testování rozptylů (funkce var.test), H0: rozptyly jsou shodné, H1: rozptyly se liší.

    T-test s rovností rozptylů je silnější, ale můžeme ho použít jen v případě splnění podmínek.

Alternativní hypotézu H1 specifikujeme parametrem alternative = “two.sided” (není rovno), “greater” (větší) nebo “less” (menší).

Příklad 1

Teorie: podíl železa ve sloučenině je 12.1 %. Máme VS s \(n = 9\) vzorky. H0: \(\mu = 12.1\), H1: \(\mu \neq 12.1\). \(\alpha = 0.05\).

Index 1 2 3 4 5 6 7 8 9
Podíl (%) 11.7 12.2 10.9 11.4 11.3 12.0 11.1 10.7 11.6

Jeden výběr, předpokládáme normální data => jednovýběrový t-test. Vzhledem k H1 oboustranný test.

x <- c(11.7, 12.2, 10.9, 11.4, 11.3, 12.0, 11.1, 10.7, 11.6)
t.test(x = x, alternative = "two.sided", mu = 12.1, conf.level = 0.95)
## 
##  One Sample t-test
## 
## data:  x
## t = -4.0406, df = 8, p-value = 0.003733
## alternative hypothesis: true mean is not equal to 12.1
## 95 percent confidence interval:
##  11.05286 11.81380
## sample estimates:
## mean of x 
##  11.43333

\(p = 0.003733 < \alpha,\) zamítáme H0 a přijímáme H1.

Poznámka: t-test souvisí s klasickým intervalovým odhadem \(\mu\) z VS. Vidíme, že \(\mu_0\) neleží v tomto intervalu, proto zamítáme H0.

Příklad 2

Dvě metody analýzy podílu fosforu, každý vzorek analyzován oběma metodami. Je mezi metodami významný rozdíl? Jinými slovy, není jedna trochu vychýlena oproti druhé? H0: \(\Delta \mu = 0,\) H1: \(\Delta \mu \neq 0.\) \(\alpha = 0.05\).

Index 1 2 3 4 5 6 7 8 9 10
Metoda 1 13.3 17.6 4.1 17.2 10.1 3.7 5.1 7.9 8.7 11.6
Metoda 2 13.4 17.3 4.1 17.0 10.3 4.0 5.1 8.0 8.8 12.0

Párový test, oboustranný.

x <- c(13.3, 17.6, 4.1, 17.2, 10.1, 3.7, 5.1, 7.9, 8.7, 11.6)
y <- c(13.4, 17.3, 4.1, 17.0, 10.3, 4.0, 5.1, 8.0, 8.8, 12.0)
t.test(x = x, y = y, alternative = "two.sided", paired = TRUE, conf.level = 0.95)
## 
##  Paired t-test
## 
## data:  x and y
## t = -1.0487, df = 9, p-value = 0.3217
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.22099887  0.08099887
## sample estimates:
## mean of the differences 
##                   -0.07

\(p = 0.3217 \geqslant \alpha,\) nezamítáme H0.

Příklad 3

Běh 50 m, 10 chlapců (M) a 8 děvčat (F). Dosahují chlapci kratších časů? H0: \(\mu _M = \mu _F,\) H1: \(\mu _M < \mu _F.\) \(\alpha = 0.05\).

Index 1 2 3 4 5 6 7 8 9 10
M 9.3 9.9 9.0 8.9 9.6 10.6 9.5 10.0 9.4 10.1
F 10.7 10.0 9.2 9.9 9.3 9.8 10.0 12.0

Dvouvýběrový test, nezávislá měření (o párech ani nemusíme uvažovat vzhledem k tomu, že v každé skupině je jiný počet subjektů), jednostranný. Nutno ověřit shodnost rozptylů.

m <- c(9.3, 9.9, 9.0, 8.9, 9.6, 10.6, 9.5, 10.0, 9.4, 10.1)
f <- c(10.7, 10.0, 9.2, 9.9, 9.3, 9.8, 10.0, 12.0)
mean(m)
## [1] 9.63
mean(f)
## [1] 10.1125
var.test(x = m, y = f, conf.level = 0.95)
## 
##  F test to compare two variances
## 
## data:  m and f
## F = 0.34652, num df = 9, denom df = 7, p-value = 0.1414
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.07184355 1.45434810
## sample estimates:
## ratio of variances 
##           0.346517

\(p = 0.1414 \geqslant \alpha,\) nezamítáme H0 (shodnost rozptylů). Volíme jednostranný dvouvýběrový t-test s rovností rozptylů.

t.test(x = m, y = f, alternative = "less", paired = FALSE, var.equal = TRUE, conf.level = 0.95)
## 
##  Two Sample t-test
## 
## data:  m and f
## t = -1.4341, df = 16, p-value = 0.0854
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
##       -Inf 0.1049052
## sample estimates:
## mean of x mean of y 
##    9.6300   10.1125

\(p = 0.0854 \geqslant \alpha,\) nezamítáme H0. Přestože průměr chlapců je nižší než průměr dívek, rozdíl je vzhledem k malému počtu dat a jejich rozptylu příliš nepatrný.

Poznámka: kdybychom pro stejná data zvolili „univerzálnější“ t-test s nerovností rozptylů (jsme opatrnější, var.equal = FALSE), obdrželi bychom vyšší p-hodnotu 0.1017. Takový test má nižší sílu (\(1 - \beta\)), a tedy menší šanci objevit H1.

Neparametrické testy středních hodnot (distribution free)

Lze je použít i pro silně nenormální rozdělení, i pro řadové proměnné. Celkově však mají menší sílu (\(1 - \beta\)) než parametrické testy (např. t–test). Znamená to, že rády zůstávají opatrně u nezamítnutí H0 (buď tedy H0 skutečně platí, nebo je málo dat, aby byl prokázán opak). Pro přijetí H1 je něco musí velmi silně přesvědčit. Často se nám tak stane, že přestože H1 platí, test to není schopen prokázat.

Pracují obecně s nesymetrickými rozděleními. Místo průměru se proto počítá s mediánem, protože leží skutečně uprostřed (značíme \(x_{0.5}\) – je to 50% kvantil). U symetrických rozdělení je medián roven střední hodnotě – průměru.

Znaménkový test

Nejjednodušší, ale také velmi slabý test, proto je malá šance prokázat H1. Vnímá totiž jen, zda naměřená hodnota je větší (+) či menší (-), nebere v potaz velikost změny. Položky, kde nedošlo ke změně, nepočítáme.

Oboustranný znaménkový test je exaktní (tedy přesnější) obdobou McNemarova testu popsaného v kapitole 15. Testy četností.

H0: počet „+“ a „-“ je stejný.

  • Oboustranný test: nejčastější varianta. H1: počty „+“ a „-“ nejsou stejné.

    Z = menší z počtu obou znamének.

    p <- min(1, 2*pbinom(Z, n, 0.5))

    Příklad. Celkem 16 dětí se zúčastnilo intenzivní výuky výslovnosti angličtiny. Ve srovnání se stavem „před“ dopadlo v následném testování 12 dětí lépe a 4 se naopak zhoršily. H0: metoda nemá vliv, H1: metoda má vliv, \(\alpha = 0.05\).

znamenkaPlus <- 12;
znamenkaMinus <- 4;
n <- znamenkaPlus + znamenkaMinus
Z <- min(znamenkaPlus, znamenkaMinus)  # menší z obou počtů

p <- min(1, 2*pbinom(Z, n, 0.5))
p
## [1] 0.07681274

\(0.0768 \geqslant \alpha,\) nezamítáme H0. Rozhodnutí testu je možná překvapivé, ale je skutečně slabý, takže tento poměr pro něj stále ještě není dostatečně přesvědčivý, aby přijal H1.

  • Jednostranný test. H1: počty „-“ jsou častější (Z = počet znamének „+“), případně opačná H1: počty „+“ jsou častější (Z = počet znamének „-“).

    p <- pbinom(Z, n, 0.5)

    Příklad. Viz předchozí zadání, H1: metoda má kladný vliv (počty „+“ jsou častější).

Z <- znamenkaMinus
p <- pbinom(Z, n, 0.5)
p
## [1] 0.03840637

\(0.0384 < \alpha,\) přijímáme H1.

Jednostranné testy mají obecně větší sílu objevit H1 (což vidíme zde, p-hodnota je skutečně nižší ve srovnání s oboustranným testem).

Wilcoxonův test jednovýběrový či dvouvýběrový párový

Používá se v podobných situacích jako znaménkový test, ale zohledňuje i velikost rozdílu. Proto je test silnější, neboli máme větší šanci odhalit malé rozdíly mezi měřeními a rozhodnout se pro H1. Rozdíly srovná podle velikosti – uvažuje tedy „menší“ a „větší“ rozdíl, nikoliv skutečnou velikost rozdílu (jak počítá t-test). Pokud je možné na data použít t-test, je vhodnější jej preferovat, protože ten má ještě větší sílu než Wilcoxonův test.

Obdobně jako u t-testů, R má jednu společnou funkci wilcox.test() pro jednovýběrový test (zadáme parametry x a referenční hodnotu mediánu mu) i pro porovnání dvou výběrů (zadáme x i y a parametr paired = TRUE).

Příklad Mějme situaci z předchozího příkladu, ale k dispozici máme u každého žáka body z páru (před a po). H1: metoda má vliv, \(\alpha = 0.05\).

pred <- c(11,  8, 8, 14, 13, 7, 12, 15, 11, 16, 9,  5, 18, 14, 15, 9)
po <-   c(15, 11, 9, 16, 11, 4, 13, 19, 15, 15, 13, 12, 20, 13, 18, 15)
wilcox.test(x = pred, y = po, alternative = "two.sided", paired = TRUE,
            exact = FALSE, correct = TRUE, conf.level = 0.95)
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  pred and po
## V = 20, p-value = 0.0136
## alternative hypothesis: true location shift is not equal to 0

\(0.0136 < \alpha,\) přijímáme H1.

A co jednostranná hypotéza? H1: metoda má pozitivní vliv (neboli medián bodů před je nižší než po).

wilcox.test(x = pred, y = po, alternative = "less", paired = TRUE,
            exact = FALSE, correct = TRUE, conf.level = 0.95)
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  pred and po
## V = 20, p-value = 0.006798
## alternative hypothesis: true location shift is less than 0

\(0.0068 < \alpha,\) přijímáme H1.

Je vidět, že se znalostí velikosti rozdílů je Wilcoxonův test silnější než znaménkový test.

Příklad Vzpomeňme si na příklad odhadu intervalu spolehlivosti pro medián z kapitoly 12. Intervalové odhady.

Řešme hypotézu, zda medián ZS je roven 11 na hladině významnosti \(\alpha = 0.05\), k dispozici máme následující vzorek. H0: medián = 11, H1: medián není 11.

cisla <- c(11,  8, 8, 14, 13, 7, 12, 15, 11, 16, 9,  5, 18, 14, 15, 9)

wilcox.test(x = cisla, mu = 11, alternative = "two.sided", paired = FALSE,
            exact = FALSE, correct = FALSE, conf.level = 0.95)
## 
##  Wilcoxon signed rank test
## 
## data:  cisla
## V = 63, p-value = 0.5079
## alternative hypothesis: true location is not equal to 11

Nemůžeme zamítnout H0. Jak tento test souvisí s intervalem spolehlivosti?

library(asbio)     # pokud není nainstalována, tak jednorázově nainstalovat install.packages("asbio")
## Warning: package 'asbio' was built under R version 3.2.5
## Loading required package: tcltk
ci.median(cisla, conf = 0.95)
## 
## 95% Confidence interval for population median 
## Estimate     2.5%    97.5% 
##     11.5      8.0     15.0

Číslo 11 v intervalu leží, a proto není možné hypotézu H0 zamítnout.

Dvouvýběrový Wilcoxonův test pro nezávislé skupiny (Mann-Whitney U-test)

Tento test má spoustu různých názvů, např. také Mann–Whitney–Wilcoxon, případně Wilcoxon–Mann–Whitney test nebo Wilcoxon rank-sum test. Důležité však je, že jde o alternativu předchozího testu, kdy ale dvě skupiny nejsou spárovány, jedná se o nezávislá měření. V R se test volá stejnou funkcí wilcox.test(), jen s parametrem paired = FALSE.

Test opět zjišťuje, zda dva výběry mají shodný medián. Nicméně způsob výpočtu není zcela jednoduchý, a tak se v literatuře objevují varování, že by se neměl používat pro porovnávání dvou výběrů s výrazně odlišným rozdělením, kdy výsledky testu nemusejí být zcela správné. Každopádně tento test je velice robustní na nejrůznější nenormální rozdělení, takže je to pro nenormální data stále výrazně lepší varianta než použít nesprávně t-test.

Příklad Ve dvou skupinách žáků byly v testu získány následující počty bodů. H0: mezi výsledky obou skupin nejsou významné rozdíly, H1: významné rozdíly jsou.

skupinaA <- c(5, 7, 12, 15, 13)
skupinaB <- c(9, 6, 16, 11, 14, 17)
wilcox.test(x = skupinaA, y = skupinaB, alternative = "two.sided", paired = FALSE,
            exact = TRUE, conf.level = 0.95)
## 
##  Wilcoxon rank sum test
## 
## data:  skupinaA and skupinaB
## W = 11, p-value = 0.5368
## alternative hypothesis: true location shift is not equal to 0

\(0.5368 \geqslant \alpha,\) nezamítáme H0.

Příklad To samé, ale jiná data.

skupinaA <- c(8, 7, 2, 1, 3, 3, 2, 7, 4, 2, 5, 5, 5, 2, 2)
skupinaB <- c(7, 7, 7, 3, 3, 1, 10, 10, 10, 7, 8, 8, 9, 9, 9, 5, 5, 5, 8, 19, 17, 6, 6)
wilcox.test(x = skupinaA, y = skupinaB, alternative = "two.sided", paired = FALSE,
            exact = FALSE, correct = TRUE, conf.level = 0.95)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  skupinaA and skupinaB
## W = 59.5, p-value = 0.0007196
## alternative hypothesis: true location shift is not equal to 0

\(0.0007 < \alpha,\) přijímáme H1.

Úlohy

  1. Viz domácí úkol.

© 29. 4. 2015 Tomáš Bořil, borilt@gmail.com