#################################
# RANKING ANAYLSIS
# w/ nested marginal means
#
#
# Part of replication files for:
# Dill, Janina, Howlett, Marnie, & Müller-Crepon, Carl (2025).
# Do Ukrainians Still Prefer Self-Defense Against Russia At Any Cost? 
# Journal of Peace Research, forthcoming.
#
#############################

# library(devtools)
# devtools::install_github(repo = "carl-mc/cjRank")

# LOCALS ########################
library(cjRank)

## Categorical attribute variables
attr.vars <- paste0("attr", 1:5, "_level")

## Binary outcome variable
outcome <- "choice"

## Unique respondent ID
respid <- "id"

## Unique task x respondent ID
taskid <- "resp_pair"

# Meta table
rank.meta.df <- data.frame(expand.grid(experiment = c(0:2)))
rank.meta.df$id <- 1:nrow(rank.meta.df)

# cjRank objects
cjr.ls <- lapply(c(as.list(0:2)), function(e){
  print(e)
  cjRank$new(data = data.long[data.long$experiment %in% e &
                                data.long$in_2024, ], 
             outcome.var = outcome,
             task_id.var = taskid,
             attr_cat.var = attr.vars,
             resp_id.var = respid,
             min_obs = 100,
             cluster_var = respid)
})

## Rankings
ranks.ls <- lapply(cjr.ls, function(x){x$get_ranks()})

## Bootstrap
bs.ranks.ls <- lapply(cjr.ls, function(x){x$get_ranks_bs(iter = 200)})

## Nested Marginal Means
nest.mm.ls <- lapply(cjr.ls, function(x){x$get_nested_margmean(adj_coocurrence = TRUE)})

## Within-Rank Marginal Means
win.mm.ls <- lapply(cjr.ls, function(x){x$get_within_rank_margmean(adj_coocurrence = TRUE)})




# FEATRUE RANKING ###############
for(i in rank.meta.df$id){
  ## Put into table
  sum.boot <- bs.ranks.ls[[i]]
  sum.boot$type <- rep(names(attr.labs), each = 3)
  sum.boot$label <- unlist(attr.labs)
  sum.boot <- sum.boot[order(sum.boot$rank, sum.boot$mean),]
  
  # Prepare table
  
  ## Copy
  tab.df <- sum.boot
  
  ## Suppress linebreaks in types and labels
  tab.df <- tab.df[, c("rank", "type", "label", "mean", "median", "lower_bound", "upper_bound")]
  tab.df$type <- gsub("\n", " ", tab.df$type, fixed = T)
  tab.df$label <- gsub("\n", " ", tab.df$label, fixed = T)
  
  ## Feed to stargazer
  tab <- stargazer(tab.df, 
                   summary = FALSE,
                   column.sep.width = "0pt",
                   label = paste0("tab:exp",rank.meta.df$experiment[i], "_ranking"),
                   rownames = FALSE,
                   digits = 2, align = T,
                   title = paste("Ranking of Strategy Features, Experiment", 
                                 rank.meta.df$experiment[i]),
                   notes.append = F, 
                   notes = "Note: Mean, median, and quantiles q2.5 and q97.5 from bootstrapped confidence intervals.")
  
  ## Fix stargazer output
  tab <- gsub("Donbas & Crimea", "Donbas \\& Crimea", tab, fixed = T)
  tab <- c(tab[1:10], 
           paste(paste(paste0("\\multicolumn{1}{c}{",
                              c("Rank", "Attribute", "Level", "Mean", "Median", "q2.5", "q97.5"), 
                              "}"), collapse = " & "), " \\\\ \\\\[-1.8ex]"),
           tab[12:length(tab)])
  tab[12:length(tab)] <- gsub("multicolumn{1}{c}", "multicolumn{1}{l}", 
                              tab[12:length(tab)], fixed = T)
  
  # Write Table
  writeLines(tab, 
             file.path(tab.path, paste0("exp", rank.meta.df$experiment[i],
                                        "_ranking.tex")))
  
  
}




