rm(list = setdiff(ls(), lsf.str())) ## To erase memory

## Loading required packages
if(!require("pacman")) install.packages('pacman')
library(pacman)
p_load(tidyverse)
p_load(igraph)
p_load(qgraph)
p_load(fastDummies)
p_load(Matrix)
p_load(knitr)
p_load(rgexf)
p_load(plotly)
p_load(htmlwidgets)
p_load(visNetwork)

## 2020 Wave of the ANES as example data
anes_timeseries_2020 <- read.csv("anes_timeseries_2020.csv")

## For reproducibility
set.seed(42)

## Selecting Res-IN node variables: 
anes_nodes <- dplyr::select(anes_timeseries_2020, c("V201336", "V202257", 
                                                    "V201417", "V201312", 
                                                    "V201416", "V201262", 
                                                    "V202337", "V201258"))

## Recoding missing values (check ANES codebook!)
anes_nodes$V201336[anes_nodes$V201336 %in% c(-9, -8, -7, -6, -5, 99)] <- NA
anes_nodes$V202257[anes_nodes$V202257 %in% c(-9, -8, -7, -6, -5, 99)] <- NA
anes_nodes$V201417[anes_nodes$V201417 %in% c(-9, -8, -7, -6, -5, 99)] <- NA
anes_nodes$V201312[anes_nodes$V201312 %in% c(-9, -8, -7, -6, -5, 99)] <- NA
anes_nodes$V201416[anes_nodes$V201416 %in% c(-9, -8, -7, -6, -5, 99)] <- NA
anes_nodes$V201262[anes_nodes$V201262 %in% c(-9, -8, -7, -6, -5, 99)] <- NA
anes_nodes$V202337[anes_nodes$V202337 %in% c(-9, -8, -7, -6, -5, 99)] <- NA
anes_nodes$V201258[anes_nodes$V201258 %in% c(-9, -8, -7, -6, -5, 99)] <- NA

## Give more recognizable names
colnames(anes_nodes) <- c("abort", "income", "immigr", "welfare", "gay_mar",
                          "environm", "gun_own", "aid_black")

## Re-code response labels to make them more recognizable 
anes_nodes <- anes_nodes %>% mutate(
                  abort = recode(abort, `1` = "never",
                                        `2` = "spec_cases",
                                        `3` = "in_need",
                                        `4` = "always",
                                        `5` = "other"),
                  
  income = recode(income, `1` = "favor",`2` = "oppose", `3` = "neither"),
  immigr = recode(immigr, `1` = "send_back", `2` = "guest_prog",`3` = "neither"),
  welfare = recode(welfare, `1` = "increase",`2` = "decrease", `3` = "same",),
  gay_mar = recode(gay_mar, `1` = "recogn",`2` = "civ_union", `3` = "no_recogn"),
  environm = recode(environm, `1` = "regul+++",`2` = "regul++", `3` = "regul+", `4` = "neut",`5` = "regul-", `6` = "regul--", `7` = "regul---"),
  gun_own = recode(gun_own, `1` = "more_diffic",`2` = "easier", `3` = "same"),
  aid_black = recode(aid_black, `1` = "yes++",`2` = "yes+", `3` = "yes", `4` = "neut",`5` = "no", `6` = "no+", `7` = "no++"))

## Using the "fastDummies" package by Jacob Kaplan & Benjamin Schlegel  
df_dummies <- fastDummies::dummy_cols(anes_nodes, ignore_na = TRUE, remove_selected_columns = TRUE)

## Number of collumns in the binarized, item-response data-frame:
dim(df_dummies)[2]

## Generate a symmetric correlation matrix & apply item response lables
res_in_cor <- cor(df_dummies, method = "spearman", use = "pairwise.complete.obs")
colnames(res_in_cor) <- colnames(df_dummies)
rownames(res_in_cor) <- colnames(df_dummies)

## Looping over the item-names to zero-out correlations pertaining to the same item
j <- 1
i <- 1
while(i <= ncol(anes_nodes)) {
  res_in_cor[j:((j+length(levels(factor(anes_nodes[, i]))))-1), j:((j+length(levels(factor(anes_nodes[, i]))))-1)] <- 0
  j <- j+length(levels(factor(anes_nodes[, i])))
  i <- i+1}

## Force all negative correlations to zero; remove any NA correlations
res_in_cor[res_in_cor<0] <- 0
res_in_cor[is.na(res_in_cor)] <- 0

