### 10 NEPARAMETRICKE TESTY
# 10.1 (Wilcoxonuv jednovyberovy test)
data <- read.delim('15-anova-means-skull.txt')
upface.HN <- data[data$pop == 'nem', 'upface.H']
upface.HN <- na.omit(upface.HN)
n <- length(upface.HN) # 19
(tab <- data.frame(n = n, min = min(upface.HN), max = max(upface.HN))) # 19,  62-76


## test predpokladu normality
shapiro.test(upface.HN)$p.val # 0.04190113


## Miaove test symetrie
# install.packages("lawstat")
lawstat::symmetry.test(upface.HN, boot = F, option = 'MGG')$p.val # 0.06815532


## graficke posouzeni predpokladu
par(mar = c(4, 4, 1, 1))
hist(upface.HN, prob = T, breaks = seq(61, 76, by = 3), 
     density = 40, col = 'yellow', border = 'orange3',
     xlab = '',
     ylab = 'relativni cetnost', main = '', axes = F)
box(bty = 'o')
axis(1, seq(62.5, 74.5, by = 3))
axis(2, las = 1)
xfit <- seq(55, 80, length = 512)
yfit <- dnorm(xfit, mean(upface.HN), sd(upface.HN))
lines(xfit, yfit, lwd = 2, col = 'orange4')
mtext('vyska horni casti tvare (v mm)', side = 1, line = 2.1)

qqnorm(upface.HN, main = '', xlab = '', pch = 21, col = 'orange4', bg = 'yellow',
       ylab = 'vyberovy kvantil', las = 1)
mtext('teoreticky kvantil', side = 1, line = 2.1)
qqline(upface.HN, lwd = 2, col = 'orange4')

hist(upface.HN, prob = T, breaks = seq(61, 76, by = 3), main = '', axes = F,
     density = 40, col = 'yellow', border = 'orange3',
     xlab = '', ylab = 'relativni cetnost')
box(bty = 'o')
axis(1, seq(62.5, 74.5, by = 3))
axis(2, las = 1)
lines(density(upface.HN), lwd = 2, col = 'orange4') 
abline(v = median(upface.HN, type = 2), lty = 2, col = 'red', lwd = 3)

mtext('vyska horni casti tvare (v mm)', side = 1, line = 2.1)
as.p <- lawstat::symmetry.test(upface.HN, boot = F)$p.val
legend('topleft', col = 'red', lty = 2, lwd = 2, legend = 'median', bty = 'n')


## Wilcoxonuv jednovyberovy test
x0 <- 70.00
m <- sum(upface.HN - x0 != 0) # 17
alpha <- 0.05
wilcox.test(upface.HN, mu = x0, conf.int = T, conf.level = 1 - alpha, 
            alt = 'greater', correct = F)
q <- qsignrank(1 - alpha, m) # 111
(tab <- data.frame(q)) 


## graficka vizualizace vysledku
par(mar = c(1, 4, 1, 1))
boxplot(upface.HN, col = rgb(1, 1, 0, 0.3), ylim = c(62, 85),
        xlab = '', ylab = 'vyska horni casti tvare (v mm)', las = 1, 
        medcol = 'orange3', border = 'orange4')
points(70.00, bg = 'red', pch = 21, cex = 1.3, col = 'darkred')
legend('topleft', col = c('chartreuse4', 'red'), bty = 'n', pch = c(NA, 20), lty = c(1, NA), lwd = 2,
       legend = c('nemecka p.', 'cernjachovska k.'))



# 10.2 (Znamenkovy jednovyberovy test)
data <- read.delim('21-goldman-shells.csv', sep = ';')
tibia.LRF <- data[data$pop == 'yos' & data$sex == 'f', 'tibia.LR']
tibia.LRF <- na.omit(tibia.LRF)
n <- length(tibia.LRF) # 8
(tab <- data.frame(n = n, min = min(tibia.LRF), max = max(tibia.LRF))) # 8, 299-331.5


## test predpokladu normality
shapiro.test(tibia.LRF)$p.val # 0.03395534


## Miaove test symetrie
lawstat::symmetry.test(tibia.LRF, boot = F, option = 'MGG')$p.val # 0.03567783 


## graficke posouzeni predpokladu
par(mar = c(4, 4, 1, 1))
hist(tibia.LRF, prob = T, breaks = seq(297, 333, by = 9), main = '', axes = F,
     density = 40, col = 'coral', ylim = c(0, 0.06),
     xlab = '',
     ylab = 'relativni cetnost')
