###----------------------------------------------###
### (c) 2016 | Hynek Cígler                      ###
### Masaryk University, Department of Psychology ###
### Simulation of the polychoric correlation...  ###
### ... bias caused by differently skewed items  ###
###------ Thu Dec 06 11:08:12 2018 --------------###

library(mvtnorm)
library(polycor)
library(parallel)
options("mc.cores"=detectCores()-1)

n <- 150
rho <- seq(0,1,length=500)
pc1 <- pc2 <- tc <- numeric(500)

for (i in 1:500) {
  data <- rmvnorm(n, c(0, 0), matrix(c(1, rho[i], rho[i], 1), 2, 2))
  x <- data[,1]; y <- data[,2]
  xb <- ifelse(x>=1, 1, 0); yb <- ifelse(y>=-1, 1, 0) ## číslo hned za >= jsou cut-offy (tau) pro korelaci; šikmost. Je-li šikmost opačná, u silných vztahů je polychorická korelace výrazně nadhodnocená (blízká 1)
  pc1[i] <- cor(x, y)
  pc2[i] <- cor(xb, yb)
  tc[i] <- polychor(xb, yb)
}


# Graf s korelací ---------------------------------------------------------

plot(pc1, pc2, cex=.6, col="red", xlab="True linear relationship",
     ylab="Observed correlation", xlim=c(0,1), ylim=c(0,1))
lines(lowess(pc1, pc2), col="red", lwd=2)
points(pc1, tc, cex=.6, col="blue")
lines(lowess(pc1, tc), col="blue", lwd=2)
legend("topleft", c("Pearson (0/1)","Tetrachoric"), col=c(2,4), lty=1, bty="n")
abline(0, 1, col="lightgray")
abline(1, 0, col="lightgray")


# Graf s vysvětleným rozptylem --------------------------------------------

plot(pc1, pc2**2, cex=.6, col="red", xlab="True correlation (r)",
     ylab="Observed shared variance (r^2)", xlim=c(0,1), ylim=c(0,1))
lines(lowess(pc1, pc2**2), col="red", lwd=2)
points(pc1, tc**2, cex=.6, col="blue")
lines(lowess(pc1, tc**2), col="blue", lwd=2)
legend("topleft", c("Pearson (0/1)","Tetrachoric","true relation"), col=c(2,4,"lightgray"), lty=1, bty="n", lwd=2,
       pch=c(1,1,NA))
lines(rho, rho**2, col="lightgray", lwd=2)
abline(1, 0, col="lightgray")


hist(tc) ## histogram tetrachorických korelací, mezi 0-1 by měl být monotónní
