library(tidyr)
library(dplyr)
library(ggplot2)
library(scales)
# install.packages("patchwork")
# install.packages("cowplot")
library(grid)
library(patchwork)
library(cowplot)


# ======================================================
# Figure 1.Mid: shock U.S. only and North America
# ======================================================


# ======================================================
# 1. Sector information
# ======================================================

Sector_list <- c(
  "Cereals",
  "Fruit & Veg",
  "Fish & Seafood",
  "Meat & Meat products",
  "Veg Oil & Fats",
  "Dairy",
  "Sugar & Sugar products",
  "Bev & Tobacco",
  "Processed food"
)


sector_size <- data.frame(
  Sector = Sector_list,
  Share = c(
    0.007422,
    0.008407,
    0.002751,
    0.008285,
    0.005419,
    0.004928,
    0.002099,
    0.008368,
    0.021580
  )
)


# Order sectors from smallest to largest sector share
sector_order_key <- sector_size %>%
  arrange(Share) %>%
  pull(Sector)


# Sector colors
sector_colors <- setNames(
  scales::hue_pal()(nrow(sector_size)),
  sector_size$Sector
)



# ======================================================
# 2. Data input
# ======================================================


# ------------------------------------------------------
# Current clinical coverage - U.S. only
# ------------------------------------------------------

Mid_qo_us <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    -1.67, -2.05, -1.97, -3.28, -1.30,
    -3.60, -4.50, -4.42, -3.68
  )
)


Mid_po_us <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    -0.28, -0.50, -0.97, -0.22, -0.26,
    -0.13, -0.16, -0.12, -0.16
  )
)


Mid_qmw_us <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    -3.66, -3.96, -3.30, -4.65, -3.47,
    -4.60, -4.85, -4.63, -4.32
  )
)


Mid_qxw_us <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    1.17, 1.22, 1.70, 1.74, 1.30,
    1.04, 0.85, 0.29, 0.53
  )
)



# ------------------------------------------------------
# Current clinical coverage - North America
# ------------------------------------------------------

Mid_qo_na <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    -1.72, -2.14, -2.03, -3.31, -1.34,
    -3.62, -4.52, -4.45, -3.73
  )
)


Mid_po_na <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    -0.28, -0.52, -0.99, -0.23, -0.27,
    -0.13, -0.16, -0.12, -0.17
  )
)


Mid_qmw_na <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    -3.65, -3.96, -3.28, -4.61, -3.49,
    -4.61, -4.81, -4.63, -4.32
  )
)


Mid_qxw_na <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    1.04, 0.95, 1.44, 1.49, 1.24,
    0.76, 0.45, -0.09, 0.10
  )
)



# ------------------------------------------------------
# Expanded weight-management access - U.S. only
# ------------------------------------------------------

Over_qo_us <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    -2.97, -3.67, -3.54, -5.84, -2.31,
    -6.42, -8.02, -7.88, -6.56
  )
)


Over_po_us <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    -0.51, -0.90, -1.74, -0.41, -0.47,
    -0.24, -0.29, -0.22, -0.30
  )
)


Over_qmw_us <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    -6.49, -7.06, -5.87, -8.23, -6.14,
    -8.17, -8.61, -8.24, -7.68
  )
)


Over_qxw_us <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    2.09, 2.19, 3.07, 3.15, 2.33,
    1.90, 1.53, 0.52, 0.96
  )
)



# ------------------------------------------------------
# Expanded weight-management access - North America
# ------------------------------------------------------

Over_qo_na <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    -3.05, -3.81, -3.62, -5.88, -2.38,
    -6.44, -8.04, -7.91, -6.62
  )
)


Over_po_na <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    -0.52, -0.93, -1.77, -0.41, -0.48,
    -0.24, -0.30, -0.23, -0.31
  )
)


Over_qmw_na <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    -6.48, -7.06, -5.84, -8.17, -6.17,
    -8.18, -8.56, -8.24, -7.68
  )
)


Over_qxw_na <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    1.89, 1.79, 2.69, 2.77, 2.24,
    1.46, 0.94, -0.03, 0.34
  )
)



# ======================================================
# 3. Scenario definitions
# ======================================================

scenario_levels <- c(
  "Current: U.S. only",
  "Current: North America",
  "Expanded: U.S. only",
  "Expanded: North America"
)


# IMPORTANT:
#
# Use exactly the SAME geometry for filled/hollow pairs.
#
# Shape 21 = circle that accepts both color and fill.
# Shape 23 = diamond that accepts both color and fill.
#
# Therefore:
# filled/hollow circles are exactly the same size
# filled/hollow diamonds are exactly the same size.

scenario_shapes <- c(
  "Current: U.S. only"       = 21,
  "Current: North America"   = 21,
  "Expanded: U.S. only"      = 23,
  "Expanded: North America"  = 23
)


# Slight horizontal offsets so all four scenarios are visible
scenario_offsets <- c(
  "Current: U.S. only"       = -0.24,
  "Current: North America"   = -0.08,
  "Expanded: U.S. only"      =  0.08,
  "Expanded: North America"  =  0.24
)



# ======================================================
# 4. Combine data and convert to long format
# ======================================================

make_panel_df <- function(
    current_us,
    current_na,
    expanded_us,
    expanded_na
) {
  
  current_us %>%
    
    rename(
      Current_US = `2030`
    ) %>%
    
    left_join(
      current_na %>%
        rename(
          Current_NA = `2030`
        ),
      by = "Sector"
    ) %>%
    
    left_join(
      expanded_us %>%
        rename(
          Expanded_US = `2030`
        ),
      by = "Sector"
    ) %>%
    
    left_join(
      expanded_na %>%
        rename(
          Expanded_NA = `2030`
        ),
      by = "Sector"
    ) %>%
    
    # Attach ONE sector-share value to each sector
    left_join(
      sector_size,
      by = "Sector"
    ) %>%
    
    mutate(
      
      Sector = factor(
        Sector,
        levels = sector_order_key
      ),
      
      x_id = as.numeric(Sector)
      
    ) %>%
    
    pivot_longer(
      
      cols = c(
        Current_US,
        Current_NA,
        Expanded_US,
        Expanded_NA
      ),
      
      names_to = "Scenario_code",
      
      values_to = "Impact"
      
    ) %>%
    
    mutate(
      
      # ------------------------------------------------
      # Scenario
      # ------------------------------------------------
      
      Scenario = case_when(
        
        Scenario_code == "Current_US" ~
          "Current: U.S. only",
        
        Scenario_code == "Current_NA" ~
          "Current: North America",
        
        Scenario_code == "Expanded_US" ~
          "Expanded: U.S. only",
        
        Scenario_code == "Expanded_NA" ~
          "Expanded: North America"
        
      ),
      
      
      Scenario = factor(
        Scenario,
        levels = scenario_levels
      ),
      
      
      # ------------------------------------------------
      # Horizontal offset
      # ------------------------------------------------
      
      x_plot =
        x_id +
        unname(
          scenario_offsets[
            as.character(Scenario)
          ]
        ),
      
      
      # ------------------------------------------------
      # Fill
      #
      # U.S. only = filled with sector color
      # North America = hollow/white
      # ------------------------------------------------
      
      PointFill = case_when(
        
        Scenario_code %in% c(
          "Current_US",
          "Expanded_US"
        ) ~ as.character(Sector),
        
        Scenario_code %in% c(
          "Current_NA",
          "Expanded_NA"
        ) ~ "Hollow"
        
      )
      
    )
}