box(bty = 'o')
axis(1, seq(301.5, 328.5, by = 9))
axis(2, las = 1)
xfit <- seq(280, 340, length = 512)
yfit <- dnorm(xfit, mean(tibia.LRF), sd(tibia.LRF))
lines(xfit, yfit, lwd = 2, col = 'firebrick')
mtext('delka lytkove kosti (v mm)', side = 1, line = 2.1)

qqnorm(tibia.LRF, main = '', xlab = '', pch = 21, col = 'firebrick3', bg = 'bisque',
       ylab = 'vyberovy kvantil', las = 1)
mtext('teoreticky kvantil', side = 1, line = 2.1)
qqline(tibia.LRF, lwd = 2, col = 'firebrick3')

hist(tibia.LRF, prob = T, breaks = seq(297, 333, by = 9), main = '', axes = F,
     density = 35, col = 'coral', ylim = c(0, 0.065), xlab = '',
     ylab = 'relativni cetnost')
box(bty = 'o')
axis(1, seq(301.5, 328.5, by = 9))
axis(2, las = 1)
lines(density(tibia.LRF), lwd = 2, col = 'firebrick3') # nesymetricke
abline(v = median(tibia.LRF, type = 2), lty = 2, col = 'blue', lwd = 2)

mtext('delka lytkove kosti (v mm)', side = 1, line = 2)
as.p <- lawstat::symmetry.test(tibia.LRF, boot = F, option = 'MGG')$p.val
legend('topleft', col = 'blue', lty = 2, lwd = 2, legend = 'median', bty = 'n')


## znamenkovy jednovyberovy test
x0 <- 329.40
m <- sum(tibia.LRF - x0 != 0) # 8
alpha <- 0.05
# install.packages("BSDA")
BSDA::SIGN.test(tibia.LRF, md = x0, conf.level = 1 - alpha, 
                alt = 'two.sided') # interpolovany IS

q1 <- qbinom(alpha / 2, m, 1 / 2) - 1 # 0
q2 <- qbinom(1 - alpha / 2, m, 1 / 2) # 7
(tab <- data.frame(q1, q2)) # 1, 0, 7


## graficka vizualizace vysledku
par(mar = c(1, 4, 1, 1))
boxplot(tibia.LRF, col = 'bisque', medcol = 'firebrick4', 
        border = 'firebrick3', ylim = c(298, 350),
        xlab = '', ylab = 'delka lytkove kosti (v mm)', las = 1)
points(329.40, bg = 'red', pch = 21, cex = 1.3, col = 'black')
legend('topright', col = c('firebrick4', 'black'), pt.bg = c(NA, 'red'), 
       bty = 'n', pch = c(NA, 21), 
       lty = c(1, NA), lwd = c(2, NA),
       legend = c('neoliticka p.', 'soucasna p.'))



# 10.3 (Wilcoxonuv parovy test)
data <- read.delim('21-goldman-tigara.csv', sep = ';')
data.F <- data[data$pop == 'Ipituaq' & data$sex == 'f', c('femur.HDR', 'femur.HDL')]
data.F <- na.omit(data.F)
femur.HDL <- data.F$femur.HDL
femur.HDR <- data.F$femur.HDR
diff <- femur.HDL - femur.HDR
n <- length(diff)
(tab <- data.frame(n = n, min.L = min(femur.HDL), max.L = max(femur.HDL),
                   min.R = min(femur.HDR), max.R = max(femur.HDR))) # 12,  38-44.02 37.85-44.09


## test predpokladu normality
shapiro.test(diff)$p.val # 0.02463559


## Miaove test symetrie
lawstat::symmetry.test(diff, boot = F, option = 'MGG')$p.val # 0.475678 


## graficke posouzeni predpokladu
par(mar = c(4, 4, 1, 1))
hist(diff, prob = T, breaks = seq(-2.9, 1.3, by = 0.84), main = '', axes = F,
     density = 40, col = 'deepskyblue', ylim = c(0, 0.45),
     xlab = '', ylab = 'relativni cetnost')
box(bty = 'o')
axis(1, seq(-2.48, 0.88, by = 0.84))
axis(2, las = 1)
xfit <- seq(-3.5, 2, length = 512)
yfit <- dnorm(xfit, mean(diff), sd(diff))
lines(xfit, yfit, lwd = 2, col = 'deepskyblue4')
mtext('rozdil prumeru L - P (v mm)', side = 1, line = 2.1)

qqnorm(diff, main = '', xlab = '', pch = 21, col = 'deepskyblue3', bg = 'aliceblue',
       ylab = 'vyberovy kvantil', las = 1)
mtext('teoreticky kvantil', side = 1, line = 2.1)
qqline(diff, lwd = 2, col = 'deepskyblue4')

