##############################
# ROBUSTNESS CHECKS
#
#
# 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.
#
#############################


# MAIN LINEAR RESULTS #########

# Estimation
coef.df <- do.call(rbind, lapply(outvars, function(o){
  # Formula
  f <- as.formula(paste(o, " ~ ", 
                        paste(paste0("attr",rep(1:5, 1), "_level"), ## As linear terms, not dummies
                              collapse = " + ")))
  # Estimate
  m <- fixest::feols(f,
                     data = data.long[data.long$experiment %in% c(1:2) & data.long$in_2024,],
                     cluster = clustvar)
  
  # Extract results
  data.frame(outcome = names(outvars)[outvars == o],
             treat = paste0("attr",rep(1:5, 1), "_level"),
             treat_type = rep(names(attr.labs),1),
             coef = c(rep(0, 0), m$coefficients[2:6]),
             se = c(rep(NA, 0), m$se[2:6]),
             p = c(rep(NA, 0), summary(m)$coeftable[2:6,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 = rev(names(attr.labs)),
                             ordered = T)

# Plot
g <- ggplot(coef.df, 
            aes(x = coef, y = treat_type)) +
  geom_vline(xintercept = 0, lty = 2, col = "grey") +
  geom_point() +
  geom_segment(aes(yend = treat_type, 
                   x = coef - 1.96*se, 
                   xend = coef + 1.96*se)) +
  facet_grid( ~ outcome,
              scales  = "free") +
  theme_minimal() +
  labs(y = NULL, x = "Average Linear Attribute Effect") +
  theme(panel.spacing = unit(1, "lines"),
        panel.grid.minor.y = element_blank())

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


# AVERAGE FEATURE CHOICE PROBABILITIES ####
#  A la Abramson et al., 2020
# Figure A1

# Contrasts:
contrasts <- expand.grid(a = 1:3, b = 1:3)
contrasts <- contrasts[contrasts$a < contrasts$b,]

# Estimation:
coef.df <- do.call(rbind, lapply(1:nrow(contrasts), function(i){
  do.call(rbind, lapply(1:5, function(a){
    # Formula
    f <- as.formula(paste("choice", " ~ 0 +", 
                          paste(paste0("attr",a, "_d", contrasts[i,]),
                                collapse = " + ")))
    # Data Selection
    # attrX_amb indicates ambiguity / variance on an attribute
    this.df <- data.long[data.long[,paste0("attr",a, "_level")] %in% 
                           contrasts[i,],]
    this.df <- this.df[this.df$resp_pair %in% this.df$resp_pair[duplicated(this.df$resp_pair)] &
                         this.df[,paste0("attr",a, "_amb")] != 1 &
                         this.df$experiment %in% c(1:2) &
                         this.df$in_2024,]
    
    # Estimate
    m <- fixest::feols(f,
                       data = this.df,
                       cluster = clustvar)
    
    # Extract results
    data.frame(outcome = names(outvars)[outvars == "choice"],
               contrast = paste(contrasts[i,1], "vs.", contrasts[i,2]),
               level = contrasts[i,],
               treat = paste0("attr",a, "_level"),
               treat_lab = attr.labs[[a]][unlist(contrasts[i,])],
               treat_type = names(attr.labs)[a],
               coef = m$coefficients,
               se = m$se,
               p = summary(m)$coeftable[,4])
  }))
}))

# Prepare plot
coef.df$outcome <- factor(coef.df$outcome, 
                          levels = names(outvars),
                          ordered = T)
coef.df$contrast <- factor(coef.df$contrast, 
                           levels = unique(coef.df$contrast),
                           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(coef.df$treat_lab)),
                            ordered = T)
# Plot
g <- ggplot(coef.df, 
            aes(x = coef, y = treat_lab)) +
  geom_vline(xintercept = 0.5, lty = 2, col = "grey") +
  geom_point() +
  geom_segment(aes(yend = treat_lab, 
                   x = coef - 1.96*se, 
                   xend = coef + 1.96*se)) +
  facet_grid(treat_type ~ contrast,
             scales  = "free_y", switch = "y") +
  theme_minimal() +
  scale_x_continuous(breaks = seq(.25, .75, by = .25)) +
  labs(y = NULL, x = "Average Feature Choice Probability") +
  scale_y_discrete(position = "right") +
  theme(panel.spacing = unit(.75, "lines"),
        panel.grid.minor.y = element_blank(),
        strip.text.y.left = element_text(angle = 0, hjust = 0))

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



# AMCIEs -- INTERACTION EFFECTS ###########
# Figure A25 and A26

# Estimation
for(o in c("score", "choice")){
  intvars <- c("attr1_level", "attr5_level")
  coef.df <- do.call(rbind, lapply(intvars, function(i){
    
    # Locals
    if((i == "attr1_level")){
      attr.vars <- paste0("attr",rep(2:5, 2), "_d", rep(2:3, each = 4))
      all.treat.vars <- paste0("attr",rep(2:5, 3), "_d", rep(1:3, each = 4))
      treat_lab = as.vector(do.call(rbind, attr.labs[2:5]))
      treat_type = rep(names(attr.labs)[2:5],3)
      int_lab <- as.vector(do.call(rbind, attr.labs[1]))
    } else {
      attr.vars <- paste0("attr",rep(1:4, 2), "_d", rep(2:3, each = 4))
      all.treat.vars <- paste0("attr",rep(1:4, 3), "_d", rep(1:3, each = 4))
      treat_lab = as.vector(do.call(rbind, attr.labs[1:4]))
      treat_type = rep(names(attr.labs)[1:4],3)
      int_lab <- as.vector(do.call(rbind, attr.labs[5]))
    }
    int.dummies <- sapply(c("d1", "d2", "d3"), function(s){gsub("level", s, i)})
    
    # Interaction terms
    int.terms <- paste0("I(",
                        rep(attr.vars, 2), "*",
                        rep(paste0("(",i, "==", 2:3, ")"), each = length(attr.vars)), ")")
    
    # Formula
    f <-  as.formula(paste(o, " ~ ",  "factor(", i, ") + ", 
                           paste(c(attr.vars, 
                                   int.terms),
                                 collapse = " + ")))
    
    # Estimate model
    m <- fixest::feols(f,
                       data = data.long[data.long$experiment %in% c(1:2) & data.long$in_2024,],
                       cluster = clustvar)
    
    # F-Test for joint significance
    restr.vars <- paste0("I(",
                         rep(attr.vars[grepl("attr2|attr3|attr4",attr.vars)], each = 2), "*",
                         rep(paste0("(",i, "==", 2:3, ")"), 2), ")")
    mf <-  as.formula(paste(o, " ~ ",  "factor(", i, ") + ", 
                            paste(c(attr.vars[grepl("attr2|attr3|attr4",attr.vars)], 
                                    restr.vars),
                                  collapse = " + ")))
    mf <- fixest::feols(mf,
                        data = data.long[data.long$experiment %in% c(1:2),],
                        cluster = clustvar)
    test <- wald(mf, 
                 keep = paste0(i, " =="), 
                 cluster = clustvar,
                 print = F)
    
    # Extract results
    res.df <- rbind(data.frame(outcome = names(outvars)[outvars == o],
                               intvar = i,
                               intvar_val = rep(1:3),
                               mvar = c(NA, names(m$coefficients[2:3])),
                               level = rep(1:3),
                               treat = int.dummies,
                               treat_lab = int_lab,
                               treat_type = names(attr.labs)[!names(attr.labs) %in% treat_type],
                               coef = c(0, m$coefficients[2:3]),
                               se = c(NA, m$se[2:3]),
                               p = c(NA, summary(m)$coeftable[2:3,4]),
                               wald_stat = test$stat, wald_pval = test$p),
                    data.frame(outcome = names(outvars)[outvars == o],
                               intvar = i,
                               intvar_val = 1,
                               mvar = NA,
                               level = rep(1:3, each = 4),
                               treat = rep(all.treat.vars[1:4], 3),
                               treat_lab = rep(treat_lab[1:4], 3),
                               treat_type = rep(treat_type[1:4], 3),
                               coef = c(rep(0, 4*3)),
                               se = c(rep(NA, 4*3)),
                               p = c(rep(NA, 4*3)),
                               wald_stat = test$stat, wald_pval = test$p),
                    data.frame(outcome = names(outvars)[outvars == o],
                               intvar = i,
                               intvar_val = rep(1:3, each = 4*2),
                               mvar = names(m$coefficients[-c(1:3)]),
                               level = rep(rep(2:3, each = 4), 3),
                               treat = rep(all.treat.vars[-c(1:4)], 3),
                               treat_lab = rep(treat_lab[-c(1:4)], 3),
                               treat_type = rep(treat_type[-c(1:4)], 3),
                               coef = m$coefficients[-c(1:3)],
                               se = m$se[-c(1:3)],
                               p = summary(m)$coeftable[-c(1:3),4],
                               wald_stat = test$stat, wald_pval = test$p)
    )
    # Return
    return(res.df)
  }))
  
  # Prepare plot // labels and positions
  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$int_type <- ifelse(coef.df$intvar == "attr1_level",
                             "Territorial Integrity", "Political Autonomy")
  
  coef.df$int_type <- factor(coef.df$int_type, levels = unique(coef.df$int_type),
                             ordered = T)
  
  coef.df$yval <- as.numeric(coef.df$treat_lab) +
    -1*(coef.df$intvar_val - 2)/4
  coef.df$yval <- ifelse(gsub("\n", " ", as.character(coef.df$treat_type)) == coef.df$int_type,
                         as.numeric(coef.df$treat_lab),
                         coef.df$yval)
  
  coef.df$intvar_lab <- c("Baseline", "x (M == 2)", "x (M == 3)")[coef.df$intvar_val]
  coef.df$intvar_lab <- factor(coef.df$intvar_lab, 
                               levels = unique(coef.df$intvar_lab),
                               ordered = T)
  # Plot
  g <- ggplot(coef.df, 
              aes(x = coef, y = yval, 
                  col = intvar_lab, 
                  pch = intvar_lab)) +
    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 ~ int_type,
               scales  = "free",
               switch = "y") +
    scale_y_continuous(breaks = as.numeric(coef.df$treat_lab) ,
                       labels = as.character(coef.df$treat_lab),
                       position = "right") +
    theme_minimal() +
    labs(y = NULL,
         x = paste0("Average Marginal Component Interaction Effect on ", names(outvars)[outvars == o]),
         col = "Level of interaction term:", pch = "Level of interaction term:",
         caption = paste0("Column names show name of moderator variable.\n",
                          "Wald-test of joint nullity of interaction terms with Attributes 2 to 4:\n",
                          paste(sapply(unique(as.character(coef.df$int_type)), function(t){
                            paste0(t, " (F: ", 
                                   signif(coef.df$wald_stat[coef.df$int_type == t][1], 2), "; p: ",
                                   signif(coef.df$wald_pval[coef.df$int_type == t][1], 2), ")")
                          }), collapse = ", "))) +
    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)
  
  # Save
  png(file.path(fig.path, paste0("exp12_amcie_",o,".png")),
      width = 7, height = 5, units = "in", res = 400)
  print(g)
  dev.off()
  
}


# LINEAR INTERACTIONS ###########
# Figure A27 and A28

# Estimation
for(o in c("score", "choice")){
  intvars <- c("attr1_level", "attr5_level")
  coef.df <- do.call(rbind, lapply(intvars, function(i){
    
    # Locals
    if((i == "attr1_level")){
      attr.vars <- paste0("attr",rep(2:5, 1), "_level")
      all.treat.vars <- paste0("attr",rep(2:5, 1), "_level")
      treat_type = rep(names(attr.labs)[2:5],1)
      
    } else {
      attr.vars <- paste0("attr",rep(1:4, 1), "_level")
      all.treat.vars <- paste0("attr",rep(1:4, 1), "_level")
      treat_type = rep(names(attr.labs)[1:4],1)
    }
    
    # Formula
    f <- as.formula(paste(o, " ~ ",  "factor(", i, ") + ", 
                          paste(paste0("I(",
                                       rep(attr.vars, each = 3), "*",
                                       rep(paste0("(",i, "==", 1:3, ")"), 2), ")"),
                                collapse = " + ")))
    # Estimate
    m <- fixest::feols(f,
                       data = data.long[data.long$experiment %in% c(1:2) & data.long$in_2024,],
                       cluster = clustvar)
    
    # Extract results
    data.frame(outcome = names(outvars)[outvars == o],
               intvar = i,
               intvar_val = rep(1:3, length(attr.vars)),
               treat = rep(all.treat.vars, each = 3),
               treat_type = rep(treat_type, each = 3),
               coef = c(m$coefficients[-c(1:3)]),
               se = c(m$se[-c(1:3)]),
               p = c(summary(m)$coeftable[-c(1:3),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 = rev(names(attr.labs)),
                               ordered = T)
  coef.df$int_type <- ifelse(coef.df$intvar == "attr1_level",
                             "Territorial Integrity", "Political Autonomy")
  
  coef.df$int_type <- factor(coef.df$int_type, levels = unique(coef.df$int_type),
                             ordered = T)
  
  coef.df$yval <- as.numeric(coef.df$treat_type) +
    -1*(coef.df$intvar_val - 2)/4
  
  # Plot
  g <- ggplot(coef.df, 
              aes(x = coef, y = yval, 
                  col = factor(intvar_val), 
                  pch = factor(intvar_val))) +
    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( ~ int_type,
                scales  = "free") +
    scale_y_continuous(breaks = as.numeric(coef.df$treat_type) ,
                       labels = as.character(coef.df$treat_type) ) +
    theme_minimal() +
    labs(y = NULL,
         x = paste0("Average Linear Attribute Effect on ", names(outvars)[outvars == o]),
         col = "Level of moderator variable:", pch = "Level of moderator variable:",
         caption = "Column names indicate moderator variables. ") +
    theme(panel.spacing = unit(1, "lines"),
          legend.position = "top",
          panel.grid.major.y = element_blank(),
          panel.grid.minor.y = element_line(size = .5)) +
    scale_color_manual(values = custom.col)
  
  # Save
  png(file.path(fig.path, paste0("exp12_moder_amce_lin_",o,".png")),
      width = 7, height = 3, units = "in", res = 400)
  print(g)
  dev.off()
  
}


# LOGIT: MAIN RESULTS #########
# Figure A4

# Estimation
coef.df <- do.call(rbind, lapply(outvars[2], function(o){
  # Formula
  f <- as.formula(paste(o, " ~ ", 
                        paste(paste0("attr",rep(1:5, 2), "_d", rep(2:3, each = 5)),
                              collapse = " + ")))
  # Estimate
  m <- fixest::feglm(f,
                     data = data.long[data.long$experiment %in% c(1:2) & data.long$in_2024,],
                     family = "binomial",
                     cluster = clustvar)
  
  # Extract results
  data.frame(outcome = names(outvars)[outvars == o],
             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)

# Plot
g <- ggplot(coef.df, 
            aes(x = coef, y = treat_lab)) +
  geom_vline(xintercept = 0, lty = 2, col = "grey") +
  geom_point() +
  geom_segment(aes(yend = treat_lab, 
                   x = coef - 1.96*se, 
                   xend = coef + 1.96*se)) +
  facet_grid(treat_type ~ outcome,
             scales  = "free", switch = "y") +
  theme_minimal() +
  labs(y = NULL, x = "Average Marginal Component Effect") +
  scale_y_discrete(position = "right") +
  theme(panel.spacing = unit(.75, "lines"),
        panel.grid.minor.y = element_blank(),
        strip.text.y.left = element_text(angle = 0, hjust = 0)) +
  scale_color_manual(values = custom.col)

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



# Weighted AMCE #############
# Figure A2 



# Estimation
coef.df <- do.call(rbind, lapply(outvars, function(o){
  # Formula
  f <- as.formula(paste(o, " ~ ", 
                        paste(paste0("attr",rep(1:5, 2), "_d", rep(2:3, each = 5)),
                              collapse = " + ")))
  # Estimate
  m <- fixest::feols(f,
                     data = data.long[data.long$experiment %in% c(1:2) & data.long$in_2024,],
                     cluster = clustvar,
                     weights = data.long[data.long$experiment %in% c(1:2) & data.long$in_2024,"hh_size"])
  
  # Extract results
  data.frame(outcome = names(outvars)[outvars == o],
             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(coef.df$treat_lab)),
                            ordered = T)


# Plot
g <- ggplot(coef.df, 
            aes(x = coef, y = treat_lab)) +
  geom_vline(xintercept = 0, lty = 2, col = "grey") +
  geom_point() +
  geom_segment(aes(yend = treat_lab, 
                   x = coef - 1.96*se, 
                   xend = coef + 1.96*se)) +
  facet_grid(treat_type ~ outcome,
             scales  = "free", switch = "y") +
  theme_minimal() +
  labs(y = NULL, x = "Average Marginal Component Effect") +
  scale_y_discrete(position = "right") +
  theme(panel.spacing = unit(.75, "lines"),
        panel.grid.minor.y = element_blank(),
        strip.text.y.left = element_text(angle = 0, hjust = 0)) +
  scale_color_manual(values = custom.col)

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



# ORDER EFFECTS #######
# Figure A6

# Values
hvals <- 1:3
h <- "Attribute position"

# Estimate
coef.df <- do.call(rbind, lapply(outvars, function(o){
  do.call(rbind, lapply(1:5, function(a){
    data.sub <- data.long[data.long$experiment %in% c(1:2) & data.long$in_2024,]
    
    # Formula
    if(a %in% 2:4){
      order <- apply(do.call(rbind, strsplit(as.character(data.sub$attr_order), "-", fixed = T)),
                     2, as.numeric)
      data.sub$pos <- apply(order, 1, function(x){which(x == a)[1]})
      f <- as.formula(paste(o, " ~ ", "factor(", "pos", ") + ", 
                            paste(paste0("I(",
                                         rep(paste0("attr",a, "_d", 2:3), each = length(hvals)), "*",
                                         rep(paste0("(","pos", "==", hvals, ")"), 2), ")"),
                                  collapse = " + ")))
    } else {
      f <- as.formula(paste(o, " ~ ", paste(paste0("attr",a, "_d", 2:3), collapse = "+")))
    }
    
    # Estimate
    m <- fixest::feols(f,
                       data =  data.sub,
                       cluster = clustvar)
    data.sub$pos <- NULL
    
    # Extract results
    if(a %in% 2:4){
      start <- length(hvals) + 1
      end <- start + length(hvals)*2 -1
      res.df <- data.frame(outcome = names(outvars)[outvars == o],
                 hetvar = h,
                 hetvar_val = 1+as.numeric(rep(hvals, 3)),
                 hetvar_lab = 1+rep(hvals, 3),
                 level = rep(1:3, each = length(hvals)),
                 treat = rep(paste0("attr",a, "_d", 1:3), each = length(hvals)),
                 treat_lab = rep(attr.labs[[a]], each = length(hvals)),
                 treat_type = names(attr.labs)[a],
                 coef = c(rep(0,length(hvals)), m$coefficients[start:end]),
                 se = c(rep(NA,length(hvals)), m$se[start:end]),
                 p = c(rep(NA,length(hvals)), summary(m)$coeftable[start:end,4]))
    } else {
      start <- 1+ 1
      end <- start + 1*2 -1
      res.df <- data.frame(outcome = names(outvars)[outvars == o],
                 hetvar = h,
                 hetvar_val = a,
                 hetvar_lab = a,
                 level = rep(1:3, each = 1),
                 treat = rep(paste0("attr",a, "_d", 1:3), each = 1),
                 treat_lab = rep(attr.labs[[a]], each = 1),
                 treat_type = names(attr.labs)[a],
                 coef = c(rep(0,1), m$coefficients[start:end]),
                 se = c(rep(NA,1), m$se[start:end]),
                 p = c(rep(NA,1), summary(m)$coeftable[start:end,4]))
    }
    return(res.df)
  }))
}))

# 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(coef.df$treat_lab)),
                            ordered = T)
coef.df$het_type <- coef.df$hetvar

coef.df$yval <- ifelse(coef.df$hetvar_val %in% 2:4,
                       as.numeric(coef.df$treat_lab) +
                         -1*((coef.df$hetvar_val - 3) /
                               3) / 2,
                       as.numeric(coef.df$treat_lab))

# Plot
g <- ggplot(coef.df, 
            aes(x = coef, y = yval, 
                col = factor(hetvar_lab), 
                pch = factor(hetvar_lab))) +
  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), 
                     position = "right") +
  theme_minimal() +
  labs(y = NULL,
       x = paste0("Average Marginal Component Effect "),
       col = "Attribute position:", pch = "Attribute position:" ,
       caption = paste("Position of attributes 2-4 varied randomly between respondents. ")
       ) +
  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)

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



# SE CLUSTERING #######
# Figure A5

# Values
hvals <- c("None" = "resp_pair_s",
           "Pair" = "resp_pair", 
           "Respondent" = "id", 
           "PSU" = "PSU")
h <- "SE Clustering"

# Estimate
coef.df <- do.call(rbind, lapply(outvars, function(o){
  do.call(rbind, lapply(hvals, function(c){
    # Formula
    f <- as.formula(paste(o, " ~ ", 
                          paste(paste0("attr",rep(1:5, 2), "_d", rep(2:3, each = 5)),
                                collapse = " + ")))
    # Estimate
    m <- fixest::feols(f,
                       data = data.long[data.long$experiment %in% c(1:2) & data.long$in_2024,],
                       cluster = c)
    
    # Extract results
    data.frame(outcome = names(outvars)[outvars == o],
               hetvar = h,
               hetvar_val = c,
               hetvar_lab = names(hvals)[hvals == 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$hetvar_lab <- factor(coef.df$hetvar_lab,
                             levels = names(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, 
                col = hetvar_lab, 
                pch = hetvar_lab)) +
  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) ,
                     minor_breaks = seq(.5, 100.5, by = 1),
                     labels = as.character(coef.df$treat_lab),
                     position = "right") +
  theme_minimal() +
  labs(y = NULL,
       x = paste0("Average Marginal Component Effect "),
       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)

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