# ======================================================
# 5. Construct the four panel data sets
# ======================================================

Po_rng <- make_panel_df(
  Mid_po_us,
  Mid_po_na,
  Over_po_us,
  Over_po_na
)


Qo_rng <- make_panel_df(
  Mid_qo_us,
  Mid_qo_na,
  Over_qo_us,
  Over_qo_na
)


Qmw_rng <- make_panel_df(
  Mid_qmw_us,
  Mid_qmw_na,
  Over_qmw_us,
  Over_qmw_na
)


Qxw_rng <- make_panel_df(
  Mid_qxw_us,
  Mid_qxw_na,
  Over_qxw_us,
  Over_qxw_na
)



# ======================================================
# 5A. Optional check
#
# All four scenarios within each sector should have
# exactly the SAME Share.
# ======================================================

check_sector_size <- Po_rng %>%
  
  group_by(Sector) %>%
  
  summarise(
    
    Number_scenarios = n(),
    
    Number_unique_sizes =
      n_distinct(Share),
    
    Share =
      first(Share),
    
    .groups = "drop"
    
  )


print(check_sector_size)

# For every sector:
#
# Number_scenarios   = 4
# Number_unique_sizes = 1
#
# Thus:
# circle U.S.
# circle North America
# diamond U.S.
# diamond North America
#
# all use exactly the same Share for that sector.



# ======================================================
# 6. Plot function
# ======================================================

make_plot_fourtype <- function(
    df_rng,
    y_min,
    y_max
) {
  
  # ----------------------------------------------------
  # Wide data for connector lines
  # ----------------------------------------------------
  
  connector_df <- df_rng %>%
    
    select(
      Sector,
      x_id,
      Scenario_code,
      Impact
    ) %>%
    
    pivot_wider(
      names_from = Scenario_code,
      values_from = Impact
    )
  
  
  ggplot() +
    
    
    # --------------------------------------------------
  # Zero reference line
  # --------------------------------------------------
  
  geom_hline(
    
    yintercept = 0,
    
    linewidth = 0.9,
    
    color = "black"
    
  ) +
    
    
    # --------------------------------------------------
  # Current:
  # U.S. only -> North America
  # --------------------------------------------------
  
  geom_segment(
    
    data = connector_df,
    
    aes(
      x = x_id - 0.24,
      xend = x_id - 0.08,
      y = Current_US,
      yend = Current_NA
    ),
    
    color = "grey50",
    
    linewidth = 0.6,
    
    alpha = 0.80
    
  ) +
    
    
    # --------------------------------------------------
  # Expanded:
  # U.S. only -> North America
  # --------------------------------------------------
  
  geom_segment(
    
    data = connector_df,
    
    aes(
      x = x_id + 0.08,
      xend = x_id + 0.24,
      y = Expanded_US,
      yend = Expanded_NA
    ),
    
    color = "grey50",
    
    linewidth = 0.6,
    
    alpha = 0.80
    
  ) +
    
    
    # --------------------------------------------------
  # Four scenario estimates
  #
  # IMPORTANT:
  # size = Share DIRECTLY.
  #
  # Therefore all four scenarios for the same
  # sector receive exactly the same nominal size.
  # --------------------------------------------------
  
  geom_point(
    
    data = df_rng,
    
    aes(
      x = x_plot,
      y = Impact,
      size = Share,
      color = Sector,
      fill = PointFill,
      shape = Scenario
    ),
    
    alpha = 0.90,
    
    # Slightly thicker outline
    stroke = 1.7
    
  ) +
    
    
    # --------------------------------------------------
  # X axis
  # Sector names hidden
  # --------------------------------------------------
  
  scale_x_continuous(
    
    breaks =
      seq_along(
        sector_order_key
      ),
    
    labels = NULL,
    
    expand =
      expansion(
        mult = c(
          0.04,
          0.04
        )
      )
    
  ) +
    
    
    # --------------------------------------------------
  # Y limits
  # --------------------------------------------------
  
  coord_cartesian(
    
    ylim = c(
      y_min,
      y_max
    ),
    
    clip = "off"
    
  ) +
    
    
    # --------------------------------------------------
  # Shapes
  #
  # Circle pair uses shape 21.
  # Diamond pair uses shape 23.
  #
  # Fill distinguishes U.S. vs North America.
  # --------------------------------------------------
  
  scale_shape_manual(
    
    values =
      scenario_shapes,
    
    breaks =
      scenario_levels,
    
    name =
      "GLP-1 adoption scenario",
    
    guide =
      guide_legend(
        
        order = 1,
        
        nrow = 2,
        
        byrow = TRUE,
        
        override.aes = list(
          
          # Same legend size for all four
          size = 8,
          
          color = "black",
          
          fill = c(
            "black",
            "NA",
            "black",
            "NA"
          ),
          
          stroke = 1.7,
          
          alpha = 0.90
          
        )
        
      )
    
  ) +
    
    
    # --------------------------------------------------
  # Sector colors
  # --------------------------------------------------
  
  scale_color_manual(
    
    values =
      sector_colors,
    
    breaks =
      sector_order_key,
    
    name =
      "Sector",
    
    guide =
      guide_legend(
        
        order = 2,
        
        nrow = 2,
        
        byrow = TRUE,
        
        override.aes = list(
          
          # Bigger colored sector examples
          size = 10,
          
          shape = 21,
          
          fill =
            unname(
              sector_colors[
                sector_order_key
              ]
            ),
          
          alpha = 0.90,
          
          stroke = 1
          
        )
        
      )
    
  ) +
    
    
    # --------------------------------------------------
  # Fill colors
  #
  # U.S. = sector color
  # North America = white
  # --------------------------------------------------
  
  scale_fill_manual(
    
    values = c(
      sector_colors,
      "Hollow" = "NA"
    ),
    
    guide = "none"
    
  ) +
    
    
    # --------------------------------------------------
  # Sector-share marker size
  #
  # All four scenario markers use the SAME Share.
  #
  # max_size = 23 makes markers larger in plots.
  #
  # 2.00% example removed from legend.
  # --------------------------------------------------
  
  scale_size_area(
    
    max_size = 23,
    
    breaks = c(
      0.0025,
      0.0050,
      0.0100,
      0.0150
    ),
    
    labels = c(
      "0.25%",
      "0.50%",
      "1.00%",
      "1.50%"
    ),
    
    name =
      "Sector share (2030)",
    
    guide =
      guide_legend(
        
        order = 3,
        
        nrow = 1,
        
        override.aes = list(
          
          shape = 21,
          
          color = "grey30",
          
          fill = "grey30",
          
          alpha = 0.90,
          
          stroke = 1
          
        )
        
      )
    
  ) +
    
    
    # --------------------------------------------------
  # Labels
  # --------------------------------------------------
  
  labs(
    x = NULL,
    y = NULL
  ) +
    
    
    # --------------------------------------------------
  # Theme
  # --------------------------------------------------
  
  theme_minimal(
    base_size = 46
  ) +
    
    
    theme(
      
      aspect.ratio = 1,
      
      
      # ------------------------------------------------
      # No sector names below panels
      # ------------------------------------------------
      
      axis.text.x =
        element_blank(),
      
      axis.ticks.x =
        element_blank(),
      
      
      # ------------------------------------------------
      # Y-axis numbers
      # ------------------------------------------------
      
      axis.text.y =
        element_text(
          
          size = 32,
          
          color = "black"
          
        ),
      
      
      # ------------------------------------------------
      # Grid
      # ------------------------------------------------
      
      panel.grid.major.x =
        element_line(
          
          linewidth = 0.5,
          
          color = "grey88"
          
        ),
      
      panel.grid.minor.x =
        element_blank(),
      
      
      # ------------------------------------------------
      # Legend
      # ------------------------------------------------
      
      legend.key.width =
        unit(
          1.7,
          "lines"
        ),
      
      legend.key.height =
        unit(
          1.7,
          "lines"
        ),
      
      legend.text =
        element_text(
          
          size = 32,
          
          color = "black"
          
        ),
      
      legend.title =
        element_text(
          
          size = 35,
          
          face = "bold",
          
          hjust = 0,
          
          color = "black"
          
        ),
      
      legend.title.align = 0,
      
      legend.spacing.y =
        unit(
          0.25,
          "cm"
        ),
      
      
      # ------------------------------------------------
      # Margin
      # ------------------------------------------------
      
      plot.margin =
        margin(
          10,
          10,
          10,
          10
        )
      
    )
  
}



