# ...........................Upravené výukové materiály :..............................


...........................................................
# .............. Cvičení 10. Úvod do testování hypotéz, jednovýběrové testy.............
# ......................... Michal Béreš, 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



## Inicializace potřebných balíčků.
library(readxl)   # načtení xlsx souboru
library(dplyr)    # pro efektivní preprocessing (práci s datovým rámcem)
library(ggplot2)  # pro hezčí grafiku
library(tibble)   # umožňuje nastavení počtu des. míst ve výstupech balíčku dplyr
library(EnvStats) # funkce varTest() test o rozptylu
library(moments)  # funkce skewness(), kurtosis()
library(BSDA)     # funkce SGN.test() pro test o mediánu
library(rstatix)  # funkce identify_outliers()






# * Co je to statistický test hypotéz? ####
# Mějme následující:
# - náhodná veličina X (například výška mužů)
# - výběr z náhodné veličiny (měření výšky 30 mužů)
# 
# Statistické testování hypotéz rozhoduje na základě získaných dat z náhodného výběru o
# platnosti:
# - $H_0$ - nulové hypotézy
# - $H_A$ - alternativní hypotézy
# 
# Například:
# $H_0$: $\mu_X = 175$
# $H_A$: $\mu_X > 175$
# Jelikož se jedná o statistické rozhodnutí, vždy bude vázáno k nějaké hladnině
# významnosti $\alpha$. Vždy můžeme dospět pouze k 2 různým rozhodnutím:

# - Zamítám $H_0$ ve prospěch $H_A$
#     - to znamená, že tvrdím, že $H_0$ neplatí
#     - toto rozhodnutí je s maximální chybou $\alpha$ (hladina významnosti, chyba I.
# druhu) - to znamená, že velikost této chyby jsme schopni ovlivnit

# - Nezamítám $H_0$
#     - to znamená, že tvrdím, že vzhledem k získaným datům (výběr) nelze vyvrátit $H_0$
#     - toto rozhodnutí je s chybou $\beta$ (chyba II. druhu), tato chyba není přímo
# ovlivnitelná a záleží na typu použitého testu
#     



# * Přehled testů ####


#načteme data: 
data = read_excel("10cvPASTJahoda.xlsx",
                    sheet = "data")


# ** Míry polohy ####
# Mírami polohy rozumíme údaj určující polohu dat, nehledě na tom jak jsou rozptýlená.
# Pro data z normálního rozdělení můžeme odhadovat střední hodnotu, pro ostatní medián.


# *** a) studentův t-test #####################################################
# - testujeme střední hodnotu
# - data musejí pocházet z normálního rozdělení
#     - exploračně: šikmost a špičatost leží v (-2,2)
#     - exploračně: QQ graf má body přibližně na čáře
#     - exaktně: pomocí statistického testu, např. Shapiro-Wilk test
# (shapiro.test(data))


# H_0: mu = 100
# H_A: mu <> 100
t.test(data$data, mu = 100, alternative = 'two.sided')$p.value

# H_0: mu = 100
# H_A: mu > 100
t.test(data$data, mu = 100, alternative = 'greater')$p.value

# H_0: mu = 100
# H_A: mu < 100
t.test(data$data, mu = 100, alternative = 'less')$p.value


# *** b) Wilcoxnův test #########################################################
# - testujeme medián
# - data musejí pocházet ze symetrického rozdělení
#     - exploračně: šikmost leží v (-2,2)
#     - exploračně: histogram vypadá přibližně symetricky
#     - exaktně: pomocí statistického testu, např. balíček "lawstat", funkce
# "symmetry.test(data,boot=FALSE)"


# H_0: X_0.5 = 100
# H_A: X_0.5 <> 100
wilcox.test(data$data, mu = 100, alternative = 'two.sided')$p.value

# H_0: X_0.5 = 100
# H_A: X_0.5 > 100
wilcox.test(data$data, mu = 100, alternative = 'greater')$p.value

# H_0: X_0.5 = 100
# H_A: X_0.5 < 100
wilcox.test(data$data, mu = 100, alternative = 'less')$p.value


# *** c) znaménkový test test #################################################
# - testujeme medián
# - výběr většího rozsahu (>10)
# - vyžaduje knihovnu "BSDA"
# - jakožto nejrobustnější test, se dá použít i na nespojitá data - např. pořadí v
# nějakém seznamu


# H_0: X_0.5 = 100
# H_A: X_0.5 <> 100
BSDA::SIGN.test(data$data, md = 100, alternative = 'two.sided')$p.value

# H_0: X_0.5 = 100
# H_A: X_0.5 > 100
BSDA::SIGN.test(data$data, md = 100, alternative = 'greater')$p.value

# H_0: X_0.5 = 100
# H_A: X_0.5 < 100
BSDA::SIGN.test(data$data, md = 100, alternative = 'less')$p.value


# ** Míry variability ####
# Mírami variability rozumíme údaj určující rozptýlenost/variabilitu dat, nehledě na
# celkových hodnotách. Pro data z normálního rozdělení můžeme odhadovat směrodatnou
# odchylku.

# *** test směrodatné odchylky ################################################
# - testujeme směrodatnou odchylku
# - data musejí pocházet z normálního rozdělení
#     - exploračně: šikmost a špičatost leží v (-2,2)
#     - exploračně: QQ graf má body přibližně na čáře
#     - exaktně: pomocí statistického testu, např. Shapiro-Wilk test
# (shapiro.test(data))
# - vyžaduje balíček "EnvStats"
# - funkce v Rku, porovnává rozptyl!!!


# H_0: sigma = 10
# H_A: sigma <> 10
EnvStats::varTest(data$data, sigma.squared = 10*10, 
                  alternative = 'two.sided')$p.value

# H_0: sigma = 10
# H_A: sigma > 10
EnvStats::varTest(data$data, sigma.squared = 10*10, 
                  alternative = 'greater')$p.value

# H_0: sigma = 10
# H_A: sigma < 10
EnvStats::varTest(data$data, sigma.squared = 10*10, 
                  alternative = 'less')$p.value


