#######################################################################################
########################## Intervalové odhady #########################################
############## Martina Litschmannová, Veronika Kubíčková ##############################
#######################################################################################

################## Přiklad 1. (Žárovky) ###############################################
#######################################################################################

# Odhadujeme stř.hodnotu a směr.odchylku životnosti žárovek
# Součástí zadání je informace o normalitě dat

n = 16  		# rozsah souboru 
x.bar = 3000 	# hodin.... průměr (bodový odhad střední hodnoty)
s = 20 	# hodin.... výběrová směrodatná odchylka (bodový odhad sm. odchylky)
alpha = 0.1 # hladina významnosti (spolehlivost 1-alpha = 0.9)

## Oboustranný intervalový odhad střední hodnoty

x.bar-qt(1-alpha/2,n-1)*s/sqrt(n)   # dolní mez IO
x.bar-qt(1-alpha/2,n-1)*s/sqrt(n)   # horní mez IO


## Oboustranný intervalový odhad směrodatné odchylky

sqrt( ((n-1)*s^2) / qchisq(1-alpha/2,n-1) )    # dolní mez IO
sqrt( ((n-1)*s^2) / qchisq(alpha/2,n-1) )      # horní mez IO


################## Přiklad 2. (Hloubka moře) ##########################################
#######################################################################################

# Určujeme odhad potřebného rozsahu výběru (počtu potřebnych měření)

# Předpokládáme normalitu dat, se známým rozpylem (dle zadání)

sigma = 20 # metrů .... známá směrodatná odchylka
alpha = 0.05 # hladina významnosti (spolehlivost 1-alpha = 0.95)
delta = 10 # metrů ... přípustná chyba měření 

# Odhad rozsahu výběru
(qnorm(1-alpha/2)*sigma/delta)^2 


################## Přiklad 3. (Cholesterol 1) #########################################
#######################################################################################

# Odhadujeme střední hladinu cholesterolu v séru 

# Předpokládáme normalitu dat (dle zadání)

n = 25  		# rozsah souboru 
x.bar = 6.3 	# mmol/l .... průměr (bodový odhad střední hodnoty)
s = 1.3 	# mmol/l .... výběrová směrodatná odchylka (bodový odhad sm. odchylky)
alpha = 0.05 # hladina významnosti (spolehlivost 1-alpha = 0.95)

## Oboustranný intervalový odhad střední hodnoty

x.bar-qt(1-alpha/2,n-1)*s/sqrt(n)   # dolní mez IO
x.bar-qt(1-alpha/2,n-1)*s/sqrt(n)   # horní mez IO


################## Přiklad 4. (Cholesterol 2) #########################################
#######################################################################################

# Odhadujeme podíl mužů s vyšší hladinou cholesterolu v celé populaci, 
# tj. pravděpodobnost,že náhodně vybraný muž bude mít vyšší hladinu cholesterolu

n = 200  		# rozsah souboru 
x = 120 	# počet "úspěchů"
p = x/n 	# relativní četnost (bodový odhad pravděpodobnosti)
p
alpha = 0.05 # hladina významnosti (spolehlivost 1-alpha = 0.95)

# Ověření předpokladů
n>30 & n > 9/(p*(1-p))
# Dále předpokládáme  n < 0.05*N, tj. že N > n/0.05 (N > 20n), 
# tj. že daná populace (mladých mužů) má rozsah alespoň 20*200 = 4000 mužů, 
# což je asi vcelku reálný předpoklad :-)

## Oboustranný Clopperův - Pearsonův (exaktní) intervalový odhad parametru binomického rozdělení
binom.test(x,n,alternative="two.sided",conf.level=0.95)

## Waldův (asymptotický) odhad (z-statistika) - aproximace normálním rozdělením dle CLV
p-qnorm(1-alpha/2)*sqrt( (p*(1-p)) / n )   # dolní mez IO
p+qnorm(0.975)*sqrt( (p*(1-p)) / n )       # horní mez IO

## Výpočet 11 nejčastěji používaných intervalů spolehlivosti parametru bin. rozdělení
## pomocí balíčku binom

install.packages("binom")
library(binom)
binom.confint(x,n)

################## Přiklad 5. (Hemoglobin) ############################################
#######################################################################################

## Odhadujeme střední hodnotu a směrodatnou odchylku hemoglobinu v séru

## Načtení dat z xlsx souboru (pomoci balíčku readxl)
library(readxl)
## adresu v následujícím příkazu je nutno upravit dle vlastního nastavení
hem = read_excel("C:/Martina/Výuka/Pravděpodobnost a statistika/DATA/intervalove_odhady.xlsx",
                  sheet = "Hemoglobin")
colnames(hem)="hodnoty"

## Explorační analýza
boxplot(hem$hodnoty)
# Data neobsahují odlehlá pozorování.
summary(hem$hodnoty)
sd(hem$hodnoty)