# ======================================================
# 7. Panel label helper
# ======================================================

panel <- function(
    p,
    title
) {
  
  label_strip <- ggplot() +
    
    
    annotate(
      
      "rect",
      
      xmin = -Inf,
      
      xmax = Inf,
      
      ymin = -Inf,
      
      ymax = Inf,
      
      fill = "grey90",
      
      color = NA
      
    ) +
    
    
    annotate(
      
      "text",
      
      x = 0.5,
      
      y = 0.5,
      
      label = title,
      
      # Smaller Panel A-D title
      size = 12,
      
      fontface = "bold"
      
    ) +
    
    
    scale_x_continuous(
      
      limits = c(
        0,
        1
      ),
      
      expand = c(
        0,
        0
      )
      
    ) +
    
    
    scale_y_continuous(
      
      limits = c(
        0,
        1
      ),
      
      expand = c(
        0,
        0
      )
      
    ) +
    
    
    theme_void() +
    
    
    theme(
      
      plot.margin =
        margin(
          0,
          0,
          0,
          0
        )
      
    )
  
  
  (p / label_strip) +
    
    plot_layout(
      
      heights = c(
        1,
        0.11
      )
      
    )
  
}



# ======================================================
# 8. Build four panels
# ======================================================

pA <- panel(
  
  make_plot_fourtype(
    Po_rng,
    -10,
    5
  ),
  
  "Panel A: Farmgate price"
  
)


pB <- panel(
  
  make_plot_fourtype(
    Qo_rng,
    -10,
    5
  ),
  
  "Panel B: Domestic production volume"
  
)


pC <- panel(
  
  make_plot_fourtype(
    Qmw_rng,
    -10,
    5
  ),
  
  "Panel C: Import volume"
  
)


pD <- panel(
  
  make_plot_fourtype(
    Qxw_rng,
    -10,
    5
  ),
  
  "Panel D: Export volume"
  
)



# ======================================================
# 9. Final figure
# ======================================================

make_final_figure_four_panels <- function(
    pA,
    pB,
    pC,
    pD
) {
  
  panel_grid <- (
    
    (pA | pB) /
      
      (pC | pD)
    
  ) +
    
    
    plot_layout(
      guides = "collect"
    ) &
    
    
    theme(
      
      legend.position =
        "bottom",
      
      legend.box =
        "vertical",
      
      legend.box.spacing =
        unit(
          0.6,
          "cm"
        )
      
    )
  
  
  # ----------------------------------------------------
  # Common Y-axis title
  # ----------------------------------------------------
  
  y_label <- ggplot() +
    
    
    annotate(
      
      "text",
      
      x = 0.5,
      
      y = 0.5,
      
      label =
        "Sectoral Impact (%)",
      
      angle = 90,
      
      size = 16,
      
      fontface = "bold"
      
    ) +
    
    
    scale_x_continuous(
      
      limits = c(
        0,
        1
      ),
      
      expand = c(
        0,
        0
      )
      
    ) +
    
    
    scale_y_continuous(
      
      limits = c(
        0,
        1
      ),
      
      expand = c(
        0,
        0
      )
      
    ) +
    
    
    theme_void()
  
  
  # ----------------------------------------------------
  # Combine
  # ----------------------------------------------------
  
  plot_grid(
    
    y_label,
    
    panel_grid,
    
    ncol = 2,
    
    rel_widths = c(
      0.05,
      0.95
    )
    
  )
  
}



# ======================================================
# 10. Show final figure
# ======================================================

Figure_current_expanded <-
  make_final_figure_four_panels(
    pA,
    pB,
    pC,
    pD
  )


Figure_current_expanded









#############################################
# Figure: Land-use change
#
# Bubble center = Current clinical coverage, Shock US
# Bubble size = same sector share as Figure 1
# Open circle = Current clinical coverage, Shock NA
# Black short bar = Expanded weight-management, Shock US
# Gray short bar = Expanded weight-management, Shock NA
# Vertical line links the two expanded weight-management bars
#############################################


# ======================================================
# 1. Land-use-change data
# Uses the same nine food sectors as Figure 1
# ======================================================

