################ KORELACE ################################################################
DATA <- read.csv("11_data.csv", sep=";", header = TRUE)
attach(DATA)
str(DATA)
DATA$Sex <- as.factor(DATA$Sex)
DATA$Side <- as.factor(DATA$Side)
str(DATA)

### PRIPRAVA DAT
## Pouze kompletni pripady
data <- DATA[complete.cases(DATA),]
head(data)
attach(data)
## Mame antropologicka data
library(reshape)
leva <- na.omit(data[Side == "sin",])
prava <- na.omit(data[Side == "dx",])
SINDX <- merge(leva, prava, by = "No_sub")
attach(SINDX)
# Chceme zjistit korelaci mezi Length_cl - Shoulder_w
# mame normalni rozlozeni dat nebo ne
# install.packages("ggpubr")
library(ggplot2)
library(ggpubr)
# moznost vykresleni grafu
ggqqplot(SINDX[Sex.x == "female","Length_cl.x"])
ggqqplot(SINDX[Sex.x == "male","Length_cl.x"])
ggqqplot(SINDX[Sex.x == "female","Length_cl.y"])
ggqqplot(SINDX[Sex.x == "male","Length_cl.y"])

ggqqplot(SINDX[Sex.x == "female","Shoulder_w.x"])
ggqqplot(SINDX[Sex.x == "male","Shoulder_w.x"])
ggqqplot(SINDX[Sex.x == "female","Shoulder_w.y"])
ggqqplot(SINDX[Sex.x == "male","Shoulder_w.y"])

# Nebo S-W test normality napr. dva pripady
shapiro.test(SINDX[Sex.x == "male","Length_cl.x"])
shapiro.test(SINDX[Sex.x == "female","Shoulder_w.y"])

# p hodnota je vyssi nez  0.05 tzn. nezamitame nulovou hypotezu
# data maji normalni rozlozeni.
# Pozor na velikost vzorku. S-W je citlivy na effect size, tzn. Q-Q plot
# je silne doporucen.

ggscatter(SINDX[Sex.x == "male",], x = "Length_cl.x", y = "Shoulder_w.x", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "Length of clavicle", ylab = "Shoulder Width")
ggscatter(SINDX[Sex.x == "female",], x = "Length_cl.x", y = "Shoulder_w.x", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "Length of clavicle", ylab = "Shoulder Width")
ggscatter(SINDX[Sex.x == "male",], x = "Length_cl.y", y = "Shoulder_w.y", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "Length of clavicle", ylab = "Shoulder Width")
ggscatter(SINDX[Sex.x == "female",], x = "Length_cl.y", y = "Shoulder_w.y", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "Length of clavicle", ylab = "Shoulder Width")


# vypocet korelacniho koeficientu
res1 <- cor.test(SINDX[Sex.x == "male","Length_cl.x"], 
                 SINDX[Sex.x == "male","Shoulder_w.x"], 
                 method = "pearson")
res1

# t - testová statistika (t-test statistic value)
# df - stupne volnosti (degree of freedom)
# p - hodnota významnosti testu (significance level of the t-test)
# 95 interval spolehlivosti
# cor - korelacni koeficient (correlation coefficient)

# moznost ziskani prislusnych hodnot z objektu res1
res1$p.value
res1$estimate

# p-hodnota je 0.002941 coz je mene nez 0.05 a muzeme tak odvodit, ze pro muze jsou hodnoty
# prave strany pro promenne Length_cl a Shoulder_W zavisle a koeficient korelace je 0.5698631
# Pozor na to, co je vasi nulovou hypotezou.

### mame vice promenych
### KORELACNI MATICE
## priprava dat
# Pr. muzi leva strana (chceme pouze rozmery!) - funkce cor neumoznuje korelovat faktory a zaroven korelovat napr.
# No_Sub nebo age zde nema vyznam
male_sin <- na.omit(data[Side == "sin" & Sex == "male",7:28])
res1 <- cor(male_sin)
round(res1, 2)

# velka matice vysledku avsak nemame p-hodnoty (tzn. nevime ktere z vysledku jsou signifikanti)
# install.packages("Hmisc")
library(Hmisc)
res2 <- rcorr(as.matrix(male_sin))
res2

# vysledek neni prilis prehledny (mame moc promennych)
# install.packages("corrplot")
library(corrplot)
windows()
corrplot(res1, type = "upper", order = "hclust", tl.col = "blaCK", tl.srt = 45)
corrplot(res1, type = "upper", tl.col = "black", tl.srt = 45)
# tl.col - barva popisku
# tl.srt - uhel popisku
# order = "hclust" - hierarchical clustering order (optional)
# Note: funkcí corrplot() zobrazujeme matici hodnot (-1 az 1), tzn. muzeme tak
# zobrazit i jinou, nez korelacni matici

# vyuzijeme res2 kde mame ulozene jak r - koeficienty korelace tak p - hodnoty
# nastavime sig. level 0.01 (insig = "pch" hodnoty pod 0.01 zakrizkuje)
corrplot(res2$r, type="upper", p.mat = res2$P, sig.level = 0.01, insig = "pch")
# nebo muzeme pouzit insig = "blank", kdy jsou jednoduse nevyznamne korelace (na hladine 0.01) prazde
corrplot(res2$r, type="upper", p.mat = res2$P, sig.level = 0.01, insig = "blank")

