##########################################################################################
################################## SOFTWARE R - CAST 2 ###################################
##########################################################################################

# Nastaveni pracovniho adresare - napr.:
setwd("c:/Users/brozova/Desktop/vyuka/vicerozmerky_podzim_2016_bimat/1. cviceni")
getwd()

###########################################################################
###############         Nacteni dat a jejich kontrola       ###############
###########################################################################

kardio <- read.table(file="Data_pro_R.txt", sep="\t", dec=".", header=T) # nacteni dat

dim(kardio)                       # pocet radku a sloupcu datove matice
head(kardio)                      # prvnich 10 radku tabulky
str(kardio)                       # sumarni informace o jednotlivych promennych
is(kardio)				    # typ objektu

###########################################################################
########################         Normalita      ###########################
###########################################################################

########## Graficke overeni

##### 1. histogram:
hist(kardio$sys_tlak,xlab="Systolicky tlak",main="Histogram")

# Histogram s krivkou hustoty normalniho rozdeleni
windows()
x <- kardio$sys_tlak
h<-hist(x, col="yellow", xlab="vyska",
        main="Histogram s krivkou hustoty normalniho rozdeleni")
xfit<-seq(min(x),max(x),length=40)
yfit<-dnorm(xfit,mean=mean(x),sd=sd(x))
yfit <- yfit*diff(h$mids[1:2])*length(x)
lines(xfit, yfit, col="blue", lwd=2) 

# Nastaveni vykresleni vice grafu do 1 okna
windows()
par(mfrow=c(1,2))                
hist(kardio$sys_tlak[kardio$pohlavi_rek=="Z"],breaks=5,cex.axis=0.8,xlab="zeny",main="Systolicky tlak")
hist(kardio$sys_tlak[kardio$pohlavi_rek=="M"],breaks=5,cex.axis=0.8,xlab="muzi",main="Systolicky tlak")

##### 2. boxplot:
boxplot(kardio$sys_tlak,ylab="Systolicky tlak")
boxplot(kardio$sys_tlak~kardio$pohlavi_rek,names=c("Muzi","Zeny"),ylab="Systolicky tlak")

##### 3. Q-Q graf:
qqnorm(kardio$sys_tlak)
par(mfrow=c(1,2))                 
qqnorm(kardio$sys_tlak[kardio$pohlavi_rek=="Z"])
qqnorm(kardio$sys_tlak[kardio$pohlavi_rek=="M"])

########## Overeni pomoci testu

##### 1. Shapiruv-Wilkuv test - doporucovan pro mensi vzorky, muze byt moc prisny pro velke vzorky:
shapiro.test(kardio$sys_tlak)
shapiro.test(kardio$sys_tlak[kardio$pohlavi_rek=="Z"])
shapiro.test(kardio$sys_tlak[kardio$pohlavi_rek=="M"])

##### 2. Kolmogorovuv-Smirnovuv test (Liliefosova varianta)

# Lilieforsovu variantu testu pocitame vzdy, kdyz nezname stredni hodnotu a rozptyl dat
install.packages("nortest")
library("nortest")
lillie.test(kardio$sys_tlak)
lillie.test(kardio$sys_tlak[kardio$pohlavi_rek=="Z"])
lillie.test(kardio$sys_tlak[kardio$pohlavi_rek=="M"])

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

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

## Priklad 1.
# 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

# pred vypoctem je nutne overit normalitu dat
hist(kardio$sys_tlak,xlab="Systolicky tlak",main="Histogram")

alpha <- 0.05                     # hladina vyznamnosti
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
test.stat <- (prumer - prumer.populace)/(s/sqrt(n)) # realizace testove statistiky
test.stat
p.hodnota <- 2 * (1 - pt(test.stat,n-1))
p.hodnota 				    #p-hodnota > 0.05 H0 nezamitame 

# interval spolehlivosti naseho vyberoveho prumeru z dat
IS.dolni <- prumer - (s/sqrt(n))*kvantil
IS.horni <- prumer + (s/sqrt(n))*kvantil
IS.dolni;prumer;IS.horni 
#prumer.populace = 131 se nachazi v 95% intervalu spolehlivosti vyberoveho prumeru -> H0 nezamitame

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

## Priklad 2.
# 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 # vstupem jsou diference hodnot

# pred vypoctem je nutne overit normalitu diferenci
hist(rozdil.tlaku,xlab="Po-pred lecbou",main="Zmena systolickeho tlaku",breaks=5)
# data se neridi normalnim rozdelenim, zvolili bychom neparametrickou variantu testu, ukazeme si ale, jak by se spocitalo parametricky 

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 3.
# 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

# pred vypoctem je nutne overit normalitu pro jednotlive skupiny zvlast
windows()
par(mfrow=c(1,2))                
hist(kardio$sys_tlak[kardio$pohlavi_rek=="Z"],breaks=5,cex.axis=0.8,xlab="zeny",main="Systolicky tlak")
hist(kardio$sys_tlak[kardio$pohlavi_rek=="M"],breaks=5,cex.axis=0.8,xlab="muzi",main="Systolicky tlak")# pro demonstrativni ucely je dostatacne 
boxplot(sys_tlak~pohlavi_rek,data=kardio)

