# .....................................................................................
# ............................    Intervalové odhady     ..............................
# ............ ..............    Martina Litschmannová    .............................
# .....................................................................................

# Nezobrazuje-li se vám text korektně, nastavte File \ Reopen with Encoding... na UTF-8
# Pro zobrazení obsahu skriptu použijte CTRL+SHIFT+O
# Pro spouštění příkazů v jednotlivých řádcích použijte CTRL+ENTER

# .....................................................................................
## Příprava prostředí ####

# Instalace knihoven (lze také provést přes okno Packages -> Install)
# Provádíme pouze jednou na daném počítači
install.packages("readxl")  # načtení xlsx souborů
install.packages("moments") # pro výpočet šikmosti a špičatosti
install.packages("dplyr")   # pro efektivní práci s datovým souborem
install.packages("tidyr")   # pro efektivní práci s datovým souborem (pivot_longer)
install.packages("ggplot2") # pro hezčí grafiku
install.packages("ggpubr")  # pro kombinování grafů z ggplot2
install.packages("rstatix") # pro identifikaci odlehlých pozorování
install.packages("forcats") # obsahuje funkci fct_infreq(), která seřadí úrovně faktorů podle četností
install.packages("EnvStats") # obsahuje funkci varTest(), která určí IO pro rozptyl
install.packages("binom")   # obsahuje funkci binom.confint(), která určí IO pro p-st

# Aktivace knihoven (nutno opakovat při každém novém spuštění Rka, vhodné mít na začátku skriptu)
library(readxl)
library(moments)
library(dplyr)
library(tidyr)
library(ggplot2)
library(ggpubr)
library(rstatix)
library(forcats)
library(EnvStats)
library(binom)

# Výpis pracovního adresáře (také je vidět v záhlaví Console)
getwd()

# Nastavení pracovního adresáře dle potřeby (lze také přes Session -> Set Working Directory -> Choose Directory)
# POZOR! Cesta musí obsahovat lomítka "/", nikoliv obrácené lomítko (backslash). 
setwd("C:/Users/Superman/SuperStatistika")

# Připravujeme-li publikaci v češtině, lze spustit příkaz, který ve VÝSTUPECH změní des. tečku na des.čárku.
# Stále ale píšeme v příkazech desetinné tečky, příkaz změní jen VÝSTUP.
# Zároveň je vhodné nastavit požadovaný počet des. míst ve výstupech generovaných základním Rkem
# i ve výstupech z balíčku dplyr.

options(OutDec= ",",           # nastavení des. čárky v grafických výstupech
        digits = 7,            # nastavení počtu des. míst ve výstupech základního Rka
        pillar.sigfig = 10)    # nastavení počtu platných cifer ve výstupech balíčku dplyr





#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------
# 1.Příklad 
# Předpokládejme, že hmotnost kaprů je náhodná veličina s normálním rozdělením pravděpodobnosti se 
# směrodatnou odchylkou σ = 0,5kg.


#-----------------------------------------------------------------------------------
# a) Měřením hmotnosti 30-ti kaprů jsme zjistili jejich průměrnou hmotnost 2,3kg. 
#    Určete 95% intervalový odhad (= interval spolehlivosti) střední hodnoty hmotnosti kaprů.

# Tzn. předpokládáme, že hmotnost kaprů ...X -> N(μ,σ^2) a σ známe.

# Odhadovaný parametr: μ

alfa = 0.05   # hladina významnosti
xprumer = 2.3 # naměřený průměr
sigma = 0.5   # směrodatná odchylka
n = 30
z = qnorm(0.975,0,1) # (1-α/2).100% kvantil rozdělení N(0,1)

# M_D ...dolní mez intervalového odhadu:

MD = xprumer - sigma*z/sqrt(n) # = 2,121081
MD

# M_H ...dolní mez intervalového odhadu:

MH = xprumer + sigma*z/sqrt(n) # = 2,478919
MH


