library("xlsx")
library("rstatix") # je třeba pro Levenův test
library("car")     # je třeba pro Durbin-Watsonův test

setwd("C:/Jirka/Environment/Výuka/VLBS0621 Biostatistika pro P-PooL/Jaro 2023")

## První příklad
data1.0<-read.xlsx("11_anova_zadani.xlsx","ANOVA",startRow=9,encoding="UTF-8")[,2:5]

colnames(data1.0)<-tolower(make.names(colnames(data1.0))) # trošku upravíme názvy sloupců
colnames(data1.0)[1]<-"datum"

data1.0$den<-weekdays(data1.0$datum)

data1.1<-reshape(data1.0,varying=list(c("čd","rj","leo")),times=c("čd","rj","leo"),direction="long") # je třeba data reformátovat do „nudle“

data1<-data.frame("den"=data1.1$den,
                  "dopravce"=data1.1$time,
                  "počet"=data1.1$čd)

model<-aov(počet~dopravce+den,data=data1) # anova hlavních efektů
model<-aov(počet~dopravce*den,data=data1) # anova s interakcemi

summary(model)

boxplot(počet~dopravce*den,data=data1,las=2,xlab="")

# Otestování předpokladů
hist(residuals(model))                  # normalita a nulový průměr reziduí
levene_test(model)                      # homoskedasticita
durbinWatsonTest(model)                 # nezávislost reziduí na hodnotě

# post-hoc Tukey HSD test - je vidět, že nízké p-hodnoty má jediná skupina
TukeyHSD(model)



## Druhý příklad
data2.0<-read.xlsx("11_anova_zadani.xlsx","Kruskal-Wallis test",startRow=7,encoding="UTF-8")[,2:4]

data2<-data2.0 # data 2 jsou použitelná hned

kruskal.test(mm.Hg~pohlaví,data=data2)
kruskal.test(mm.Hg~věk,data=data2)

data2$pohlaví_a_věk<-paste(data2$pohlaví,data2$věk) # spojení sloupců pohlaví a věk do jednoho (protože Kruskal-Wallisův test neumí víc faktorů)

kruskal.test(mm.Hg~pohlaví_a_věk,data=data2)