##########################################################################################
################################## SOFTWARE R - CAST 3 ###################################
##########################################################################################

# Nastaveni pracovniho adresare - napr.:
setwd("f:/Prace/Vyuka/FSTA_2015")


############### Nacteni dat ###############

rm(list=ls())
kardio <- read.table(file="Data_pro_R.txt", sep="\t", dec=".", header=T) # nacteni dat
dim(kardio)                       # rozmery datove tabulky
head(kardio)                      # vypsani prvnich 10 radku tabulky


###########################################################################################
################################### PARAMETRICKE TESTY ####################################
###########################################################################################

############### Test o prumeru pri neznamem rozptylu - jednovyberovy t-test ###############

# Chceme srovnat na hladine vyznamnosti alpha=0.05 prumerny systolicky tlak v nasem souboru 
# s prumernym systolickym tlakem 131 mm Hg zjistenym pri populacnim epidemiologickem pruzkumu.
#
# Tzn. H0: prumer = 131; H1: prumer != 131

alpha <- 0.05                     # hladina spolehlivosti
prumer.populace <- 131            # prumerny systolicky tlak v cele populaci

?t.test
x <- kardio$sys_tlak
test <- t.test(x,mu=prumer.populace)
test

# overeni s manualnim vypoctem:
prumer <- mean(x)                 # vyberovy prumer
s <- sd(x)                        # vyberova smerodatna odchylka
n <- length(x)                    # velikost vzorku
kvantil <-  qt(1-alpha/2,n-1)     # kvantil Studentova rozdeleni
IS.dolni <- prumer - (s/sqrt(n))*kvantil
IS.horni <- prumer + (s/sqrt(n))*kvantil
IS.dolni
prumer
IS.horni
test.stat <- (prumer - prumer.populace)/(s/sqrt(n)) # realizace testove statistiky
test.stat
p.hodnota <- 2 * (1 - pt(test.stat,n-1))
p.hodnota



###########################################################################################

############### Parovy t-test ###############

# Chceme overit na hladine vyznamnosti alpha=0.05, ze lecba vedla ke zmene systolickeho tlaku.
# 
# Tzn. H0: (prumer_po_lecbe - prumer_pred_lecbou) = 0; 
#      H1: (prumer_po_lecbe - prumer_pred_lecbou) != 0

rozdil.tlaku <- kardio$sys_tlak_po_lecbe - kardio$sys_tlak
t.test(rozdil.tlaku,mu=0)         # vypocet paroveho t-testu
t.test(kardio$sys_tlak_po_lecbe,kardio$sys_tlak,paired=TRUE) # jiny zpusob vypoctu paroveho t-testu

# Pokud bychom chteli overit, ze lecba vedla ke snizeni systolickeho tlaku, byly by hypotezy:
# H0: (prumer_po_lecbe - prumer_pred_lecbou) = 0;
# H1: (prumer_po_lecbe - prumer_pred_lecbou) < 0.
t.test(rozdil.tlaku,alternative="less",mu=0) # vsimnete si, jak vypada interval spolehlivosti;
                                             # p-hodnota je bohuzel tak mala, ze na ni neni videt rozdil oproti oboustranne alternativni hypoteze


############### Dvouvyberovy t-test ############### 

## Priklad 1.
# Chceme overit na hladine vyznamnosti alpha=0.05, zda se systolicky tlak lisi podle pohlavi.
#
# Tzn. H0: (prumer_zeny - prumer_muzi) = 0; 
#      H1: (prumer_zeny - prumer_muzi) != 0
# Lze take testovat, ze se prumery lisi o nejakou konstantu c, pak by H0 vypadala takto: (prumer_zeny - prumer_muzi) = c;

t.test(kardio$sys_tlak[kardio$pohlavi_rek=="M"],kardio$sys_tlak[kardio$pohlavi_rek=="Z"])
# V R je automaticky u dvouvyberoveho t-testu provadena Welchova korekce pro nestejne rozptyly.
# Pocet stupnu volnosti tedy neni roven (n1+n2-2), protoze pri Welchove korekci se pouziva slozitejsi vzorec.
# Interval spolehlivosti neobsahuje zvolenou konstantu c (v tomto pripade nulu), proto zamitame H0. 

# jine zpusoby vypoctu:
with(kardio,t.test(sys_tlak[pohlavi_rek=="M"],sys_tlak[pohlavi_rek=="Z"]))
t.test(kardio$sys_tlak~kardio$pohlavi_rek)
t.test(sys_tlak~pohlavi_rek,data=kardio)

boxplot(sys_tlak~pohlavi_rek,data=kardio)

# Ukol 1: Overte na hladine vyznamnosti alpha=0.05, zda se lisi vaha u muzu a zen.
# Napoveda: Overte graficky (pomoci histogramu), zda se vaha u muzu a zen ridi normalnim rozdelenim. Pokud ne, vahu
# zlogaritmujte.
hist(kardio$vaha)
hist(log(kardio$vaha)) # neni to idealni, ale pro demonstracni ucely nam to bude stacit
t.test(log(vaha)~pohlavi_rek,data=kardio)