# * Pravděpodobnost výskytu u jednoho výběru ##################################
# *** IO pravděpodobnosti ####
# - testujeme pravděpodobnost
# - vyžadujeme dostatečný počet dat: $n>\frac{9}{p(1-p)}$
# - Clopperův - Pearsonův odhad (binom.test)
#     - jako parametr nebere data, ale počet úspěchů a počet pozorování 


pi = 0.3
data_bin = runif(n = 100, min = 0, max = 1) < pi

n = length(data_bin)
x = sum(data_bin)
n
x

# H_0: pi = 0.2
# H_A: pi <> 0.2
binom.test(x = x, n = n, p = 0.2, alternative = 'two.sided')$p.value

# H_0: pi = 0.2
# H_A: pi > 0.2
binom.test(x = x, n = n, p = 0.2, alternative = 'greater')$p.value

# H_0: pi = 0.2
# H_A: pi < 0.2
binom.test(x = x, n = n, p = 0.2, alternative = 'less')$p.value
















#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------
# 1.Příklad 
# Měřili jsme velikost tíhového zrychlení na Zemi. Provedli jsme celkem 10 měření,
# průměrná hodnota naměřených hodnot vyšla 9.71 m.s^{-2} a výběrová směrodatná odchylka 
# vyšla 0,05 m.s^{-2}. 
# Předpokládejme, že data pocházejí z normálního rozdělení (tj. že chyby měření mají 
# normální rozdělení pravděpodobnosti).
# Rozhodněte, zda je naměřená hodnota 9.71 m.s^{-2} statisticky významně menší, než 
# hodnota tíhového zrychlení 8,81 m.s^{-2} udávaná pro tuto zeměpisnou šířku.
#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------

# Použijeme jednovýběrový t-test

# Předpoklady testu: Normalita dat - podle zadání splněno

# 1.) NULOVÁ A ALTERNATIVNÍ HYPOTÉZA:

#     H_0:  μ = 9,81
#     H_A:  μ < 9,81

# 2.) U KLASICKÉHO TESTU VOLÍME HLADINU VÝZNAMNOSTI 
#    (U ČISTÉHO TESTU VÝZNAMNOSTI NEURČUJEME PŘEDEM, ROZHODUJEME SE PODLE p-hodnoty):

#     α = 0,05

alfa=0.05

# 3.) TESTOVÁ STATISTIKA:  

#     T = (\bar{X} -  μ)sqrt(n)/S --> t_{n-1}



# 4.) URČÍME OBOR PŘIJETÍ H_0 A JEHO DOPLNĚK KRITICKÝ OBOR 
#     (U ČISTÉHO TESTU VÝZNAMNOSTI NEURČUJEME, ROZHODUJEME SE PODLE p-hodnoty)
#     

n = 10
#     obor přijetí: u této alternativy (5% kvantil rozdělení t_{n-1} až nekonečno):

t1=qt(alfa,n-1) # = -1.833113 ...dolní mez oboru přijetí
                #                horní mez oboru přijetí = nekonečno

# 5.) VÝPOČET POZOROVANÉ HODNOTY TESTOVÉ STATISTIKY
#    

#     T_{obs} = (\bar{X} -  μ)sqrt(n)/S = 
prumer =  9.71
s = 0.05 
n = 10
mi = 9.81

Tobs = (prumer - mi)*sqrt(n)/s # = -6.324555

#Můžeme si nakreslit:----------------------------------------------------------
T = seq(from = -8, to = 8, by = 0.01)
f = dt(T, n-1)
t1=qt(alfa,n-1) #dolní mez oboru přijetí
y1 = dt(t1,n-1)
y3 = dt(Tobs,n-1)
plot(T, f, ylab='f_T',ylim=c(0,0.5),xlim=c(-8,8) ,main="t test", col=c("black"), type='s')
grid()
par(new=TRUE) # že chceme kreslit do jednoho grafu
lines(c(t1,t1),c(0,2*y1), pch=4,  ylab='',ylim=c(0,0.5),xlim=c(-8,8), main="", col=c("blue"), type='l')
lines(c(Tobs,Tobs),c(0,0.2), pch=4,  ylab='',ylim=c(0,0.5),xlim=c(-8,8), main="", col=c("orange"), type='l')
#-------------------------------------------------------------------------------


# 6.) U KLASICKÉHO TESTU ROZHODNEME O VÝSLEDKU TESTU:

# Tobs = -0.6324555 leží v oboru přijetí (-1.833113, ∞), proto


# NEZAMÍTÁME NULOVOU HYPOTÉZU H_0, ŽE STŘEDNÍ HODNOTA MĚŘENÍ JE ROVNA 9,81m.s^{-1}.
# =========================================================================



# 7.) U ČISTÉHO TESTU VÝZNAMNOSTI URČÍME p-hodnotu (p-value) A ROZHODNEME

# p-hodnota je u této alternativy rovna F(Tobs)

FTobs= pt(Tobs,n-1) # = 6.846828e-05


# phodnota = F(Tobs) = 6.846828e-05 < 0.05 => 


# Na hladině významnosti alfa=0.05 ZAMÍTÁME H_0, 
# ===============================================
# naměřená hodnota JE statisticky významně menší než udávaná.
# ============================================================





#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------
# 2.Příklad 
# Výrobce tvrdí, že jeho motory typu A dosahují v průměru maximálního výkonu 100 kW. 
#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------

#-----------------------------------------------------------------------------------
# a) Testujte hypotézu, že výrobce říká pravdu, to jest, že střední hodnota 
#    výkonů motorů typu A je rovna 100kW. Naměřené hodnoty výkonů jsou v souboru:
#    "10cvPASTJahoda.xlsx" list "motory". 


#Studentův t-test ########################################################
# - testujeme střední hodnotu
# - data musejí pocházet z normálního rozdělení
#     - exploračně: šikmost a špičatost leží v (-2,2)
#     - exploračně: QQ graf má body přibližně na čáře
#     - exaktně: pomocí statistického testu, např. Shapiro-Wilk test
# (shapiro.test(data))
###########################################################################