Land_mid_us <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(-1.12,-1.49,-1.65,-0.96,-0.62,-1.12,-1.78,-1.52,-1.24)
)
Land_mid_na <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(-1.15,-1.56,-1.70,-0.96,-0.64,-1.10,-1.77,-1.51,-1.24)
)
Land_over_us <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(-1.99,-2.68,-2.96,-1.73,-1.10,-2.01,-3.22,-2.74,-2.23)
)
Land_over_na <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(-2.05,-2.79,-3.03,-1.72,-1.13,-2.00,-3.20,-2.73,-2.24)
)
scenario_levels <- c("Current: U.S. only","Current: North America","Expanded: U.S. only","Expanded: North America")
scenario_shapes <- c(
  "Current: U.S. only" = 21,
  "Current: North America" = 21,
  "Expanded: U.S. only" = 23,
  "Expanded: North America" = 23
)
scenario_offsets <- c(
  "Current: U.S. only" = -0.24,
  "Current: North America" = -0.08,
  "Expanded: U.S. only" = 0.08,
  "Expanded: North America" = 0.24
)
Land_rng <- Land_mid_us %>%
  rename(Current_US = `2030`) %>%
  left_join(Land_mid_na %>% rename(Current_NA = `2030`),by = "Sector") %>%
  left_join(Land_over_us %>% rename(Expanded_US = `2030`),by = "Sector") %>%
  left_join(Land_over_na %>% rename(Expanded_NA = `2030`),by = "Sector") %>%
  left_join(sector_size,by = "Sector") %>%
  mutate(
    Sector = factor(Sector,levels = sector_order_key),
    x_id = as.numeric(Sector)
  ) %>%
  pivot_longer(
    cols = c(Current_US,Current_NA,Expanded_US,Expanded_NA),
    names_to = "Scenario_code",
    values_to = "Impact"
  ) %>%
  mutate(
    Scenario = case_when(
      Scenario_code == "Current_US" ~ "Current: U.S. only",
      Scenario_code == "Current_NA" ~ "Current: North America",
      Scenario_code == "Expanded_US" ~ "Expanded: U.S. only",
      Scenario_code == "Expanded_NA" ~ "Expanded: North America"
    ),
    Scenario = factor(Scenario,levels = scenario_levels),
    x_plot = x_id + unname(scenario_offsets[as.character(Scenario)]),
    PointFill = case_when(
      Scenario_code %in% c("Current_US","Expanded_US") ~ as.character(Sector),
      Scenario_code %in% c("Current_NA","Expanded_NA") ~ "Hollow"
    )
  )
check_land_sector_size <- Land_rng %>%
  group_by(Sector) %>%
  summarise(
    Number_scenarios = n(),
    Number_unique_sizes = n_distinct(Share),
    Sector_share = first(Share),
    .groups = "drop"
  )
print(check_land_sector_size)
make_land_plot <- function(df_rng,y_min,y_max) {
  connector_df <- df_rng %>%
    select(Sector,x_id,Scenario_code,Impact) %>%
    pivot_wider(names_from = Scenario_code,values_from = Impact)
  ggplot() +
    geom_hline(yintercept = 0,linewidth = 0.9,color = "black") +
    geom_segment(
      data = connector_df,
      aes(x = x_id - 0.24,xend = x_id - 0.08,y = Current_US,yend = Current_NA),
      color = "grey50",linewidth = 0.6,alpha = 0.80
    ) +
    geom_segment(
      data = connector_df,
      aes(x = x_id + 0.08,xend = x_id + 0.24,y = Expanded_US,yend = Expanded_NA),
      color = "grey50",linewidth = 0.6,alpha = 0.80
    ) +
    geom_point(
      data = df_rng,
      aes(x = x_plot,y = Impact,size = Share,color = Sector,fill = PointFill,shape = Scenario),
      alpha = 0.90,stroke = 1.7
    ) +
    scale_x_continuous(
      breaks = seq_along(sector_order_key),
      labels = NULL,
      expand = expansion(mult = c(0.04,0.04))
    ) +
    coord_cartesian(ylim = c(y_min,y_max),clip = "off") +
    scale_shape_manual(
      values = scenario_shapes,
      breaks = scenario_levels,
      name = "GLP-1 adoption scenario",
      guide = guide_legend(
        order = 1,
        nrow = 2,
        byrow = TRUE,
        override.aes = list(
          size = 8,
          color = "black",
          fill = c("black",NA,"black",NA),
          stroke = 1.7,
          alpha = 0.90
        )
      )
    ) +
    scale_color_manual(
      values = sector_colors,
      breaks = sector_order_key,
      name = "Sector",
      guide = guide_legend(
        order = 2,
        nrow = 2,
        byrow = TRUE,
        override.aes = list(
          size = 10,
          shape = 21,
          fill = unname(sector_colors[sector_order_key]),
          alpha = 0.90,
          stroke = 1
        )
      )
    ) +
    scale_fill_manual(
      values = c(sector_colors,"Hollow" = NA),
      guide = "none"
    ) +
    scale_size_area(
      max_size = 23,
      breaks = c(0.0025,0.0050,0.0100,0.0150),
      labels = c("0.25%","0.50%","1.00%","1.50%"),
      name = "Sector share (2030)",
      guide = guide_legend(
        order = 3,
        nrow = 1,
        override.aes = list(
          shape = 21,
          color = "grey30",
          fill = "grey30",
          alpha = 0.90,
          stroke = 1
        )
      )
    ) +
    labs(x = NULL,y = NULL) +
    theme_minimal(base_size = 46) +
    theme(
      aspect.ratio = 1,
      axis.text.x = element_blank(),
      axis.ticks.x = element_blank(),
      axis.title.x = element_blank(),
      axis.text.y = element_text(size = 32,color = "black"),
      panel.grid.major.x = element_line(linewidth = 0.5,color = "grey88"),
      panel.grid.minor.x = element_blank(),
      legend.key.width = unit(1.7,"lines"),
      legend.key.height = unit(1.7,"lines"),
      legend.text = element_text(size = 32,color = "black"),
      legend.title = element_text(size = 35,face = "bold",hjust = 0,color = "black"),
      legend.title.align = 0,
      legend.spacing.y = unit(0.25,"cm"),
      plot.margin = margin(10,10,10,10)
    )
}
p_land <- make_land_plot(
  Land_rng,
  y_min = -4,
  y_max = 1
)
land_panel_grid <- p_land +
  theme(
    legend.position = "bottom",
    legend.box = "vertical",
    legend.box.spacing = unit(0.6,"cm")
  )
land_y_label <- ggplot() +
  annotate(
    "text",
    x = 0.5,
    y = 0.5,
    label = "Land-use Change (%)",
    angle = 90,
    size = 16,
    fontface = "bold"
  ) +
  scale_x_continuous(limits = c(0,1),expand = c(0,0)) +
  scale_y_continuous(limits = c(0,1),expand = c(0,0)) +
  theme_void()