## Generate an igrpah object from our adjacency matrix:
res_in_graph_igraph <- igraph::graph_from_adjacency_matrix(res_in_cor,
                                                           weighted=TRUE, 
                                                           mode="undirected", 
                                                           diag=FALSE)

## Here we deliberately model the graph layout using the F-R algorithm 
res_in_layout  <- layout_with_fr(res_in_graph_igraph) %>% as.data.frame()
res_in_layout$node_names <- colnames(res_in_cor)
colnames(res_in_layout) <- c("x", "y", "node_names")

## Performing PCA with the R-native "prcomp" function
res_in_layout_pca <- prcomp(res_in_layout[,1:2])$x
res_in_layout[1:2] <- res_in_layout_pca

## qgraph package (Epskamp and Costantini et.al.)
res_in_qgraph <- qgraph::qgraph(res_in_cor, layout = "spring", maximum = 1, 
                              vsize = 6, DoNotPlot = TRUE, 
                              sampleSize = nrow(anes_nodes),
                              title = "Plot generated with qgraph",
                              mar = c(100,100,100,100), normalize = FALSE)

## Extracting edgelist to manually run F-R force directed algorithim
temp_edges <- cbind(res_in_qgraph$Edgelist$from, res_in_qgraph$Edgelist$to)

## Run FR using qgraph-supplied function
layout_fr <-  qgraph.layout.fruchtermanreingold(
  temp_edges, weights = res_in_qgraph$Edgelist$weight, vcount = nrow(temp_edges))

## Orthogonal rotation and re-scaling using PCA
layout_fr_pca <- princomp(layout_fr[1:nrow(res_in_cor), ])
layout_fr_pca$scores <- layout_fr_pca$scores/max(layout_fr_pca$scores)*10
layout_fr_pca$scores[, 1] <- -(layout_fr_pca$scores[, 1])

## Replace the orig. layout with PCA rotated & rescaled layout
res_in_qgraph$layout <- layout_fr_pca$scores 
res_in_qgraph$layout.orig <- layout_fr_pca$scores 

## First define an igraph object from the ResIn adjecency matrix
res_in_igraph <- igraph::graph_from_adjacency_matrix(
                          res_in_cor, weighted=TRUE, 
                          mode="undirected", diag=FALSE)

## Create a ggplot-friendly data-frame using the igraph F-R algorithm:
res_in_layout  <- layout_with_fr(res_in_igraph) 
res_in_layout <- princomp(res_in_layout)$scores %>% as.data.frame()

res_in_layout$node_names <- colnames(res_in_cor)
colnames(res_in_layout) <- c("x", "y", "from")

## Use igraph's basic plotting structure as scaffolding for edgelist:
g_temp <- igraph::as_data_frame(res_in_graph_igraph)

## We'll replace the igraph edgelist values with PCA-rotated solution
g_temp$from.x <- res_in_layout$x[match(g_temp$from, res_in_layout$from)]
g_temp$from.y <- res_in_layout$y[match(g_temp$from, res_in_layout$from)]
g_temp$to.x <- res_in_layout$x[match(g_temp$to, res_in_layout$from)]
g_temp$to.y <- res_in_layout$y[match(g_temp$to, res_in_layout$from)]

## Final plotting frame:
plotting_frame <- left_join(g_temp, res_in_layout, by = "from")

ANES_net_fig <- ggplot() +
    geom_curve(data = plotting_frame, aes(x = from.x, xend = to.x, y = from.y,
                yend = to.y), curvature = 0, color = "grey") +
        ## geom_curve requires the x, xend, y and yend aesthetics to be able to plot the network edges
    geom_text(data = plotting_frame, aes(x = from.x, y = from.y, label = from), size = 3) +
    geom_text(data = plotting_frame, aes(x = to.x, y = to.y, label = to), size = 3) +
        ## Make sure to plot node labels both based on the to and from columns in the plotting frame!
    ggtitle("Belief Structure among ANES 2020 Respondents")+
    theme_bw()+
    expand_limits(x = c(min(plotting_frame$x-2), max(plotting_frame$x+2)),
                  y = c(min(plotting_frame$y-2), max(plotting_frame$y+2))) +
    theme(axis.text.x = element_blank(), axis.title.x = element_blank(),
          axis.text.y = element_blank(), axis.title.y = element_blank(),
          axis.ticks = element_blank(), panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(), legend.position = "bottom",
          legend.text = element_blank(), plot.title = element_text(hjust = 0.5))

