variacni.rada <- function(X, row.names){
  X  <- as.numeric(X)
  U  <- sort(unique(X))
  nu <- length(U)
  nj <- rep(0, nu)
  for(j in 1:nu){
    nj[j] <- sum(X == U[j])
  }
  
  n  <- sum(nj)
  pj <- nj / n
  Nj <- cumsum(nj)
  Fj <- cumsum(pj)
  
  variacni.rada <- data.frame(nj = nj, pj = pj, Nj = Nj, Fj = Fj)
  row.names(variacni.rada) <- row.names
  variacni.rada
}

rel.barplot <- function(N, col = 1:length(N), border = 'black', names = 1:length(N), main = '', xlab = '', 
                        ylab = 'relativni cetnost', xlim = c(0.2, 2), ylim = c(-0.03, 1.03), 
                        density = 60, cex = 1, mtext = '', a = 0, axes = axes){
  n       <- sum(N)
  l       <- length(N)
  
  barplot(matrix(N / n, l, 1), col = col, border = border, density = density,
          main = main, xlim = xlim, ylim = ylim, ylab = ylab, axes = T, las = 1)
  legend('topright',legend = rev(names), fill = rev(col), density = density, bty = 'n')
  mtext(xlab, 1, line = 1)
  
  stred <- 0.7
  cn    <- cumsum(N) / n
  cn2   <- (N / n) / 2
  vyska <- c(cn2[1], cn[1:(length(cn) - 1)] + cn2[2:length(cn2)])
  vyska <- vyska + a
  text(stred, vyska, paste(N, '; ',round(N / n * 100, 2),'%',sep = ''), cex = cex)
}


dotplot <- function(x, y, xlim = c(min(x), max(x)), ylim = c(min(y), max(y)), 
                    col = 'black', pch = 21, bg = 'white', lwd = 1, cex = 1,
                    xlab = 'x', ylab = 'y', main = '', las = 1){
  rand <- rnorm(length(x), 0, 0.03)
  x2 <- x + rand
  plot(x2, y, type = 'p', xlim = xlim, ylim = ylim, pch = pch, col = col, bg = bg, lwd = lwd,
       main = main, xlab = xlab, ylab = ylab, las = las)
}

norm2 <- function(x, y, mu1, mu2, sigma1, sigma2){
  rho <- 0
  Sigma <- matrix(c(sigma1 ^ 2, sigma2 * sigma1 * rho, sigma1 * sigma2 * rho, sigma2 ^ 2), 2, 2, byrow = T)
  xy <- c(x[1] - mu1, y[1] - mu2)
  konstanta <- 1 / (2 * pi * sqrt(sigma1 * sigma2 * (1 - rho ^ 2)))
  hustota <- konstanta * exp(- 1 / 2 * t(xy) %*% solve(Sigma) %*% xy)
  return(hustota)
}



rel.barplot.two <- function(n, col = 'grey40', col.number = 'white', density = NULL, border = 'black',
                            xlab = '', ylab = 'Relative Frequency', main = '', names = rep('', dim(t(t(n)))[2]), 
                            legend = 1:dim(t(t(n)))[1], cex.main = 1.2, cex = 1, las = 1){
  
  n <- t(t(n))
  s <- apply(n, 2, sum)
  r <- t(t(n) / s)
  d <- dim(n)[1]
  b <- dim(n)[2]
  
  n <- n[d:1, ]
  s <- s[b:1]
  r <- r[d:1, ]
  
  v <- (lower.tri(matrix(1, d, d)) + 1 / 2 * diag(d)) %*% r
  
  x <- rep(0.5:(b - 0.5), rep(d, b), main = '', xlab = xlab)
  
  barplot(t(t(r)), width = 1, space = 0, density = density, col = col, border = border,
          ylim = c(0, 1), xlim = c(0, b + 1), xlab = '', ylab = ylab, names = names, main = '', 
          cex.names = cex, axes = F )
  legend('topright', fill = rev(col), density = density, legend = legend, bty = 'n', cex = cex)
  text(x = x, y = v, paste(round(r, 4) * 100, '%'), col = col.number, cex = cex)
  axis(2, las = T)
  mtext(main, side = 3, font = 2, line = 1.7, cex = cex.main)
  mtext(xlab, side = 1, font = 1, line = 2.5, cex = cex)
}

corZ.test <- function(X, Y, rho0 = 0, conf.level = 0.95, alternative){
  n <- length(X)
  alpha <- 1 - conf.int
  r <- cor(X, Y, method = 'pearson') 
  z <- 1 / 2 * log((1 + r) / (1 - r))       
  ksi0 <- 1 / 2 * log((1 + rho0) / (1 - rho0)) 
  t0   <- (z - ksi0) * sqrt(n - 3) 
  
  if (alternative == 'two.sided'){
    dh <- tanh(z - qnorm(1 - alpha / 2) / sqrt(n - 3)) 
    hh <- tanh(z - qnorm(alpha / 2) / sqrt(n - 3)) 
    p.hodnota <- 2 * min(pnorm(t0), 1 - pnorm(t0))}
  
  if (alternative == 'greater'){
    dh <- tanh(z - qnorm(1 - alpha) / sqrt(n - 3)) 
    hh <- 1
    p.hodnota <- 1 - pnorm(t0)
  }
  
  
  if (alternative == 'less'){
    dh <- -1
    hh <- tanh(z - qnorm(alpha) / sqrt(n - 3)) 
    p.hodnota <- pnorm(t0)
  }
  
  return(data.frame(rho = r, rho0 = rho0, zW = t0, dh = dh, hh = hh, p.val = p.hodnota))
}