## Priklad 2. 
# Chceme overit na hladine vyznamnosti testu a) alpha=0.05, b) alpha=0.01, zda se systolicky tlak lisi u muzu a zen o 21 mm Hg.
#
# Tzn. H0: (prumer_zeny - prumer_muzi) = 21; 
#      H1: (prumer_zeny - prumer_muzi) != 21

tlak_zeny <- kardio$sys_tlak[kardio$pohlavi_rek=="Z"]
tlak_muzi <- kardio$sys_tlak[kardio$pohlavi_rek=="M"]

# a) alpha = 0.05:
t.test(tlak_muzi,tlak_zeny,mu=21) # H0 zamitame (tzn. statisticky vyznamny vysledek), protoze p-value=0.0324<0.5 a CI neobsahuje 21

# b) alpha = 0.01:
t.test(tlak_muzi,tlak_zeny,mu=21,conf.level=0.99) # H0 nezamitame, protoze p-value=0.0324>0.01 a CI obsahuje 21


#### Dvouvyberovy t-test pro jednostrannou alternativu:

# Chceme overit na hladine vyznamnosti alpha=0.05, zda je systolicky tlak u zen nizsi nez u muzu.
#
# Tzn. H0: (prumer_zeny - prumer_muzi) = 0; 
#      H1: (prumer_zeny - prumer_muzi) < 0.

t.test(kardio$sys_tlak[kardio$pohlavi_rek=="Z"],kardio$sys_tlak[kardio$pohlavi_rek=="M"],alternative="less")
t.test(sys_tlak~pohlavi_rek,data=kardio,alternative="less") # tohle neni spravne, protoze kategorie "M" je v abecede driv nez "Z", takze se pocita, zda je systolicky tlak u muzu nizsi nez u zen




###########################################################################################
################################## NEPARAMETRICKE TESTY ###################################
###########################################################################################

############### Wilcoxonuv test - neparametricky test pro jeden vyber ###############

## Ukol 4: Chceme overit na hladine vyznamnosti alpha=0.05, zda lide v nasem souboru trpi nadvahou (tedy
# zda maji BMI vetsi nez 25).
# a) Stanovte nulovou a alternativni hypotezu.
# b) Overte graficky normalitu BMI.
# c) Vzhledem k tomu, ze BMI se neridi normalnim rozlozenim (nelze pouzit t.test) a rozdeleni BMI neni ani symetricke 
#    (nelze pouzit Wilcoxonuv test), zlogaritmujte BMI (pouzijte dekadicky logaritmus) a overte graficky normalitu zlogaritmovaneho BMI. 
# d) Spocitejte median BMI, najdete funkci pro jednovyberovy Wilcoxonuv test a overte pomoci Wilcoxonova testu na hladine vyznamnosti alpha=0.05, zda lide trpi nadvahou.
# e) Spocitejte i geometricky prumer BMI a overte pomoci t-testu na hladine vyznamnosti alpha=0.05, zda lide trpi nadvahou.

x <- kardio$BMI[!is.na(kardio$BMI)] # vybereme validni hodnoty (tzn. odstraneni hodnot NA)

## Reseni:
# a) Nulova a alternativni hypoteza:
# H0: median < 25; H1: median > 25.

# b) Graficke overeni normality BMI:
hist(x)

# c) Zlogaritmovani BMI.
x <- log10(x)
hist(x) # neni to idealni, ale pro demonstracni ucely nam to bude stacit.

# d) Vypocteni medianu, nalezeni funkce pro jednovyberovy Wilcoxonuv test a overeni pomoci Wilcoxonova testu na hladine vyznamnosti alpha=0.05, zda lide trpi nadvahou.
median.log <- median(x)
median <- 10^median.log
median
median(kardio$BMI[!is.na(kardio$BMI)]) # puvodni hodnota medianu - neni divu, ze je to stejna hodnota, kdyz pri logaritmicke transformaci se nemeni poradi hodnot!
??wilcoxon
?wilcox.test
wilcox.test(x,alternative="greater",mu=log10(25),conf.level=0.95)

# e) Vypocteni geometrickeho prumeru a overeni pomoci jednovyberoveho t-testu na hladine vyznamnosti alpha=0.05, zda lide trpi nadvahou.
prumer.log <- mean(x)
geom.prumer <- 10^prumer.log
geom.prumer
t.test(x,alternative="greater",mu=log10(25),conf.level=0.95)

# Zaver: H0 nezamitame. Tzn. soubor lidi netrpi nadvahou (tzn. vetsina lidi nema nadvahu (nekolik lidi ale nadvahu ma)).



############### Mann-Whitneyho test - neparametricky test pro dva vybery ###############