## Calculating summary statistics based on select co-variates & stats; here we use Dem. & Rep. feelings thermometers
node_covars <- c("V201157", "V201156") 

## We would like to extract the mean of each thermometer variable for the subset of respondents who selected each item response
node_costats <- c("mean", "mean") 

covars_frame <- dplyr::select(anes_timeseries_2020, node_covars)
cov_stats <- as.data.frame(matrix(NA, nrow(res_in_layout), ncol(covars_frame)))
    
for(i in 1:nrow(res_in_layout)) {
  for(j in 1:ncol(covars_frame)) {
        cov_stats[i, j] <- do.call(node_costats[j], c(list(
          x = covars_frame[, j][df_dummies[, i] == 1], na.rm = TRUE)))
  }
}

colnames(cov_stats) <- paste(node_covars, node_costats, sep = "_")
cov_stats$node_label <- colnames(res_in_cor)

## We also calculate the level of affective polarization based on the Dem & Rep feelings thermometer scores
cov_stats$aff_pola <- abs(cov_stats$V201157_mean - cov_stats$V201156_mean)
res_in_layout_stats <- cbind(res_in_layout, cov_stats)

## Here, we integrate our node-level covariate statistics back into the original plotting frame
plotting_frame <- left_join(g_temp, res_in_layout_stats, by = "from")

p_load(ggdark)

ANES_net_fig <- ggplot() +
    geom_curve(data = plotting_frame, aes(x = from.x, xend = to.x, y = from.y,
                yend = to.y, size = weight^10), curvature = 0.15, color = "darkgrey", alpha = 0.4) +
                  ## Note that we use the size aesthetic to control the thickness of the graph edges
    geom_text(data = plotting_frame, aes(x = from.x, y = from.y, label = from, color = aff_pola), size = 6) +
    geom_text(data = plotting_frame, aes(x = to.x, y = to.y, label = to, color = aff_pola), size = 6) +
                  ## Make sure to plot node labels both based on the to and from columns in the plotting frame
    ggtitle("Belief Structure among ANES 2020 Respondents")+
    expand_limits(x = c(min(plotting_frame$x-1.5), max(plotting_frame$x+1.5)),
                  y = c(min(plotting_frame$y-1.5), max(plotting_frame$y+1.5))) +
    scale_color_continuous(name = "Intensity of Partisan Affective Polarization", high = "yellow", low = "blue")+
    scale_size_continuous(guide = "none")+
    dark_theme_gray()+
    theme(axis.text.x = element_blank(), axis.title.x = element_blank(),
          axis.text.y = element_blank(), axis.title.y = element_blank(),
          axis.ticks = element_blank(), panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(), legend.position = "bottom",
          legend.text = element_blank(), plot.title = element_text(hjust = 0.5, size=18))

test_run_1 <- ResIN(anes_nodes, node_vars = c("abort", "income", "immigr",
                    "welfare", "gay_mar", "environm", "gun_own", "aid_black"), 
                   output = c("igraph"), plot_graph = TRUE,
                   plot_title = "Example plot generated with igraph package")

kable(round(test_run_1$adj_matrix[1:10, 1:6], 2), format = "html", booktabs = TRUE)

test_run_2 <- ResIN(anes_nodes, 
                    node_vars = c("abort", "income", "welfare", "gay_mar",
                                  "environm", "gun_own", "aid_black"), 
                    output = c("qgraph"), plot_graph = FALSE, same_item_groups = TRUE,
                    qgraph_arglist = list(layout = "spring", maximum = 1, vsize = 6,
                                          title = "Example plot generated with qgraph package",
                                          mar = c(2,2,2,1), DoNotPlot = TRUE, normalize = FALSE,
                                          theme = "Borkulo", shape = "diamond", pastel = TRUE,
                                          label.scale = FALSE, label.cex = 0.8, label.prop = 0.5,
                                          legend.cex = 0.6))

## Importing ESS 2018 dataset and selecting all relevant variables
ESS_9 <- read.csv("ESS_9.csv")

ESS_9_n <- ESS_9 %>% dplyr::select(cntry, lrscale, gincdif, freehms, imsmetn, imdfetn, impcntr,
                                   ipstrgv, ipeqopt, impenv, pspwght)

