##############################
# WITHIN LOCATION ANALYSIS
#
#
# 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.
#
#############################

# Affectedness mean by round and location
data.long$exp1 <- ifelse(data.long$experiment == 0, 1, 0)
data.long <- data.long %>%
  dplyr::group_by(settle_id, exp1) %>%
  dplyr::mutate(selfaff_mean = mean(as.numeric(selfaffected)-1, na.rm = T),
         famaff_mean = mean(as.numeric(famaffected)-1, na.rm = T)) %>%
  dplyr::ungroup()
data.long <- as.data.frame(data.long)

# Values
hvals <- c("Shelling" = "shelling_mon24_d")
h <- "Effect of violence"



# Estimate
for(c in hvals){
  print(c)
  coef.df <- do.call(rbind, lapply(outvars, function(o){
    data.long$mod <- data.long[,c]
    data.long$mod <- as.numeric(unlist(data.long$mod))
    if(grepl("mon", c)){
      data.long[data.long$experiment == 0,"mod"] <- 1
    }
    
    # Code max per settlement ID
    data.long <- data.long %>%
      dplyr::group_by(settle_id, experiment) %>%
      dplyr::mutate(mod = max(mod)) %>% 
      dplyr::ungroup() 
    
    # Formula
    f <- as.formula(paste(o, " ~ mod + ", 
                          paste(paste0("I(mod==2)*attr",rep(1:5, 2), "_d", rep(2:3, each = 5)),
                                collapse = " + "),
                          "| ", paste(paste0("settle_id[", paste0("attr",rep(1:5, 2), "_d", rep(2:3, each = 5)), "]"),
                                      collapse = " +"),
                          "+ ", paste(paste0("experiment[", paste0("attr",rep(1:5, 2), "_d", rep(2:3, each = 5)), "]"),
                                      collapse = " +")
    ))
    # Estimate
    m <- fixest::feols(f,
                       data = data.long[data.long$in_2022, ],
                       cluster = "settle_id")
    
    # Extract results
    data.frame(outcome = names(outvars)[outvars == o],
               hetvar = h,
               hetvar_val = c,
               level = rep(1:3, each = 5),
               treat = paste0("attr",rep(1:5, 3), "_d", rep(1:3, each = 5)),
               treat_lab = as.vector(do.call(rbind, attr.labs)),
               treat_type = rep(names(attr.labs),3),
               coef = c(rep(0, 5), m$coefficients[2:11]),
               se = c(rep(NA, 5), m$se[2:11]),
               p = c(rep(NA, 5), summary(m)$coeftable[2:11,4]))
    # }))
  }))
  
  # Prepare plot
  coef.df$outcome <- factor(coef.df$outcome, 
                            levels = names(outvars),
                            ordered = T)
  coef.df$treat_type <- factor(coef.df$treat_type, 
                               levels = names(attr.labs),
                               ordered = T)
  coef.df$treat_lab <- factor(coef.df$treat_lab, 
                              levels = rev(unique(unlist(attr.labs))),
                              ordered = T)
  
  coef.df$het_type <- coef.df$hetvar
  coef.df$hetvar_val <- as.numeric(factor(coef.df$hetvar_val,
                                          levels = hvals, ordered = T))

  coef.df$yval <- as.numeric(coef.df$treat_lab) +
    -1*((coef.df$hetvar_val - mean(coef.df$hetvar_val)) /
          max(coef.df$hetvar_val)) / 1.5
  
  # Plot
  g <- ggplot(coef.df, 
              aes(x = coef, y = yval)) +
    geom_vline(xintercept = 0, lty = 2, col = "grey") +
    geom_point() +
    geom_segment(aes(yend = yval, 
                     x = coef - 1.96*se, 
                     xend = coef + 1.96*se)) +
    facet_grid(treat_type ~ outcome,
               scales  = "free", switch = "y") +
    scale_y_continuous(breaks = as.numeric(coef.df$treat_lab) ,
                       labels = as.character(coef.df$treat_lab),
                       minor_breaks = seq(.5, 100.5, by = 1),
                       position = "right", expand = expansion(add = .4)) +
    theme_minimal() +
    labs(y = NULL,
      x = paste0("Differential Change in AMCE since July 2022 "),
      col = h, pch = h,
      caption = paste("Column names indicates outcome variable. ")) +
    theme(panel.spacing = unit(.75, "lines"),
          panel.grid.major.y = element_blank(),
          strip.text.y.left = element_text(angle = 0, hjust = 0),
          panel.grid.minor.y = element_line(size = .5),
          legend.position = "top",
          legend.box.margin=margin(-10,-10,-10,-10),
          plot.caption = element_text(hjust = 0, colour = "grey20")) +
    scale_color_manual(values = custom.col)
  
  
  png(file.path(fig.path, paste0("winloc_", c,".png")),
      width = 7, height = 6, units = "in", res = 400)
  print(g)
  dev.off()
}

