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).
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)
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
t.test(VS1, conf.level = 1-alfa)$conf.int
## [1] 98.95998 104.30664
## attr(,"conf.level")
## [1] 0.95
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.
Čí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.
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
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
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
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
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
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.
© 25. 3. 2015 Tomáš Bořil, borilt@gmail.com