## Sub-setting only Austrian respondents
ESS_9_at <- ESS_9_n[ESS_9_n$cntry=="AT", ]

## Using the advanced function to generate a plotable data-frame
ResIN_gen_AT <- ResIN_gen(ESS_9_at, node_vars =
                            c("lrscale", "gincdif", "freehms", "impcntr", 
                              "ipstrgv", "ipeqopt", "impenv"), 
                          weights = "pspwght", EBICglasso = TRUE, network_stats = TRUE)

## ResIn_gen stores the ggplot-ready dataframe in an object called "ggplot_frame"
plotting_frame <- ResIN_gen_AT$ggplot_frame

## Here we order the plotting frame by the statistics we are using as a visual aesthetic.
  ## We'd like the most central nodes and edges to appear on top of the more peripheral nodes which is why we sort in ascending order.
plotting_frame <- plotting_frame[order(plotting_frame$Strength, decreasing = FALSE), ]

## Now we can plot the network using geom_curve!
ESS_net_fig <- ggplot() +
    geom_curve(data = plotting_frame, aes(x = from.x, xend = to.x, y = from.y,
                yend = to.y, size = weight^10, color = Strength), curvature = 0.2) +
                  ## Here, we specified edge size according to the correlation weight and color edges based on node strength centrality
    geom_text(data = plotting_frame, aes(x = from.x, y = from.y, label = from), size = 4.2, color = "white") +
    geom_text(data = plotting_frame, aes(x = to.x, y = to.y, label = to), size = 4.2, color = "white") +
    ggtitle("Belief System Structure among Austrian ESS 2018 Respondents")+
    theme_dark()+
    scale_color_continuous(name = "Node Strength Centrality (lowest to highest): ", high = "#132B43", low = "#56B1F7")+
    scale_size_continuous(guide = "none")+
    theme(axis.text.x = element_blank(), axis.title.x = element_blank(),
          axis.text.y = element_blank(), axis.title.y = element_blank(),
          axis.ticks = element_blank(), panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(), legend.position = "bottom",
          legend.text = element_blank())

## Very simple example employing the ResIn_gen function to ANES data: 
ResIN_gen_ANES <- ResIN_gen(anes_nodes, EBICglasso = TRUE, network_stats = TRUE)

## Here we need to do just a little bit of prep-work for the visNetwork package
    ## At the node level, visNetwork minimally needs a unique id column.
      ## Additional columns such as "label", "group", and "title" are optional but very useful for interactive plotting
ResIN_gen_ANES$outcome_frame$id <- ResIN_gen_ANES$outcome_frame$node_names
ResIN_gen_ANES$outcome_frame$label <- ResIN_gen_ANES$outcome_frame$node_names
ResIN_gen_ANES$outcome_frame$group <- ResIN_gen_ANES$same_items

    ## The "title" column allows for a flexible HTML pop-up window as part of the integration. 
      ## Here, we customize this window with network descriptive statistics:
ResIN_gen_ANES$outcome_frame$title <- paste0("<center> <strong>", ResIN_gen_ANES$outcome_frame$label, "</strong> </center> <br>", 
                                          "Strength centrality: ", round(ResIN_gen_ANES$outcome_frame$Strength, 2), "<br>",
                                          "Betweeness centrality: ", ResIN_gen_ANES$outcome_frame$Betweenness, "<br>",
                                          "Closeness centrality: ", round(ResIN_gen_ANES$outcome_frame$Closeness, 3) , "<br>",
                                          "Expected influence: ", round(ResIN_gen_ANES$outcome_frame$ExpectedInfluence, 2), "<br>")

## At the edge-level, visNetwork minimally requires a "from" and "to" vector. 
    ## Here we additionally specify the edge weights as "value" collumn.
ResIN_gen_ANES$ggplot_frame$value <- ResIN_gen_ANES$ggplot_frame$weight

## Note how we keep the visualization as sparse as possible by specifying a 0-degree nearest neighbor highlighting option:
ResIN_gen_ANES_plot <- visNetwork(ResIN_gen_ANES$outcome_frame, ResIN_gen_ANES$ggplot_frame,
                                  main = "Interactive ANES 2020 Belief System Visualization") %>% 
  visOptions(highlightNearest = list(enabled = TRUE, degree = 0, hover = TRUE)) %>%
  visNodes(size = 24) %>% visInteraction(hideEdgesOnDrag = TRUE)