#  OVĚŘÍME, ŽE JSOU SPLNĚNY PODMÍNKY PRO POUŽITÍ TESTU:

#     X...výkon motoru typu A musí mít normální rozdělení pravděpodobnosti => OVĚŘÍME:

#     I.) Načteme data:

vykony = read_excel("10cvPASTJahoda.xlsx",
                    sheet = "motory")

#     II.) Identifikace odlehlých pozorování

boxplot(vykony$vykon)  # Data neobsahují odlehlá pozorování => 
# OK, jinak bychom museli zvážit jeich odstranění.

#     III.) Ověření normality - exploračně

moments::skewness(vykony$vykon)    # šikmost -0.3286819 => OK (musí být mezi -2 a 2)
moments::kurtosis(vykony$vykon)-3  # špičatost -0.4337902 => OK (musí být mezi -2 a 2)

qqnorm(vykony$vykon)
qqline(vykony$vykon) # asi OK body jsou cca na linii

hist(vykony$vykon) # asi OK, připomíná gaussovu křivku

#     IV.) Ověření normality - testem 

#Pro konečné rozhodnutí o normalitě dat použijeme test normality.
# Předpoklad normality ověříme Shapirovovým - Wilkovovým testem.

# H0: Data jsou výběrem z normálního rozdělení.
# Ha: Data nejsou výběrem z normálního rozdělení.

shapiro.test(vykony$vykon)
# p-value = 0.5374 > 0.05  => Na hl. významnosti 0,05 nelze předpoklad normality zamítnout.

# NEZAMÍTÁME NORMALITU n.v. X...výkon motoru typu A, můžeme proto pokračovat:




# 1.) NULOVÁ A ALTERNATIVNÍ HYPOTÉZA:

#     H_0:  μ = 100
#     H_A:  μ <> 100

# 2.) U KLASICKÉHO TESTU VOLÍME HLADINU VÝZNAMNOSTI 
#    (U ČISTÉHO TESTU VÝZNAMNOSTI NEURČUJEME PŘEDEM, ROZHODUJEME SE PODLE p-hodnoty):

#     α = 0,05

alfa=0.05

# 3.) TESTOVÁ STATISTIKA:  

#     T = (\bar{X} -  μ)sqrt(n)/S --> t_{n-1}


             
# 4.) URČÍME OBOR PŘIJETÍ H_0 A JEHO DOPLNĚK KRITICKÝ OBOR 
#     (U ČISTÉHO TESTU VÝZNAMNOSTI NEURČUJEME, ROZHODUJEME SE PODLE p-hodnoty)
#     
             
#     obor přijetí: (2,5% kvantil rozdělení t_{n-1} až 97,5% kvantil rozdělení t_{n-1}):
      qt(0.025,29) # = -2.04523
      qt(0.975,29) # =  2.04523
      
# 5.) VÝPOČET POZOROVANÉ HODNOTY TESTOVÉ STATISTIKY
#    
      
#     T_{obs} = (\bar{X} -  μ)sqrt(n)/S = 
prumer =  mean(vykony$vykon) # = 99.18333 
smerodch = sd(vykony$vykon) # = 11.30856
n = length(vykony$vykon) # = 30
mi = 100

Tobs = (prumer - mi)*sqrt(n)/smerodch # = -0.3955472

#Můžeme si nakreslit:----------------------------------------------------------
T = seq(from = -5, to = 5, by = 0.01)
f = dt(T, n-1)
t1=qt(alfa/2,n-1) #dolní mez oboru přijetí
t2=qt(1-alfa/2,n-1) #horní mez oboru přijetí
y1 = dt(t1,n-1)
y2 = dt(t2,n-1)
y3 = dt(Tobs,n-1)
plot(T, f, ylab='f_T',ylim=c(0,0.5),xlim=c(-5,5) ,main="t test", col=c("black"), type='s')
grid()
par(new=TRUE) # že chceme kreslit do jednoho grafu
lines(c(t1,t1),c(0,2*y1), pch=4,  ylab='',ylim=c(0,0.5),xlim=c(-5,5), main="", col=c("blue"), type='l')
lines(c(t2,t2),c(0,2*y2), pch=4,  ylab='',ylim=c(0,0.5),xlim=c(-5,5), main="", col=c("blue"), type='l')
lines(c(Tobs,Tobs),c(0,2*y3), pch=4,  ylab='',ylim=c(0,0.5),xlim=c(-5,5), main="", col=c("orange"), type='l')
#-------------------------------------------------------------------------------


# 6.) U KLASICKÉHO TESTU ROZHODNEME O VÝSLEDKU TESTU:

# Tobs = -0.3955472 leží v oboru přijetí (-2.04523, 2.04523 ), proto


# NEZAMÍTÁME NULOVOU HYPOTÉZU H_0, ŽE STŘEDNÍ HODNOTA VÝKONŮ JE ROVNA 100KW.
# =========================================================================



# 7.) U ČISTÉHO TESTU VÝZNAMNOSTI URČÍME p-hodnotu (p-value) A ROZHODNEME

# p-hodnota je u oboustranného testu rovna 2*minimum{F(Tobs), 1-F(Tobs)}

FTobs= pt(Tobs,n-1) # = 0.347667
1-FTobs             # = 0.652333

phodnota = 2*FTobs  # = 0.695334 > 0.05 => NEZAMÍTÁME H_0
#                                          ==============



# NEBO TO CELÉ NECHÁME SPOČÍTAT R-ko:

t.test(vykony$vykon, mu=100, alternative = "two.sided")


# p-value = 0.6953 > 0.05 -> Na hl. významnosti 0,05 NEZAMÍTÁME nulovou hypotézu 
# ve prospěch hypotézy alternativní. 
# Střední hodnota výkonů se statisticky významně neliší od 100 kW.



#-----------------------------------------------------------------------------------
# b) Pomocí čistého testu významnosti ověřte, zda je střední hodnota výkonů 
#    motorů typu A statisticky významně menší než 105kW.

# V a) jsme ověřili podmínky pro použití t-testu => můžeme pokračovat