# NESTED MARGINAL MEANS ##############
all.coef.df <- NULL
for(i in rank.meta.df$id){
  # Get results
  coef.df <- nest.mm.ls[[i]]
  
  # F-tests
  ftest.df <- cjr.ls[[i]]$get_rank_fstats()
  
  # Bootstrap results
  sum.boot <- bs.ranks.ls[[i]]
  sum.boot$type <- rep(names(attr.labs), each = 3)
  sum.boot$label <- unlist(attr.labs)
  sum.boot <- sum.boot[order(sum.boot$rank, sum.boot$mean),]
  
  # Prepare Plot
  coef.df$treat_type <- factor(coef.df$attribute ,
                               levels = unique(coef.df$attribute),
                               labels = names(attr.labs),
                               ordered = T)
  coef.df$rank_num <- coef.df$rank
  rank_labs <- c()
  coef.df$rank <- c("Baseline", 
                    paste0("", 
                           c("Russ.-contr. gov", "Full integrity",
                             "Neg. neutrality", rep("", 10))))[coef.df$rank_num + 1] ## Plot only first 4
  
  coef.df$rank_lab <- paste0(coef.df$rank,
                             ifelse(coef.df$rank == "Baseline", 
                                    paste0("\n \n N: ", coef.df$N_full, "\n "),
                                    paste0(c(NA, paste0("\n Rank ", sum.boot$rank, " [", 
                                                        sum.boot$lower_bound, ";", sum.boot$upper_bound, "]")[
                                                          !duplicated(sum.boot$rank)])[coef.df$rank_num + 1],
                                           "\n N: ", coef.df$N_full,"\n ",
                                           "F-stat: ", c(NA, signif(ftest.df$fstat, 2))[coef.df$rank_num + 1])))
  coef.df$rank_lab <- factor(coef.df$rank_lab, levels = unique(coef.df$rank_lab), ordered = T)
  coef.df$experiment <- rank.meta.df$experiment[i]
  
  # Save
  all.coef.df <<- rbind(all.coef.df, coef.df)
  
  # # Plot // Only plot combined graph below
  # g <- ggplot(coef.df[as.numeric(coef.df$rank_lab) <= 4, ], 
  #             aes(x = factor(feature_key), y = coef)) +
  #   geom_hline(yintercept = .5, lty = 2, col = "grey") +
  #   geom_point() +
  #   geom_segment(aes(xend =  factor(feature_key), 
  #                    y = coef - 1.96*se, 
  #                    yend = coef + 1.96*se)) +
  #   scale_y_continuous(breaks = seq(.25, .75, by = .25)) +
  #   facet_grid(treat_type ~ rank_lab) +
  #   theme_minimal() +
  #   labs(x = "Attribute Level", y = "Marginal Means") +
  #   theme(panel.spacing = unit(1, "lines"),
  #         panel.grid.minor.x = element_blank())
  # 
  # # Save
  # png(file.path(fig.path, paste0("exp", rank.meta.df$experiment[i],"_ranks_mm.png")),
  #     width = 6.5, height = 6.5, units = "in", res = 400)
  # print(g)
  # dev.off()
}

# All coefficients in one
all.coef.df$experiment_lab <- c("0 (2022-07)", "I (2024-12/2025-01)", "II (2024-12/2025-01)")[all.coef.df$experiment + 1]
all.coef.df$xpos <- as.numeric(all.coef.df$feature_key) - .2 + all.coef.df$experiment/5
all.coef.df$rank <- factor(all.coef.df$rank, levels = unique(all.coef.df$rank), 
                           ordered = T)
# Plot
g <- ggplot(all.coef.df[all.coef.df$rank != "", ], 
            aes(x = xpos, y = coef, col = experiment_lab)) +
  geom_hline(yintercept = .5, lty = 2, col = "grey") +
  geom_point() +
  geom_segment(aes(xend =  xpos, 
                   y = coef - 1.96*se, 
                   yend = coef + 1.96*se)) +
  scale_y_continuous(breaks = seq(.25, .75, by = .25)) +
  scale_x_continuous(breaks = 1:3, expand = expansion(add = .25)) +
  facet_grid(treat_type ~ rank) +
  theme_minimal() +
  labs(x = "Attribute Level", y = "Marginal Means",
       col = "Experiment") +
  theme(panel.spacing = unit(1, "lines"),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_line(size = .5),
        legend.position = "top")

