# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # 
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # 
#               R code for running Empirical Application 2 in paper                   #
#   "Generalized Latent Variable Models for Location, Scale, and Shape parameters"    #
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # 

rm(list = ls())
set.seed(1234)

# Set working directory
# Save file "glvmlss_0.0.1.tar.gz" in working directory
# Uncomment below for installing the glvmlss package
# if(!"glvmlss" %in% installed.packages()) install.packages("glvmlss_0.0.1.tar.gz", type = "source", repos = NULL)

library(glvmlss)
data <- ANES2020[,1:8]
dorg <- ANES2020[,9:11]

# Uncomment to load estimation results and avoid running the code
# Save file "Blind Submission - Auxiliary files/bsub_Results_ANES2020.Rdata" in working directory
# load("bsub_Results_ANES2020.Rdata")

# Summary statistics (Table A1 in Online Supplementary Material)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
vtable::st(data, summ = c('notNA(x)','mean(x, na.rm = T)','sd(x, na.rm = T)','PerformanceAnalytics::skewness(x, na.rm = T)','PerformanceAnalytics::kurtosis(x, na.rm = T,method = "excess")'),
           summ.names = c('N','Mean','Std. Dev.','Skewness','Kurtosis'))

# Estimation
# ~~~~~~~~~~
p = ncol(data)
famt <- rep("beta",p)

mod1 <- glvmlss::glvmlss(data = data, family = famt, mu.eq = ~ Z1, sg.eq = ~ 1,
                         EM.maxit = 30, DM.maxit = 300, nQP = 100, "Fisher" = T, "verbose" = T)

mod2 <- glvmlss::glvmlss(data = data, family = famt, mu.eq = ~ Z1, sg.eq = ~ Z1,
                          EM.maxit = 30, DM.maxit = 300, nQP = 100, "Fisher" = T, "verbose" = T)

# Save results for future use:
# save.image("bsub_Results_ANES2020.Rdata")

# Table 1 (Section 3.1)
# ~~~~~~~~~~~~~~~~~~~~~
xtable::xtable(
  rbind(c(mod1$AIC, mod1$BIC, mod1$K),
        c(mod2$AIC, mod2$BIC, mod2$K)))

# Table 2 (Section 3.1)
# ~~~~~~~~~~~~~~~~~~~~~
xtable::xtable(
  cbind(mod2$b[,1,1], mod2$seb[,1,1],
        mod2$b[,2,1], mod2$seb[,2,1],
        mod2$b[,1,2], mod2$seb[,1,2],
        mod2$b[,2,2], mod2$seb[,2,2]))

# Figures 1 and 2 (Section 3.1)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Save file "Blind Submission - Auxiliary files/plotAuxE2.R" in working directory
source("plotAuxE1.R")
datatot <- as.data.frame(cbind(data,dorg))
datatot$zA1 <- mod1$Z$Z1
datatot$zA2 <- mod2$Z$Z1

cor.test(-1*datatot$zA1, datatot$LRscale, use = "pairwise.complete.obs")
cor.test(-1*datatot$zA1, datatot$LCscale, use = "pairwise.complete.obs")

cor.test(-1*datatot$zA2, datatot$LRscale, use = "pairwise.complete.obs")
cor.test(-1*datatot$zA2, datatot$LCscale, use = "pairwise.complete.obs")

