# ......................................................................................
# ...........................11. CV PAVEL JAHODA........................................
# ..............................UPRAVENÉ MATERIÁLY:.....................................
# ......................... Cvičení 11 - Vybrané dvouvýběrové testy ....................
# ...........................        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

# Aktivace všech potřebných balíčků (v případě potřeby nutno nainstalovat)

## Nastavte si pracovní adresář a uložte si soubor testy_dvouvyberove.xlsx do pracovního adresáře.
#setwd("C:/Users/Martina/OneDrive - VŠB-TU Ostrava/Výuka/Pravděpodobnost a statistika/DATA/aktualni/")

## 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(moments)  # funkce skewness(), kurtosis()
library(rstatix)  # funkce identify_outliers()
library(tidyr)   # funkce pivot_longer()






# Pokud je třeba, nastavíme požadovaný formát vypočtených číselných charakteristik
options(scipen = 100, # 100 (desetinné číslo), -100 (vědecká notace)
        digits = 7,  # max. počet platných cifer
        pillar.sigfig = 8) # max. počet platných cifer ve výstupech z balíčku dplyr





#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
# 1. Příklad
# Zvážili jsme hmotnosti 30 brambor odrůdy A a 20 brambor odrůdy B. Výběrová
# směrodatná odchylka hmotnosti u brambor A je s1 = 0,3 kg a u odrůdy B je to
# s2 = 0,41 kg. Předpokládejme, že hmotnosti brambor odrůdy A i B jsou náhodné
# veličiny s normálním rozdělením pravděpodobnosti.
# 
#Testujte hypotézu, že odrůda B má větší variabilitu hmotnosti.
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------

# Chceme použít test o shodě rozptylů (dvouvýběrový F-test)

# PŘEDPOKLADY: předpokladem testu je, že výběry jsou nezávislé a
#
# X_1 ... hmotnost brambor odrůdy A -> N(μ1, σ1^2)
#
# X_2 ... hmotnost brambor odrůdy B -> N(μ1, σ2^2)
#
# PODLE ZADÁNÍ JSOU PŘEDPOKLADY SPLNĚNY

# H_0: σ1 = σ2
# H_A: σ1 < σ2

#Testová statistika:

# F = S1^2/S2^2 -> F_{n1-1,n2-1} ... Fisherovo snedecorovo rozdělení

# Pozorovaná hodnota:

Fobs = 0.3^2/0.41^2

# p-hodnota = F(Fobs) =
pf(Fobs,29,19) # = 0.06292607

# ZÁVĚR: Na hladině významnosti 0.05 nezamítáme nulovou hypotézu, že rozptyl hmotností
# je u odrůdy A stejný jako u odrůdy B.



#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
# 2. Příklad
# Dva stroje vyrábějí podložky o daném poloměru. Byly změřeny poloměry 20-ti podložek, které
# vyrobil první stroj a 22 podložek, které vyrobil druhý stroj. Testujte hypotézu, 
# že střední hodnoty poloměrů vyrobených prvním a druhým strojem jsou stejné.
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------

polomery = read_excel("testy_dvouvyberove.xlsx",
                  sheet = "podlozky") 

## Explorační analýza za účelem identifikace odlehlých pozorování
boxplot(polomery$x1)$out #zobrazi boxplot a vypíše odlehlá poz. pokud jsou
boxplot(polomery$x2)$out #zobrazi boxplot a vypíše odlehlá poz. pokud jsou

# ŽÁDNÁ ODLEHLÁ POZOROVÁNÍ - MŮŽEME POKRAČOVAT

# CHCEME POUŽÍT DVOUVÝBĚROVÝ t-TEST = > ověříme předpoklady testu:

# Normalita X1:

boxplot(polomery$x1) # OK
qqnorm(polomery$x1)  
qqline(polomery$x1) # OK
šikmost = moments::skewness(polomery$x1,na.rm = T) # =cca -0.229 => OK
špičatost = moments::kurtosis(polomery$x1,na.rm = T)-3 # =cca -0.664 => OK

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

shapiro.test(polomery$x1) # p-value = 0.9375 => NEZAMÍTÁME NORMALITU X1


# Normalita X2:

boxplot(polomery$x2) # OK
qqnorm(polomery$x2)  
qqline(polomery$x2) # OK?
šikmost = moments::skewness(polomery$x2,na.rm = T) # =cca -0.448 => OK
špičatost = moments::kurtosis(polomery$x2,na.rm = T)-3 # =cca -0.181 => OK