Figure_land_use <- plot_grid(
  land_y_label,
  land_panel_grid,
  ncol = 2,
  rel_widths = c(0.05,0.95)
)
Figure_land_use









#############################################
# Figure 2. Decomposed U.S. Trade Difference goes to Appendix
# Difference = Shock NA - Shock US
# Obesity: bubble center
# Overweight: open circle
# Mid pathway only
#############################################

# ======================================================
# 0. Required packages
# ======================================================
library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)
library(patchwork)
library(cowplot)
library(grid)

# ======================================================
# 1. Sector information
# ======================================================

Sector_list <- c("Cereals","Fruit & Veg","Fish & Seafood","Meat & Meat products","Veg Oil & Fats","Dairy","Sugar & Sugar products","Bev & Tobacco","Processed food")
sector_size <- data.frame(Sector = Sector_list,Share = c(0.007422,0.008407,0.002751,0.008285,0.005419,0.004928,0.002099,0.008368,0.021580))
sector_order_key <- sector_size %>% arrange(Share) %>% pull(Sector)
sector_colors <- setNames(scales::hue_pal()(nrow(sector_size)),sector_size$Sector)

# ======================================================
# 2. Obesity results
# Difference = Shock NA - Shock US
# ======================================================

# Export from U.S. to CA, MX, ROW
Obesity_export_diff <- data.frame(
  Sector = Sector_list,
  CA  = c(-0.66, -0.89, -0.74, -0.81, -0.54, -1.28, -1.32, -1.33, -1.23),
  MX  = c(-0.70, -0.58, -0.73, -0.95, -0.66, -0.79, -0.88, -0.77, -0.72),
  ROW = c( 0.03,  0.05,  0.06,  0.03,  0.03,  0.02,  0.02,  0.01,  0.02),
  Type = "Export",
  check.names = FALSE
)

# Import to U.S. from CA, MX, ROW
Obesity_import_diff <- data.frame(
  Sector = Sector_list,
  CA  = c( 0.09, -0.01,  0.20,  0.18,  0.07,  0.17,  0.06,  0.04,  0.07),
  MX  = c( 0.48,  0.10,  0.30,  0.40,  0.12,  0.20,  0.27,  0.05,  0.12),
  ROW = c(-0.08, -0.11, -0.14, -0.11, -0.08, -0.03, -0.07, -0.02, -0.04),
  Type = "Import",
  check.names = FALSE
)


# ======================================================
# 3. Overweight results
# Difference = Shock NA - Shock US
# ======================================================

# Export from U.S. to CA, MX, ROW
Overweight_export_diff <- data.frame(
  Sector = Sector_list,
  CA  = c(-0.94, -1.28, -1.07, -1.17, -0.77, -1.84, -1.90, -1.91, -1.77),
  MX  = c(-1.06, -0.91, -1.16, -1.47, -1.01, -1.23, -1.36, -1.17, -1.10),
  ROW = c( 0.04,  0.08,  0.08,  0.05,  0.04,  0.04,  0.03,  0.01,  0.02),
  Type = "Export",
  check.names = FALSE
)

# Import to U.S. from CA, MX, ROW
Overweight_import_diff <- data.frame(
  Sector = Sector_list,
  CA  = c( 0.12, -0.01,  0.28,  0.25,  0.09,  0.24,  0.07,  0.05,  0.09),
  MX  = c( 0.73,  0.15,  0.46,  0.60,  0.18,  0.29,  0.41,  0.08,  0.18),
  ROW = c(-0.12, -0.15, -0.21, -0.15, -0.12, -0.05, -0.11, -0.02, -0.05),
  Type = "Import",
  check.names = FALSE
)

# ======================================================
# 4. Long format and combine obesity and overweight
# ======================================================

to_long_obesity <- function(df) {df %>% pivot_longer(cols = c(CA,MX,ROW),names_to = "Partner",values_to = "Current") %>% select(Sector,Type,Partner,Current)}
to_long_overweight <- function(df) {df %>% pivot_longer(cols = c(CA,MX,ROW),names_to = "Partner",values_to = "Expanded") %>% select(Sector,Type,Partner,Expanded)}
obesity_long <- bind_rows(to_long_obesity(Obesity_export_diff),to_long_obesity(Obesity_import_diff))
overweight_long <- bind_rows(to_long_overweight(Overweight_export_diff),to_long_overweight(Overweight_import_diff))
trade_mid <- obesity_long %>% left_join(overweight_long,by = c("Sector","Type","Partner")) %>% left_join(sector_size,by = "Sector") %>% mutate(Sector = factor(Sector,levels = sector_order_key),Partner = factor(Partner,levels = c("CA","MX","ROW")),Type = factor(Type,levels = c("Import","Export"))) %>% arrange(Type,Partner,Sector)


# ======================================================
# 5. Plot function
# Bubble center = obesity result
# Bubble size = relative sector size
# Open circle = overweight result
# ======================================================