text_A <- textGrob("(Higher)", gp=gpar(fontsize=10, fontface="italic"))
text_B <- textGrob("(Lower)", gp=gpar(fontsize=10, fontface="italic"))
gg_q_data(data$feminists) +
  theme_classic() +
  stat_ecdf(aes(data$metoo), geom = "step", linewidth = 0.7, na.rm = T,alpha = 0.3, colour = "gray70", pad = F) +
  stat_ecdf(aes(data$gaylesbian), geom = "step", linewidth = 1, na.rm = T,alpha = 1, colour = "black", pad = F, linetype="longdash") + 
  stat_ecdf(aes(data$transgend), geom = "step", linewidth = 0.7, na.rm = T,alpha = 0.3, colour = "gray70", pad = F) +
  stat_ecdf(aes(data$BLM), geom = "step", linewidth = 1, na.rm = T,alpha = 1, colour = "black", pad = F, linetype="dotted") +
  stat_ecdf(aes(data$journalist), geom = "step", linewidth = 0.7, na.rm = T,alpha = 0.3, colour = "gray70", pad = F) + 
  stat_ecdf(aes(data$scientists), geom = "step", linewidth = 1, na.rm = T,alpha = 1, colour = "black", pad = F, linetype="dotdash") +
  theme(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 = 15, r = 0, b = 0, l = 0), size = 13),
        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_A,xmin=1, xmax=1, ymin= -0.15, ymax= -0.15) + 
  annotation_custom(text_B,xmin=0, xmax= 0,ymin= -0.15, ymax= -0.15) +
  coord_cartesian(clip="off") # save as svg, 600 x 500

text_high <- textGrob("(Progressive)", gp=gpar(fontsize=10, fontface="italic"))
text_low <- textGrob("(Conservative)", gp=gpar(fontsize=10, fontface="italic"))

# mod <- mod1
mod <- mod2
itm <- "metoo"
# itm <- "scientists"

ggplot(datatot, aes(x = zA2, y = feminists)) + 
  xlim(-2.5,2.5) + 
  ylim(0,1) +
  theme_classic() +
  xlab(TeX(r"(Latent variable (Z))")) + 
  ylab(TeX(r"(Thermometer rating (scaled))")) + 
  geom_function(fun = ~ plogis(mod$b[itm,1,1] + mod$b[itm,2,1]*.x), colour = "black", linewidth = 1) +
  geom_function(fun = q1den, args = list(quant = 0.10, item = itm, model = mod), colour = "black", linetype="dotted", linewidth = 0.7) +
  geom_function(fun = q1den, args = list(quant = 0.25, item = itm, model = mod), colour = "black", linetype="dotted", linewidth = 0.7) + 
  geom_function(fun = q1den, args = list(quant = 0.50, item = itm, model = mod), colour = "black", linetype="longdash", linewidth = 0.7) + 
  geom_function(fun = q1den, args = list(quant = 0.75, item = itm, model = mod), colour = "black", linetype="dotted", linewidth = 0.7) + 
  geom_function(fun = q1den, args = list(quant = 0.90, item = itm, model = mod), colour = "black", linetype="dotted", linewidth = 0.7) + 
  theme(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=2, xmax=2, ymin= -0.175, ymax= -0.175) + 
  annotation_custom(text_low,xmin=-2, xmax=-2,ymin= -0.175, ymax= -0.175) +
  coord_cartesian(clip="off") # save as svg, 550 x 500

gg_qq_empirical(-1*datatot$zA2, datatot$LRscale) +
  theme_light() +
  xlim(-2.25,2.25) + 
  ylim(-2.25,2.25) + 
  xlab(TeX(r"(Empirical quantiles for Empirical Bayes factor scores)")) + 
  ylab(TeX(r"(Empirical quantiles for the  \textit{Left-Right} (standardised) scale)")) + 
  theme(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 = 15, r = 0, b = 0, l = 0), size = 13),
        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)) # save as svg, 600 x 500

gg_qq_empirical(-1*datatot$zA2, datatot$LCscale) +
  theme_light() +
  xlim(-2.25,2.25) + 
  ylim(-2.25,2.25) + 
  xlab(TeX(r"(Empirical quantiles for Empirical Bayes factor scores)")) + 
  ylab(TeX(r"(Empirical quantiles for the  \textit{Lib.-Cons.} (standardised) scale)")) + 
  theme(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 = 15, r = 0, b = 0, l = 0), size = 13),
        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)) # save as svg, 600 x 500
