Intervalové odhady

Průměr je sám o sobě náhodnou veličinou. Počítáme ho totiž z konkrétního jednoho výběrového souboru (VS), což je náhodný výběr ze základního souboru (ZS), jehož střední hodnotu chceme pomocí průměru VS odhadnout. Mluvíme zde o tzv. bodovém odhadu, protože výsledkem je jedno číslo.

Bylo by dobré vyjádřit, jak moc můžeme konkrétnímu průměru jakožto odhadu střední hodnoty ZS věřit. Statistika nám pro to naštěstí nabízí nástroj, a dokonce si sami můžeme zvolit, jak přesný odhad chceme vytvořit. Výstupem není jedno číslo (průměr), ale rovnou celý interval, ve kterém se s určitou námi zvolenou pravděpodobností skutečná střední hodnota ZS nachází.

Je však zřejmé, že čím větší jistotu chceme mít, že střední hodnota bude v našem intervalovém odhadu ležet, tím širší interval musíme vytvořit. Stoprocentní jistotu budeme mít pouze v případě, kdy bude interval nekonečně široký, což samozřejmě není ideální. Optimálním řešením se jeví zvolit rozumný kompromis mezi šířkou intervalu tak, aby byly získané hodnoty k něčemu použitelné, a zároveň byla dostatečně vysoká pravděpodobnost, že skutečná střední hodnota v tomto intervalu bude ležet. Mluvíme tak například o 95% konfidenčním intervalu (tzv. confidence level), což znamená, že když ze ZS uděláme mnoho jednotlivých náhodných VS, pro každý vypočteme interval, tak 95 % intervalů bude zahrnovat správnou střední hodnotu. Ve statistice se často používá termín hladina významnosti (significance level) \(\alpha\) udávající, kolik procent případů bude naopak mimo, v našem případě bychom se tedy bavili o hladině významnosti \(\alpha = 0.05\) (5 % intervalů nebude obsahovat správnou střední hodnotu).

Interval spolehlivosti střední hodnoty normálního rozdělení

Abychom mohli dobře sledovat, jak intervalový odhad funguje, vymysleme si ZS, u kterého budeme znát skutečné parametry.

Mějme např. ZS s normálním rozdělením, střední hodnotou \(\mu = 100\) a směrodatnou odchylkou \(\sigma = 15.\)

Vybereme náhodný výběr, ze kterého se budeme snažit odhadnout skutečnou střední hodnotu.

mi <- 100
sigma <- 15

VS1 <- rnorm(100, mean = mi, sd = sigma)   # VS o velikosti 100 prvků
hist(VS1, breaks = 20)

Tradiční „učebnicový“ výpočet krok po kroku

alfa <- 0.05    # hladina významnosti 0.05 odpovídá 95% konfidenčnímu intervalu

n <- length(VS1)    # automatické zjištění počtu prvků VS
prumer <- mean(VS1) # průměr výběru
s <- sd(VS1)        # výběrová směrodatná odchylka

SE <- s/sqrt(n)     # tzv. standardní chyba (standard error)
kvantil <- qt(1-alfa/2, df = n-1)  # kvantil t-rozdělení s n-1 stupni volnosti

c(prumer - kvantil*SE, prumer + kvantil*SE)  # interval spolehlivosti střední hodnoty
## [1]  98.95998 104.30664

Výhodný výpočet využívající vedlejší efekt t-testu v R

t.test(VS1, conf.level = 1-alfa)$conf.int
## [1]  98.95998 104.30664
## attr(,"conf.level")
## [1] 0.95

Poznámka pro zvídavé

když zavoláme samotný t-test, dostaneme spoustu zajímavých údajů

t.test(VS1, conf.level = 1-alfa)
## 
##  One Sample t-test
## 
## data:  VS1
## t = 75.435, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##   98.95998 104.30664
## sample estimates:
## mean of x 
##  101.6333

Je to však jen textová forma různých hodnot, ke kterým se dá dostat, když známe jejich jména. Ta zjistíme funkcí names().

names(t.test(VS1, conf.level = 1-alfa))
##  [1] "statistic"   "parameter"   "p.value"     "conf.int"    "estimate"   
##  [6] "null.value"  "stderr"      "alternative" "method"      "data.name"

A právě takto jsme zjistili, že za t.test() stačí napsat $conf.int a dostaneme se tak k našemu intervalu, což je zvláště výhodné, pokud bychom chtěli s hodnotami dále někde automaticky počítat.

Co se stane, když bude mít VS menší počet prvků?

