###########################################
# SAMPLE CHARACTERISTICS
#
#
# 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.
#
#############################


# VARIABLES ###############################

# Variables along which to compute heterogeneous effects
these.hetvars <- c("Gender" =  "sex",  
             "Age" = "age_5gr", 
             "Children" = "children",  
             "Education" = "educ",  
             "Econ.\ndepriv." = "econdepr",  
             "Rural\nUrban" = "urbrur", 
             "Interview\nlanguage" = "lang",  
             "Native\nlanguage" = "natlang",
             "Ethnic\nidentity" = "ethnic", 
             
             "Affectedness score" = "aff_score",
             "East vs. West" = "region",  
             "Oblast first attacked" = "diraffected",
             "Self war-affected" = "selfaffected",
             "Family war-affected" = "famaffected",
             "Any oneside violence" = "onesided_all_d",
             "Any battles" = "battle_all_d" ,
             "Any shelling" = "shelling_all_d",
             
             "Ukr. nation at stake" = "survnat" ,
             "Importance of victory" = "victimp",  
             "Trust in president" = "trustpres")


# SUMMARY STATS FOR VARIABLES #################
# Tables A2 -- A4

## IDs
stats <- list("Demographics" = c(1:9), 
              "Affectedness" = c(10:17), 
              "Political attitudes" = c(18:20))

## Summary stats tables via stargazer

### This data
data.sum <- data.long[data.long$experiment %in% 1:2 & data.long$in_2024, ]

### Estimate Aff Score
pca <- prcomp(data.frame(lapply(na.omit(data.sum[,these.hetvars[c(11:17)]]), function(x){as.numeric(as.factor(x))})),
              center = TRUE, 
              scale. = TRUE) 
print(summary(pca))
print(pca)

#### Save to data
data.sum[apply(data.sum[,these.hetvars[c(11:17)]], 1, function(x){!any(is.na(x))}), "aff_score"] <- 
  pca$x[,1]

#### Cut into 3 values
data.sum$aff_score <- cut(data.sum$aff_score, c(-99, quantile(data.sum$aff_score, 1:3/3, na.rm = T)),
                          labels = c("low", "medium", "high"))


for(s in seq_along(stats)){
  stat.tab <- NULL
  stat.ls <- lapply(stats[[s]], function(i){
    var <- these.hetvars[i]
    vals <- sort(na.omit(unique(data.sum[,var])))
    df <- data.frame(Statistic = names( table(data.sum[!duplicated(data.sum$id),var])), 
                     N =  as.vector(table(data.sum[!duplicated(data.sum$id),var])), 
                     Mean = as.vector(table(data.sum[!duplicated(data.sum$id),var])/
                                        sum(table(data.sum[!duplicated(data.sum$id),var]))))
    df <- df[df$N > 0,]
    df
  })
  tab.df <- do.call(rbind, stat.ls)
  
  # Stargazer
  tab <- stargazer(tab.df,
                   title = paste0("Respondent-level summary statistics: ", names(stats)[s]),
                   summary = FALSE,
                   rownames = FALSE,
                   notes.align = "l", align = T,
                   digits = 2, 
                   digit.separate = 0, 
                   column.sep.width ="-10pt",
                   type  = "latex")
  
  # Add subheadings and lines
  tab.fin <- tab[1:11]
  running <- 13
  for(i in seq_along(stat.ls)){
    tab.fin <- c(tab.fin, 
                 paste0("\\hline \\\\[-1.8ex] " , "\\multicolumn{1}{c}{\\textbf{",
                        names(these.hetvars)[stats[[s]][i]], "}} && \\\\"),
                 tab[running:(running + nrow(stat.ls[[i]]) - 1)])
    running <- running + nrow(stat.ls[[i]]) 
  }
  tab.fin <- c(tab.fin, 
               tab[running:length(tab)])
  tab.fin[12:length(tab.fin)] <- gsub("multicolumn{1}{c}", "multicolumn{1}{l}", 
                                      tab.fin[12:length(tab.fin)], fixed = T)
  
  # Write to disc
  writeLines(tab.fin, 
             file.path(tab.path, paste0("sumstats_", s,".tex")))
}


