### 9 DVOUVYBEROVE PARAMETRICKE TESTY
# 9.1 Klasicky test o rozdilu strednich hodnot
data <- 'TO-DO'
data <- 'TO-DO' # odstraneni NA
skull.BM <- 'TO-DO' # vyber muzi, promenna skull.B
skull.BF <- 'TO-DO' # vyber zeny, promenna skull.B
n1 <- 'TO-DO' # pocet muzi
n2 <- 'TO-DO' # pocet zeny
'TO-DO' # rozsah promennych

## test predpokladu normality
'TO-DO' # normalita pro zeny
'TO-DP' # normalita pro muze

## graficke posouzeni normality
# muzi
r <- round(1 + 3.3*log10(n1)) 
'TO-DO' # histogram a qqplot

# zeny
r <- round(1 + 3.3*log10(n2)) 
'TO-DO' # histogram a qqplot

## pro vyber testu musime otestovat, zda plati shodnost rozptylu
alpha <- 'TO-DO'
(var.test(skull.BM, skull.BF, alt = 'TO-DO', conf.level = 1 - alpha))

# pomoci vzorce pro kriticky obor
s1 <- sd(skull.BM)
s2 <- sd(skull.BF)
tw <- s1^2 / s2^2
q1 <- qf(alpha / 2, n1 - 1, n2 - 1)
q2 <- qf(1 - alpha / 2, n1 - 1, n2 - 1)
(tab <- data.frame(q1, q2))

# pokud chceme vypsat pouze interval spolehlivosti z funkce var.test, muzeme
# to udelat pomoci $conf.int[1:2] (vypise samotny interval)
# nebo $conf.int (vypise krome intervalu i spolehlivost)
var.test('TO-DO')$conf.int[1:2]

# obdobne si muzeme vypsat pouze samotnou p-hodnotu pomoci $p.value
var.test('TO-DO')$p.value


## klasicky test o rozdilu strednich hodnot
t.test(skull.BM, skull.BF, var.equal = 'TO-DO', conf.level = 1 - alpha)

# pomoci vzorcu
m1 <- mean(skull.BM)
m2 <- mean(skull.BF)
s1 <- sd(skull.BM)
s2 <- sd(skull.BF)
s_h2 <- ((n1 - 1) * s1^2 + (n2 - 1) * s2^2) / (n1 + n2 - 2)
s_h <- sqrt(s_h2)
c <- 0

tw <- (m1 - m2 - c) / (s_h * sqrt(1 / n1 + 1 / n2)) # test. statistika
q1 <- qt( alpha / 2, n1 + n2 - 2)
q2 <- qt( 1 - alpha / 2, n1 + n2 - 2)

(tab <- data.frame(tw, q1, q2)) # test. statistika a hranice kritickeho oboru

dh <- m1 - m2 - s_h * sqrt(1 / n1 + 1 / n2) * qt(1 - alpha / 2, n1 + n2 - 2)
hh <- m1 - m2 - s_h * sqrt(1 / n1 + 1 / n2) * qt(alpha / 2, n1 + n2 - 2)
(tab <- data.frame(dh, hh)) # interval spolehlivosti

(p.hodnota <- 2 * min(pt(tw, n1 + n2 - 2), 1 - pt(tw, n1 + n2 - 2))) # p-hodnota


## graficka interpretace vysledku
par(mar = c(2, 4, 1, 1))
boxplot(skull.BM, skull.BF, col = c('mintcream', 'blanchedalmond'),
        medcol = c('deepskyblue4', 'brown4'),  border = c('darkblue', 'darkred'),
        ylim = c(120, 155),  xlab = '', ylab = 'sirka mozkovny (v mm)', 
        names = c('muzi', 'zeny'), las = 1)
mtext('pohlavi', line = 2.5, side = 1)
points(c(m1, m2), pch = 19, col = c('blue', 'red'), cex = 1.2)
legend('topright', pch = c(NA, 19), lty = c(1, NA), lwd = c(2, NA), 
       legend = c('median', 'prumer'), bty = 'n', horiz = T)



# 9.2 Welchuv test o rozdilu strednich hodnot
data <- 'TO-DO'
data <- 'TO-DO' # odstraneni NA
nose.BC <- 'TO-DO' # vyber populace 'cin', promenna nose.B
nose.BB <- 'TO-DO' # vyber populace 'ban', promenna nose.B
n1 <- 'TO-DO' # pocet pro poulaci cin
n2 <- 'TO-DO' # pocet pro poulaci ban
'TO-DO' # rozsah promennych pro obe populace

## test predpokladu normality
'TO-DO'
'TO-DO' # zvlast pro kazdou populaci

## graficke posouzeni normality
# cinska populace
'TO-DO' # histogram a qqplot

# bantuska populace
'TO-DO' # histogram a qqplot

## test o podilu rozptylu
alpha <- 'TO-DO'
'TO-DO'

# pomoci vzorcu
s1 <- sd(nose.BC)
s2 <- sd(nose.BB)
tw <- s1^2 / s2^2
q1 <- qf(alpha / 2, n1 - 1, n2 - 1)
q2 <- qf(1 - alpha / 2, n1 - 1, n2 - 1)
(tab <- data.frame(q1, q2))


## Welchuv test o rovnosti strednich hodnot
alpha <- 'TO-DO'
t.test('TO-DO', var.equal = 'TO-DO', alt = 'TO-DO', conf.level = 'TO-DO')