hist(polomery$x2) # asi OK?

shapiro.test(polomery$x2) # p-value =  0.3114 => NEZAMÍTÁME NORMALITU X2


# Rovnost rozptylů X1 a X2 (homoskedasticita):
var.test(polomery$x1,polomery$x2,ratio=1,conf.level=0.95) # p-value = 0.7604

# => nezamítáme rovnost rozptylů

# => PROVEDEME t-TEST:

# H_0: μ1 = μ2
# H_A: μ1 <> μ2

t.test(polomery$x1,polomery$x2,mu=0,alternative="two.sided",
       var.equal=TRUE,conf.level=0.95)

# t = 16.511, df = 40, p-value < 0.00000000000000022 (1-pt(16.511,40) vychází 1)

# t = 16.511, df = 40, p-value < 0.00000000000000022
# => ZAMÍTÁME NULOVOU HYPOTÉZU VE PROSPĚCH ALTERNATIVY, STŘEDNÍ HODNOTY SE STATISTICKY 
# VÝZNAMNĚ LIŠÍ


# => PROVEDEME JEDNOSTRANNÝ t-TEST:

# H_0: μ1 = μ2
# H_A: μ1 > μ2

t.test(polomery$x1,polomery$x2,mu=0,alternative="greater",
       var.equal=TRUE,conf.level=0.95)


# => ZAMÍTÁME NULOVOU HYPOTÉZU VE PROSPĚCH ALTERNATIVY, STŘEDNÍ HODNOTA 
# POLOMĚRŮ PODLOŽEK, KTERÉ VYROBIL PRVNÍ STROJ JE STATISTICKY 
# VÝZNAMNĚ VĚTŠÍ NEŽ TĚCH, KTERÉ VYROBIL DRUHÝ STROJ.











#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
# 3. Příklad (Cholesterol - minimalistická verze)
# Data v souboru "testy_dvouvyberove.xlsx", sheet = "cholesterol2" udávají 
# hladinu cholesterolu v krvi mužů dvou různých věkových skupin (20-30 letých a
# 40-50 letých). Ověřte na hladině významnosti 0,05 hypotézu, zda sehladina 
# cholesterolu v krvi starších mužů neliší od hladiny cholesterolu v krvi mladších mužů.
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------

# Test o shodě středních hodnot / mediánů

## Načtení dat z xlsx souboru (pomoci balíčku readxl)

chol = read_excel("testy_dvouvyberove.xlsx",
                  sheet = "cholesterol2",
                  skip = 1)   
colnames(chol)=c("mladsi","starsi")


mladsi = chol[ , c(1)]
starsi = chol[ , c(2)]


#_________________________________________________________
## Explorační analýza za účelem identifikace odlehlých pozorování u MLADŠÍCH
boxplot(mladsi)$out #zobrazi boxplot a vypíše odlehlá poz. pokud jsou

mladsi$ID = 1:length(mladsi$mladsi) # přidáme sloupec s identifikátorem řádku

# POKUD se rozhodneme odstranit odlehlá pozorování =>
# seřadíme podle velikosti a vidíme, že je třeba odstranit 44. a 73. řádek

x1 = mladsi[-c(44,73),] # ... Hodnoty cholesterolu mladších mužů


#_________________________________________________________
## Explorační analýza za účelem identifikace odlehlých pozorování u STARŠÍCH
boxplot(starsi)$out #zobrazi boxplot a vypíše odlehlá poz. pokud jsou

starsi$ID = 1:length(starsi$starsi) # přidáme sloupec s identifikátorem řádku

# POKUD se rozhodneme odstranit odlehlá pozorování =>
# seřadíme podle velikosti a vidíme, že je třeba odstranit 52. řádek

x2 = starsi[-c(52),] # ... Hodnoty cholesterolu starších mužů

x2 = na.omit(x2) # odstranilijsme NA



#_________________________________________________________
## OVĚŘENÍ PŘEDPOKLADŮ

# CHCEME POUŽÍT DVOUVÝBĚROVÝ t-TEST = > ověříme předpoklady testu:

# Normalita X1 (mladsi):

