#######################################################################################
########################### Vícevýběrové testy ########################################
######################### Martina Litschmannová, Adéla Vrtková ########################
#######################################################################################

################## Přiklad 1. (Doplnění tabulky ANOVA) ################################
#######################################################################################

n=c(40,40,42)          # rozsahy výběrů
prum=c(300,290,310)    # průměry
s=c(33,34,31)          # směrodatné odchylky
n.total=sum(n)         # celkový rozsah výběrů
k=3                    # počet tříd
df.b=k-1               # počet stupňů volnosti - skupinový
df.e=n.total-k         # počet stupňů volnosti - reziduální

# celkový průměr
prum.total=weighted.mean(prum,n)
prum.total

# součet čtverců mezi třídami
ss.b=sum(n*(prum-prum.total)^2)
ss.b

# součet čtverců uvnitř tříd
ss.e=sum((n-1)*s^2)
ss.e

# rozptyl mezi třídami
ms.b=ss.b/df.b
ms.b

# rozptyl uvnitř tříd
ms.e=ss.e/df.e
ms.e

# F-poměr
F=ms.b/ms.e
F

# p-hodnota
p=1-pf(F,df.b,df.e)
p

# Na hladině významnosti 0,05 zamítáme hypotézu o shodě středních hodnot (p-hodnota=0,024, ANOVA), 
# tj. střední hodnoty alespoň jedné dvojice skupin se statisticky významně liší.

# odhady skupinových efektů
efekt=prum-prum.total
efekt

# Oproti celkovému průměru vykazuje nejvíce podprůměrné výsledky skupina 2 
# (o cca 10 jednotek nižší než celkový průměr). Naopak průměr skupiny 3 je 
# o cca 10 jednotek vyšší než celkový průměr. Průměrné výsledky skupiny 1 
# odpovídají celkovému průměru.

################## Přiklad 2. (Kyselina listová) ######################################
#######################################################################################
library(openxlsx)
kysel = readWorkbook("C:/Users/Desktop/vicevyberove_testy.xlsx",
                     sheet=1,
                     colNames=TRUE,
                     startRow = 1)

colnames(kysel)=c("sk1","sk2","sk3")   # přejmenování sloupců

# převod do standardního datového formátu
kysel.s=stack(kysel)
colnames(kysel.s)=c("hodnoty","skupina")
kysel.s=na.omit(kysel.s)

boxplot(kysel)

# Ověření normality
par(mfrow=c(1,3))
for (i in 1:3){
  qqnorm(kysel[,i])
  qqline(kysel[,i])
}

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

tapply(kysel.s$hodnoty,kysel.s$skupina,shapiro.test)

# Na hladině významnosti 0,05 nezamítáme předpoklad normality.

# Ověření shody rozptylů
s=tapply(kysel.s$hodnoty,kysel.s$skupina,sd,na.rm=TRUE)
s

s2=s^2
s2

max(s2)/min(s2)
# -> Dle krabicového grafu a informace o poměru největšího a nejmenšího rozptylů (<2) nepředpokládáme,
# že se rozptyly statisticky významně liší

# předpoklad normality nebyl zamítnut -> Bartlettův test
bartlett.test(kysel) 

# Na hladině významnosti 0,05 nelze zamítnout předpoklad o shodě rozptylů 
# (p-hodnota = 0,647, Bartlettův test).

# Chceme srovnávat stř. hodnoty nezávislých výběrů z normálních rozdělení se stejnými rozptyly -> ANOVA
# příkaz aov() vyžaduje data ve standardním datovém formátu

vysledky=aov(kysel.s$hodnoty~kysel.s$skupina)
summary(vysledky)  
# Na hladině významnosti 0,05 zamítáme hypotézu o shodě středních hodnot (p-hodnota<<0,001, ANOVA) -> mnohonásobné porovnávání
TukeyHSD(vysledky)