hist(diff, prob = T, breaks = seq(-2.9, 1.3, by = 0.84), main = '', axes = F,
     density = 40, col = 'deepskyblue', ylim = c(0, 0.55),
     xlab = '', ylab = 'relativni cetnost')
box(bty = 'o')
axis(1, seq(-2.48, 0.88, by = 0.84))
axis(2, las = 1)

lines(density(diff), lwd = 2, col = 'deepskyblue4') # nesymetricke
abline(v = median(diff, type = 2), lty = 2, col = 'red', lwd = 2)
mtext('rozdil prumeru L - P (v mm)', side = 1, line = 2.1)
legend('topleft', col = 'red', lty = 2, lwd = 2, legend = 'median', bty = 'n')


## Wilcoxonuv parovy test
z0 <- 0
m <- sum(diff - z0 != 0)
alpha <- 0.10
wilcox.test(diff, mu = z0, conf.int = T, conf.level = 1 - alpha, 
            alternative = 'less', correct = F)
wilcox.test(femur.HDL, femur.HDR, paired = T, conf.int = T, conf.level = 1 - alpha, 
            alternative = 'less', correct = F)
q <- qsignrank(alpha , m) - 1 # 17
(tab <- data.frame(q))


## graficka vizualizace vysledku
par(mar = c(1, 4, 1, 1))
boxplot(diff, col = 'lightcyan1', ylim = c(-5, 1.5),
        xlab = '', ylab = 'rozdily a-p. prumeru (v mm)', las = 1, 
        medcol = 'royalblue4', border = 'deepskyblue4')
segments(0.8,0,1.2,0, lty = 2, col = 'red', lwd = 3)
legend('bottomleft', col = c('royalblue4', 'red'), bty = 'n',
       lty = c(1, 2), lwd = c(2, 2), 
       legend = c('namerene rozdily', expression(paste(z[0] == 0))))

par(mar = c(2, 4, 1, 1))
boxplot(femur.HDL, femur.HDR, col = c('lightcyan1', 'mintcream'), 
        names = c('prava strana', 'leva strana'),
        xlab = '', ylab = 'a-p prumer (v mm)', las = 1, 
        medcol = c('deepskyblue4', 'deepskyblue4'), border = c('deepskyblue3', 'deepskyblue3'))



# 10.4 (Wilcoxonuv dvouvyberovy test)
data <- read.delim('19-more-samples-correlations-skull.txt')
intorb.BB <- data[data$pop == 'ban', 'intorb.B']
intorb.BB <- na.omit(intorb.BB)
intorb.BP <- data[data$pop == 'per', 'intorb.B']
intorb.BP <- na.omit(intorb.BP)
n1 <- length(intorb.BB) # 14
n2 <- length(intorb.BP) # 46
(tab <- data.frame(n1, n2))#, min1 = min(intorb.BB), max1 = max(intorb.BB),
#min2 = min(intorb.BP), max2 = max(intorb.BP))) # 14, 46, 21-30, 19-28


## test predpokladu normality
shapiro.test(intorb.BB)$p.val # 0.4537105
shapiro.test(intorb.BP)$p.val # 0.08481863


## graficke posouzeni predpokladu
# bantuska
par(mar = c(3, 4, 0, 1))
hist(intorb.BB, prob = T, breaks = seq(20, 30, by = 2), main = '', axes = F,
     col = 'green3', density = 40, ylim = c(0, 0.20),
     xlab = '', ylab = 'relativni cetnost')
box(bty = 'o')
axis(1, seq(21, 29, by = 2))
axis(2, las = 1)
xfit <- seq(15, 35, length = 512)
yfit <- dnorm(xfit, mean(intorb.BB), sd(intorb.BB))
lines(xfit, yfit, lwd = 2, col = 'green4') 
mtext('interorb. sirka - bantuska p.', side = 1, line = 2.1)

qqnorm(intorb.BB, main = '', xlab = '', pch = 21, col = 'green4', bg = 'green1',
       ylab = 'vyberovy kvantil', las = 1)
mtext('teoreticky kvantil', side = 1, line = 2.1)
qqline(intorb.BB, lwd = 2, col = 'green4')

# peruanska
hist(intorb.BP, prob = T, breaks = seq(17, 29, by = 2), main = '', axes = F,
     col = 'orange2', density = 40, ylim = c(0, 0.25), 
     xlab = '', ylab = 'relativni cetnost')