# Alternativou ke klasickému testu (kde vytváříme IO - v terminologii klasických testů
# tzv. obor přijetí a jeho doplněk do R kritický obor) je tzv. čistý test významnosti.
# Výsledkem čistého testu významnosti je p-hodnota. Na jejím základě rozhodujeme o
# zamítnutí či nezamítnutí $H_0$.
# p-hodnota se dá chápat jako nejvyšší možná hladina významnosti, taková aby naše
# rozhodnutí bylo - nezamítám. Tedy IO/obor přijetí by obsahoval zkoumanou hodnotu 105.

# 1.) NULOVÁ A ALTERNATIVNÍ HYPOTÉZA:

#     H_0:  μ = 105
#     H_A:  μ < 100

# 2.) TESTOVÁ STATISTIKA:  

#     T = (\bar{X} -  μ)sqrt(n)/S --> t_{n-1}


# 4.) VÝPOČET POZOROVANÉ HODNOTY TESTOVÉ STATISTIKY

# ve skutečnosti nemusíme tutro hodnotu počítat - udělá to za nás R-ko v 5.)...

#      T_{obs} = (\bar{X} -  μ)sqrt(n)/S = 
prumer =  mean(vykony$vykon) # = 99.18333 
smerodch = sd(vykony$vykon) # = 11.30856
n = length(vykony$vykon) # = 30
mi = 105

Tobs = (prumer - mi)*sqrt(n)/smerodch # = -2.817265


#Můžeme si nakreslit:----------------------------------------------------------
T = seq(from = -5, to = 5, by = 0.01)
f = dt(T, n-1)
t1=qt(alfa,n-1) #dolní mez oboru přijetí
#horní mez oboru přijetí je nekonečno
y1 = dt(t1,n-1)
y2 = dt(t2,n-1)
y3 = dt(Tobs,n-1)
plot(T, f, ylab='f_T',ylim=c(0,0.5),xlim=c(-5,5) ,main="t test", col=c("black"), type='s')
grid()
# že chceme kreslit do jednoho grafu
lines(c(t1,t1),c(0,2*y1), pch=4,  ylab='',ylim=c(0,0.5),xlim=c(-5,5), main="", col=c("blue"), type='l')
lines(c(Tobs,Tobs),c(0,0.2), pch=4,  ylab='',ylim=c(0,0.5),xlim=c(-5,5), main="", col=c("orange"), type='l')



# 5.) VÝPOČET HODNOTY p-value A ROZHODNUTÍ

# p-hodnota je u jednostranného testu (H_A:  μ < 100) rovna F(Tobs)

FTobs= pt(Tobs,n-1) # = 0.004314271

phodnota = FTobs  # = 0.004314271 < 0.05 => ZAMÍTÁME H_0 ve prospěch H_A,


# Střední hodnota motorů typu A je ( na hladině významnosti vyšší než 0.004314271) 
# ========================================================================
# statisticky významně menší, než 105kW.
# ======================================


#Tzn. Nulovou hypotzézu můžeme zamítnout na jakékoli hladině významnosti vyšší než 0.004314271

# NEBO TO CELÉ NECHÁME SPOČÍTAT R-ko:

 t.test(vykony$vykon, mu=105, alternative = "less")

# p-value = 0.004314 < 0.05 -> Na hl. významnosti 0,05 zamítáme nulovou hypotézu 
# ve prospěch hypotézy alternativní 
# Střední hodnota motorů typu A je statisticky významně menší, než 105kW
# ========================================================================

 


#-----------------------------------------------------------------------------------
# c) Pomocí čistého testu významnosti ověřte, zda je střední hodnota výkonů 
#    motorů typu A je vyšší než 105kW.
 
# V a) jsme ověřili podmínky pro použití t-testu => můžeme pokračovat
 
# výpočet pozorované statistiky Tobs by byl stejný jako v b) -> stejná hodnota
# 
 

# NEBO TO CELÉ NECHÁME SPOČÍTAT R-ko:

t.test(vykony$vykon, mu=105, alternative = "greater")

# p-value = 0.9957 > 0.05 -> Na hl. významnosti 0,05 nezamítáme nulovou hypotézu 
# ve prospěch hypotézy alternativní.
# ========================================================================


# Všimněte si, že tato jednostranná alternativa vedla k "nezamítnutí" $H_0$, že μ = 105.
# Je to z důvodu porovnávání nepravděpodobné $H_0$ s ještě méně pravděpodobnou $H_A$.





#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------
#  3. Příklad 
# Máme výběr 218 pacientů a změřili jsme jejich bílkovinné sérum (soubor
# testy_jednovyberove.xlsx list bilk_serum). Ověřte, zda se průměrné bílkovinné sérum
# (Albumin) všech pacientů tohoto typu (populační průměr µ) statisticky významně liší od
# hodnoty 35 g/l.
#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------



#Studentův t-test ########################################################
# - testujeme střední hodnotu
# - data musejí pocházet z normálního rozdělení
#     - exploračně: šikmost a špičatost leží v (-2,2)
#     - exploračně: QQ graf má body přibližně na čáře
#     - exaktně: pomocí statistického testu, např. Shapiro-Wilk test
# (shapiro.test(data))
###########################################################################

#  OVĚŘÍME, ŽE JSOU SPLNĚNY PODMÍNKY PRO POUŽITÍ TESTU:

#     X...obsah albuminu v krvi musí mít normální rozdělení pravděpodobnosti => OVĚŘÍME:

#     I.) Načteme data:

albumin = read_excel("10cvPASTJahoda.xlsx",
                     sheet = "albumin")
head(albumin)

colnames(albumin)="hodnoty"

#     II.) Identifikace odlehlých pozorování

boxplot(albumin$hodnoty)  # Data neobsahují odlehlá pozorování => 
# OK, jinak bychom museli zvážit jeich odstranění.

#     III.) Ověření normality - exploračně

moments::skewness(albumin$hodnoty)    # šikmost -0.3286819 => OK (musí být mezi -2 a 2)
moments::kurtosis(albumin$hodnoty)-3  # špičatost -0.4337902 => OK (musí být mezi -2 a 2)