# corrplot umoznuje mnoho zpusobu zobrazeni
# priklad: ponechame jak horni, tak dolni trojuhelnik (vuci diagonale) - note: parametr type
corrplot(res1, tl.col = "black", tl.srt = 45)
# misto velikosti kruhu lze pouzit pie charts
corrplot(res1, tl.col = "black", tl.srt = 45, method = "pie")
# nebo muzeme pouzit primo barevne ctverce
corrplot(res1, tl.col = "black", tl.srt = 45, method = "color")
# nebo primo barevne hodnoty
corrplot(res1, tl.col = "black", tl.srt = 45, method = "number")
# chceme ctverce (barevne) a hodnoty koeficientu korelace (cerne)
windows()
corrplot(res2$r, method="color", type="upper", order= "hclust", 
         addCoef.col = "black", tl.col="black", tl.srt=45, diag=FALSE )
# addCoef.col - barva pridaneho koeficientu korelace
# diag = FALSE - nechceme diagonalu

# moznost jiz zmineneho order = "hclust" a type = "lower
corrplot(res1, type = "lower", order = "hclust", tl.col = "black", tl.srt = 45)
# vlevo dole spise zaporne korelace...
corrplot(res1, type = "upper", order = "hclust", tl.col = "black", tl.srt = 45)

# moznost vyuzit vlastni barevne schema
col <- colorRampPalette(c("darkmagenta", "white", "darkgreen"))(20)
corrplot(res1, type="upper", order="hclust", col=col, tl.col = "black")
# nebo barevne schema pomoci RcolorBrewer
library(RColorBrewer)
corrplot(res1, type="upper", order="hclust", col=brewer.pal(n = 11, name="BrBG"))
# u RColorBrewer pozor na mnozstvi barev, ne vsechny palety umoznuji neomezene mnozstvi (spise mnohem mene)

# muzeme nastavit i pozadi
corrplot(res1, type="upper", order="hclust", col=col, bg = "dimgrey", tl.col = "black")

# Bonus: zajimava moznost: balicek performanceAnalytics
# install.packages("PerformanceAnalytics", dependencies=TRUE, type="win.binary")
library(PerformanceAnalytics)

chart.Correlation(male_sin[,1:6], histogram=TRUE, pch=19)
# note muzeme zobrazit i pro vsechny hodnoty z objektu male_sin, avsak pozor, jedna se o vetsi objekt
windows()
chart.Correlation(male_sin, histogram=TRUE, pch=19)
# levy trojuhelnik ukazuje bivariantni rozlozeni (bodove grafy) s prolozenou krivkou
# pravy trojuhelnik ukazuje hodnotu korelace a zaroven jeji vyznam pomoci hvezdicek "***","**",...
# p-values(0, 0.001, 0.01, 0.05, 0.1, 1) = symbols(“***”, “**”, “*”, “.”, " “)


######################## DATUM V R #################################
# install.packages("chron")
library(chron)

Sys.time() # systemovy cas
Sys.Date() # systemove datum
date()     # datum

cas <- Sys.time()
cas
typeof(cas)

cas2 <- as.POSIXlt(cas)
cas2
typeof(cas2)
unclass(cas2)

cas2$`sec` # sekundy
cas2$min   # minuty
cas2$hour  # hodiny
cas2$mday  # den mesice (1 - 31)
cas2$mon   # mesic (0 - 11, 11 = prosinec)
cas2$year  # rok od 1900 (tzn. pto 2020 je to 120ty rok)
cas2$wday  # den tydne (0 - 6, 0 = nedele)
cas2$yday  # den roku (0 - 365)
cas2$isdst # daylight savings indicator
cas2$zone  # casova zona

cas3 <- as.POSIXct(cas)
cas3
typeof(cas3)

cas4 <- as.Date(cas)
cas4
typeof(cas4)

as.Date("13/3/1992", format = "%d/%m/%Y")
as.Date("Březen 13, 1992", format = "%B %d, %Y")
as.Date("13-Březen-1992", format= "%d-%B-%Y")

letapane <- c("10/7/1856", "3/1/1892", "21/6/1948")
letapane
# note: 
# 10. 7. 1856 *Nikola Tesla
# 3. 1. 1892 *J. R. R. Tolkien
# 21. 6. 1948 *A. Sapkowski

rok_narozeni <- as.Date(letapane, format = "%d/%m/%Y")
rok_narozeni

### Import dat numerickeho formatu
## Excel
data <- c(17431, 15349)
data2 <- as.Date(data, origin = "1899-12-30")
data2
# note:
# 21. 9. 1947 * S. King
# 8. 1 .1942 * S. Hawking

## Excel on Mac
data <- c(17426, 38540)
data3 <- as.Date(data, origin = "1904-01-01")
data3

### zmena formatu data
rok_narozeni
format(rok_narozeni, "%a %b %d")
format(rok_narozeni, "%A %B %d %Y")

min(rok_narozeni)
max(rok_narozeni)

rok_narozeni
plus <- rok_narozeni + 5
plus