# Save
png(file.path(fig.path, paste0("all_ranks_mm.png")),
    width = 6.5, height = 6.5, units = "in", res = 400)
print(g)
dev.off()


# WITHIN RANK MARGINAL MEANS ###########

all.coef.df <- NULL
for(i in rank.meta.df$id){
  # Get results
  coef.df <- win.mm.ls[[i]]
  
  # Prepare Plot
  coef.df$treat_type <- factor(coef.df$attribute ,
                               levels = unique(coef.df$attribute),
                               labels = (names(attr.labs)),
                               ordered = T)
  coef.df$rank_num <- coef.df$rank
  coef.df$rank <- c("Russ.-contr. gov", "Full integrity",
                    "Neg. neutrality", "Crimea+Luhansk+Donetsk", rep("", 5))[coef.df$rank_num]
  coef.df$rank <- paste0(coef.df$rank,
                         paste0("\n N: ", coef.df$N_full))
  coef.df$rank <- factor(coef.df$rank, levels = unique(coef.df$rank), ordered = T)
  coef.df$experiment <- rank.meta.df$experiment[i]
  
  # Save
  all.coef.df <- rbind(all.coef.df, coef.df)
  
  # # Plot // Only combined plot below
  # g <- ggplot(coef.df[coef.df$rank_num <= 3, ], 
  #             aes(x = factor(feature_key), y = coef)) +
  #   geom_hline(yintercept = .5, lty = 2, col = "grey") +
  #   geom_point() +
  #   geom_segment(aes(xend =  factor(feature_key), 
  #                    y = coef - 1.96*se, 
  #                    yend = coef + 1.96*se)) +
  #   scale_y_continuous(breaks = seq(.25, .75, by = .25)) +
  #   facet_grid(treat_type ~ rank) +
  #   theme_minimal() +
  #   labs(x = "Attribute Level", y = "Within-Rank Marginal Means") +
  #   theme(panel.spacing = unit(1, "lines"),
  #         panel.grid.minor.x = element_blank())
  # 
  # 
  # # Save
  # png(file.path(fig.path, paste0("exp", rank.meta.df$experiment[i],"_win_rank_mm.png")),
  #     width = 6.5, height = 6.5, units = "in", res = 400)
  # print(g)
  # dev.off()
}


# All coefficients in one
all.coef.df$experiment_lab <- c("0 (2022-07)", "I (2024-12/2025-01)", "II (2024-12/2025-01)")[all.coef.df$experiment + 1]
all.coef.df$xpos <- as.numeric(all.coef.df$feature_key) - .2 + all.coef.df$experiment/5
all.coef.df$rank <- gsub("\\n N.*","",all.coef.df$rank)
all.coef.df$rank <- factor(all.coef.df$rank, levels = unique(all.coef.df$rank), 
                           ordered = T)
# Plot
g <- ggplot(all.coef.df[all.coef.df$rank_num <= 3, ], 
            aes(x = xpos, y = coef, col = experiment_lab)) +
  geom_hline(yintercept = .5, lty = 2, col = "grey") +
  geom_point() +
  geom_segment(aes(xend =  xpos, 
                   y = coef - 1.96*se, 
                   yend = coef + 1.96*se)) +
  scale_y_continuous(breaks = seq(.25, .75, by = .25)) +
  scale_x_continuous(breaks = 1:3, expand = expansion(add = .25)) +
  facet_grid(treat_type ~ rank) +
  theme_minimal() +
  labs(x = "Attribute Level", y = "Within-Rank Marginal Means",
       col = "Experiment") +
  theme(panel.spacing = unit(1, "lines"),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_line(size = .5),
        legend.position = "top")

# Save
png(file.path(fig.path, paste0("all_win_rank_mm.png")),
    width = 6.5, height = 8.5, units = "in", res = 400)
print(g)
dev.off()