## Ukol 5: Chceme overit na hladine vyznamnosti alpha=0.05, zda se lisi BMI podle pohlavi.
# a) Stanovte nulovou a alternativni hypotezu.
# b) Overte graficky normalitu BMI a stejneho rozdeleni dat u muzu a u zen.
# c) Vzhledem k tomu, ze BMI se neridi normalnim rozlozenim, pouzijte Mann-Whitneyho test (ten je schovan v R pod funkci 
#    wilcox.test, protoze Mann-Whitneyho test lze nazvat tez dvouvyberovy Wilcoxonuv test).

x <- kardio$BMI
pohlavi <- kardio$pohlavi_rek

## Reseni:
# a) Nulova a alternativni hypoteza:
# H0: F(x)=F(y); H1: F(x)!=F(y).

# b) Graficke overeni normality BMI u muzu a u zen:
hist(x)
hist(x[pohlavi=="M"])
hist(x[pohlavi=="Z"])
boxplot(x~pohlavi)

# c) Overeni na hladine vyznamnosti alpha=0.05, zda se lisi BMI podle pohlavi.
wilcox.test(x~pohlavi,conf.level=0.95)

# Zaver: H0 zamitame. Tzn. BMI se lisi podle pohlavi.



############### ANOVA ################

### Rekategorizace veku:
kardio$vek.3kat <- 1+(kardio$vek>=30)+(kardio$vek>=50) # vytvoreni trech kategorii: 1...<30, 2...30-50, 3...>=50 
table(kardio$vek.3kat,useNA="ifany") # frekvencni tabulka 


### 1. Overeni normality:

# 1.1 histogram:
windows()
par(mfrow=c(1,3))                 # nastaveni vykresleni 3 grafu do 1 okna
hist(kardio$sys_tlak[kardio$vek.3kat==1])
hist(kardio$sys_tlak[kardio$vek.3kat==2])
hist(kardio$sys_tlak[kardio$vek.3kat==3])
# Normalita neni splnena! ANOVA by se nemela pocitat - tyhle skripty jsou jen na ukazku, jak se funkce volaji!!!

# 1.2 boxplot:
par(mfrow=c(1,1))
boxplot(kardio$sys_tlak~kardio$vek.3kat)

# 1.3 Q-Q graf:
par(mfrow=c(1,3))                 # nastaveni vykresleni 3 grafu do 1 okna
qqnorm(kardio$sys_tlak[kardio$vek.3kat==1])
qqnorm(kardio$sys_tlak[kardio$vek.3kat==2])
qqnorm(kardio$sys_tlak[kardio$vek.3kat==3])

# 1.4 Shapiruv-Wilkuv test - doporucovan pro mensi vzorky, muze byt moc prisny pro velke vzorky:
shapiro.test(kardio$sys_tlak[kardio$vek.3kat==1])
shapiro.test(kardio$sys_tlak[kardio$vek.3kat==2])
shapiro.test(kardio$sys_tlak[kardio$vek.3kat==3])


### 2. Overeni homogenity rozptylu:

# 2.1 Histogram, boxplot - viz overeni normality

# 2.2 Bartlettuv test (srovnani homogenity rozptylu vice nez dvou vyberu z normalniho rozdeleni):
bartlett.test(sys_tlak~vek.3kat,data=kardio) 
# Bartlettuv test predpoklada normalitu puvodnich hodnot!!!!

# 2.3 Levenuv test:
Levene <- function(y,skupina){
  skupina <- as.factor(skupina) # preventivni prevedeni na faktor
  meds <- tapply(y, skupina, median)
  resp <- abs(y - meds[skupina])
  anova(lm(resp ~ skupina))[1, 4:5]
}
Levene(kardio$sys_tlak,kardio$vek.3kat) # zavolani vytvorene funkce

# Ani rozptyl neni homogenni! ANOVA by se nemela pocitat - tyhle skripty jsou jen na ukazku, jak se funkce volaji!!!


### ANOVA:
anova.test <- anova(lm(sys_tlak~vek.3kat,data=kardio)) # je potreba nejdriv vytvorit linearni model pomoci lm
anova.test
anova.tabulka <- as.matrix(anova.test)
anova.tabulka
anova.tabulka[1,5]                # vypsani p-hodnoty  

# 2. zpusob: 
anova.test <- aov(sys_tlak~vek.3kat,data=kardio)
summary(anova.test)

# zjisteni, ktere skupiny se od sebe lisi:
pairwise.t.test(kardio$sys_tlak,kardio$vek.3kat,p.adj="bonferroni")

# Ukol: zjistete, ktere dalsi metody korekce p-hodnoty existuji:

# Reseni:
? pairwise.t.test
? p.adjust

tukey_hsd <- TukeyHSD(aov(kardio$sys_tlak~as.factor(kardio$vek.3kat)))
tukey_hsd


################ Kruskaluv - Wallisuv test ###################

kruskal.test(sys_tlak~vek.3kat,data=kardio)