# pomoci vzorcu
m1 <- mean(nose.BC)
m2 <- mean(nose.BB)
s1 <- sd(nose.BC)
s2 <- sd(nose.BB)
c <- 0

tw <- (m1 - m2 - c) / sqrt(s1 ^ 2 / n1 + s2 ^ 2 / n2)
df <- (s1 ^ 2 / n1 + s2 ^ 2 / n2) ^ 2 / ((s1 ^ 2 / n1) ^ 2 / (n1 - 1) + (s2 ^ 2 / n2) ^ 2 / (n2 - 1))
q <- qt(alpha, df)
(tab <- data.frame(tw, q)) # test. statistika a kriticky obor

hh <- m1 - m2 - sqrt(s1 ^ 2 / n1 + s2 ^ 2 / n2) * qt(alpha, df)
(tab <- data.frame(hh)) # interval spolehlivosti

(p.hodnota <- pt(tw, df)) # p-hodnota


## graficka interpretace vysledku
par(mar = c(2, 4, 1, 1))
boxplot('TO-DO')
points('TO-DO') # pridani prumeru do boxplotu
legend('TO-DO') # pridani legendy



# 9.3 Test o rozdilu korelacnich koeficientu
data <- 'TO-DO'
'TO-DO' # odstraneni NA
data.M <- 'TO-DO' # vyber muzi
data.F <- 'TO-DO' # vyber zeny
lowex.LM <- 'TO-DO' # promenna lowex.L pro muze
tru.LM  <- 'TO-DO' # promenna tru.L pro muze
lowex.LF <- 'TO-DO' # promenna lowex.L pro zeny
tru.LF <- 'TO-DO' # promenna tru.L pro zeny

n1 <- 'TO-DO' # rozsah prvniho nah. vyberu
n2 <- 'TO-DO' # rozsah druheho nah. vyberu

rho1 <- 'TO-DO' # korelace prvniho nah. vyberu
rho2 <- 'TO-DO' # korelace pdruheho nah. vyberu

(tab <- data.frame(n1, rho1 = rho1, n2, rho2 = rho2))

## test predpokladu dvourozmerne normality
'TO-DO' # test pro muze
'TO-DO' # test pro zeny

# graficke posouzeni dvourozmerne normality
alpha <- 'TO-DO'

# muzi
Z <- MASS::kde2d(lowex.LM, tru.LM, n = 50, lims = c(860, 1180, 350, 580))
par(mar = c(0, 2, 0, 2))
GA::persp3D(Z$x, Z$y, Z$z, ticktype = 'simple', xlab = 'delka dolni koncetiny - muzi',
            ylab = 'delka trupu', zlab = 'relativni cetnost', 
            col.palette = hcl.colors, border = 'black', theta = 15, phi = 30)

par(mar = c(3, 4, 2, 1))
car::dataEllipse(lowex.LM, tru.LM, level = 1 - alpha, xlim = c(880, 1130), 
                 ylim = c(380, 560), xlab = '', ylab = 'delka trupu (v mm)',
                 main = '', pch = 21,  col = 'darkblue', bg = 'yellow', lwd = 2, las = 1)
mtext('delka dolni koncetiny - muzi (v mm)', side = 1, line = 2.1)

# zeny
'TO-DO'


## Test o rozdilu korelacnich koeficientu

alpha <- 'TO-DO'
ksi0  <- 0
r1 <- 'TO-DO' # korelace muzi
r2 <- 'TO-DO' # korelace zeny
zR1 <- 1 / 2 * log((1 + r1) / (1 - r1)) 
zR2 <- 1 / 2 * log((1 + r2) / (1 - r2)) 
zw <- (zR1 - zR2 - ksi0) / sqrt(1 / (n1 - 3) + 1 / (n2 - 3)) # test. statistika
q <- qnorm(alpha) # kriticka hodnota

(tab <- data.frame(zR1, zR2, zw, q))

(hh <- tanh(zR1 - zR2 - sqrt(1 / (n1 - 3) + 1 / (n2 - 3)) * qnorm(alpha))) # interval spolehlivosti

p.hodnota <- pnorm(zw) # p-hodnota


## graficka vizualizace vysledku
par(mar = c(4, 4, 1, 1))
plot(lowex.LM, tru.LM, xlab = '', ylab = 'delka trupu (v mm)', 
     xlim = c(820, 1100), ylim = c(350, 550), las = 1,
     pch = 21, col = 'blue', bg = 'cornflowerblue')
k <- lm(tru.LM ~ lowex.LM)$coef
x <- seq(min(lowex.LM), max(lowex.LM), length = 1000)
y <- k[1] + x * k[2]
lines(x, y, col = 'darkblue', lwd = 2)

points(lowex.LF, tru.LF, pch = 21, col = 'red', bg = 'salmon')
k <- lm(tru.LF ~ lowex.LF)$coef
x <- seq(min(lowex.LF), max(lowex.LF), length = 1000)
y <- k[1] + x * k[2]
lines(x, y, col = 'darkred', lwd = 2)

r1 <- round(r1, digit = 4)
r2 <- round(r2, digit = 4)
mtext('delka dolni koncetiny (v mm)', side = 1, line = 2.3)
legend('topleft', pch = 21, pt.bg = c('cornflowerblue', 'salmon'), 
       col = c('blue', 'red'), legend = c('muzi', 'zeny'), 
       bty = 'n')

