Marginal Effects Plots

R Code for Marginal Effects Plots

library(readxl)
library(tidyverse)
library(patchwork)

# --- Load data ---
file_path <- "C:/Users/Michael/Desktop/Temp One Drive for R/ME Values.xlsx"
principal_data <- read_excel(file_path, sheet = 1)

# --- Clean and prep ---
principal_data$Variable <- trimws(principal_data$Variable)
principal_data$Category <- trimws(principal_data$Category)

# --- Extract base category (e.g., "Crash Characteristics") ---
principal_data <- principal_data %>%
  mutate(CategoryGroup = str_remove(Category, "\\s*\\d+$"))  # remove trailing numbers

# --- Pivot longer ---
plot_long <- principal_data %>%
  pivot_longer(cols = c(KA, BC, O), names_to = "Severity", values_to = "Effect") %>%
  mutate(
    Severity = factor(Severity, levels = c("O", "BC", "KA")),
    Category = factor(Category),  # e.g., "Crash Characteristics 1"
    CategoryGroup = factor(CategoryGroup)  # e.g., "Crash Characteristics"
  )

# --- Colors ---
severity_colors <- c("KA"="#C93312", "BC"="#FDDDA0", "O"="#74A089")

# --- Layout ---
ncol <- 3
categories <- unique(plot_long$Category)
nplots <- length(categories)
nrows <- ceiling(nplots / ncol)
bottom_row_start <- (nrows - 1) * ncol + 1
bottom_row_indices <- bottom_row_start:nplots
bottom_row_indices <- bottom_row_indices[bottom_row_indices <= nplots]

# --- Global y-limits with buffer ---
y_min <- min(plot_long$Effect, na.rm = TRUE)
y_max <- max(plot_long$Effect, na.rm = TRUE)
buffer <- 0.0001 * (y_max - y_min)
y_limits <- c(y_min - buffer, y_max + buffer)

# --- Create plots ---
plot_list <- lapply(seq_along(categories), function(i) {
  cat <- categories[i]
  df <- plot_long %>% filter(Category == cat)
  df$Variable <- factor(df$Variable, levels = rev(unique(df$Variable)))
  
  is_bottom_row <- i %in% bottom_row_indices
  x_lab <- NULL
  y_lab <- if(is_bottom_row) "Marginal Effect" else NULL
  group_title <- unique(df$CategoryGroup)
  
  ggplot(df, aes(x = Variable, y = Effect, fill = Severity)) +
    geom_col(position = position_dodge(width = 0.8), width = 0.8,
             color = "black", size = 0.25) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
    scale_fill_manual(values = severity_colors) +
    coord_flip() +
    scale_y_continuous(limits = y_limits, expand = c(0, 0)) +
    labs(title = group_title, x = x_lab, y = y_lab) +
    theme_bw(base_size = 10) +
    theme(
      plot.title = element_text(size=18, face="bold", hjust = 0.5),
      axis.title = element_text(size=18, color = "black", face = "bold"),
      axis.text.x = if(is_bottom_row) element_text(size = 16, color = "black") else element_blank(),
      axis.ticks.x = if(is_bottom_row) element_line() else element_blank(),
      axis.text.y = element_text(size = 13, color = "black"),
      legend.position = "none",
      panel.grid.major.y = element_line(color = "gray30"),
      panel.grid.minor.y = element_blank(),
      panel.border = element_blank()
    )
})

# --- Arrange plots dynamically ---
final_plot <- wrap_plots(plot_list, ncol = ncol)
final_plot

# --- Save ---
ggsave("marginal_effect_plots.png", plot = final_plot,
       path = "C:/Users/Michael/OneDrive - Texas State University/NLOGIT6/CRIS/New_Studies/Embankment_Michael/Plots",
       width = 18, height = 12, dpi = 300, bg = "white")