# Skript slouzi k:
# 1. Nacteni souboru s daty.
# 2. Uprave hodnot kategorialnich promennych.
# 3. Vypoctu ANOVy pro vicerozmerna data.
# Vytvoreno 20.10.2020; EK


# Nastaveni cesty 
setwd("c:/Users/koritakova/Documents/Vyuka/Vicerozmerky a AKD/Vicerozmerky cviceni/2020/cv03_ANOVA") 

# Nacteni (a pripadne naistalovani) knihoven
#install.packages("readxl")
library("readxl") # nacteni knihovny s funkci read_excel
library("plyr")   # nacteni knihovny s funkci mapvalues
library("dplyr")  # nacteni knihovny pro popisnou sumarizaci
library("car")    # nacteni knihovny s funkci leveneTest

# Nacteni dat a vyber promennych pro analyzu
data <- read_excel("Data_neuro.xlsx",sheet = "data",col_names=T) # nacteni dat
names(data) # nazvy jednotlivych promennych
data <- data[,c("Group_3kat","Gender_rek","Hippocampus_volume (mm3)")] # vyber promennych pro analyzu
names(data)[3] <- "Hippocampus" # uprava nazvu 3. promenne

# Uprava popisku promenne Group_3kat
data$Group_3kat <- mapvalues(data$Group_3kat,from=levels(as.factor(data$Group_3kat)),to=c("CN","MCI","AD")) # vytvoreni nove promenne s textovymi hodnotami kategorii
data$Group_3kat <- factor(data$Group_3kat, levels=c("CN","MCI","AD")) # prerazeni poradi kategorii

# pristup k promennym podle nazvu (tzn. na promenne bude mozno se odkazovat primo nazvem promenne misto data$nazev_promenne)
data <- as.data.frame(data) # vytvoreni data frame
attach(data) 


# Popisna sumarizace
popis <- data %>% 
  group_by(Group_3kat,Gender_rek) %>% 
  summarize(N = n_distinct(Hippocampus),
            mean = mean(Hippocampus),
            sd = sd(Hippocampus),
            median = median(Hippocampus),
            min = min(Hippocampus),
            max = max(Hippocampus))
popis_group_celkem <- data %>% 
  group_by(Group_3kat) %>% 
  summarize(N = n_distinct(Hippocampus),
            mean = mean(Hippocampus),
            sd = sd(Hippocampus),
            median = median(Hippocampus),
            min = min(Hippocampus),
            max = max(Hippocampus))
popis_gender_celkem <- data %>% 
  group_by(Gender_rek) %>% 
  summarize(N = n_distinct(Hippocampus),
            mean = mean(Hippocampus),
            sd = sd(Hippocampus),
            median = median(Hippocampus),
            min = min(Hippocampus),
            max = max(Hippocampus))
popis # popisna sumarizace dle kombinaci skupiny a pohlavi
popis_group_celkem # popisna sumarizace dle skupiny
popis_gender_celkem # popisna sumarizace dle pohlavi


# Krabicovy graf
boxplot(Hippocampus ~(Gender_rek*Group_3kat))


# Overeni normality v kazde kombinaci skupin (+ p-hodnota dle Shapirova-Wilkova testu)
a = length(levels(as.factor(Group_3kat))) # pocet skupin dle jedne kategorialni promenne
b = length(levels(as.factor(Gender_rek))) # pocet skupin dle druhe kategorialni promenne
windows()
par(mfrow=c(a,b))
for (i in 1:a) {
  for (j in 1:b) {
    Hip_podskupina <- Hippocampus[which(Gender_rek==unique(Gender_rek)[j] & Group_3kat==unique(Group_3kat)[i])]
    hist(Hip_podskupina,
         main=paste0(unique(Gender_rek)[j],"-",unique(Group_3kat)[i],", p=",round(shapiro.test(Hip_podskupina)$p.value,3)),
         freq = F, xlab="Hippocampus")
    lines(seq(min(Hip_podskupina), max(Hip_podskupina), by=.1), 
          dnorm(seq(min(Hip_podskupina), max(Hip_podskupina), by=.1),mean(Hip_podskupina), sd(Hip_podskupina)), col="blue")
}}

# Overeni homogenity rozptylu
leveneTest(Hippocampus ~ (Gender_rek*Group_3kat),center=mean)

# ANOVA s interakci
anova_s_interakci <- aov(Hippocampus ~ Group_3kat*Gender_rek)
#anova_s_interakci <- aov(Hippocampus ~ Group_3kat + Gender_rek + Group_3kat:Gender_rek) # jinak zapsano
summary(anova_s_interakci)

# ANOVA bez interakce
anova_bez_interakce <- aov(Hippocampus ~ Group_3kat + Gender_rek)
summary(anova_bez_interakce)
#anova_bez_interakce <- anova(lm(Hippocampus ~ Group_3kat + Gender_rek)) # 2. zpusob
#anova_bez_interakce

# Post-hoc testovani
TukeyHSD(anova_bez_interakce, which = "Group_3kat") # Post-hoc testy pro model bez interakce
TukeyHSD(anova_s_interakci, which = "Group_3kat:Gender_rek") # Post-hoc testy pro model s interakci