qqnorm(albumin$hodnoty)
qqline(albumin$hodnoty) # asi OK body jsou cca na linii

hist(albumin$hodnoty) # asi OK, připomíná gaussovu křivku

#     IV.) Ověření normality - testem 

#Pro konečné rozhodnutí o normalitě dat použijeme test normality.
# Předpoklad normality ověříme Shapirovovým - Wilkovovým testem.

# H0: Data jsou výběrem z normálního rozdělení.
# Ha: Data nejsou výběrem z normálního rozdělení.

shapiro.test(albumin$hodnoty)
# p-value = 0.2358 > 0.05  => Na hl. významnosti 0,05 nelze předpoklad normality zamítnout.

#======================
# NEZAMÍTÁME NORMALITU n.v. X...obsah albuminu v krvi, můžeme proto pokračovat:
#======================


#___________________________________________________________________
#POZNÁMKA K ZAOKROUHLOVÁNÍ:
# 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
length(albumin$hodnoty) # = 218 => sd zaokrouhlujeme na 3 platné cifry
sd(albumin$hodnoty)     # = 0.393683  sd a míry polohy zaokrouhlujeme na tisíciny
#_____________________________________________________________________

mean(albumin$hodnoty)

# normalita OK -> t.test

# H0: mu = 35 g/l
# Ha: mu <> 35 g/l

xprum = mean(albumin$hodnoty) # = 34.48677
mi = 35
s = sd(albumin$hodnoty) # = 0.393683
n = length(albumin$hodnoty) # = 218



Tobs = (xprum - mi)*sqrt(n)/s # = -2.817265


#Můžeme si nakreslit:----------------------------------------------------------
alfa = 0.05
T = seq(from = -25, to = 25, by = 0.01)
f = dt(T, n-1)
t1=qt(alfa/2,n-1) #dolní mez oboru přijetí
t2=qt(1-alfa/2,n-1)#horní mez oboru přijetí je nekonečno
y1 = dt(t1,n-1)
y2 = dt(t2,n-1)
y3 = dt(Tobs,n-1)
plot(T, f, ylab='f_T',ylim=c(0,0.5),xlim=c(-25,25) ,main="t test", col=c("black"), type='s')
grid()
# že chceme kreslit do jednoho grafu
lines(c(t1,t1),c(0,2*y1), pch=4,  ylab='',ylim=c(0,0.5),xlim=c(-25,25), main="", col=c("blue"), type='l')
lines(c(t2,t2),c(0,2*y2), pch=4,  ylab='',ylim=c(0,0.5),xlim=c(-25,25), main="", col=c("blue"), type='l')
lines(c(Tobs,Tobs),c(0,0.2), pch=4,  ylab='',ylim=c(0,0.5),xlim=c(-25,25), main="", col=c("orange"), type='l')
#Můžeme si nakreslit:----------------------------------------------------------

phodnota = 2*pt(Tobs,n-1) # = 7.948785e-49 << 0.05 => zamítáme nulovou hypotézu ve prospěch alternativy


t.test(albumin$hodnoty, mu=35, alternative = "two.sided")$p.value # => zamítáme Ho
t.test(albumin$hodnoty, mu=35, alternative = "greater")# => nezamítáme Ho
t.test(albumin$hodnoty, mu=35, alternative = "less") # => zamítáme Ho

# u prostředního testu nebyla zamítnuta Ho proto, že alternativa je ještě méně přijatelná!


# VÝSLEDEK OBOUSTRANNÉHO TESTU:

# p-value << 0.05 -> Na hl. významnosti 0,05 zamítáme nulovou hypotézu 
# ve prospěch hypotézy alternativní:
# Střední hodnota albuminu se statisticky významně liší od 35 g/l.


# VÝSLEDEK JEDNOSTRANNÝCH TESTŮ:

# p-value << 0.05 -> Na hl. významnosti 0,05 zamítáme nulovou hypotézu 
# ve prospěch hypotézy alternativní:
# Střední hodnota albuminu je statisticky významně MENŠÍ od 35 g/l.








#-------------------------Jednovýběrový t-test/ Znaménkový test---------------------
#-----------------------------------------------------------------------------------
# 4.Příklad 
# V souboru testy_jednovyberove.xlsx list preziti jsou uvedeny doby přežití pro 100
# pacientů s rakovinou plic léčených novým lékem. Z předchozích studií je známo, že
# průměrné přežití takových pacientů bez podávání nového léku je 22,2 měsíce. Lze na
# základě těchto dat usoudit, že nový lék prodlužuje přežití?


# Načtení dat z xlsx souboru (pomoci balíčku readxl)
preziti = readxl::read_excel("10cvPASTJahoda.xlsx",
                             sheet = "preziti")   


head(preziti)

colnames(preziti)="hodnoty"

## Explorační analýza
par(mfrow = c(1, 2))          # matice grafů 1x2

boxplot(preziti$hodnoty)
hist(preziti$hodnoty)

# **Data obsahují OP -> můžeme je odstranit. Nebo si také všimnout, že se pravděpdobně
# jedná o exponenciální rozdělení a OP tam ve skutečnosti nejsou (rozdělení se tak
# prostě chová).**


# Data obsahují odlehlá pozorování. Pomoci f-ce boxplot je umíme vypsat.
pom=boxplot(preziti$hodnoty, plot = FALSE)
pom$out
# rozhodli-li jsme se pro odstranění odlehlých hodnot, pak 
preziti$hodnoty.bez=preziti$hodnoty # doporučujeme nepřepisovat původní data
preziti$hodnoty.bez[preziti$hodnoty %in% pom$out]=NA

## Explorační analýza pro data bez odlehlých pozorování
boxplot(preziti$hodnoty.bez)
summary(preziti$hodnoty.bez,na.rm=TRUE)

length(na.omit(preziti$hodnoty.bez))   # sd zaokrouhlujeme na 3 platné cifry
sd(preziti$hodnoty.bez,na.rm=TRUE)     # sd a míry polohy zaokr. na desetiny

# **Test o míře polohy (střední hodnotě / mediánu)**