boxplot(x1$mladsi) # OK
qqnorm(x1$mladsi)  
qqline(x1$mladsi) # OK
šikmost = moments::skewness(x1$mladsi,na.rm = T) # =cca  0.01591997 => OK
špičatost = moments::kurtosis(x1$mladsi,na.rm = T)-3 # =cca -0.3434796 => OK

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

shapiro.test(x1$mladsi) # p-value = 0.4639 => NEZAMÍTÁME NORMALITU X1


# Normalita X2:

boxplot(x2$starsi) # OK
qqnorm(x2$starsi)  
qqline(x2$starsi) # OK
šikmost = moments::skewness(x2$starsi,na.rm = T) # =cca 0.057 => OK
špičatost = moments::kurtosis(x2$starsi,na.rm = T)-3 # =cca -0.406 => OK

hist(x2$starsi) # OK

shapiro.test(x2$starsi) # p-value =  0.9395 => NEZAMÍTÁME NORMALITU X2


# Rovnost rozptylů X1 a X2 (homoskedasticita):
var.test(x1$mladsi,x2$starsi,ratio=1,conf.level=0.95) # p-value = cca 0

# => ZAMÍTÁME ROVNOST ROZPTYLŮ

# => NEPROVEDEME t-TEST, ale:


#_________________________________________________________
# Aspinové-Welchův test shody středních hodnot 
# (PŘEDPOKLADY:  NORMALITA X1 i X2 - ověřeno výše) 

# X1 ...cholesterol mladších  -> N(μ1, σ1^2)
# X2 ...cholesterol starších  -> N(μ2, σ2^2)

# H_0: μ1 - μ2 = 0
# H_A: μ1 - μ2 <> 0

t.test(x1$mladsi,x2$starsi,mu=0,alternative="two.sided",
       var.equal=FALSE,conf.level=0.95)

# t = -13.106, p-value < 0.00000000000000022
# 95 percent confidence interval:
# -0.6034840 -0.4447047

0.5*(-0.6034840 + (-0.4447047)) # = -0.5240944


# => Na hladině významnosti 0,05 zamítáme H_0, střední hodnoty se statisticky významně liší.
#==========================================================================================
# Očekáváme, že u mladších mužů je hodnota choelsterolu 
# cca o 0.524 mmol/l menší než u starších.



# jednostranný test

# H_0: μ1 - μ2 = 0
# H_A: μ1 - μ2 < 0

t.test(x1$mladsi,x2$starsi,mu=0,alternative="less",
       var.equal=FALSE,conf.level=0.95)

# t = -13.936, p-value < 0.00000000000000022
# 95 percent confidence interval:
#  -Inf -0.4576698


# => Na hladině významnosti 0,05 zamítáme H_0, střední hodnota je u mladších mužů statisticky
#===========================================================================================

# významně menší, než u starších  (minimálně o 0,457 mmol/l).
#===========================================================



# Dle výsledků výběrového šetření očekáváme, že střední obsah cholesterolu v krvi straších mužů
# bude cca o 0,524 mmol/l vyšší než střední obsah cholesterolu u mladších mužů.
# Dle 95% levostranný intervalového odhadu daného rozdílu očekáváme střední obsah cholesterolu
# u starších mužů minimálně o 0,457 mmol/l větší než stř. hodnota cholesterolu u mladších mužů. 

# Na hladině významnosti 0,05 zamítáme předpoklad o shodě středních hodnot cholesterolu ve skupinách 
# mladších a starších mužů ve prospěch alternativy, že starší muži mají vyšší střední hladinu 
# cholesterolu než muži mladší (t-test, t = 13,1, df = 95, p-hodnota < 0,001).





#----------------------čačí postup :)------------------------------------------
#------------------------------------------------------------------------------
#  Příklad 3. (Cholesterol) ####
#
#
#
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------

# Test o shodě středních hodnot / mediánů

## Načtení dat z xlsx souboru (pomoci balíčku readxl)

chol = read_excel("testy_dvouvyberove.xlsx",
                  sheet = "cholesterol2",
                  skip = 1)   
colnames(chol)=c("mladsi","starsi")

## Převod do standardního datového formátu
chol.s = pivot_longer(chol,cols = c("mladsi","starsi"),
                      names_to = "skupina",
                      values_to = "hodnoty")
chol.s = na.omit(chol.s)   # Funkci na.omit() používejte velmi ooutřetně!


## Explorační analýza
boxplot(chol.s$hodnoty~chol.s$skupina)