# Values
hvals <- c("Self-affectedness" = "selfaff_mean",  "Fam-affectedness" = "famaff_mean")
h <- "Effect of violence"


# Estimate
for(c in hvals){
  print(c)
  coef.df <- do.call(rbind, lapply(outvars, function(o){

    data.long$mod <- data.long[,c]
    data.long$mod <- as.numeric(unlist(data.long$mod))

    # Code max per settlement ID
    data.long <- data.long %>%
      dplyr::group_by(settle_id, experiment) %>%
      dplyr::mutate(mod = mean(mod)) %>% 
      dplyr::ungroup() 
    
    # Formula
    f <- as.formula(paste(o, " ~ mod + ", 
                          paste(paste0("I(mod)*attr",rep(1:5, 2), "_d", rep(2:3, each = 5)),
                                collapse = " + "),
                          "| ", paste(paste0("settle_id[", paste0("attr",rep(1:5, 2), "_d", rep(2:3, each = 5)), "]"),
                                      collapse = " +"),
                          "+ ", paste(paste0("experiment[", paste0("attr",rep(1:5, 2), "_d", rep(2:3, each = 5)), "]"),
                                      collapse = " +")
    ))
    # Estimate
    m <- fixest::feols(f,
                       data = data.long[data.long$in_2022, ],
                       cluster = "settle_id")
    
    # Extract results
    data.frame(outcome = names(outvars)[outvars == o],
               hetvar = h,
               hetvar_val = c,
               level = rep(1:3, each = 5),
               treat = paste0("attr",rep(1:5, 3), "_d", rep(1:3, each = 5)),
               treat_lab = as.vector(do.call(rbind, attr.labs)),
               treat_type = rep(names(attr.labs),3),
               coef = c(rep(0, 5), m$coefficients[2:11]),
               se = c(rep(NA, 5), m$se[2:11]),
               p = c(rep(NA, 5), summary(m)$coeftable[2:11,4]))
    # }))
  }))
  
  # Prepare plot
  coef.df$outcome <- factor(coef.df$outcome, 
                            levels = names(outvars),
                            ordered = T)
  coef.df$treat_type <- factor(coef.df$treat_type, 
                               levels = names(attr.labs),
                               ordered = T)
  coef.df$treat_lab <- factor(coef.df$treat_lab, 
                              levels = rev(unique(unlist(attr.labs))),
                              ordered = T)
  
  coef.df$het_type <- coef.df$hetvar
  coef.df$hetvar_val <- as.numeric(factor(coef.df$hetvar_val,
                                          levels = hvals, ordered = T))

  coef.df$yval <- as.numeric(coef.df$treat_lab) +
    -1*((coef.df$hetvar_val - mean(coef.df$hetvar_val)) /
          max(coef.df$hetvar_val)) / 1.5
  
  # Plot
  g <- ggplot(coef.df, 
              aes(x = coef, y = yval
              )) +
    geom_vline(xintercept = 0, lty = 2, col = "grey") +
    geom_point() +
    geom_segment(aes(yend = yval, 
                     x = coef - 1.96*se, 
                     xend = coef + 1.96*se)) +
    facet_grid(treat_type ~ outcome,
               scales  = "free", switch = "y") +
    scale_y_continuous(breaks = as.numeric(coef.df$treat_lab) ,
                       labels = as.character(coef.df$treat_lab),
                       minor_breaks = seq(.5, 100.5, by = 1),
                       position = "right", expand = expansion(add = .4)) +
    theme_minimal() +
    labs(y = NULL,
      x = paste0("Differential Change in AMCE since July 2022 "),
      col = h, pch = h,
      caption = paste("Column names indicates outcome variable. ")) +
    theme(panel.spacing = unit(.75, "lines"),
          panel.grid.major.y = element_blank(),
          strip.text.y.left = element_text(angle = 0, hjust = 0),
          panel.grid.minor.y = element_line(size = .5),
          legend.position = "top",
          legend.box.margin=margin(-10,-10,-10,-10),
          plot.caption = element_text(hjust = 0, colour = "grey20")) +
    scale_color_manual(values = custom.col)
  
  
  png(file.path(fig.path, paste0("winloc_", c,".png")),
      width = 7, height = 6, units = "in", res = 400)
  print(g)
  dev.off()
}