corZ.two.test <- function(data1, data2, conf.level = 0.95, alternative = 'two.sided'){
  alpha <- 1 - conf.level
  X1 <- data1[, 1]
  Y1 <- data1[, 2]
  X2 <- data2[, 1]
  Y2 <- data2[, 2]
  n1 <- length(X1) # 25
  n2 <- length(X2) # 24
  r1 <- cor(X1, Y1, method = 'pearson') 
  r2 <- cor(X2, Y2, method = 'pearson') 
  z1 <- 1 / 2 * log((1 + r1) / (1 - r1)) 
  z2 <- 1 / 2 * log((1 + r2) / (1 - r2)) 
  u0 <- (z1 - z2) / sqrt(1 / (n1 - 3) + 1 / (n2 - 3)) 
  sg <- sqrt(1 / (n1 - 3) + 1 / (n2 - 3))
  
  if (alternative == 'two.sided'){
    dh <- tanh(z1 - z2 - sg * qnorm(1 - alpha / 2))
    hh <- tanh(z1 - z2 - sg * qnorm(alpha / 2))
    p.hodnota <- 2 * min(pnorm(u0), 1 - pnorm(u0))}
  
  if (alternative == 'greater'){
    dh <- tanh(z1 - z2 - sg * qnorm(1 - alpha))
    hh <- 2
    p.hodnota <- 1 - pnorm(u0)
  }
  
  if (alternative == 'less'){
    dh <- -2
    hh <- tanh(z1 - z2 - sg * qnorm(alpha))
    p.hodnota <- pnorm(u0)
  }
  
  return(data.frame(R1 = r1, R2 = r2, u0 = u0, dh = dh, hh = hh, p.val = p.hodnota))
}
7

cor.plot <- function(data1, data2, col = c('blue', 'red'), bg = c('cornflowerblue', 'salmon'), xlab = '',
                     line.col = c('darkblue', 'darkred'), lwd = c(2, 2),
                     xlim = c(min(X1, X2) - 10, max(X1, X2) + 10),
                     ylim = c(min(Y1, Y2) - 10, max(Y1, Y2) + 10)) {
  X1 <- data1[ , 1]
  Y1 <- data1[ , 2]
  X2 <- data2[ , 1]
  Y2 <- data2[ , 2] 
plot(X1, Y1, xlab = '', ylab = 'delka trupu (v mm)', 
     xlim = xlim, ylim = ylim, las = 1,
     pch = 21, col = col[1], bg = bg[1])
points(X2, Y2, pch = 21, col = col[2], bg = bg[2])

k <- lm(Y1 ~ X1)$coef
x1 <- seq(min(X1), max(X1), length = 1000)
y1 <- k[1] + x1 * k[2]
lines(x1, y1, col = line.col[1], lwd = lwd[1])

k <- lm(Y2 ~ X2)$coef
x2 <- seq(min(X2), max(X2), length = 1000)
y2 <- k[1] + x2 * k[2]
lines(x2, y2, col = line.col[2], lwd = lwd[2])
}

Spearman.test <- function(X, Y, alternative = 'two.sided', conf.level = 0.95, exact = T){
  alpha <- 1 - conf.level
  n <- length(X)
  sE <- cor(X, Y, method = 'spearman')
  sA <- sE * sqrt(n - 1)
  
  if (alternative == 'two.sided'){
    p.ex <- 2 * min(SuppDists::pSpearman(sE, n), 1 - SuppDists::pSpearman(sE, n))
    p.as <- 2 * min(pnorm(sA), 1 - pnorm(sA))}
  
  if (alternative == 'greater'){
    p.ex <- 1 - SuppDists::pSpearman(sE, n)
    p.as <- 1 - pnorm(sA)
  }
  
  if (alternative == 'less'){
    p.ex <- SuppDists::pSpearman(sE, n)
    p.as <- pnorm(sA)
  }
  if(exact == T) tab <- data.frame(rS = sE, sE, p.value = p.ex)
  if(exact == F) tab <- data.frame(rS = sE, sA, p.value = p.as)
  return(tab)
}

odds.ratio.test <- function(x, conf.level = 0.95){
  a <- x[1, 1]
  b <- x[1, 2]
  c <- x[2, 1]
  d <- x[2, 2]
  alpha <- 1 - conf.level
  OR <- (a / b) / (c / d)
  lnOR <- log(OR)
  s <- sqrt(1 / a + 1 / b + 1 / c + 1 / d)
  t0 <- lnOR / s
  dh <- lnOR - s * qnorm(1 - alpha / 2)
  hh <- lnOR + s * qnorm(1 - alpha / 2)
  df.exp <- exp(dh)
  hh.exp <- exp(hh)
  p <- 2 * min(pnorm(t0), 1 - pnorm(t0))
  tab <- data.frame(OR, lnOR, t0, dh, hh, p)
  return(tab)
}