# Data obsahují odlehlá pozorování.
chol.s$ID = 1:length(chol.s$skupina)

outliers = chol.s %>% 
  group_by(skupina) %>% 
  identify_outliers(hodnoty)

chol.s$hodnoty.out = ifelse(chol.s$ID %in% outliers$ID,NA,chol.s$hodnoty)

boxplot(chol.s$hodnoty.out~chol.s$skupina)

# Samostatné proměnné
mladsi = chol.s$hodnoty %>% filter(skupina == "mladsi") %>% droplevels()
starsi = chol.s$hodnoty %>% filter(skupina == "starsi") %>% droplevels()


## Explorační analýza pro data po odstranění odlehlých pozorování
boxplot(chol.s$hodnoty.out~chol.s$skupina)

chol.s %>% 
  group_by(skupina) %>% 
  summarise(rozsah = length(na.omit(hodnoty.out)),
            sm.odch = sd(hodnoty.out,na.rm = T))

# sd budeme zaokrouhlovat na 3 platné cifry, tj. sd i míry polohy budeme zaokrouhlovat na tisíciny

# Ověření normality
ggplot(chol.s,aes(x = hodnoty.out)) +
  geom_histogram() +
  facet_wrap(~skupina,nrow = 2)

ggplot(chol.s,aes(sample = hodnoty.out)) +
  stat_qq() +
  stat_qq_line() +
  facet_wrap(~skupina,scales = "free")

chol.s %>% 
  group_by(skupina) %>% 
    summarise(
    šikmost = moments::skewness(hodnoty.out,na.rm = T),
    špičatost = moments::kurtosis(hodnoty.out,na.rm = T)-3)

chol.s %>% 
  group_by(skupina) %>% 
  shapiro.test(hodnoty.out)$p.value
  
# Šikmost i špičatost odpovídá norm. rozdělení. Pro konečné rozhodnutí o normalitě dat použijeme
# test normality.

# Na hl. významnosti 0.05 nelze předpoklad normality zamítnout 
# (p-hodnota=0.464 (mladší), p-hodnota=0.940 (starší), Shapirův-Wilkův test).

# Ověření shody rozptylů

# Chceme-li ověření shody rozptylů řešit raději pomocí intervalového odhadu, 
# je pro snazší interpretaci lepší zajistit, aby byl nalezen odhad pro poměr
# max(rozptyl 1, rozptyl 2)/min(rozptyl 1, rozptyl 2), tj. větší rozptyl v čitateli.

boxplot(chol.s$hodnoty.out~chol.s$skupina) 

# V našem případě volíme poměr rozptyl.starsi/rozptyl.mladsi.

# H0: sigma.starsi = sigma.mladsi
# Ha: sigma.starsi <> sigma.mladsi
 
var.test(starsi$hodnoty.out,mladsi$hodnoty.out,ratio=1,conf.level=0.95)
# Na hl. významnosti 0.05 zamítáme předpoklad o shodě rozptylů (p-hodnota<<0.001,test o shodě rozptylů).
# Pozorovanou neshodu mezi rozptyly lze na hladině významnosti 0,05 označit 
# za statisticky významnou (F test, F = 12,2, df_num = 83, df_denom = 96, p-hodnota < 0,001).

# Doplnění:
# Rozptyl hodnot cholesterolu ve skupině starších mužů je cca 12x větší než 
# rozptyl hodnot cholesterolu ve skupině mladších mužů. 95% intervalový odhad  
# tohota poměru je (8;19). 

# Ověření shody středních hodnot (Aspinové-Welchův test)

# Pro snadnou interpretaci intervalového odhadu je vhodné pořadí výběru stanovit tak,
# aby rozdíl jejich průměrů byl kladný.

# H0: mu.starsi = mu.mladsi (mu.starsi - mu.mladsi = 0)
# Ha: mu.starsi > mu.mladsi (mu.starsi - mu.mladsi > 0)
t.test(starsi$hodnoty.out,mladsi$hodnoty.out,mu=0,alternative="greater",
       var.equal=FALSE,conf.level=0.95)

# Dle výsledků výběrového šetření očekáváme, že střední obsah cholesterolu v krvi straších mužů
# bude cca o 0,524 mmol/l vyšší než střední obsah cholesterolu u mladších mužů.
# Dle 95% levostranný intervalového odhadu daného rozdílu očekáváme střední obsah cholesterolu
# u starších mužů minimálně o 0,457 mmol/l větší než stř. hodnota cholesterolu u mladších mužů. 