# skupinové efekty
prum=tapply(kysel.s$hodnoty,kysel.s$skupina,mean,na.rm=TRUE)
n=tapply(kysel.s$hodnoty,kysel.s$skupina,length)
prum.total=weighted.mean(prum,n)
efekt=prum-prum.total
efekt
# Považujeme-li vysoký obsah kyseliny listové za pozitivní, pak statisticky významně nejlepších
# výsledků dosáhli pacienti ze skupiny 1 (průměrný obsah kys. listové o cca 27 jednotek vyšší než
# průměrný obsah kys. listové v krvi všech testovaných pacientů) a statisticky významně nejhorších výsledků
# dosáhli pacienti ze skupiny 2 (průměrný obsah kys. listové o cca 26 jednotek nižší než
# průměrný obsah kys. listové v krvi všech testovaných pacientů). Obsah kys. listové v krvi pacientů ze skupiny 3
# odpovídá celkovému průměru. Všechny tři skupiny pacientů jsou navzájem dle obsahu kys. listové
# v krvi statisticky významně odlišné.

# Další část skriptu již neobsahuje tak podrobné komentáře =)

################## Přiklad 3. (Králíci) ###############################################
#######################################################################################
kralici = readWorkbook("C:/Users/Desktop/vicevyberove_testy.xlsx",
                       sheet=2,
                       colNames=TRUE,
                       startRow = 1)

colnames(kralici)=c("viden","cesky","kalif")   # přejmenování sloupců

# převod do standardního datového formátu
kralici.s=stack(kralici)
colnames(kralici.s)=c("hodnoty","skupina")
kralici.s=na.omit(kralici.s)

par(mfrow=c(1,1))
boxplot(kralici)

# Odstranění odlehlého pozorování
pom=boxplot(kralici)
pom
kralici.s$hodnoty.bez=kralici.s$hodnoty
kralici.s$hodnoty.bez[kralici.s$hodnoty.bez %in% pom$out]=NA

# Převod dat do tabulky
viden=kralici.s$hodnoty.bez[kralici.s$skupina=="viden"]
cesky=kralici.s$hodnoty.bez[kralici.s$skupina=="cesky"]
kalif=kralici.s$hodnoty.bez[kralici.s$skupina=="kalif"]

n.max=max(length(viden),length(cesky),length(kalif))
length(viden)=n.max
length(cesky)=n.max
length(kalif)=n.max

kralici.bez=as.data.frame(cbind(viden,cesky,kalif))

# Krabicový graf
boxplot(kralici.bez)

# Ověření normality
par(mfrow=c(1,3))
for (i in 1:3){
  qqnorm(kralici.bez[,i])
  qqline(kralici.bez[,i])
}

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

tapply(kralici.s$hodnoty.bez,kralici.s$skupina,shapiro.test)

# Na hladině významnosti 0,05 nezamítáme předpoklad normality.

# Ověření shody rozptylů
s=tapply(kralici.s$hodnoty.bez,kralici.s$skupina,sd,na.rm=TRUE)
s

s2=s^2
s2

max(s2)/min(s2)
# -> Dle krabicového grafu a informace o poměru největšího a nejmenšího rozptylů (blízký 2, avšak rozsah výběrů < 30)
# je těžší odhadnout, zda lze předpokládat shodu rozptylů. Rozhodnout nám pomůže test.

# předpoklad normality nebyl zamítnut -> Bartlettův test
bartlett.test(kralici.bez) 

# Na hladině významnosti 0,05 nelze zamítnout předpoklad o shodě rozptylů 
# (p-hodnota = 0,217, Bartlettův test).

# Chceme srovnávat stř. hodnoty nezávislých výběrů z normálních rozdělení se stejnými rozptyly -> ANOVA
# příkaz aov() vyžaduje data ve standardním datovém formátu

vysledky=aov(kralici.s$hodnoty.bez~kralici.s$skupina)
summary(vysledky)  
# Na hladině významnosti 0,05 zamítáme hypotézu o shodě středních hodnot (p-hodnota<<0,001, ANOVA) -> mnohonásobné porovnávání
TukeyHSD(vysledky)

# skupinové efekty
prum=tapply(kralici.s$hodnoty.bez,kralici.s$skupina,mean,na.rm=TRUE)
n=tapply(kralici.s$hodnoty.bez,kralici.s$skupina,length)
prum.total=weighted.mean(prum,n)
efekt=prum-prum.total
efekt

################## Přiklad 4. (Jakost) ################################################
#######################################################################################
jakost.s = readWorkbook("C:/Users/Desktop/vicevyberove_testy.xlsx",
                        sheet=3,
                        colNames=TRUE,
                        startRow = 1)

colnames(jakost.s)=c("poradi","skupina")   # přejmenování sloupců

par(mfrow=c(1,1))
boxplot(jakost.s$poradi~jakost.s$skupina)