Čím více hodnot, tím přesnější odhad, tím pádem užší interval. Šířka intervalu obecně závisí na počtu prvků i na rozptylu samotných hodnot. Pokud je hodnot málo nebo jsou velmi rozptýleny, statistika je chytrá a pozná, že jim nemůže věřit tak dobře, jako kdyby jich bylo hodně nebo byly málo rozptýleny.

VS2 <- rnorm(5, mean = mi, sd = sigma)   # VS o velikosti 5 prvků
t.test(VS2, conf.level = 1-alfa)$conf.int
## [1]  85.30079 102.78484
## attr(,"conf.level")
## [1] 0.95
t.test(VS2, conf.level = 1 - 0.01)$conf.int  # chceme mít větší jistotu, změníme hladinu významnosti
## [1]  79.54619 108.53944
## attr(,"conf.level")
## [1] 0.99

Je vidět, že pro VS s menším počtem prvků vychází širší interval, protože si zkrátka není tak jistý, a proto je opatrnější. Pokud chceme zvýšit pravděpodobnost úspěchu a mít 99% konfidenční interval, hranice budou ještě více do šířky.

Experiment: úspěšnost odhadu při velkém počtu VS

Nr <- 50  # Počet realizací VS

mi <- 100    # Parametry ZS
sigma <- 15

alfa <- 0.05     # námi zvolená hladina významnosti -> ovlivní šířku intervalů, ale i spolehlivost
pocetSpatne <-0  # počítadlo neúspěchů

plot(c(1, Nr), c(mi, mi), type = 'l', lty = 'dashed', xlab = 'realizace', ylab = '')  # vodorovná čára skutečné střední hodnoty

for (i in 1:Nr) {
    VS <- rnorm(5, mean = mi, sd = sigma)  # jeden VS má 5 hodnot
    interval <- t.test(VS, conf.level = 1-alfa)$conf.int
    
    if (mi >= interval[1]  &  mi <= interval[2]) {  # mi leží uvnitř intervalu
        segments(x0 = i, y0 = interval[1], x1 = i, y1 = interval[2], col = 'blue')
    } else {  # mi leží mimo meze intervalu
        pocetSpatne <- pocetSpatne + 1
        segments(x0 = i, y0 = interval[1], x1 = i, y1 = interval[2], col = 'orange')
    }
}

# zastoupení intervalů, které neobsahují správnou střední hodnotu, pro nekonečné Nr bude přesně rovno alfa
pocetSpatne/Nr
## [1] 0.04

Interval spolehlivosti rozptylu a směrodatné odchylky normálního rozdělení

Výše jsme se bavili jen o tom, jak určit odhad střední hodnoty ZS. Často ale chceme vědět, jak jsou data rozptýlena. Rozptyl je však opět jen bodový odhad zatížený náhodnou chybou, a proto je dobré umět určit interval spolehlivosti pro odhad rozptylu ZS.

Zatím se mi nepodařilo najít jednoduchou funkci v R, která by toto přímo dělala, a proto definujme svou vlastní. Její použití je pak jednoduché.

var.interval <- function(data, conf.level = 0.95) {
    a <- 1 - conf.level
    df <- length(data) - 1
    chilower <- qchisq((1 - conf.level)/2, df)
    chiupper <- qchisq((1 - conf.level)/2, df, lower.tail = FALSE)
    v <- var(data)
    c(df * v/chiupper, df * v/chilower)
}
mi <- 100
sigma <- 15
VS1 <- rnorm(100, mean = mi, sd = sigma)   # VS o velikosti 100 prvků

var.interval(VS1)       # interval spolehlivosti rozptylu
## [1] 139.9342 244.9613
sqrt(var.interval(VS1)) # interval spolehlivosti směrodatné odchylky
## [1] 11.82938 15.65124

Experiment: úspěšnost při velkém počtu opakování

Nr <- 10000  # Počet realizací VS

mi <- 100    # Parametry ZS
sigma <- 15

alfa <- 0.05     # námi zvolená hladina významnosti -> ovlivní šířku intervalů, ale i spolehlivost
pocetSpatne <- 0  # počítadlo neúspěchů

for (i in 1:Nr) {
    VS <- rnorm(5, mean = mi, sd = sigma)  # jeden VS má 5 hodnot
    interval <- sqrt(var.interval(VS, conf.level = 1-alfa))
    
    if (sigma >= interval[1]  &  sigma <= interval[2]) {  # sigma leží uvnitř intervalu
        
    } else {  # sigma leží mimo meze intervalu
        pocetSpatne <- pocetSpatne + 1
    }
}

# zastoupení intervalů, které neobsahují správnou směrodatnou odchylku, pro nekonečné Nr bude přesně rovno alfa
pocetSpatne/Nr
## [1] 0.0484