# Na hladině významnosti 0,05 zamítáme předpoklad o shodě středních hodnot cholesterolu ve skupinách 
# mladších a starších mužů ve prospěch alternativy, že starší muži mají vyšší střední hladinu 
# cholesterolu než muži mladší (t-test, t = 13,1, df = 95, p-hodnota < 0,001).












#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
#  Příklad 4. (Deprese)
# Údaje v souboru ???? představují délku remise (období bez známek nemoci) ve dnech z prostého 
# náhodného výběru ze dvou různých skupin pacientů - pacienti s endogenní depresi (fyziologické příčiny)
# a pacienti s neurotickou depresí (psychologické příčiny).
# Ověřte, zda je pozorovaný rozdíl mezi průměrnou délkou remise u těchto dvou skupin pacientů
# statisticky významný.
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------

# Test o shodě středních hodnot / mediánů

## Načtení dat z xlsx souboru (pomoci balíčku readxl)
deprese = read_excel("testy_dvouvyberove.xlsx",
                     sheet = "deprese")   
colnames(deprese)=c("endo","neuro")

## Převod do standardního datového formátu
deprese.s = pivot_longer(deprese,cols = c("endo","neuro"),
                         names_to = "skupina",
                         values_to = "hodnoty")
deprese.s = na.omit(deprese.s)

## Explorační analýza
boxplot(deprese.s$hodnoty~deprese.s$skupina)

# Data neobsahují odlehlá pozorování.

# Samostatné proměnné
endo = deprese.s %>% filter(skupina == "endo") %>% droplevels()
neuro = deprese.s %>% filter(skupina == "neuro") %>% droplevels()

deprese.s %>% group_by(skupina) %>% 
  summarise(rozsah = length(na.omit(hodnoty)),
            sm.odch = sd(hodnoty,na.rm = T))
#___________________________________________________________________
#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
# sd budeme zaokrouhlovat na 3 platné cifry, sd i míry polohy budeme zaokrouhlovat na jednotky

# Ověření normality
ggplot(deprese.s,aes(x = hodnoty)) +
  geom_histogram() +
  facet_wrap(~skupina,nrow = 2)

# HISTOGRAM => NEVYPADÁ NA NORMÁLNÍ ROZDĚLENÍ!

ggplot(deprese.s,aes(sample = hodnoty)) +
  stat_qq() +
  stat_qq_line() +
  facet_wrap(~skupina,scales = "free")

# QQ GRAF => NEVYPADÁ NA NORMÁLNÍ ROZDĚLENÍ!

deprese.s %>% 
  group_by(skupina) %>% 
  summarise(
    šikmost = moments::skewness(hodnoty,na.rm = T),
    špičatost = moments::kurtosis(hodnoty,na.rm = T)-3)

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

shapiro.test(endo$hodnoty)
shapiro.test(neuro$hodnoty)
# (Shapirův-Wilkův test, p-hodnota_endo<<0,001, p-hodnota_neuro<<0,001).

# Na hl. významnosti 0,05 ZAMÍTÁME NORMALITU DAT! 




# => NEPROVEDEME ANI DVOUVÝBĚROVÝ t-TEST, ANI ASPINOVÉ-WELSCHŮV TEST, ALE:


#_______________________________________________________
# Ověření shody mediánů pomocí Mannova - Whitneyho testu:

# POZOR! V nulové hypotéze testujeme, zda se rozdělení populací, z nichž
# byly výběry pořízeny liší pouze svou polohou (konkrétně mediánem).

# Dle histogramů předpokládáme, že data mají stejný typ rozdělení.


# H0: medián.neuro = medián.endo (med.neuro - med.endo = 0)
# Ha: medián.neuro > medián.endo (med.neuro - med.endo > 0)



#  Mann - Whitneyho test je rozšířením Wilcoxonova testu, proto:
wilcox.test(neuro$hodnoty,endo$hodnoty,
            mu=0,
            alternative="greater",
            conf.level=0.95,
            conf.int = T) # POZOR! Přednastavená hodnota conf.int je FALSE.

# VÝSLEDEK TESTU:
#  W = 16366, p-value < 0.00000000000000022
# alternative hypothesis: true location shift is greater than 0
# 95 percent confidence interval:
#  168 Inf
# sample estimates:
#   difference in location 
# 191.0001 