# MAKE AFFECTEDNESS SCORE PCA #################
sum.df <- NULL
for(sample in c("in_2022", "in_2024")){
  for(e in 0:2){
    # Subset 
    data.sub <- data.long[data.long$experiment == e &
                            data.long[,sample] == 1, ]
    
    
    # Estimate
    pca <- prcomp(data.frame(lapply(na.omit(data.sub[,these.hetvars[c(11:17)]]), function(x){as.numeric(as.factor(x))})),
                  center = TRUE, 
                  scale. = TRUE) 
    print(summary(pca))
    print(pca)
    
    # Save to data
    data.sub[apply(data.sub[,these.hetvars[c(11:17)]], 1, function(x){!any(is.na(x))}), "aff_score"] <- 
      pca$x[,1]
    
    # Cut into 3 values
    data.sub$aff_score <- cut(data.sub$aff_score, c(-99, quantile(data.sub$aff_score, 1:3/3, na.rm = T)),
                              labels = c("low", "medium", "high"))
    
    
    # SUMMARY STATS FOR VARIABLES #################
    # Tables A2 -- A4
    
    ## IDs
    stats <- list("Demographics" = c(1:9), 
                  "Affectedness" = c(10:17), 
                  "Political attitudes" = c(18:20))
    
    ## Summary stats tables via stargazer
    for(s in seq_along(stats)){
      stat.tab <- NULL
      stat.ls <- lapply(stats[[s]], function(i){
        var <- these.hetvars[i]
        vals <- sort(na.omit(unique(data.sub[,var])))
        df <- data.frame(Statistic = names( table(data.sub[!duplicated(data.sub$id),var])), 
                         N =  as.vector(table(data.sub[!duplicated(data.sub$id),var])), 
                         Mean = as.vector(table(data.sub[!duplicated(data.sub$id),var])/
                                            sum(table(data.sub[!duplicated(data.sub$id),var]))),
                         var = names(these.hetvars)[i])
        df <- df[df$N > 0,]
        df
      })
      tab.df <- do.call(rbind, stat.ls)
      
      # Save
      sum.df <- rbind(sum.df, 
                      cbind(experiment = rep(e, nrow(tab.df)),
                            sample = sample,
                            tab.df))
      
      # Don't Print individual summary stats for replication data
      # # Stargazer
      # tab <- stargazer(tab.df[, !colnames(tab.df) %in% c("var")],
      #                  title = paste0("Experiment ", e, " Respondent-level summary statistics: ", names(stats)[s]),
      #                  summary = FALSE,
      #                  rownames = FALSE,
      #                  notes.align = "l", align = T,
      #                  digits = 2, 
      #                  digit.separate = 0, 
      #                  column.sep.width ="-10pt",
      #                  type  = "latex")
      # 
      # # Add subheadings and lines
      # tab.fin <- tab[1:11]
      # running <- 13
      # for(i in seq_along(stat.ls)){
      #   tab.fin <- c(tab.fin, 
      #                paste0("\\hline \\\\[-1.8ex] " , "\\multicolumn{1}{c}{\\textbf{",
      #                       names(these.hetvars)[stats[[s]][i]], "}} && \\\\"),
      #                tab[running:(running + nrow(stat.ls[[i]]) - 1)])
      #   running <- running + nrow(stat.ls[[i]]) 
      # }
      # tab.fin <- c(tab.fin, 
      #              tab[running:length(tab)])
      # tab.fin[12:length(tab.fin)] <- gsub("multicolumn{1}{c}", "multicolumn{1}{l}", 
      #                                     tab.fin[12:length(tab.fin)], fixed = T)
      # 
      # # Write to disc
      # writeLines(tab.fin, 
      #            file.path(tab.path, paste0("exp", e, sample, "_sumstats_", s,".tex")))
    }
  }
}

# Summary figure
sum.df$Statistic <- gsub("\\s*\\([^\\)]+\\)","",sum.df$Statistic)
sum.df$Statistic[grepl("Vocational", sum.df$Statistic)] <- "1 Vocational"
sum.df$Statistic[grepl("Secondary special", sum.df$Statistic)] <- "4 Special Secondary"
sum.df$Statistic[grepl("Higher", sum.df$Statistic)] <- "5 Higher"
sum.df$Statistic[grepl("Complete general", sum.df$Statistic)] <- "3 Complete Secondary"
sum.df$Statistic[grepl("Basic general", sum.df$Statistic)] <- "2 Basic Secondary"
sum.df$Statistic <- factor(sum.df$Statistic, levels = rev(sort(unique(sum.df$Statistic))), ordered = T)
sum.df$var <- factor(sum.df$var, levels = unique(sum.df$var), ordered = T)

# Drop Exp 0 in 2024
sum.df <- sum.df[!(sum.df$experiment == 0 & sum.df$sample == "in_2024"),]
sum.df$sample <- c(in_2022 = "Sample A", in_2024 = "Sample B")[sum.df$sample]
sum.df$experiment_lab <- c("2022-07", "2024-12.1", "2024-12.2")[sum.df$experiment + 1]


g <- ggplot(sum.df[sum.df$var %in% names(these.hetvars[1:9]),], 
            aes(x = Mean, y = Statistic, col = factor(experiment_lab), pch = factor(sample))) +
  geom_point(alpha = .8, size = 3) +
  theme_minimal() +
  # scale_y_discrete(position = "right") +
  facet_grid(var ~ 0, scales = "free_y", space = "free_y", switch = "y") +
  labs(color = "Experiment", pch = "Sample", y = NULL) +
  theme(strip.placement = "outside",
        strip.background = element_blank(),
        strip.text.x = element_blank(),
        legend.position = "top",
        legend.box="vertical", legend.margin=margin())

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


