packs <- c("haven","ggplot2", "gridExtra", "grid", "gridtext", "latex2exp")
packs <- packs[!packs %in% installed.packages()]
if(length(packs) != 0) install.packages(packs)
rm(packs)

library(haven)
library(ggplot2)
library(gridExtra)
library(grid)
library(gridtext)
library(latex2exp)

q1D <- function(x, quant, item, model){
  mu <- (model$b[item,1,1] + model$b[item,2,1]*x)
  sg <- exp(model$b[item,1,2] + model$b[item,2,2]*x)
  nu <- plogis(model$b[item,1,3] + model$b[item,2,3]*x)
  
  p2nuc <- function(p){
    g1max <- sqrt(2)*(4-pi)/(pi-2)^(3/2)
    nuc <- (p - g1max/(2*g1max))*(2*g1max)
    nuc
  }
  mc2md <- function(muc,sgc,nuc){
    b <- sqrt(2/pi)
    nud <- nc2nd(nuc)
    d <- nud/sqrt(1+nud^2)
    sgd <- sc2sd(sgc,nuc)
    mud <- muc - b*sgd*d
    mud
  }
  sc2sd <- function(sgc,nuc){
    b <- sqrt(2/pi)
    nud <- nc2nd(nuc)
    d <- nud/sqrt(1+nud^2)
    sgd <- sgc/sqrt(1-b^2*d^2)
    sgd
  }
  nc2nd <- function(nuc){
    b <- sqrt(2/pi)
    R <- ((2*abs(nuc))/(4-pi))^(1/3) * sign(nuc)
    nud <- R/sqrt(b^2-(1-b^2)*R^2)
    nud 
  }
  
  nuc = p2nuc(nu)
  nud = nc2nd(nuc = nuc)
  sgd = sc2sd(sgc = sg, nuc = nuc)
  mud = mc2md(muc = mu, sgc = sg, nuc = nuc)
  
  sapply(1:length(x), function(i) sn::qsn(quant, xi = mud[i], omega = sgd[i], alpha = nud[i]))
}

text_high <- textGrob("(Faster)", gp=gpar(fontsize=10, fontface="italic"))
text_low <- textGrob("(Slower)", gp=gpar(fontsize=10, fontface="italic"))
ylablim <- c(-1.65,-1.7, -1.65,-2.45,-2.75,-3,-1.85,-2.15,-1.85)

gg_Dist <- function(item, mod, data){
  texttitle <- paste0("Item ",item, " $(log(t_",item,") )$")
  p1 <- ggplot(data, aes(y = get(paste0("T",item)))) + 
    # ylim(min(data[,paste0("T",item)]), max(data[,paste0("T",item)])) +
    scale_x_continuous(limit = c(-3,3), breaks = round(seq(-3,3,by = 1),1)) +
    theme_classic() +
    xlab(TeX(r"(Latent speed factor ($z_2$))")) + 
    ylab(TeX(r"(Response time (log-minutes))")) +
    ggtitle(TeX(texttitle)) +
    geom_function(fun = ~ mod$b[paste0("T",item),1,1] + mod$b[paste0("T",item),2,1]*.x, colour = "black", linewidth = 1) +
    geom_function(fun = q1D, args = list(quant = 0.025, item = paste0("T",item), model = mod), colour = "black", linetype="dotted", linewidth = 0.7) +
    geom_function(fun = q1D, args = list(quant = 0.10, item = paste0("T",item), model = mod), colour = "black", linetype="dotted", linewidth = 0.7) +
    geom_function(fun = q1D, args = list(quant = 0.25, item = paste0("T",item), model = mod), colour = "black", linetype="dotted", linewidth = 0.7) + 
    geom_function(fun = q1D, args = list(quant = 0.50, item = paste0("T",item), model = mod), colour = "black", linetype="longdash", linewidth = 0.7) + 
    geom_function(fun = q1D, args = list(quant = 0.75, item = paste0("T",item), model = mod), colour = "black", linetype="dotted", linewidth = 0.7) + 
    geom_function(fun = q1D, args = list(quant = 0.90, item = paste0("T",item), model = mod), colour = "black", linetype="dotted", linewidth = 0.7) + 
    geom_function(fun = q1D, args = list(quant = 0.975, item = paste0("T",item), model = mod), colour = "black", linetype="dotted", linewidth = 0.7) + 
    theme(plot.title = element_text(margin = margin(t = 0, r = 0, b = 10, l = 0), size = 15),
          axis.title.y = element_text(margin = margin(t = 0, r = 15, b = 0, l = 0), size = 13),
          axis.title.x = element_text(margin = margin(t = 17, r = 0, b = 0, l = 0), size = 13),
          axis.title =   element_text(margin = margin(t = 0, r = 0, b = 10, l = 0)),
          axis.text.x = element_text(size = 12),
          axis.text.y = element_text(size = 12),
          plot.margin = margin(t = 15, r = 15, b = 15, l = 15)) + 
    annotation_custom(text_high,xmin=3, xmax=3, ymin= ylablim[item], ymax=ylablim[item]) +
    annotation_custom(text_low,xmin=-3, xmax=-3,ymin= ylablim[item], ymax= ylablim[item]) +
    coord_cartesian(clip="off")
  p1
}

gg_Prob <- function(item, mod, data){
  texttitle <- paste0("Item ",item, " $(y_",item,")$")
  axistitle <- paste0("$P(y_",item," = 1 \\,|\\, z_1)$")
  p1 <- ggplot(data, aes(y = get(paste0("Y",item)))) + 
    ylim(0,1) +
    scale_x_continuous(limit = c(-3,3), breaks = round(seq(-3,3,by = 1),1)) +
    theme_classic() +
    xlab(TeX(r"(Latent ability factor ($z_1$))")) + 
    ylab(TeX(axistitle)) +
    ggtitle(TeX(texttitle)) +
    geom_function(fun = ~ plogis(mod$b[paste0("Y",item),1,1] + mod$b[paste0("Y",item),3,1]*.x), colour = "black", linewidth = 1) +
    theme(plot.title = element_text(margin = margin(t = 0, r = 0, b = 10, l = 0), size = 15),
          axis.title.y = element_text(margin = margin(t = 0, r = 15, b = 0, l = 0), size = 13),
          axis.title.x = element_text(margin = margin(t = 17, r = 0, b = 0, l = 0), size = 13),
          axis.title =   element_text(margin = margin(t = 0, r = 0, b = 10, l = 0)),
          axis.text.x = element_text(size = 12),
          axis.text.y = element_text(size = 12),
          plot.margin = margin(t = 15, r = 15, b = 15, l = 15))
  p1
}