# => ZÁVĚR:
# Doba remise pacientů s neurotickou depresí je cca o 191 dnů delší než u pacientů 
# s endogenní depresí. Dle 95% levostranného intervalového odhadu očekáváme,
# že pacienti s neurotickou depresí mají minimálně o 168 dní delší dobu remise
# než pacienti s neurotickou depresí.

# Na hladině významnosti 0,05 zamítáme hypotézu o shodě 
# mediánů dob do remise onemocnění pro obě skupiny pacientů ve prospěch alternativy 
# (Wilcoxonův test s korekcí na spojitost, W = 16 366, p-hodnota << 0,001).
# Medián doby remise je u pacientů s neurotickou depresí statisticky významně
# větší než u pacientů s endogenní depresí.





#--------------------------Minimalistická verze---------------------------------
#------------------------------------------------------------------------------
# 5. Příklad (Osmolalita - párová data) 
# Sledujeme osmolalitu (koncentraci aktivních částic) moči na lůžkové stanici 
# v 08:00 hodin a v 11:00 hodin u 16 mužů. Na základě výsledků uvedených v 
# souboru ? ověřte, zda se osmolalita statisticky významně
# zvýšila.
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------

#____________________________________________
# Načtení dat 
osmolalita = read_excel("testy_dvouvyberove.xlsx",
                        sheet = "osmolalita", skip = 1) 

colnames(osmolalita)=c("ID","o8","o11")

#____________________________________________
# Výpočet nárůstu osmolality
osmolalita$rozdil = osmolalita$o8 -  osmolalita$o11


#____________________________________________
# IDENTIFIKACE ODLEHLÝCH POZOROVÁNÍ 
boxplot(osmolalita$rozdil)$out

# Data obsahují odlehlá pozorování ID = 3 a ID = 10.

# Odstranění odlehlých hodnot
osmolalita = osmolalita[-c(3,10),]


#____________________________________________
# OVĚŘENÍ NORMALITY ROZDÍLU bez odlehlých pozorování
boxplot(osmolalita$rozdil) # asi OK
hist(osmolalita$rozdil)    # asi OK
qqnorm(osmolalita$rozdil)  
qqline(osmolalita$rozdil) # ?OK
šikmost = moments::skewness(osmolalita$rozdil,na.rm = T) # =cca -0.436 => OK
špičatost = moments::kurtosis(osmolalita$rozdil,na.rm = T)-3 # =cca -0.467 => OK

shapiro.test(osmolalita$rozdil) # p-value =  0.5452 => NEZAMÍTÁME NORMALITU ROZDÍLŮ


#____________________________________________
# => provedeme oboustranný Párový t-test
# H0: μ = 0 mm.....střední hodnota rozdílu: "v 8" - "v 11"
# Ha: μ <> 0 mm

t.test(osmolalita$rozdil,mu=0,alternative="two.sided")

# t = -3.0514, df = 13, p-value = 0.009276 < 0.05
# 95 percent confidence interval:
# -41.480046  -7.091383
# st5ed tohoto intervalu je:
0.5*(-41.480046  -7.091383) # = -24.28571

# Na hladině významnosti 0,05 zamítáme nulovou hypotézu ve prospěch alternativy,
# že střední hodnota osmolality je v 8:00 menší (asi o 24mmol/l) než v 11:00.


#____________________________________________
# => provedeme jednostranný Párový t-test
# H0: μ = 0 mm.....střední hodnota rozdílu: "v 8" - "v 11"
# Ha: μ < 0 mm

t.test(osmolalita$rozdil,mu=0,alternative="less")

# t = -3.0514, df = 13, p-value = 0.004638
# alternative hypothesis: true mean is less than 0
# 95 percent confidence interval:
#  -Inf -10.19089
# sample estimates:
#  mean of x 
# -24.28571 


# => 
# Dle výběrového šetření lze očekávat, že osmolalita moči se
# mezi 8 a 11 hodinou zvýší o cca 24 mmol/kg. Dle 95% intervalového odhadu
# lze očekávat, že dojde k navýšení osmolality minimálně o 10 mmol/kg).
# Na hladině významnosti 0,05 lze tento nárůst označit za statisticky
# významný (párový t-test, t = 3,1, df = 13, p-hodnota = 0,005).