t.test(kardio$sys_tlak~kardio$pohlavi_rek)
# 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=="M"],kardio$sys_tlak[kardio$pohlavi_rek=="Z"])
t.test(sys_tlak~pohlavi_rek,data=kardio)

# Priklad 4: 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.

windows()
par(mfrow=c(1,2)) 
hist(kardio$vaha[kardio$pohlavi_rek=="Z"],cex.axis=0.8,xlab="zeny",main="Vaha")
hist(kardio$vaha[kardio$pohlavi_rek=="M"],cex.axis=0.8,xlab="muzi",main="Vaha")

windows()
par(mfrow=c(1,2)) 
hist(log(kardio$vaha[kardio$pohlavi_rek=="Z"]),cex.axis=0.8,xlab="zeny",main="log(Vaha)")
hist(log(kardio$vaha[kardio$pohlavi_rek=="M"]),cex.axis=0.8,xlab="muzi",main="log(Vaha)")
 # neni to idealni, ale pro demonstracni ucely nam to bude stacit

t.test(log(vaha)~pohlavi_rek,data=kardio)

## Priklad 5. 
# 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: ###########################

## Priklad 6. 
# 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
t.test(sys_tlak~pohlavi_rek,data=kardio,alternative="greater") # abychom dostali odpovidajici p-hodnotu musime "less" nahradit za "greater"

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

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

## Samostatny ukol 1. 
# 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) Spocitejte median BMI, aritmeticky a geometricky prumer.
# c) Overte graficky normalitu BMI.
# d) Vzhledem k tomu, ze BMI se neridi normalnim rozlozenim overte pomoci Wilcoxonova testu na hladine 
#    vyznamnosti alpha=0.05, zda lide trpi nadvahou.
# e) Interpretujte zaver testovani.

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

## Reseni:
# a) Nulova a alternativni hypoteza:
# b) Median BMI, aritmeticky a geometricky prumer:
# c) Graficke overeni normality BMI:
# d) Testovani:
# e) Zaver


############### Mannuv-Whitneyuv test - neparametricky test pro dva vybery ###############

## Samostatny ukol 2: 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 u muzu a zen zvlast.
# c) Vzhledem k tomu, ze BMI se neridi normalnim rozlozenim, pouzijte Mannuv-Whitneyuv test (ten je schovan v R pod funkci 
#    wilcox.test, protoze Mannuv-Whitneyuv test lze nazvat tez dvouvyberovy Wilcoxonuv test).
# d) Zaver testovani.

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

## Reseni:
# a) Nulova a alternativni hypoteza:
# b) Graficke overeni normality BMI u muzu a u zen:
# c) Overeni na hladine vyznamnosti alpha=0.05, zda se lisi BMI podle pohlavi.
# d) Zaver: 


###########################################################################################
############################## Testovani pro 3 a vice skupin ##############################
###########################################################################################

############### parametricky - ANOVA ################

### Rekategorizace veku:
kardio$vek.3kat <- 1+(kardio$vek>=30)+(kardio$vek>=50) # vytvoreni trech kategorii: 1...<30, 2...30-49, 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],main="Vek: <30",xlab="Systolicky tlak")
hist(kardio$sys_tlak[kardio$vek.3kat==2],main="Vek: 30-49",xlab="Systolicky tlak")
hist(kardio$sys_tlak[kardio$vek.3kat==3],main="Vek: >=50",xlab="Systolicky tlak")
# 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,ylab="Systolicky tlak",names=c("<30","30-49",">=50"),xlab="Vek")

# 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

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

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

###########################################################################################
######################################### Korelace  #######################################
###########################################################################################

# Priklad: Overte vztah celkoveho cholesterolu a LDL cholesterolu.

# Bodovy graf:
plot(kardio$ldl_cholesterol,kardio$cel_cholesterol)

# Pearsonuv korelacni koeficient:
cor.test(kardio$cel_cholesterol,kardio$ldl_cholesterol)

# Spearmanuv korelacni koeficient:
cor.test(kardio$cel_cholesterol,kardio$ldl_cholesterol,method="spearman")

# Korelacni matice:
cor(cbind(kardio$cel_cholesterol,kardio$ldl_cholesterol))
cor(kardio[,10:16]) # je patrne, ze jsou tam nejake chybejici hodnoty
nevalidni <- (rowSums(is.na(kardio[,10:16]))>0) # nalezeni radku s chybejicimi hodnotami
kardio.valid <- kardio[nevalidni==F,] # vytvoreni validniho souboru
cor(kardio.valid[,10:16]) 

###########################################################################################
################################### Linearni regrese ######################################
###########################################################################################

y=kardio$cel_cholesterol
x=kardio$ldl_cholesterol
lin.model <- lm(y~x)
plot(x,y)
lines(x,as.numeric(lin.model$coefficients[1])+as.numeric(lin.model$coefficients[2])*x,col="red")