Interval spolehlivosti pro binomická rozdělení

Př. V náhodně vybraném vzorku ze 100 lidí hlasovalo 58 pro zavření velkých obchodů během svátků. Dá se očekávat, že to je skutečně většinový názor celé populace?

binom.test(58, 100)$conf.int
## [1] 0.4771192 0.6780145
## attr(,"conf.level")
## [1] 0.95

Zdá se, že to tak jisté není. Klidně to také může být názor menší části společnosti. Z takto malého počtu lidí jednoznačný závěr neuděláme. Pokud bychom se ale zeptali 1000 lidí a vyšel stejný poměr (tedy 580 pro), je to už pro statistiku dostatečně velký vzorek pro signifikantní potvrzení nadpoloviční většiny:

binom.test(580, 1000)$conf.int
## [1] 0.5487112 0.6108164
## attr(,"conf.level")
## [1] 0.95

Př. V percepčním testu bylo jedno slovo rozpoznáno správně 12 lidmi z celkového počtu 30. Určete interval spolehlivosti pro úspěšnost rozpoznání tohoto slova v rámci celé populace, ze které byli tito lidé náhodně vybráni.

binom.test(12, 30)$conf.int
## [1] 0.2265576 0.5939651
## attr(,"conf.level")
## [1] 0.95

Skutečný poměr v populaci je tedy něco mezi 0.23 a 0.59, což je v tento moment dost nepřesně určená hodnota, opět by to chtělo větší vzorek pro zpřesnění (zúžení) intervalu.

Př. Řešíme hypotézu, zda lidé skutečně slyší znělost hlásek. V percepčním testu byly páry s/z a š/ž. Každá hláska z páru byla namluvena 5 různými mluvčími. Testu se zúčastnilo 30 hodnotitelů. U každé ze čtyř hlásek chceme vyhodnotit, jaký je poměr úspěšnosti jejího rozpoznání a doplnit o konfidenční interval pro hladinu významnosti 0.05.

Příklad pro hlásku /z/, která byla tedy celkově hodnocena 150krát (5 mluvčích krát 30 hodnotících). Když si za každé úspěšné rozpoznání započítáme jeden bod, můžeme tedy v součtu dostat hodnotu mezi 0 a 150. V našem případě to dopadlo 131 úspěšnými rozpoznáními (tedy ve 37 případech byla označena mylně za neznělou /s/).

Jednoduše by člověk vypočetl poměr 131 / 150 = 0.8733 a byl by hotov s tím, že položka byla rozpoznána cca v 87 % případů. To je ale jen poměr pro tento konkrétní náhodný výběr, správné by bylo přidat statistiku určující odhad poměru pro základní soubor, tedy konfidenční interval.

Nabízí se využít konfidenční interval binomického rozdělení, ale jsou tu dvě záležitosti, které je nutné poladit.

První jsou replikace, kde jeden hodnotitel slyšel hlásku od 5 mluvčích. Ve výsledném součtu může být ten samý hodnotitel započítán vícekrát a kdybychom vypočetli interval pro prostý poměr 131 ze 150, interval by si myslel, že máme 150 nezávislých hodnocení, čímž bychom se dopustili typické chyby (tzv. pseudoreplikace). Musíme tedy za každého náhodně a nezávisle vybraného hodnotícího vypočítat jednu průměrnou hodnotu mezi 0 a 1. Vyřešíme to jednoduše tak, že obě čísla poměru vydělíme počtem mluvčích. Místo 131 ku 150 dodáme do funkce 131/5 ku 150/5, tedy celkový počet bude roven správně 30 hodnotitelům. Přestože poměr vychází stejně, pro konfidenční interval je celkový počet 30 důležité číslo, které má zásadní dopad na šířku intervalu. Pro číslo 150 by interval vycházel přehnaně optimisticky úzký, nereflektoval by realitu pouhých 30 hodnotících posluchačů.

Druhým problémem je, že máme jednu hlavní hypotézu (lidé poznají poslechem znělost hlásek), zde máme ale čtyři hlásky, pro každou děláme samostatný interval. Pokud tento interval nebude zahrnovat hodnotu 0.5 (náhodné tipování padesát na padesát), prohlásíme to za důkaz, že se nejedná o náhodný trend. Jedná se o typický příklad násobného testování a je potřeba provést korekci hladiny významnosti tak, aby nebyla 0.05 pro každý dílčí test (každou hlásku), ale aby platila dohromady pro celou úlohu všech testovaných hlásek. Jinak by mohly přirozeně vyskakovat příliš často falešné objevy jen díky chybě 1. druhu z důvodu většího množství dílčích testů v rámci celkové hypotézy (my ale chceme, aby chyba 1. druhu mohla nastat s poměrem 0.05 pro celou hlavní hypotézu). Více viz multiple testing na xkcd. Nejbezpečnější je provedení tzv. Bonferroniho korekce, kdy 0.05 vydělíme počtem dílčích testů (zde 4 hlásky).