#--------------------------------čačí verze-------------------------------------
#------------------------------------------------------------------------------
# 5. Příklad (Osmolalita - párová data) 
# Sledujeme osmolalitu (koncentraci aktivních částic) moči na lůžkové stanici 
# v 08:00 hodin a v 11:00 hodin u 16 mužů. Na základě výsledků uvedených v 
# souboru ? ověřte, zda se osmolalita statisticky významně
# zvýšila.
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------

# Načtení dat 
osmolalita = read_excel("testy_dvouvyberove.xlsx",
                        sheet = "osmolalita", skip = 1)  
osmolalita = osmolalita[,c(2,3)]
colnames(osmolalita)=c("o8","o11")

# Výpočet nárůstu osmolality
osmolalita$narust = osmolalita$o11 -  osmolalita$o8

## Explorační analýza 
boxplot(osmolalita$narust)

# Data obsahují odlehlá pozorování.

# Odstranění odlehlých hodnot
osmolalita$ID = 1:length(osmolalita$narust)

outliers = osmolalita %>% 
  identify_outliers(narust)

osmolalita$narust.out = ifelse(osmolalita$ID %in% outliers$ID,NA,osmolalita$narust)

## Explorační analýza pro data bez odlehlých pozorování
boxplot(osmolalita$narust.out)
abline(h=0,lty=2) # přímka definující nulový efekt

# vizualizace vhodná pro publikaci
ggplot(osmolalita,
       aes(x = "",
           y = narust.out))+ # estetika
  stat_boxplot(geom = "errorbar", 
               width = 0.2) +   # ohraničení fousů (musí být před geom_boxplot)
  geom_boxplot()+
  labs(x = "",
       y = "osmolalita (Osm/kg)")+  # popisky
  theme_bw()+
  geom_point(aes(x = "",
                 y = mean(osmolalita$narust.out, na.rm=T)),  # vykreslení průměru
             color = "red",
             shape = 3) +
  geom_hline(yintercept=0, linetype="dashed", 
             color = "red", size=1)


osmolalita %>% 
  summarise(rozsah = length(na.omit(narust)),
            sm.odch = sd(narust,na.rm = T),
            prumer = mean(narust,na.rm = T),
            median = median(narust,na.rm = T))

# sd zaokrouhlujeme na 2 platné cifry, tj. sd a míry polohy zaokrouhlujeme na jednotky

# Ověření normality
ggplot(osmolalita,aes(x = narust.out)) +
  geom_histogram() 

ggplot(osmolalita,aes(sample = narust.out)) +
  stat_qq() +
  stat_qq_line() 

osmolalita %>% 
  summarise(
    šikmost = moments::skewness(narust.out,na.rm = T),
    špičatost = moments::kurtosis(narust.out,na.rm = T)-3)

osmolalita %>% 
  shapiro.test(narust.out)$p.value

# Šikmost i špičatost odpovídají normálnímu rozdělení.

# Na hl. významnosti 0.05 nelze předpoklad normality zamítnout 
# (Shapirův-Wilkův test, W = 0,949, p-hodnota=0,545).

# Párový t-test
# H0: mu.narust = 0 mm
# Ha: mu.narust > 0 mm

t.test(osmolalita$narust.out,mu=0,alternative="greater")

# Dle výběrového šetření lze očekávat, že osmolalita moči se
# mezi 8 a 11 hodinou zvýší o cca 24 mmol/kg. Dle 95% intervalového odhadu
# lze očekávat, že dojde k navýšení osmolality minimálně o 10 mmol/kg).
# Na hladině významnosti 0,05 lze tento nárůst označit za statisticky
# významný (párový t-test, t = 3,1, df = 13, p-hodnota = 0,005).




#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
# 5. Příklad (rovnost parametrů binomického rozdělení) 
# Při 150-ti laparoskopických operacích se komplikace vyskytly v 18-ti případech
# a při 146-ti otevřených operacích se vyskytly ve 21 případech. Testujte hypotézu, 
# že pravděpodobnost komplikací u laparoskopických operací je stejná jako u otevřených.
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------

# TEST O SHODĚ PARAMETRŮ DVOU BINOMICKÝCH ROZDĚLENÍ

n1 = 150
x1 = 18
p1 = 18/150 # = 0.12

n2 = 146
x2 = 21
p2 = 21/146 # = 0.1438...


# PŘEDPOKLADY TESTU: V konstrukci testu je použita CLV => 

