# Skript ke cviceni aplikovana statistika II - cviceni pro antropology

# Resene priklady (tj. pracovni listy)

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

#############################################################################
# 2. cviceni ###
#############################################################################

data <- read.table('Data/newborns.txt', header=TRUE)
summary(data)

data <- na.omit(data)
data$edu.M <- factor(data$edu.M, labels=c('ZS', 'SS', 'SSmat', 'VS'))

table(data$sex.C, data$edu.M)

boxplot(data$weight.C ~ data$sex.C + data$edu.M, varwidth=T, notch=T,
        xlab='Skupina', ylab='Porodni hmotnost [g]')
mm1 <- mean(data[data$edu.M=='ZS' & data$sex=='f','weight.C'])
mm2 <- mean(data[data$edu.M=='ZS' & data$sex=='m','weight.C'])
mm3 <- mean(data[data$edu.M=='SS' & data$sex=='f','weight.C'])
mm4 <- mean(data[data$edu.M=='SS' & data$sex=='m','weight.C'])
mm5 <- mean(data[data$edu.M=='SSmat' & data$sex=='f','weight.C'])
mm6 <- mean(data[data$edu.M=='SSmat' & data$sex=='m','weight.C'])
mm7 <- mean(data[data$edu.M=='VS' & data$sex=='f','weight.C'])
mm8 <- mean(data[data$edu.M=='VS' & data$sex=='m','weight.C'])
points(1:8, c(mm2,mm2,mm3,mm4,mm5,mm6,mm7,mm8), col='red', pch=16)

model.newborns <- aov(weight.C ~ sex.C*edu.M, data=data)
par(mfrow=c(2,2)) #nastavi zobrazeni 4 grafu najednou (na 2 radky a 2 sloupce)
plot(model.newborns)

shapiro.test(model.newborns$residuals)

# library(nortest)
# lillie.test(model.newborns$residuals)

anova(model.newborns)

# pro interakci
qf(0.95, 3, 1376)
## [1] 2.61137
# pro faktor B
qf(0.95, 3, 1376)
## [1] 2.61137
# pro faktor A
qf(0.95, 1, 1376)
## [1] 3.848226

TukeyHSD(model.newborns, which=c('sex.C','edu.M'))

a <- 2
b <- 4
n <- nrow(data)
m.boys <- mean(data[data$sex.C =='m','weight.C'])
m.girls <- mean(data[data$sex.C =='f','weight.C'])
n.sex <- as.numeric(table(data$sex.C))
m.zs <- mean(data[data$edu.M =='ZS','weight.C'])
m.ss <- mean(data[data$edu.M =='SS','weight.C'])
m.ssmat <- mean(data[data$edu.M =='SSmat','weight.C'])
m.vs <- mean(data[data$edu.M =='VS','weight.C'])
n.edu <- as.numeric(table(data$edu.M))
res.mean.sq <- anova(model.newborns)['Residuals','Mean Sq']

# chlapci-devcata
abs(m.boys - m.girls)
## [1] 106.7074
sqrt((a-1) * (1/n.sex[1] + 1/n.sex[2]) * res.mean.sq * qf(0.95, a-1, n-a*b))
## [1] 72.64113

# ZS-SS
abs(m.zs - m.ss)
## [1] 198.9039
sqrt((b-1) * (1/n.edu[1] + 1/n.edu[2]) * res.mean.sq * qf(0.95, b-1, n-a*b))
## [1] 130.9452
# ZS-SSmat
abs(m.zs - m.ssmat)
## [1] 219.0501
sqrt((b-1) * (1/n.edu[1] + 1/n.edu[3]) * res.mean.sq * qf(0.95, b-1, n-a*b))
## [1] 131.9525
# ZS-VS
abs(m.zs - m.vs)
## [1] 225.536
sqrt((b-1) * (1/n.edu[1] + 1/n.edu[4]) * res.mean.sq * qf(0.95, b-1, n-a*b))
## [1] 233.9152
# SS-SSmat
abs(m.ss - m.ssmat)
## [1] 20.14623
sqrt((b-1) * (1/n.edu[2] + 1/n.edu[3]) * res.mean.sq * qf(0.95, b-1, n-a*b))
## [1] 129.4559
# SS-VS
abs(m.ss - m.vs)
## [1] 26.6321
sqrt((b-1) * (1/n.edu[2] + 1/n.edu[4]) * res.mean.sq * qf(0.95, b-1, n-a*b))
## [1] 232.516
# SSmat-VS
abs(m.ssmat - m.vs)
## [1] 6.48587
sqrt((b-1) * (1/n.edu[3] + 1/n.edu[4]) * res.mean.sq * qf(0.95, b-1, n-a*b))
## [1] 233.0847

coefs <- aov(weight.C ~ sex.C + edu.M, data=data)$coefficients
girls <- coefs[1] + c(0, coefs[3], coefs[4], coefs[5])
boys <- coefs[2] + girls
boxplot(data$weight.C~data$edu.M, ylim=c(2800, 3300), border='white',
        xlab='Vzdelani matky', ylab='Porodni hmotnost')
points(1:4,girls, col='red', pch=16)
lines(1:4,girls, col='red', lty=2)
points(1:4, boys, col='blue', pch=17)
lines(1:4, boys, col='blue', lty=2)
legend('topleft', legend=c('divky', 'chlapci'), col=c('red', 'blue'), pch=c(16,17))