#-----------------------------------------------------------------------------------
# b) Kolik kaprů je třeba zvážit, aby 95% IO μ : (xprůměr - Δ, xprůměr + Δ), kde
#   Δ \leq 0.1?

# n \geq

(sigma*z/0.1)^2 # = 96,03647

#-----------------------------------------------------------------------------------
# Zaokrouhlení směrodatné odchylky podle velikosti výberu (n):
# [2; 10] => počet platných cifer 1
# (10; 30] => počet platných cifer 2
# (30; 2000] => počet platných cifer 3 
# (2000; ...) => počet platných cifer 4

# Např. n = 30, tj. směr. odchylku zaokrouhlíme na 2 platné cifry 

#-----------------------------------------------------------------------------------
# intervalové odhady střední hodnoty i směr. odchylky zaokrouhlujeme na stejný řád jako směr. odchylku
# dolní mez IO zaokrouhlujeme dolů
# horní mez IO zaokrouhlujeme nahoru
#-----------------------------------------------------------------------------------                                              


#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------
#  Příklad (původně 6.) (Hemoglobin - odhad rozsahu výběru) 

# Jaký musí být počet pozorování, jestliže chceme s pravděpodobností 0,95 stanovit průměrnou
# hodnotu hemoglobinu u novorozenců s chybou nejvýše 1,0 g/l. Populační rozptyl hodnot se
# odhaduje hodnotou 46,0 g2/l2.
#

# 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 n \geq
(qnorm(1-alpha/2)*sigma/delta)^2 # = 176,7071

# => Počet pozorování musí být alespoň 177.


#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------
# 2.Příklad 

# Předpokládejme, že délka trubky je náhodná veličina X s normálním rozdělením 
# pravděpodobnosti, tj. X -> N(μ,σ^2). Změřením 18-ti trubek jsme zjistili jejich
# průměrnou délku xprum = 300cm s výběrovou směrodatnou odchylkou s = 2cm. 

#-----------------------------------------------------------------------------------
# a) Určete 90% intervalový odhad (= interval spolehlivosti) střední hodnoty EX = μ. 

n = 18  		  # rozsah souboru 
xprum = 300 	# cm.... průměr (bodový odhad střední hodnoty)
s = 2 	      # cm.... výběrová směrodatná odchylka (bodový odhad sm. odchylky)
alfa = 0.1           # hladina významnosti (spolehlivost 1-alpha = 0.9)
t = qt(1-alfa/2,n-1) # (1-α/2).100% kvantil studentova rozdělení t_(n-1)

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

DM = xprum-t*s/sqrt(n)   # dolní mez IO
HM = xprum+t*s/sqrt(n)   # horní mez IO

#Poznámky:
#-----------------------------------------------------------------------------------
# Zaokrouhlení směrodatné odchylky podle velikosti výberu (n):
# [2; 10] => počet platných cifer 1
# (10; 30] => počet platných cifer 2
# (30; 2000] => počet platných cifer 3 
# (2000; ...) => počet platných cifer 4

# n = 16, tj. směr. odchylku zaokrouhlíme na 2 platné cifry (v našem případě na desetiny)
# intervalové odhady střední hodnoty i směr. odchylky zaokrouhlujeme na stejný řád jako směr. odchylku
# dolní mez IO zaokrouhlujeme dolů
# horní mez IO zaokrouhlujeme nahoru


# 90% intervalový odhad střední hodnoty délky trubky je (299,1; 300,9) cm.
#                                                       =================


#-----------------------------------------------------------------------------------
# b) Určete 90% intervalový odhad (= interval spolehlivosti) směrodatné odchylky σ. 

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

n = 18  		  # rozsah souboru 
s = 2 	      # cm.... výběrová směrodatná odchylka (bodový odhad sm. odchylky)
alfa = 0.1           # hladina významnosti (spolehlivost 1-alpha = 0.9)