scenario_levels <- c("Current clinical coverage","Expanded weight-management access")
scenario_shapes <- c("Current clinical coverage" = 21,"Expanded weight-management access" = 23)
scenario_offsets <- c("Current clinical coverage" = -0.14,"Expanded weight-management access" = 0.14)
make_trade_mid_plot <- function(df,partner,type,y_min,y_max) {
  df2 <- df %>% filter(Partner == partner,Type == type) %>% mutate(Sector = factor(Sector,levels = sector_order_key),x_id = as.numeric(Sector))
  plot_df <- df2 %>%
    pivot_longer(cols = c(Current,Expanded),names_to = "Scenario_code",values_to = "Impact") %>%
    mutate(
      Scenario = case_when(
        Scenario_code == "Current" ~ "Current clinical coverage",
        Scenario_code == "Expanded" ~ "Expanded weight-management access"
      ),
      Scenario = factor(Scenario,levels = scenario_levels),
      x_plot = x_id + unname(scenario_offsets[as.character(Scenario)]),
      PlotSize = case_when(
        Scenario_code == "Current" ~ Share,
        Scenario_code == "Expanded" ~ Share
      )
    )
  
  ggplot() +
    geom_hline(yintercept = 0,linewidth = 1,color = "black") +
    geom_segment(data = df2,aes(x = x_id - 0.14,xend = x_id + 0.14,y = Current,yend = Expanded),color = "grey50",linewidth = 0.6,alpha = 0.80) +
    geom_point(
      data = plot_df,
      aes(
        x = x_plot,
        y = Impact,
        size = PlotSize,
        color = Sector,
        shape = Scenario
      ),
      fill = NA,
      alpha = 0.90,
      stroke = 1.7
    ) +
    scale_x_continuous(breaks = seq_along(sector_order_key),labels = NULL,expand = expansion(mult = c(0.04,0.04))) +
    coord_cartesian(ylim = c(y_min,y_max),clip = "off") +
    theme(aspect.ratio = 1) +
    scale_shape_manual(
      values = scenario_shapes,
      breaks = scenario_levels,
      name = "Coverage scenario",
      guide = guide_legend(
        order = 1,
        nrow = 1,
        byrow = TRUE,
        override.aes = list(
          size = 8,
          color = "black",
          fill = NA,
          stroke = 1.7,
          alpha = 0.90
        )
      )
    ) +
    scale_color_manual(values = sector_colors,breaks = sector_order_key,name = "Sector",guide = guide_legend(order = 2,nrow = 2,byrow = TRUE,override.aes = list(size = 10,shape = 21,fill = unname(sector_colors[sector_order_key]),alpha = 0.90,stroke = 1))) +
    scale_size_area(max_size = 23,breaks = c(0.0025,0.0050,0.0100,0.0150),labels = c("0.25%","0.50%","1.00%","1.50%"),name = "Sector share (2030)",guide = guide_legend(order = 3,nrow = 1,override.aes = list(shape = 21,color = "grey30",fill = "grey30",alpha = 0.90,stroke = 1))) +
    labs(x = NULL,y = NULL) +
    theme_minimal(base_size = 40) +
    theme(axis.text.x = element_blank(),axis.ticks.x = element_blank(),axis.title.x = element_blank(),axis.text.y = element_text(size = 30,color = "black"),panel.grid.major.x = element_line(linewidth = 0.5,color = "grey88"),panel.grid.minor.x = element_blank(),legend.key.width = unit(1.7,"lines"),legend.key.height = unit(1.7,"lines"),legend.text = element_text(size = 32),legend.title = element_text(size = 34,face = "bold"),plot.margin = margin(10,10,10,10))
}
# ======================================================
# 6. Panel wrapper with grey label bar
# Keep original plot format
# ======================================================

panel_diff <- function(plot_obj,label_text) {
  label_plot <- ggplot(data.frame(x = 0.5,y = 0.5),aes(x,y)) +
    geom_rect(aes(xmin = 0,xmax = 1,ymin = 0,ymax = 1),fill = "grey90",color = NA) +
    geom_text(label = label_text,size = 12,fontface = "bold",hjust = 0.5,vjust = 0.5) +
    scale_x_continuous(limits = c(0,1),expand = c(0,0)) +
    scale_y_continuous(limits = c(0,1),expand = c(0,0)) +
    theme_void()
  plot_obj / label_plot + plot_layout(heights = c(12,2))
}
pI_CA <- panel_diff(make_trade_mid_plot(trade_mid,partner = "CA",type = "Import",y_min = -5,y_max = 5),"Import from Canada")
pE_CA <- panel_diff(make_trade_mid_plot(trade_mid,partner = "CA",type = "Export",y_min = -5,y_max = 5),"Export to Canada")
pI_MX <- panel_diff(make_trade_mid_plot(trade_mid,partner = "MX",type = "Import",y_min = -5,y_max = 5),"Import from Mexico")
pE_MX <- panel_diff(make_trade_mid_plot(trade_mid,partner = "MX",type = "Export",y_min = -5,y_max = 5),"Export to Mexico")
pI_ROW <- panel_diff(make_trade_mid_plot(trade_mid,partner = "ROW",type = "Import",y_min = -10,y_max = 5),"Import from Rest of World")
pE_ROW <- panel_diff(make_trade_mid_plot(trade_mid,partner = "ROW",type = "Export",y_min = -10,y_max = 5),"Export to Rest of World")
panel_grid_2 <- ((pI_CA | pE_CA) / (pI_MX | pE_MX) / (pI_ROW | pE_ROW)) + plot_layout(guides = "collect") & theme(legend.position = "bottom",legend.box = "vertical",legend.box.spacing = unit(0.6,"cm"))
y_label_2 <- ggplot() + annotate("text",x = 0.5,y = 0.5,label = "U.S. Trade Response Differences (percentage points)",angle = 90,size = 14,fontface = "bold") + scale_x_continuous(limits = c(0,1),expand = c(0,0)) + scale_y_continuous(limits = c(0,1),expand = c(0,0)) + theme_void()
figure2_final <- plot_grid(y_label_2,panel_grid_2,ncol = 2,rel_widths = c(0.04,0.96))
figure2_final








#############################################
# Figure: GHG and CO2 Emissions by Activities
# Bubble center = Obesity, Shock US
# Bubble size = sector share
# Open circle = Obesity, Shock NA
# Short bars = Overweight, Shock US and Shock NA
# Vertical line links the two overweight bars
#############################################

# ======================================================
# 0. Required packages
# ======================================================
library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)
library(patchwork)
library(cowplot)
library(grid)

# ======================================================
# 1. Sector information
# ======================================================

Sector_list <- c(
  "Cereals",
  "Fruit & Veg",
  "Fish & Seafood",
  "Meat & Meat products",
  "Veg Oil & Fats",
  "Dairy",
  "Sugar & Sugar products",
  "Bev & Tobacco",
  "Processed food"
)

sector_size <- data.frame(
  Sector = Sector_list,
  
  # CO2 sectoral shares from Excel table
  CO2_share = c(
    0.002572671,  # Cereal
    0.000902421,  # Veg_Fru
    0.000404898,  # Fish
    0.002338355,  # Meat
    0.002401265,  # Vegoil_F
    0.001744601,  # Dairy
    0.001547423,  # Sugar
    0.001911849,  # Beve_Tob
    0.007876525   # Oth_Fod
  ),
  
  # Total GHG sectoral shares from Excel table
  GHG_share = c(
    0.018127851,  # Cereal
    0.004815414,  # Veg_Fru
    0.000309525,  # Fish
    0.001740481,  # Meat
    0.006855848,  # Vegoil_F
    0.001310898,  # Dairy
    0.001148547,  # Sugar
    0.001473797,  # Beve_Tob
    0.008712234   # Oth_Fod
  )
)

sector_legend_order <- sector_size %>%
  arrange(GHG_share) %>%
  pull(Sector)

# Keep colors fixed by sector name across both plots
sector_colors <- setNames(
  scales::hue_pal()(length(Sector_list)),
  Sector_list
)
# ======================================================
# 2. Select year to plot
# ======================================================

Year_used <- "2030"

# ======================================================
# 3. Data input: GHG emissions
# ======================================================

# ------------------------------
# Obesity, Shock US, GHG
# ------------------------------

