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"  "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 100 lidí jich 58 hlasovalo 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í.

Př. V percepčním testu bylo jedno slovo správně rozpoznáno 12 lidmi z celkových 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

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")
## 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, borilt@gmail.com