# Ověření normality nemá smysl provádět - z povahy jde o diskrétní data (pořadí)

# Ověření shody rozptylů
s=tapply(jakost.s$poradi,jakost.s$skupina,sd,na.rm=TRUE)
s

s2=s^2
s2

max(s2)/min(s2)
# -> Dle krabicového grafu a informace o poměru největšího a nejmenšího rozptylů (<2)
# ze předpokládat shodu rozptylů. (Kruskalův - Wallisův test má větší sílu testu, jsou-li data homoskedasticitní.)

# jde o "pořadová" data, nemá smysl uvažovat o předpokladu normality -> Leveneův test
install.packages("car")
library(car)
leveneTest(jakost.s$poradi~jakost.s$skupina) 

jakost.s$skupina=as.factor(jakost.s$skupina)
leveneTest(jakost.s$poradi~jakost.s$skupina)
# Na hladině významnosti 0,05 nelze zamítnout předpoklad o shodě rozptylů 
# (p-hodnota = 0,750, Leveneův test).

# Chceme srovnávat mediány nezávislých výběrů (data nelze považovat za výběry z norm. rozdělení) -> Kruskalův-Wallisův test

kruskal.test(jakost.s$poradi,jakost.s$skupina)
# Na hladině významnosti 0,05 nelze zamítnout hypotézu o shodě mediánů (p-hodnota=0,295, Kruskalův-Wallisův test). Tj. statisticky významné rozdíly 
# mezi výrobci (z hlediska pořadí výrobků v soutěži) neexistují.

################## Přiklad 5. (Trombin) ###############################################
#######################################################################################
trombin.s = readWorkbook("C:/Users/Desktop/vicevyberove_testy.xlsx",
                         sheet=4,
                         colNames=TRUE,
                         startRow = 2)

colnames(trombin.s)=c("hodnoty","skupina")   # přejmenování sloupců

par(mfrow=c(1,1))
boxplot(trombin.s$hodnoty~trombin.s$skupina)

# Převod dat do tabulky
a=trombin.s$hodnoty[trombin.s$skupina=="A"]
b=trombin.s$hodnoty[trombin.s$skupina=="B"]
c=trombin.s$hodnoty[trombin.s$skupina=="C"]


n.max=max(length(a),length(b),length(c))
length(a)=n.max
length(b)=n.max
length(c)=n.max

trombin=as.data.frame(cbind(a,b,c))

# Ověření normality
par(mfrow=c(1,3))
for (i in 1:4){
  qqnorm(trombin[,i])
  qqline(trombin[,i])
}

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

tapply(trombin.s$hodnoty,trombin.s$skupina,shapiro.test)

# Na hladině významnosti 0,05 zamítáme předpoklad normality.

# Ověření shody rozptylů
s=tapply(trombin.s$hodnoty,trombin.s$skupina,sd,na.rm=TRUE)
s

s2=s^2
s2

max(s2)/min(s2)
# -> Dle krabicového grafu a informace o poměru největšího a nejmenšího rozptylů (>>2)
# nelze předpokládat shodu rozptylů. Rozhodnout nám pomůže test.

# předpoklad normality byl zamítnut -> Leveneův test
library(car)
leveneTest(trombin.s$hodnoty~trombin.s$skupina) 

trombin.s$skupina=as.factor(trombin.s$skupina)
leveneTest(trombin.s$hodnoty~trombin.s$skupina)
# Na hladině významnosti 0,05 zamítáme předpoklad o shodě rozptylů 
# (p-hodnota << 0,001, Leveneův test).

# Chceme srovnávat mediány nezávislých výběrů, která nemají normální rozdělení a mají různé rozptyly -> Kruskalův - Wallisův test

kruskal.test(trombin.s$hodnoty,trombin.s$skupina)
# Na hladině významnosti 0,05 zamítáme hypotézu o shodě mediánů (p-hodnota<<0,001, Kruskalův-Wallisův test). Tj. trombinový čas je statisticky významně
# ovlivněn preparátem. -> mnohonásobné porovnávání
install.packages("dunn.test")
library(dunn.test)
dunn.test(trombin.s$hodnoty,trombin.s$skupina,method="bonferroni")

# skupinové efekty
med=tapply(trombin.s$hodnoty,trombin.s$skupina,quantile,prob=0.5,na.rm=TRUE)
med.total=quantile(trombin.s$hodnoty,probs=0.5)
efekt=med-med.total
efekt