GHG_obesity_us <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2026` = c(-1.07, -1.31, -1.36, -1.95, -0.83, -2.13, -2.15, -2.61, -2.09),
  `2027` = c(-1.22, -1.49, -1.56, -2.25, -0.95, -2.46, -2.50, -3.01, -2.42),
  `2028` = c(-1.40, -1.73, -1.82, -2.62, -1.09, -2.87, -2.93, -3.52, -2.82),
  `2029` = c(-1.60, -1.99, -2.10, -3.04, -1.25, -3.33, -3.42, -4.09, -3.28),
  `2030` = c(-1.71, -2.13, -2.26, -3.28, -1.34, -3.61, -3.71, -4.43, -3.55)
)

# ------------------------------
# Obesity, Shock NA, GHG
# ------------------------------

GHG_obesity_na <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2026` = c(-1.09, -1.34, -1.39, -1.96, -0.85, -2.13, -2.15, -2.62, -2.11),
  `2027` = c(-1.25, -1.54, -1.60, -2.26, -0.97, -2.46, -2.50, -3.02, -2.44),
  `2028` = c(-1.44, -1.80, -1.86, -2.64, -1.12, -2.88, -2.93, -3.54, -2.85),
  `2029` = c(-1.65, -2.07, -2.15, -3.06, -1.29, -3.34, -3.42, -4.11, -3.31),
  `2030` = c(-1.76, -2.23, -2.33, -3.31, -1.38, -3.62, -3.72, -4.46, -3.59)
)

# ------------------------------
# Overweight, Shock US, GHG
# ------------------------------

GHG_overweight_us <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2026` = c(-1.18, -1.43, -1.50, -2.15, -0.91, -2.34, -2.37, -2.86, -2.30),
  `2027` = c(-1.53, -1.87, -1.96, -2.81, -1.19, -3.07, -3.12, -3.77, -3.02),
  `2028` = c(-1.96, -2.43, -2.55, -3.67, -1.53, -4.02, -4.10, -4.92, -3.95),
  `2029` = c(-2.50, -3.12, -3.29, -4.75, -1.96, -5.21, -5.34, -6.39, -5.13),
  `2030` = c(-3.04, -3.82, -4.04, -5.84, -2.38, -6.43, -6.63, -7.89, -6.33)
)

# ------------------------------
# Overweight, Shock NA, GHG
# ------------------------------

GHG_overweight_na <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2026` = c(-1.20, -1.48, -1.53, -2.16, -0.93, -2.35, -2.37, -2.88, -2.32),
  `2027` = c(-1.56, -1.94, -2.00, -2.83, -1.21, -3.08, -3.13, -3.78, -3.05),
  `2028` = c(-2.01, -2.51, -2.61, -3.69, -1.56, -4.03, -4.11, -4.94, -3.99),
  `2029` = c(-2.57, -3.24, -3.37, -4.78, -2.00, -5.23, -5.35, -6.42, -5.17),
  `2030` = c(-3.12, -3.97, -4.13, -5.89, -2.44, -6.45, -6.64, -7.93, -6.39)
)

# ======================================================
# 4. Data input: CO2 emissions
# ======================================================

# ------------------------------
# Obesity, Shock US, CO2
# ------------------------------

CO2_obesity_us <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2026` = c(-1.14, -1.38, -1.37, -1.96, -0.91, -2.14, -2.72, -2.61, -2.21),
  `2027` = c(-1.29, -1.58, -1.57, -2.26, -1.04, -2.47, -3.13, -3.02, -2.55),
  `2028` = c(-1.49, -1.83, -1.83, -2.63, -1.20, -2.88, -3.66, -3.53, -2.98),
  `2029` = c(-1.70, -2.10, -2.11, -3.05, -1.37, -3.34, -4.25, -4.10, -3.45),
  `2030` = c(-1.81, -2.25, -2.27, -3.29, -1.47, -3.62, -4.60, -4.44, -3.73)
)

# ------------------------------
# Obesity, Shock NA, CO2
# ------------------------------

CO2_obesity_na <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2026` = c(-1.16, -1.42, -1.39, -1.97, -0.93, -2.14, -2.72, -2.62, -2.23),
  `2027` = c(-1.32, -1.63, -1.60, -2.27, -1.06, -2.47, -3.14, -3.03, -2.57),
  `2028` = c(-1.52, -1.90, -1.87, -2.65, -1.22, -2.89, -3.67, -3.55, -3.00),
  `2029` = c(-1.75, -2.19, -2.16, -3.07, -1.41, -3.35, -4.26, -4.12, -3.48),
  `2030` = c(-1.87, -2.36, -2.33, -3.32, -1.51, -3.63, -4.61, -4.46, -3.77)
)

# ------------------------------
# Overweight, Shock US, CO2
# ------------------------------

CO2_overweight_us <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2026` = c(-1.25, -1.52, -1.50, -2.16, -1.00, -2.35, -2.99, -2.87, -2.43),
  `2027` = c(-1.62, -1.98, -1.96, -2.83, -1.30, -3.09, -3.92, -3.78, -3.19),
  `2028` = c(-2.08, -2.56, -2.56, -3.68, -1.67, -4.03, -5.12, -4.94, -4.16),
  `2029` = c(-2.66, -3.30, -3.30, -4.76, -2.14, -5.22, -6.64, -6.40, -5.39),
  `2030` = c(-3.22, -4.04, -4.05, -5.86, -2.60, -6.44, -8.19, -7.91, -6.65)
)

# ------------------------------
# Overweight, Shock NA, CO2
# ------------------------------

CO2_overweight_na <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2026` = c(-1.27, -1.57, -1.53, -2.17, -1.02, -2.36, -2.99, -2.88, -2.45),
  `2027` = c(-1.65, -2.05, -2.01, -2.84, -1.32, -3.09, -3.93, -3.79, -3.22),
  `2028` = c(-2.13, -2.66, -2.62, -3.71, -1.71, -4.04, -5.13, -4.96, -4.20),
  `2029` = c(-2.72, -3.42, -3.38, -4.79, -2.19, -5.24, -6.65, -6.43, -5.44),
  `2030` = c(-3.31, -4.19, -4.15, -5.90, -2.67, -6.46, -8.21, -7.94, -6.71)
)


# ======================================================
# 5. Combine data for each emission type
# ======================================================

