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) }