# Ověření normality - exploračně
moments::skewness(preziti$hodnoty.bez,na.rm=TRUE) # = 0.8330227
moments::kurtosis(preziti$hodnoty.bez,na.rm=TRUE)-3 # = -0.1748022 

par(mfrow = c(1, 2))          # matice grafů 1x2

qqnorm(preziti$hodnoty.bez)
qqline(preziti$hodnoty.bez)
hist(preziti$hodnoty.bez)

# QQ - graf i hist. ukazují, že výběr pravd. není výběrem z norm. rozdělení.
# Šikmost i špičatost odpovídá norm. rozdělení.  
# použijeme test normality.

# Předpoklad normality ověříme Shapirovovým . Wilkovovým testem.
shapiro.test(preziti$hodnoty.bez)
# p-value < 0.05 -> Na hl. významnosti 0.05 zamítáme předpoklad normality 



# KDYBY DATA POCHÁZELA ZE SYMETRICKÉHO ROZDĚLENÍ, MOHLI BYCHOM POUŽÍT WILCOXONŮV
# TEST => explorační posouzení symetrie - výše hist. a šikmost 

# Předpoklad symetrie - ověření testem
# H0: data pocházejí ze symetrického rozdělení
# HA: ~H0

lawstat::symmetry.test(preziti$hodnoty.bez,boot=FALSE)
# p-value = 0.001774 < 0.05 -> Na hl. významnosti 0.05 zamítáme předpoklad symetrie 



# normalita zamítnuta -> symetrie zamítnuta -> Sign. test
# H0: median = 22,2 měsíců
# Ha: median > 22,2 měsíců

BSDA::SIGN.test(preziti$hodnoty.bez, md=22.2,
                alternative="greater", conf.level=0.95)

# p-value = 0.9393 > 0.05 -> Na hl. významnosti 0,05 nelze zamítnout nulovou hypotézu 
# Medián doby přežití není statisticky významně větší než 22,2 měsíců.

median(preziti$hodnoty.bez, na.rm = TRUE) # namereny median = 20 => meli jsme pouzit jinou alternativu

# H0: median = 22,2 měsíců
# Ha: median < 22,2 měsíců

BSDA::SIGN.test(preziti$hodnoty.bez, md=22.2,
                alternative="less", conf.level=0.95)


# p-value = 0.08983 > 0.05 -> Na hl. významnosti 0,05 nelze zamítnout nulovou hypotézu 
# Medián doby přežití není ani statisticky významně menší než 22,2 měsíců.


BSDA::SIGN.test(preziti$hodnoty.bez, md=22.2,"less",0.95)









#----------------------------Jednovýběrový F-test-----------------------------------
#-----------------------------------------------------------------------------------
# 4.Příklad 
# Výběrová směrodatná odchylka maximálních výkonů jednotlivých motorů typu A 
# by neměla být větší než 5 kW. Testujte hypotézu, že tento požadavek byl splněn. 
#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------

#  test směrodatné odchylky  (Jednovýběrový F-test)#################################
# - testujeme směrodatnou odchylku
# - data musejí pocházet z normálního rozdělení
#     - exploračně: šikmost a špičatost leží v (-2,2)
#     - exploračně: QQ graf má body přibližně na čáře
#     - exaktně: pomocí statistického testu, např. Shapiro-Wilk test
# (shapiro.test(data))
# - vyžaduje balíček "EnvStats"
# - funkce v Rku, porovnává rozptyl!!!



# PŘEDPOKLADY TESTU: Jde o stejná data jako v 1.Příkladu => NORMALITA DAT JE JIŽ OVĚŘENA => OK 

# Načteme data:
vykony = read_excel("10cvPASTJahoda.xlsx",
                    sheet = "motory")

# 1.) NULOVÁ A ALTERNATIVNÍ HYPOTÉZA:

#     H_0:   σ  = 5
#     H_A:   σ  > 5

# 2.) TESTOVÁ STATISTIKA:  

#     CH = (n-1)*s^2/σ^2 --> χ_{n-1} ... chí kvadrát rozdělení s n-1 stupni volnosti


# 4.) VÝPOČET POZOROVANÉ HODNOTY TESTOVÉ STATISTIKY

n = length(vykony$vykon) # = 30
s = sd(vykony$vykon) # = 11.30856
sig = 5

#      CH_{obs} = (n-1)*s/σ = 
CHobs = (n-1)*(s)^2/(sig)^2 # = 148.3448


#Můžeme si nakreslit:----------------------------------------------------------
alfa = 0.05  #zvolená hladina významnosti
x = seq(from = -1, to = 150, by = 0.01)
f = dchisq(x, n-1)
t2=qchisq(1-alfa,n-1) #horní mez oboru přijetí
                      #dolní mez oboru přijetí je nula
y2 = dchisq(t2,n-1)
y3 = dchisq(CHobs,n-1)
plot(x, f, ylab='f_CH',ylim=c(0,0.2),xlim=c(-1,150) ,main="F test", col=c("black"), type='s')
grid()
# že chceme kreslit do jednoho grafu
lines(c(t2,t2),c(0,2*y2), pch=4,  ylab='',ylim=c(0,0.2),xlim=c(-1,150), main="", col=c("blue"), type='l')
lines(c(CHobs,CHobs),c(0,0.2), pch=4,  ylab='',ylim=c(0,0.2),xlim=c(-1,150), main="", col=c("orange"), type='l')
#Můžeme si nakreslit-----------------------------------------------------------


# 5.) VÝPOČET HODNOTY p-value A ROZHODNUTÍ

# p-hodnota je u jednostranného testu (H_A:   σ  > 5 ) rovna 1 - F(CHobs)

FCHobs= 1 - pchisq(CHobs,n-1) # = 0 ...velmi malé číslo zaokrouhleno na 0

phodnota = FCHobs  # = cca 0 < 0.05 => ZAMÍTÁME H_0 ve prospěch H_A,


# Směrodatná odchylka výkonů motorů typu A je 
# ===========================================
# statisticky významně větší, než 105kW.
# ======================================


# NEBO TO CELÉ NECHÁME SPOČÍTAT R-ko:
# - vyžaduje balíček "EnvStats"
# - funkce v Rku, testuje rozptyl!!!