scenario_levels <- c("Current: U.S. only","Current: North America","Expanded: U.S. only","Expanded: North America")
scenario_shapes <- c("Current: U.S. only" = 21,"Current: North America" = 21,"Expanded: U.S. only" = 23,"Expanded: North America" = 23)
scenario_offsets <- c("Current: U.S. only" = -0.24,"Current: North America" = -0.08,"Expanded: U.S. only" = 0.08,"Expanded: North America" = 0.24)
make_emission_panel_df <- function(obesity_us,obesity_na,overweight_us,overweight_na,emission_type,share_column) {
  obesity_us %>%
    select(Sector,all_of(Year_used)) %>%
    rename(Current_US = all_of(Year_used)) %>%
    left_join(obesity_na %>% select(Sector,all_of(Year_used)) %>% rename(Current_NA = all_of(Year_used)),by = "Sector") %>%
    left_join(overweight_us %>% select(Sector,all_of(Year_used)) %>% rename(Expanded_US = all_of(Year_used)),by = "Sector") %>%
    left_join(overweight_na %>% select(Sector,all_of(Year_used)) %>% rename(Expanded_NA = all_of(Year_used)),by = "Sector") %>%
    left_join(sector_size,by = "Sector") %>%
    mutate(Share = .data[[share_column]],Emission = emission_type) %>%
    arrange(Share) %>%
    mutate(Sector = factor(Sector,levels = Sector),x_id = row_number()) %>%
    pivot_longer(cols = c(Current_US,Current_NA,Expanded_US,Expanded_NA),names_to = "Scenario_code",values_to = "Impact") %>%
    mutate(
      Scenario = case_when(
        Scenario_code == "Current_US" ~ "Current: U.S. only",
        Scenario_code == "Current_NA" ~ "Current: North America",
        Scenario_code == "Expanded_US" ~ "Expanded: U.S. only",
        Scenario_code == "Expanded_NA" ~ "Expanded: North America"
      ),
      Scenario = factor(Scenario,levels = scenario_levels),
      x_plot = x_id + unname(scenario_offsets[as.character(Scenario)]),
      PointFill = case_when(
        Scenario_code %in% c("Current_US","Expanded_US") ~ as.character(Sector),
        Scenario_code %in% c("Current_NA","Expanded_NA") ~ "Hollow"
      )
    )
}
GHG_rng <- make_emission_panel_df(GHG_obesity_us,GHG_obesity_na,GHG_overweight_us,GHG_overweight_na,"Total GHG emissions","GHG_share")
CO2_rng <- make_emission_panel_df(CO2_obesity_us,CO2_obesity_na,CO2_overweight_us,CO2_overweight_na,"Total CO2 emissions","CO2_share")

# ======================================================
# 6. Plot function
# Original grid retained
# Bubble center = Obesity, Shock US
# Open circle = Obesity, Shock NA
# Short bars = Overweight Shock US and Shock NA
# Vertical line links the two overweight bars
# Bubble size = emission-type-specific sector share
# Sector order = small to large sector share within each plot
# Sector color legend order = Total GHG sector share, small to large
# ======================================================

make_emission_plot <- function(df_rng,y_min,y_max) {
  sector_order_key <- levels(df_rng$Sector)
  ggplot(df_rng) +
    geom_hline(yintercept = 0,linewidth = 1,color = "black") +
    geom_point(aes(x = x_plot,y = Impact,size = Share,color = Sector,fill = PointFill,shape = Scenario),alpha = 0.90,stroke = 1.7) +
    scale_x_continuous(breaks = seq_along(sector_order_key),labels = NULL,expand = expansion(mult = c(0.05,0.05))) +
    coord_cartesian(ylim = c(y_min,y_max),clip = "off") +
    theme(aspect.ratio = 1) +
    scale_shape_manual(
      values = scenario_shapes,
      breaks = scenario_levels,
      name = "GLP-1 adoption scenario",
      guide = guide_legend(
        order = 1,
        nrow = 2,
        byrow = TRUE,
        override.aes = list(size = 8,color = "black",fill = c("black",NA,"black",NA),stroke = 1.7,alpha = 0.90)
      )
    ) +
    scale_color_manual(
      values = sector_colors,
      breaks = sector_legend_order,
      limits = sector_legend_order,
      name = "Sector",
      guide = guide_legend(
        order = 2,
        nrow = 2,
        byrow = TRUE,
        override.aes = list(size = 10,shape = 21,fill = unname(sector_colors[sector_legend_order]),alpha = 0.90,stroke = 1)
      )
    ) +
    scale_fill_manual(values = c(sector_colors,"Hollow" = NA),guide = "none") +
    scale_size_area(
      max_size = 26,
      limits = c(0,0.020),
      breaks = c(0.0025, 0.0050,0.0100,0.0150),
      labels = c("0.25%","0.50%","1.00%","1.50%"),
      name = "Sector share (2030)",
      guide = guide_legend(order = 3,nrow = 1,override.aes = list(shape = 21,color = "grey30",fill = "grey30",alpha = 0.90,stroke = 1))
    ) +
    labs(x = NULL,y = NULL) +
    theme_minimal(base_size = 40) +
    theme(axis.text.x = element_blank(),axis.ticks.x = element_blank(),axis.title.x = element_blank(),axis.text.y = element_text(size = 32,color = "black"),panel.grid.major.x = element_line(linewidth = 0.5,color = "grey88"),panel.grid.minor.x = element_blank(),legend.key.width = unit(1.7,"lines"),legend.key.height = unit(1.7,"lines"),legend.text = element_text(size = 32),legend.title = element_text(size = 34,face = "bold"),plot.margin = margin(10,10,10,10))
}

# ======================================================
# 7. Panel wrapper with grey label bar
# ======================================================


panel_emission <- function(plot_obj,label_text) {
  label_plot <- ggplot(data.frame(x = 0.5,y = 0.5),aes(x,y)) +
    geom_rect(aes(xmin = 0,xmax = 1,ymin = 0,ymax = 1),fill = "grey90",color = NA) +
    geom_text(label = label_text,size = 12,fontface = "bold",hjust = 0.5,vjust = 0.5) +
    scale_x_continuous(limits = c(0,1),expand = c(0,0)) +
    scale_y_continuous(limits = c(0,1),expand = c(0,0)) +
    theme_void()
  plot_obj / label_plot + plot_layout(heights = c(12,1))
}
p_GHG <- panel_emission(make_emission_plot(GHG_rng,y_min = -10,y_max = 2),paste0("Total GHG Emissions, ",Year_used))
p_CO2 <- panel_emission(make_emission_plot(CO2_rng,y_min = -10,y_max = 2),paste0("Total CO2 Emissions, ",Year_used))
panel_grid <- (p_GHG | p_CO2) + plot_layout(guides = "collect") & theme(legend.position = "bottom",legend.box = "vertical",legend.box.spacing = unit(0.6,"cm"))
y_label <- ggplot() + annotate("text",x = 0.5,y = 0.5,label = "Sectoral Emission Impact (%)",angle = 90,size = 14,fontface = "bold") + scale_x_continuous(limits = c(0,1),expand = c(0,0)) + scale_y_continuous(limits = c(0,1),expand = c(0,0)) + theme_void()
figure_emissions <- plot_grid(y_label,panel_grid,ncol = 2,rel_widths = c(0.04,0.96))
figure_emissions