################ MAPA SVETA - ANTROPOLOGICKA DATA (STUDIE)
# nacteni dat
data <- read.table("11_mapy_data.txt",header = T)
str(data)
attach(data)
newmap <- getMap(resolution = "low")                      # mapa sveta - nizke rozliseni
windows()
plot(newmap, col="grey",  border="white", main = "World") # barva oznacuje barvu vykreslenych zemi, hranice chceme bile a nadpis grafu chceme World
points(data$long, data$lat, col = alpha(4,0.2), cex = data$log, pch = 16) # stejne jako u 2D grafu vykreslime data
# v rámci argumentu col je pouzita funkce alpha z balicku scales, ktera nam umoznuje zobrazit polopruhledne body
# prvni argument je barva (v tomto pripade 4) a druhy je pro nastaveni pruhlednosti

mp <- NULL                                                   # prazdny objekt
mapWorld <- borders("world", colour="gray50", fill="gray50") # vytvorime vrstvu ktera obsahuje hranice sveta
mp <- ggplot() +   mapWorld                                  # do prazdenoho objektu vlozime ggplot() vcetne vrstvy mapWorld
# mp zobrazime s body v sam. okne
windows(width = 20, height = 10)
mp + geom_point(data = data, aes(x = long, y = lat, alpha = 0.8),
                size = c(data$log*2), shape = 19, col = "darkblue") + 
  guides(fill=FALSE, alpha=FALSE, size=FALSE)
# argument alpha v aes funguje podobne jako funkce alpha z balicku scales, tzn. nastavujeme pruhlednost bodu
# note: typ bodu je zde nastavovan ne pres pch, ale skrz shape

################ Bonus - globe ################
# install.packages("maptools")
# install.packages("maps")
# install.packages("threejs")
library(threejs)
library(maptools)
library(maps)
data(world.cities, package="maps")
cities <- world.cities[order(world.cities$pop,decreasing=TRUE)[1:1000],]
value  <- 100 * cities$pop / max(cities$pop)

globejs(bg="black", lat=cities$lat, long=cities$long, value=value, rotationlat=-0.34, rotationlong=-0.38, fov=30)

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

### PRIPRAVA DAT
## Pouze kompletni pripady
data <- DATA[complete.cases(DATA),]
head(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(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.002602379 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.622
# 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[,1:6], histogram=TRUE, pch=19)
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(“***”, “**”, “*”, “.”, " “)

########
# A) Nactěte data "11_data.csv" a vytvorte correlogram z promennych:
# Length_cl, Shoulder_W, Proportion, Inner_C, Outer_C, Total_C (7 až 12)
# pouze pro zeny a pro pravou stranu.
# Nezapomente vytvorit matici p-hodnot.
# Nezobrazujte diagonalu.
# Barvy correlogramu volte libovolne, stejne tak typ zobrazeni.