# H_0: sigma = 5
# H_A: sigma > 5
EnvStats::varTest(vykony$vykon, sigma.squared = 5*5, 
                  alternative = 'greater')$p.value





#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------
# 5. Příklad  ####
# Automat vyrábí pístové kroužky o daném průměru. Výrobce udává, že směrodatná odchylka
# průměru kroužku je 0,05 mm. K ověření této informace bylo náhodně vybráno 80 kroužků a
# byl změřen jejich průměr (soubor testy_jednovyberove.xlsx list krouzky). Lze zjištěné
# výsledky považovat za statisticky významné ve smyslu zlepšení kvality produkce? Ověřte
# čistým testem významnosti.
#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------


#  test směrodatné odchylky  (Jednovýběrový F-test)#############################
# - testujeme směrodatnou odchylku
# - data musejí pocházet z normálního rozdělení
#     - exploračně: šikmost a špičatost leží v (-2,2)
#     - exploračně: QQ graf má body přibližně na čáře
#     - exaktně: pomocí statistického testu, např. Shapiro-Wilk test
# (shapiro.test(data))
# - vyžaduje balíček "EnvStats"
# - funkce v Rku, porovnává rozptyl!!!
################################################################################


# NAČTENÍ DAT z xlsx souboru 
#------------------------------------------------------
krouzky = read_excel("10cvPASTJahoda.xlsx",
                     sheet = "krouzky")

head(krouzky)

colnames(krouzky)="hodnoty"


#IDENTIFIKACE ODLEHLÝCH POZOROVÁNÍ
#------------------------------------------------------
## Explorační analýza
boxplot(krouzky$hodnoty)

# Data obsahují odlehlá pozorování. Pomoci f-ce boxplot je umíme vypsat.
pom = boxplot(krouzky$hodnoty, plot = FALSE)
pom$out


# rozhodli-li jsme se pro odstranění odlehlých hodnot, pak 
krouzky$hodnoty.bez = krouzky$hodnoty
krouzky$hodnoty.bez[krouzky$hodnoty %in% pom$out] = NA



#OVĚŘENÍ PŘEDPOKLADŮ TESTU - NORMALITA DAT
#------------------------------------------------------

# Explorační analýza pro data bez odlehlých pozorování
summary(krouzky$hodnoty.bez,na.rm=TRUE)
boxplot(krouzky$hodnoty.bez)

length(na.omit(krouzky$hodnoty.bez))# = 78 sd zaokrouhlujeme na 3 platné cifry
sd(krouzky$hodnoty.bez,na.rm=TRUE)  # = 0.02484618 sd a míry polohy zaokr. na tisíciny

# Ověření normality - exploračně
moments::skewness(krouzky$hodnoty.bez,na.rm=TRUE) # = 0.05492792
moments::kurtosis(krouzky$hodnoty.bez,na.rm=TRUE)-3 # = -0.3960399
# Šikmost i špičatost odpovídá norm. rozdělení. 

qqnorm(krouzky$hodnoty.bez)
qqline(krouzky$hodnoty.bez)
hist(krouzky$hodnoty.bez)


# Pro konečné rozhodnutí o normalitě dat použijeme test normality
# Předpoklad normality ověříme Shapirovovým . Wilkovovým testem:
shapiro.test(krouzky$hodnoty.bez)
# p-value = 0.0628 > 0.05 -> Na hl. významnosti 0,05 nelze předpoklad norm. zamítnout 



# test na míru variability -> test o rozptylu
#------------------------------------------------------

# 1.) NULOVÁ A ALTERNATIVNÍ HYPOTÉZA:

#     H_0:   σ  = 5
#     H_A:   σ  > 5

# 2.) TESTOVÁ STATISTIKA:  

#     CH = (n-1)*s^2/σ^2 --> χ_{n-1} ... chí kvadrát rozdělení s n-1 stupni volnosti


# 4.) VÝPOČET POZOROVANÉ HODNOTY TESTOVÉ STATISTIKY

n = length(na.omit(krouzky$hodnoty.bez))# = 78 sd zaokrouhlujeme na 3 platné cifry
s = sd(krouzky$hodnoty.bez,na.rm=TRUE)  # = 0.02484618
sig = 0.05

#      CH_{obs} = (n-1)*s/σ = 
CHobs = (n-1)*(s)^2/(sig)^2 # = 19.01385


#Můžeme si nakreslit:----------------------------------------------------------
alfa = 0.05  #zvolená hladina významnosti
x = seq(from = -1, to = 150, by = 0.01)
f = dchisq(x, n-1)
t1=qchisq(alfa,n-1) #horní mez oboru přijetí
#dolní mez oboru přijetí je nula
y1 = dchisq(t1,n-1)
y3 = dchisq(CHobs,n-1)
plot(x, f, ylab='f_CH',ylim=c(0,0.06),xlim=c(-1,150) ,main="F test", col=c("black"), type='s')
grid()
# že chceme kreslit do jednoho grafu
lines(c(t1,t1),c(0,2*y1), pch=4,  ylab='',ylim=c(0,0.06),xlim=c(-1,150), main="", col=c("blue"), type='l')
lines(c(CHobs,CHobs),c(0,0.03), pch=4,  ylab='',ylim=c(0,0.06),xlim=c(-1,150), main="", col=c("orange"), type='l')
#Můžeme si nakreslit-----------------------------------------------------------


# 5.) VÝPOČET HODNOTY p-value A ROZHODNUTÍ

# p-hodnota je u jednostranného testu (H_A:   σ  < 0.05 ) rovna F(CHobs)

FCHobs= pchisq(CHobs,n-1) # = 0 ...velmi malé číslo zaokrouhleno na 0

phodnota = FCHobs  # = 1.353973e-12 < 0.05 => ZAMÍTÁME H_0 ve prospěch H_A,


# Směrodatná odchylka průměrů kroužků je 
# ===========================================
# statisticky významně menší, než 0,05mm.
# ======================================