Celý kód vypadá následovně. Funkce binconf() z balíčku Hmisc umožňuje zadávat i necelá čísla (je to interpolace původní úlohy binomického rozdělení, kde byla vstupem jen přirozená čísla) a aby byly výsledky srovnatelné s hodnotami z binom.test(), je nutné vybrat přesnou metodu “exact”. Vidíme, že poměr vychází skutečně 0.8733, navíc ale vidíme, že pro základní soubor se pohybuje mezi cca 0.65 a 0.98 s korigovanou hladinou významnosti. Jelikož tento interval nezahrnuje hladinu čistě náhodného tipování 0.5, je to pro nás důkaz, že znělost /z/ lze správně rozpoznat.

library(Hmisc)
binconf(131/5, 150/5, alpha = 0.05/4, method = "exact")
##   PointEst     Lower     Upper
##  0.8733333 0.6522646 0.9778362

Interval spolehlivosti pro Poissonovo rozdělení

Př. Datový FTP server spadl 7x během 100 dní. Určete 95% konfidenční interval průměrného počtu spadnutí za 1 den.

Pouhý průměr je 7/100 = 0.07, ale konfidenční interval dává lepší odhad skutečné střední hodnoty.

poisson.test(7, T = 100)$conf.int
## [1] 0.02814363 0.14422675
## attr(,"conf.level")
## [1] 0.95

Interval spolehlivosti pro medián

Pokud naše data nemají normální rozdělení a rozhodneme se použít medián jako ukazatel střední hodnoty, opět bude užitečné určit interval spolehlivosti pro tento odhad.

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

library(asbio)     # pokud není nainstalována, tak jednorázově nainstalovat install.packages("asbio")
## Warning: package 'asbio' was built under R version 4.0.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

Vidíme, že medián na 95 % leží v poměrně širokém intervalu, který pokrývá téměř všechna naše data. Pro přesnější odhad (užší interval) by bylo nutné použít vzorek s větším počtem hodnot.

Úlohy

  1. Zkopírujte si program z podkapitoly Experiment: úspěšnost odhadu při velkém počtu VS (v sekci o střední hodnotě) do vlastního skriptu a experimentujte s nastavením parametrů.
    • Měňte hladinu významnosti \(\alpha\), např. 0.01 a 0.001. Budou se intervaly spolehlivosti rozšiřovat nebo zužovat?
    • Vraťte hladinu významnosti na \(\alpha = 0.05\) a měňte velikost VS. Lze spočítat interval spolehlivosti z VS obsahujícího jen jednu hodnotu?
    • Zkoušejte velikosti 2, 5, 10, 15, 20, 25. Rozšiřují se intervaly nebo zužují?
    • Platí přímá úměra ve smyslu „zvýšení velikosti o každých dalších 5 vzorků přináší stále stejnou míru zlepšení intervalů“? Jinými slovy, platí, že zvýšení velikosti z 5 na 10 přinese stejný poměr zúžení intervalu, jako z 20 na 25?
    • Zkuste velikosti 10, 100, 1000. Zde by se měla s každým dalším krokem šířka intervalu snížit ve stejném poměru, jelikož šířka intervalu je úměrná vztahu \(1/\sqrt n.\) 10x větší VS přinese zúžení zhruba na třetinu, neboť \(1/\sqrt {10} \approx 32\%.\) To znamená, že chceme-li 10x lepší odhad střední hodnoty, musíme pořídit 100x větší VS.
  2. Poslechového testu se zúčastnilo 20 respondentů. Jednu položku správně rozpoznalo 15 lidí, druhou položku jen 10 lidí. Lze na hladině významnosti \(\alpha = 0.05\) říci, že první položka je lépe srozumitelná? Návod: pro oba případy sestrojte interval spolehlivosti pro binomické rozdělení. Pokud se částečně překrývají, nelze vyloučit, že pravděpodobnosti rozpoznání obou položek jsou shodné. [Odpověď: intervaly se překrývají. Rozdíl obdržených počtů je příliš malý vzhledem k možné náhodě, kterou jsou čísla zatížena, a proto nemůžeme rozhodnout, že jedna položka je určitě srozumitelnější.]

© 25. 3. 2015 Tomáš Bořil,