# Skript ke cviceni aplikovana statistika II - cviceni pro antropology

# Resene priklady (tj. pracovni listy)

# Nastaveni pracovniho adresare:
setwd("...")

#############################################################################
# 1. cviceni ###
#############################################################################

# Nacteme data
data <- read.table("Data/clavicle.txt", header = TRUE)

# Zakladni prehled o datech
summary(data)

# Aby dalsi prikazy fungovaly, musi byt data$population kodovano jako faktor
is.factor(data$population)

data$population <- as.factor(data$population)

# Odstraneni NA-hodnot
data <- na.omit(data)
(t <- table(data$population))

(m1 <- mean(data[data$population=='gre','cla.L']))

(n <- sum(table(data$population)))

(m2 <- mean(data[data$population=='ind1','cla.L']))

(m3 <- mean(data[data$population=='ind2','cla.L']))

(m <- mean(data[,'cla.L']))

(s1 <- sd(data[data$population=='gre','cla.L']))

(s2 <- sd(data[data$population=='ind1','cla.L']))

(s3 <- sd(data[data$population=='ind2','cla.L']))

(s <- sd(data[,'cla.L']))

boxplot(cla.L ~ population, data=data, var.width=T, notch=T, xlab="Populace",
        ylab="delka klicni kosti")
points(1:3, c(m1,m2,m3), col="red", pch=16)


par(mfrow=c(1,3))
qqnorm(data[data$population=='gre','cla.L'],main="Populace recka", xlab="teoreticky kvantil",
       ylab='pozorovany kvantil')
qqline(data[data$population=='gre','cla.L'])
shapiro.test(data[data$population=='gre','cla.L'])

qqnorm(data[data$population=='ind1','cla.L'],main="Populace indicka z Armitsaru",
       xlab="teoreticky kvantil", ylab='pozorovany kvantil')
qqline(data[data$population=='ind1','cla.L'])
shapiro.test(data[data$population=='ind1','cla.L'])

qqnorm(data[data$population=='ind2','cla.L'],main="Populace indicka z Varansi",
       xlab="teoreticky kvantil", ylab='pozorovany kvantil')
qqline(data[data$population=='ind2','cla.L'])

shapiro.test(data[data$population=='ind2','cla.L'])

bartlett.test(cla.L ~ population, data=data)

library("car")

leveneTest(cla.L ~ population, data=data)

an1 <- aov(cla.L ~ population, data=data)
anova(an1)

an2 <- lm(cla.L ~ population, data=data)
summary(an2)

anova(an2)

TukeyHSD(aov(cla.L~population, data=data))

K <- 3
alpha <- 0.05
#gre vs ind1
abs(m1-m2)

s*sqrt((K-1)*(1/t[1] + 1/t[2]) * qf(1-alpha,K-1,n-K))

#gre vs ind2
abs(m1-m3)

s*sqrt((K-1)*(1/t[1] + 1/t[3]) * qf(1-alpha,K-1,n-K))

#ind1 vs ind2
abs(m2-m3)

s*sqrt((K-1)*(1/t[2] + 1/t[3]) * qf(1-alpha,K-1,n-K))

pairwise.t.test(data$cla.L, data$population, p.adjust.method="bonferroni", pool.sd=T)