# NEBO TO CELÉ NECHÁME SPOČÍTAT R-ko:
# - vyžaduje balíček "EnvStats"
# - funkce v Rku, testuje rozptyl!!!


#následující test je vlastně testem pro rozptyl = směrodatná odchylka na druhou
EnvStats::varTest(krouzky$hodnoty.bez, sigma.squared = 0.05^2,
                  alternative = "less")

# p-value =  1.353973e-12 << 0.05 -> 


# Na hladině významnosti 0,05 zamítáme H0 ve prospěch Ha 
# Na hladině významnosti 0,05 tvrdíme, že se produkce statisticky významně zlepšila
# =================================================================================





#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------
# 6. Příklad  
# Matematický model udává, že pravděpodobnost pooperačních komplikací je 30%.
# Testujte hypotézu, že model je správný, víte-li, že z 50 pacientů mělo 18
# pooperační komplikace.
#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------

# p = 18/50 = 0.36 naměřená relativní četnost komplikací

##Test pomocí CLV - PŘEDPOKLADY: n = 50 > 30, n= 50 > 9/(p*(1-p)) = cca 39
# 1.) NULOVÁ A ALTERNATIVNÍ HYPOTÉZA:

#     H_0:   π  = 0,3
#     H_A:   π  > 0,3

# 2.) TESTOVÁ STATISTIKA:  

#     Z = (P-π)sqrt(n)/sqrt(π(1-π)) --> N(0,1) ... normální rozdělení 


# 3.) VÝPOČET POZOROVANÉ HODNOTY TESTOVÉ STATISTIKY

pinula = 0.3
n = 50
x = 18
p= 18/50

#      Z_{obs} =  (p-π_0)sqrt(n)/sqrt(π_0(1-π_0))

Zobs = (p - pinula)*sqrt(n)/sqrt(pinula*(1- pinula)) # = 

# 4.) VÝPOČET p-HODNOTY:

phodnota = 1 - pnorm(Zobs,0,1) # = 0.177


# 5.) ZÁVĚR: Nezamítáme nulovou hypotézu, že pravděpodobnost komplkací je 0,3. 


## NEBO POMOCÍ BINOMICKÉHO ROZDĚLENÍ:

#     H_0:   π  = 0,3
#     H_A:   π  > 0,3

# X ...počet komplikací při n pokusech => X -> Bi(n,π)

#X_{obs} = 18

# pvalue = P(X >= 18) = 1 - P(X \leq 17)

pvalue = 1 - pbinom(17,50,0.3) # = 0.2178069

# v R-ku:

binom.test(x = x, n= n, p = 0.3, alternative="greater")








#_____________________________________________________________________________
# Oboustranný test počítá R-ko trochu jinak!
binom.test(x = 18, n= 50, p = 0.3, alternative="two.sided")
#95 percent confidence interval pro π: ( 0.2291571 0.5080686), p-value =  0.3568

qbinom(0.025,50,0.3) # = 9
qbinom(0.975,50,0.3) # = 22

binom.test(x = 8, n= 50, p = 0.3, alternative="two.sided") # p-value = 0.03053
binom.test(x = 9, n= 50, p = 0.3, alternative="two.sided") # p-value = 0.06532
binom.test(x = 21, n= 50, p = 0.3, alternative="two.sided") #p-value = 0.088
binom.test(x = 22, n= 50, p = 0.3, alternative="two.sided") #p-value = 0.04334

1 - (pbinom(21,34,0.15) - pbinom(9,34,0.15)) # = 0.976019
1 - (pbinom(21,34,0.15) - pbinom(8,34,0.15)) # = 0.9413413

1 - (pbinom(17,34,0.15) - pbinom(4,34,0.15)) # = 0.4076221
1 - (pbinom(17,34,0.15) - pbinom(3,34,0.15)) # = 0.04334
#_____________________________________________________________________________






#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------
# 7. Příklad  
# Mezi 34 výrobky bylo 12 vadných. Testujte hypotézu, že pravděpodobnost vady je 0,15.
#-----------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------


binom.test(x = 12, n= 34, p = 0.15, alternative="greater") 

#95 percent confidence interval pro π: (0.2178796 1.0000000), p-value = 0.002791

# p-hodnota = P(X >= Xobs) = 1 - P(X < Xobs) = 1 - P(X =< Xobs - 1) = 1 - P(X =< 12-1)

1 - pbinom(11,34,0.15) # = 0.002791454









#_____________________________________________________________________________
# Oboustranný test počítá R-ko trochu jinak
binom.test(x = 12, n= 34, p = 0.15, alternative="two.sided")
#95 percent confidence interval pro π: (0.1974586, 0.5351136), p-value = 0.002791

qbinom(0.025,34,0.15) # = 1
qbinom(0.975,34,0.15) # = 9

1 - (pbinom(9,34,0.15) - pbinom(0,34,0.15)) # = 0.02796427
#_____________________________________________________________________________





# * Příklad 5. ####
# Firma TT udává, že 1% jejich rezistorů nesplňuje požadovaná kritéria. V testované
# dodávce 1000 ks bylo nalezeno 15 nevyhovujících rezistorů. Potvrzuje tento výsledek
# tvrzení TT? Ověřte čistým testem významnosti.


n = 1000   # rozsah výběru 
x = 15     # počet "úspěchů"
p = x/n    # relativní četnost (bodový odhad pravděpodobnosti)
p 

# Ověření předpokladů
9/(p*(1-p))
# Dále předpokládáme  n/N < 0.05, tj. že daná populace (rezistorů) má rozsah 
# alespoň 1000/0.05 = 1000*20 = 20 000 rezistorů

## Clopperův - Pearsonův (exaktní) test
## H0: pi = 0.01
## Ha: pi <> 0.01

binom.test(x = x, n= n, p = 0.01, alternative="two.sided")

## Clopperův - Pearsonův (exaktní) test
## H0: pi = 0.01
## Ha: pi > 0.01

binom.test(x = x, n= n, p = 0.01, alternative="greater")

# Na hladině významnosti 0,05 nezamítáme H0 
# Nelze očekávat, že podíl vadných rezistorů ve výrobě statisticky významně
# převyšuje 1 %.