9/(p1*(1-p1)) # =  85.22727 < n = 150 OK
9/(p2*(1-p2)) # =  73.08343 < n = 146 OK


# 1. NULOVÁ HYPOTÉZA
#     H_0:   π1 - π2 = 0
#     H_A:   π1 - π2 <> 0,3


# 2. TESTOVÁ STATISTIKA

# Zobs = (p1 - p2 - 0)/sqrt(p1*(1-p1)/n1 + p2*(1-p2)/n2)


# 3. POZOROVANÁ HODNOTA TESTOVÉ STATISTIKY

Zobs = (p1 - p2 - 0)/sqrt(p1*(1-p1)/n1 + p2*(1-p2)/n2) # = -0.60591


# URČÍME p-hodnotu 2*min{F(Zobs); 1 - F(Zobs)}

phodnota = 2*pnorm(Zobs,0,1) # = 0.5445681...


# ZÁVĚR: Protože p-hodnota = 0.545 > 0.05, nezamítáme na hladině významnosti 
# nulovou hypotézu. Rozdíl mezi pravděpodobnostmi pooperačních komplikací 
# po otevřených a po laparoskopických operacích není statisticky významný.

#NEBO
## Pearsonův X2 test #######################################################
## H0: pi.PP = pi.MM
## Ha: pi.PP > pi.MM

x.MM = 18
n.MM = 150
p.MM = x.MM/n.MM

x.PP = 21
n.PP = 146
p.PP = x.PP/n.PP

prop.test(c(x.PP,x.MM),c(n.PP,n.MM),alternative="two.sided",conf.level=0.95)
prop.test(c(x.PP,x.MM),c(n.PP,n.MM),alternative="greater",conf.level=0.95)






# * Příklad 4. (Porovnání kvality) ####
x.MM = 14
n.MM = 200
p.MM = x.MM/n.MM

x.PP = 15
n.PP = 150
p.PP = x.PP/n.PP


# Ověření předpokladů
n.MM > 9/(p.MM*(1-p.MM))
n.PP > 9/(p.PP*(1-p.PP))
# Dále pro obě firmy předpokládáme, že  n/N < 0.05, tj. že daná populace (součástek) má rozsah 
# alespoň 20*n, tj. 20*200 (4 000), resp. 20*150 (3 000) součástek, což je asi vcelku reálný předpoklad.

## Waldův test #############################################################
# "Pořadí výběrů" volíme kvůli snazší interpretaci tak, aby rozdíl "p.1 - p.2" vycházel kladný.
## H0: pi.PP = pi.MM
## Ha: pi.PP > pi.MM

x.obs = (p.PP - p.MM)/sqrt((p.PP*(1-p.PP)/n.PP)+(p.MM*(1-p.MM)/n.MM))
x.obs
p.value = 1-pnorm(x.obs,0,1)
p.value

# Pravostranný intervalový odhad
alpha = 0.05
p = (x.PP + x.MM)/(n.PP + n.MM)
s.error = qnorm(1-alpha/2)*sqrt(p*(1-p)*(1/n.MM + 1/n.PP))
dolni.mez = (p.PP - p.MM) - s.error
dolni.mez

# Dle výběrového šetření očekáváme, že pravděpodobnost výroby vadné součástky
# je u firmy PP cca o 3 % větší než u firmy MM. Na základě 95% Waldova pravostranného intervalového odhadu
# (-0,028; 1,000) lze pozorovaný rozdíl v kvalitě výroby označit za statisticky nevýznamný. Ke stejným závěrům
# můžeme dojít i na základě Waldova pravostranného testu (x_obs = 0,986,
# p-hodnota = 0,162).


## Pearsonův X2 test #######################################################
## H0: pi.PP = pi.MM
## Ha: pi.PP > pi.MM

prop.test(c(x.PP,x.MM),c(n.PP,n.MM),alternative="greater",conf.level=0.95)

# Dle výběrového šetření očekáváme, že pravděpodobnost výroby vadné součástky
# je u firmy PP cca o 3 % větší než u firmy MM. Na základě 95% Clopperova - Pearsonova pravostranného intervalového odhadu
# (-0,025; 1,000) lze pozorovaný rozdíl v kvalitě výroby označit za statisticky nevýznamný. Ke stejným závěrům
# můžeme dojít i na základě Pearsonova pravostranného testu (x_obs = 0,659,
# df = 1, p-hodnota = 0,209).