# Ověření normality
qqnorm(hem$hodnoty)
qqline(hem$hodnoty)
library(moments)
skewness(hem$hodnoty)
moments::kurtosis(hem$hodnoty)-3
# Šikmost i špičatost odpovídá norm. rozdělení. Pro konečné rozhodnutí o normalitě dat použijeme
# test normality.

# Známe-li testování hypotéz, ověříme předpoklad normality Shapirovým . Wilkovým testem.
shapiro.test(hem$hodnoty)
# Na hl. významnosti 0.05 nelze předpoklad normality zamítnout (p-hodnota=0.522, Shapirův-Wilkův test).

# 95% oboustranný intervalový odhad střední hodnoty
t.test(hem$hodnoty,conf.level=0.95)

## 95% oboustranný intervalový odhad směrodatné odchylky
library(EnvStats)
varTest(hem$hodnoty,conf.level = 0.95) #95% intervalový odhad rozptylu

# Jak získat 95% intervalový odhad směrodatné odchylky?
var_hem = varTest(hem$hodnoty,conf.level = 0.95) # výsledek f-ce varTest uložíme do pomecné proměnné
names(var_hem) # zjistíme, kde najdeme informaci o intervalovém odhadu rozptylu
sqrt(var_hem$conf.int) # využijeme toho, že sm. odchylka je odmocninou z rozptylu


################## Přiklad 6. (Hemoglobin - odhad rozsahu výběru) #####################
#######################################################################################

# Určujeme odhad potřebného rozsahu výběru (počtu novorozenců, které musíme testovat)

# Předpokládáme normalitu dat, bez tohoto předpokladu je příklad neřešitelný

sigma = sqrt(46) # g/l .... známá směrodatná odchylka
alpha = 0.05 # hladina významnosti (spolehlivost 1-alpha = 0.95)
delta = 1 # g/l ... přípustná chyba měření 

# Odhad rozsahu výběru
(qnorm(1-alpha/2)*sigma/delta)^2 

################## Přiklad 7. (Poškození tkáně slinivky břišní) #######################
#######################################################################################

# Odhadujeme střední poškození tkáně pro pro dvě skupiny pacientů (podle způsobu ochlazování tkáně)
# a rozdíl středních poškození tkáně pro tyto dvě skupiny pacientů

## Načtení dat z xlsx souboru (pomoci balíčku readxl)
library(readxl)
## adresu v následujícím příkazu je nutno upravit dle vlastního nastavení
slinivka = read_excel("C:/Martina/Výuka/Pravděpodobnost a statistika/DATA/intervalove_odhady.xlsx",
                 sheet = "Slinivka")   
colnames(slinivka) = c("sk1","sk2")

## Převod dat do standardního datového formátu
slinivka.s = stack(slinivka)
colnames(slinivka.s) = c("hodnoty","skupina")

## Explorační analýza
boxplot(slinivka.s$hodnoty~slinivka.s$skupina)
# Data neobsahují odlehlá pozorování.

tapply(slinivka.s$hodnoty,slinivka.s$skupina,summary)
tapply(slinivka.s$hodnoty,slinivka.s$skupina,sd,na.rm = T)

# Ověření normality
qqnorm(slinivka$sk1)
qqline(slinivka$sk1)
# zvažte, popř. si vyzkoušejte, jak zhotovit graf přímo pomoci datové matice slinivka.s

qqnorm(slinivka$sk2)
qqline(slinivka$sk2)

library(moments)
tapply(slinivka.s$hodnoty,slinivka.s$skupina,skewness,na.rm=TRUE)
tapply(slinivka.s$hodnoty,slinivka.s$skupina,moments::kurtosis,na.rm=TRUE)-3

# Šikmost i špičatost odpovídá norm. rozdělení. Pro konečné rozhodnutí o normalitě dat použijeme
# test normality.

# Známe-li testování hypotéz, ověříme předpoklad normality Shapirovým . Wilkovým testem.
tapply(slinivka.s$hodnoty,slinivka.s$skupina,shapiro.test)
# Na hl. významnosti 0.05 nelze předpoklad normality zamítnout (p-hodnota1=0.839, p-hodnota2=0.137,
# Shapirův-Wilkův test).

# 95% oboustranné intervalové odhady středních poškození tkáně pro jednotlivé skupiny
t.test(slinivka$sk1,conf.level=0.95)
t.test(slinivka$sk2,conf.level=0.95)

# 95% oboustranný intervalový odhad rozdílu středních poškození tkáně daných skupin pacientů
t.test(slinivka.s$hodnoty~slinivka.s$skupina,conf.level=0.95) # Pozor na nutnost kontrolovat jak je rozdíl definován!
# nebo
t.test(slinivka$sk1,slinivka$sk2,conf.level=0.95)