sqrt( ((n-1)*s^2) / qchisq(1-alfa/2,n-1) )    # dolní mez IO = 1,570006
sqrt( ((n-1)*s^2) / qchisq(alfa/2,n-1) )      # horní mez IO = 2,800276

# 90% intervalový odhad směrodatné odchylky délky trubky je (1,5; 2,9) cm.
#                                                       =================



#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------
# Příklad 3. (Cholesterol) 

# Úkolem je určit průměrnou hladinu cholesterolu v séru v určité populaci mužů. V náhodném
# výběru (pocházejícím z normálního rozdělení ) 25 mužů je výběrový průměr 6,3 mmol/l a výběrová
# směrodatná odchylka 1,3 mmol/l.

#-----------------------------------------------------------------------------------
# a) Určete 95% intervalový odhad (= interval spolehlivosti) střední hladiny 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 = 5,763386
x.bar+qt(1-alpha/2,n-1)*s/sqrt(n)   # horní mez IO = 6,836614


# 95% intervalový odhad střední hladiny cholesterolu je (5,8; 6,8) mmol/l.
# (výsledek ze skript)                                  ==================


#-----------------------------------------------------------------------------------
# b) Určete 95% intervalový odhad (= interval spolehlivosti) směrodatné odchylky σ. 

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

n = 25  		  # rozsah souboru 
s = 1.3 	      # mmol/l.... výběrová směrodatná odchylka (bodový odhad sm. odchylky)
alfa = 0.05           # hladina významnosti (spolehlivost 1-alpha = 0.95)

sqrt( ((n-1)*s^2) / qchisq(1-alfa/2,n-1) )    # dolní mez IO = 1,015077
sqrt( ((n-1)*s^2) / qchisq(alfa/2,n-1) )      # horní mez IO = 1,808498

# 90% intervalový odhad směrodatné odchylky délky trubky je (1,0; 1,9) mmol/l.
#                                                        ====================


#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------
# * Příklad 4. (Cholesterol 2) ####

# Předpokládejme, že v náhodném výběru 200 mladých mužů má 120 z nich vyšší než doporučenou
# hladinu cholesterolu v séru. Určete 95% interval spolehlivosti pro procento mladých mužů
# s vyšší hladinou cholesterolu v populaci.


# 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 > 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(1-alpha/2)*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

binom.confint(x,n) #Waldův odhad (asymptotic), Clopperův - Pearsonův odhad (exact)

# [2; 10] => počet platných cifer 1
# (10; 30] => počet platných cifer 2
# (30; 2000] => počet platných cifer 3 
# (2000; ...) => počet platných cifer 4



#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------
# * Příklad 5. (Hemoglobin) ####
# V rámci výzkumné studie pracujeme s náhodným výběrem 70 žen z české populace. U každé
# z žen byl změřen hemoglobin s přesností 0,1 g/100 ml. Naměřené hodnoty jsou v uvedeny
# v souboru "intervalove_odhady.xlsx", sheet = "Hemoglobin". Nalezněte 95% intervalové odhady
# směrodatné odchylky a střední hodnoty hemoglobinu v populaci českých žen. 
# (Normalitu ověřte na základě exploračních grafů, případně testováním normality.)


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

#---------------------------------------------------
## Načtení dat z xlsx souboru (pomoci balíčku readxl)
## adresu v následujícím příkazu je nutno upravit dle vlastního nastavení
hem = read_excel("intervalove_odhady.xlsx",
                 sheet = "Hemoglobin")
colnames(hem)="hodnoty"


#---------------------------------------------------
## Explorační analýza - zjistit, zda není třeba odstranit odlehlá pozorování
boxplot(hem$hodnoty)

# Data neobsahují odlehlá pozorování.
length(hem$hodnoty)
sd(hem$hodnoty)
summary(hem$hodnoty)


#---------------------------------------------------
# Ověření normality
qqnorm(hem$hodnoty)
qqline(hem$hodnoty)
hist(hem$hodnoty)