box(bty = 'o')
axis(1, seq(18, 28, by = 2))
axis(2, las = 1)
xfit <- seq(15, 35, length = 512)
yfit <- dnorm(xfit, mean(intorb.BP), sd(intorb.BP))
lines(xfit, yfit, lwd = 2, col = 'orange3') 
mtext('interorb. sirka - peruanska p.', side = 1, line = 2.1)

qqnorm(intorb.BP, main = '', xlab = '', pch = 21, col = 'orange3', bg = 'orange1',
       ylab = 'vyberovy kvantil', las = 1)
mtext('teoreticky kvantil', side = 1, line = 2.1)
qqline(intorb.BP, lwd = 2, col = 'orange3')


## Wilcoxonuv dvouvyberovy test
alpha <- 0.10
wilcox.test(intorb.BB, intorb.BP, conf.level = 1 - alpha, alt = 'greater', 
            exact = F, conf.int = T) # interopolovany IS
q <- qwilcox(1 - alpha, n1, n2)
(tab <- data.frame(q)) # 396


## graficka vizualizace vysledku
par(mar = c(4, 4, 1, 0))
boxplot(intorb.BB, intorb.BP, names = c('bantuska', 'peruanska'), 
        col = c('green1', 'orange1'),
        xlab = '', ylab = 'interorb. sirka (v mm)', las = 1,   
        border = c('green4', 'orange4'))
mtext('populace', side = 1, line = 2)



# 10.5 (Spearmanuv poradovy exaktni test o nezavislosti)
data <- read.delim('23-3D-sken.txt', sep = '&', header = F, row.names = 1,
                   col.names = c('tym', paste('sf', 1:12, sep = '_')))
it <- as.numeric(data[1, ])
ant <- as.numeric(data[2, ])
n <- length(it) # 12

## Spearmanuv poradovy exaktni test o nezavislosti
alpha <- 0.01
source('funkce-verze-03.R')
Spearman.test(it, ant, alternative = 'greater', 
              conf.level = 1 - alpha, exact = T)
# install.packages("SuppDists")
q <- SuppDists::qSpearman(1 - alpha, n)
(tab <- data.frame(q)) # 0.6678322


## graficka vizualizace vysledku
par(mar = c(3, 4, 1, 1))
plot(it, ant, pch = 21, col = 'purple4', bg = 'plum', xlab = '', 
     ylab = 'hodnoceni antrop. tymu', las = 1)
mtext('hodnoceni IT tymu', side = 1, line = 2.1)



# 10.6 (Spearmanuv poradovy asymptoticky test o nezavislosti)
data <- read.delim('22-kralik-WHR.csv', sep = ';')
data.M <- data[data$sex == 'm' & data$age >= 10, c('age', 'WHR')]
data.M <- na.omit(data.M)
age.M <- data.M$age
WHR.M <- data.M$WHR
n <- length(age.M)
(tab <- data.frame(n))


## test predpokladu dvourozmerne normality
alpha <- 0.05
MVN::mvn(data.M, mvnTest = 'mardia')$multivariateNormality 


## vizualni posouzeni predpokladu dvourozmerne normality
Z <- MASS::kde2d(age.M, WHR.M, n = 50, lims = c(9, 17, 0.6, 1))
par(mar = c(0, 2, 0, 4))
GA::persp3D(Z$x, Z$y, Z$z, ticktype = 'simple', xlab = 'vek (v letech)',
            ylab = 'pomer obvodu pasu a boku', zlab = "relativni cetnost", 
            col.palette = terrain.colors, border = 'black', theta = 15, phi = 30)

par(mar = c(4, 4, 2, 1))
car::dataEllipse(age.M, WHR.M, level = 1 - alpha, pch = 21,  col = 'darkgreen', 
                 bg = 'bisque', lwd = 2, xlim = c(9, 17), 
                 ylim = c(0.65, 1), xlab = '', ylab = 'pomer obvodu pasu a boku',
                 main = '', las = 1)
mtext('vek (v letech)', side = 1, line = 2.1)


## Spearmanuv poradovy asymptoticky test o nezavislosti
alpha <- 0.05
source('funkce-verze-03.R')
Spearman.test(age.M, WHR.M, alternative = 'two.sided', conf.level = 1 - alpha, exact = F)
q1 <- qnorm(alpha / 2)
q2 <- qnorm(1 - alpha / 2)
(tab <- data.frame(q1, q2)) # -3.014078, -1.959964, 1.959964


## graficka vizualizace vysledku
par(mar = c(3, 4, 2, 1))
plot(age.M, WHR.M, pch = 21, col = 'darkgreen', bg = 'bisque', xlab = '', 
     ylab = 'pomer obvodu pasu a boku', las = 1)
mtext('vek (v letech)', side = 1, line = 2.1)