skewness(hem$hodnoty)             # koef šikmosti = 0 u normálního rozdělení; \in [-2,2] => OK
moments::kurtosis(hem$hodnoty)-3  # standardizovaný koef špičatosti = 0 u normálního rozdělení; \in [-2,2] => OK
# Š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).
#---------------------------------------------------



#-----------------------------------------------------------------------------------
# a) Určíme 95% oboustranný intervalový odhad střední hodnoty
t.test(hem$hodnoty,alternative = "two.sided",conf.level=0.95)

# Vyšlo:  (11,64594 12,31977)
#-------------------------------------------------------------
# Zaokrouhlení směrodatné odchylky podle velikosti výberu (n):
# [2; 10] => počet platných cifer 1
# (10; 30] => počet platných cifer 2
# (30; 2000] => počet platných cifer 3 
# (2000; ...) => počet platných cifer 4

# Např. n = 70, tj. směr. odchylku zaokrouhlíme na 3 platné cifry 

# intervalové odhady střední hodnoty i směr. odchylky zaokrouhlujeme na stejný řád jako směr. odchylku
# dolní mez IO zaokrouhlujeme dolů
# horní mez IO zaokrouhlujeme nahoru


# 95% IO μ : (11,64; 12,32)
#            =============



#-----------------------------------------------------------------------------------
# b) Určíme 95% oboustranný intervalový odhad směrodatné odchylky


varTest(hem$hodnoty,conf.level = 0.95) #je 95% intervalový odhad rozptylu

# Vyšlo: 95% intervalový odhad rozptylu = σ^2 : (LCL, UCL)

LCL = 1.467767 # dolní mez intervalu spolehlivosti pro rozptyl = σ^2
UCL = 2.874530 # horní mez intervalu spolehlivosti pro rozptyl = σ^2

# Vyšlo: 95% intervalový odhad rozptylu = σ : (sqrt(LCL), sqrt(UCL))

sqrt(LCL) # = 1,211514
sqrt(UCL) # = 1,695444


# 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

# Vyjde:     

#      LCL      UCL 
# 1,211514 1,695444 


# 95% IO σ : (1,21; 1,70)
#            ============








#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------
# * Příklad 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)
## adresu v následujícím příkazu je nutno upravit dle vlastního nastavení
slinivka = read_excel("intervalove_odhady.xlsx",
                      sheet = "Slinivka")   
colnames(slinivka) = c("sk1","sk2")

## Převod dat do standardního datového formátu
slinivka.s = pivot_longer(slinivka,
                          cols = 1:2,
                          values_to = "hodnoty",
                          names_to = "skupina")

slinivka.s = na.omit(slinivka.s)

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

slinivka.s %>% 
  group_by(skupina) %>% 
  summarise(n = length(na.omit(hodnoty)),
            sm.odch. = sd(hodnoty),
            prumer = mean(hodnoty),
            median = median(hodnoty))



# Ověření normality
qqnorm(slinivka$sk1)
qqline(slinivka$sk1)

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

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

hist(slinivka$sk2)

slinivka.s %>% 
  group_by(skupina) %>% 
  summarise(sikmost = skewness(hodnoty,na.rm = TRUE),
            spicatost = moments::kurtosis(hodnoty,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.
slinivka.s %>% 
  group_by(skupina) %>% 
  summarise(test.norm.p.value = shapiro.test(hodnoty)$p.value)


#Nebo testujeme normalitu u každé skupiny zvlášť:
shapiro.test(slinivka$sk1) # p-value = 0.8391 => nezamítáme normalitu
shapiro.test(slinivka$sk2) # p-value = p-value = 0.1367 => nezamítáme normalitu


# 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,var.equal = F) # Pozor na nutnost kontrolovat jak je rozdíl definován!
# nebo
t.test(slinivka$sk1,slinivka$sk2,var.equal = F,conf.level=0.95)